Efficient Deployment of Multi-UAVs in Massively Crowded Events †

In this paper, the efficient 3D placement of UAV as an aerial base station in providing wireless coverage for users in a small and large coverage area is investigated. In the case of providing wireless coverage for outdoor and indoor users in a small area, the Particle Swarm Optimization (PSO) and K-means with Ternary Search (KTS) algorithms are invoked to find an efficient 3D location of a single UAV with the objective of minimizing its required transmit power. It was observed that a single UAV at the 3D location found using the PSO algorithm requires less transmit power, by a factor of 1/5 compared to that when using the KTS algorithm. In the case of providing wireless coverage for users in three different shapes of a large coverage area, namely square, rectangle and circular regions, the problems of finding an efficient placement of multiple UAVs equipped with a directional antenna are formulated with the objective to maximize the coverage area and coverage density using the Circle Packing Theory (CPT). Then, the UAV efficient altitude placement is formulated with the objective of minimizing its required transmit power. It is observed that the large number of UAVs does not necessarily result in the maximum coverage density. Based on the simulation results, the deployment of 16, 19 and 26 UAVs is capable of providing the maximum coverage density of 78.5%, 82.5% and 80.3% for the case of a square region with the dimensions of 2 km × 2 km, a rectangle region with the dimensions of 6 km × 1.8 km and a circular region with the radius of 1.125 km, respectively. These observations are obtained when the UAVs are located at the optimum altitude, where the required transmit power for each UAV is reasonably small.


Introduction
Unmanned Aerial Vehicles (UAVs) have increasingly diverse and ubiquitous roles in today's society. They have been deployed in many civilian and military applications [1][2][3]. UAV can be used as an aerial base station and as a supplement to the existing ground base station when the cellular network is overloaded during a massively crowded special event or when the infrastructure of the terrestrial base stations is damaged due to natural disasters [4,5].
Al-Hajj (pilgrimage to Mecca) is an annual massively crowded event for the Muslims around the world. Every year, millions of pilgrims travel to the city of Mecca from all over the world for this event.
During the peak of the event, more than two million pilgrims perform a series of rituals in Mecca city within a one-week duration. The movement of a huge number of pilgrims occurred at the four holy sites of Mina, Arafat, Mozdalifah and the holy Masjid Al-Haram. As a result of the huge gathering of pilgrims, where large demand for communication is required, terrestrial base station overloading may happen, and sometimes, the ground base station could not provide simultaneous coverage to all pilgrims in this event [6,7]. Consequently, instead of building a cellular infrastructure network for temporary crowded events [8], UAVs can be utilized as a supplement to the terrestrial base station, in providing uninterrupted coverage [9][10][11][12][13].
Some of the key technical challenges in the deployment of UAVs as aerial base stations are the efficient 3D deployment of UAVs, power consumptions, wireless coverage optimization and interference management [3,10,14,15]. There has been an increase of research interest in the efficient 3D deployment of UAV strategies. This is because the UAV deployment problem has an impact on the power consumptions [3,12,16]. Thus, this motivates research interest in the efficient 3D deployment of UAV strategies that aim for minimum transmit power.

Related Works
The efficient 3D deployment of UAVs strategies can be categorized into two categories, namely providing wireless coverage for outdoor users [9][10][11]13,14] and providing wireless coverage for indoor users [16,17].
Two issues that must be taken into consideration in the deployment of UAV as an aerial base station are the relation between the coverage area and the altitude of UAV and the effects of urban environment on the performance of communications.
There are two effects of the UAV placement, namely the coverage performance in terms of the number of users within the coverage area, as well as the quality of the Air-To-Ground (ATG) links. ATG communication occurs in accordance to two main propagation groups, namely, Line-Of-Sight (LOS) and Non-Line-Of-Sight (NLOS) conditions. Several ATG channel models in a dense urban area have been studied in [9,18,19]. As presented in [20], the probability of LOS for ATG communication is a function of elevation angle and the average height of buildings. Furthermore, the characteristics of the ATG channel depends on the height of the aerial base stations due to the path loss and shadowing effects of obstacles, as discussed in [19].
In [9], the authors provided a generic statistical propagation radio model for predicting the ATG path loss between a Low Altitude Platform (LAP) and its corresponding ground terminal for LOS and NLOS connections based on statistical parameters for different environments of the target area. Furthermore, both issues were addressed in [9], namely the impact of the UAV placements on the number of users within the coverage area, as well as the effect of the environment on the occurrence of LOS. Many researches employ this path-loss model in addressing the UAV deployment problem [10][11][12][13]. However, the ATG path loss model [9] is suitable for outdoor ground terminals only.
The authors in [10] studied the impact of UAV's altitude position on the coverage performance in terms of coverage area and minimum transmit power. More specifically, the ATG path loss model developed by Al-Hourani et al. [9] was utilized to find the optimal altitude of UAV that maximized the UAV coverage and minimized its required transmit power.
The authors in [11] utilized the ATG model to find the optimal deployment of UAVs acting as flying base stations to minimize the total required UAV transmit power while satisfying the users' data rate requirements. In [12], the authors used this model to find the efficient UAV placement that maximized the network throughput. In [13], the study utilized this model to find the minimum number of UAVs and their 3D placement so that all the users were served. However, all previous studies assumed that all receivers are outdoors and located at 2D points. These assumptions limit the applicability of this model when we need to consider the indoor users.
Thus, a suitable path loss model needs to be studied, in the case of providing wireless coverage for indoor users. The authors in [17] utilized the outdoor-to-indoor path loss model certified by the International Telecommunication Union (ITU) [21] to provide coverage for indoor receivers. More specifically, in [17], the authors proposed an efficient 3D UAV placement that minimized the total transmit power required to provide wireless coverage for all indoor receivers, where the indoor receivers are symmetrically distributed.
In [3], the authors proposed an efficient 3D deployment of UAV using the Particle Swarm Optimization algorithm (PSO) to find an efficient 3D UAV location that minimized the total transmit power required to cover all indoor receivers where the indoor receivers were uniformly distributed inside the building. Due to the limited UAV transmit power, the authors in [16] minimized the number of UAVs required to cover the indoor receivers. However, the studies in [3,16,17] assumed that all receivers were located indoors. Furthermore, the above-mentioned strategies considered the 3D location of a single UAV and covering a small area of communication. As an example, in [3], it is assumed that all receivers are located in a building with the dimensions: horizontal building width x b is 20 m; vertical building width y b is 50 m; building heights are 200 m, 250 m and 300 m.
Thus, when considering providing wireless coverage over a larger coverage area, this leads to concerns about the number of UAVs that need to be employed. In this situation, the UAV deployment strategy that maximizes the coverage area can be used. In [14], the authors utilized Circle Packing Theory (CPT) to find an efficient placement of multiple UAVs acting as aerial base stations that maximizes the coverage area. However, this method is appropriate only for a circular coverage area.
Explicitly, this paper extends the work in [22] by proposing the UAV deployment strategy for a large coverage area. More specifically, in this paper, an efficient 3D deployment strategy is presented to find an efficient 3D placement of a single UAV that minimizes the transmit power. The performance of the 3D location of the UAV is analyzed in terms of the required transmit power when using two different algorithms, namely the PSO and K-means with Ternary Search (KTS) algorithms. Next, to provide wireless coverage over a large area, a 3D deployment strategy is proposed in a way that the total coverage over a large area is maximized using CPT. Thus, this leads to finding the number of UAVs equipped with a directional antenna needed to cover the area. The coverage performance is analyzed for three different shapes of area, namely square, rectangle and circular, in terms of coverage density.

Paper Contributions
The contributions of this paper are summarized as follows: • The existing ATG path loss model [9] and outdoor to indoor path loss model [21] are used to study the problem of a single UAV placement to provide coverage in crowded events for both outdoor and indoor receivers simultaneously, with the objective to minimize the required UAV transmit power. Due to the intractability of the formulated problem, two algorithms are developed to find an efficient 3D UAV placement using two optimization techniques, namely the PSO and KTS algorithms. The proposed algorithms consider the problem in providing wireless coverage for indoor and outdoor users, in a small area using a single UAV. • The efficient 3D placements of multiple UAVs that provide maximum wireless coverage and minimize the transmission power are found for each UAV. • The CPT is utilized to find the number of UAVs needed for providing wireless coverage for outdoor users in a large coverage area having three different shapes of coverage area, namely square, rectangle and circular. The problem is formulated with the objective to maximize the wireless coverage area using multiple UAVs. In each subarea, the UAV altitude is optimized using the algorithm to provide wireless coverage using a single UAV above.
The rest of this paper is organized as follows. In Section 2, the case of providing wireless coverage for outdoor and indoor users using a single UAV is studied, where the objective is to minimize the UAV transmit power. Section 3 presents the study of the problem in providing wireless coverage using multiple UAVs equipped with directional antennas, where the objective is to maximize the wireless coverage over a large area. Three different shapes of 2D region are considered, with each UAV employing a directional antenna, which produces a circular coverage pattern. Each circular cell is referred to as a subarea. Section 4 presents the simulation results and analysis for the deployment of the single UAV and multiple UAVs scenarios. Section 4.3 discusses the main observations obtained from the simulation results. The conclusion and future work are presented in Section 5.

System Model
Consider a coverage area that is divided into K subareas. In each subarea, there are N outdoor receivers and M indoor receivers, which are non-uniformly distributed using a beta random distribution, denoted as function f (x, y). Let [x 1 ,y 1 ] × [x 2 ,y 2 ] denote the dimension of the subarea. Then, the location of each outdoor receiver i out is denoted by (x i out , y i out , 0); the location of indoor receiver i in is denoted by (x i in , y i in , z i in ); and the 3D UAV location is represented by (x U AV , y U AV , z U AV ). In this model, it is considered that a single UAV serves as an aerial base station, as shown in Figure 1.

ATG Path Loss Model
In this system, the ATG path loss prediction between an LAP and ground terminals presented in [9,23] is utilized. More specifically, the ATG path loss is modeled by considering the probability of Line Of Sight (LOS) and Non-Line Of Sight (NLOS) links, which depend on the urban environment parameters, building heights and elevation angle. The probability of LOS and NLOS links can be given as [9]: where a and b are constant values that depend on the environment; such as urban, suburban, dense urban and rural environments; whilst the elevation angle θ out is given as θ out = sin −1 (h/d), where the distance between UAV and the ground receiver is d= h 2 + r 2 d and h is the UAV altitude. r d denotes the distance between UAV projection at the xy plane (x U AV , y U AV ) and the ground receiver coordinates located at (x i ,y i ), which is given as r d = (x U AV − x i ) 2 + (y U AV − y i ) 2 . The average path lossL (dB) for LOS and NLOS links can be formulated as [9]: where the path loss for the LOS link, L LOS , and the path loss for the NLOS link, L NLOS , are given as: where f c is the carrier frequency, c is the speed of the light, whilst η LOS and η NLOS are additional losses, which depend on the environment for LOS and NLOS links, respectively.

Outdoor-to-Indoor Path Loss Model
Furthermore, in order to provide the wireless coverage for indoor users, the path loss model for outdoor-to-indoor certified by ITU is utilized. The path loss is given as [21]: where L FSP is the free space path loss, L BP is the penetration loss for the building and L I N is the loss inside the building. The free space path loss L FSP is given as L FSP = 20log(d 3D ) + 20log( f c(GHz) ) + 32.4, where d 3D is the distance between an indoor receiver i in and the UAV, and the carrier frequency f c(GHz) is equal to 2 GHz.
The building penetration loss L BP is given as L BP = 14 + 15 (1 − cosθ in ) 2 , where the incident angle is denoted as θ in and the loss inside the building L I N is given by L I N = 0.5 d 2D in , where d 2D in is the 2D indoor distance between UAV and receiver i in located at (x i in , y i in , z i in ). More specifically, d 2D in denotes the indoor part of the 2D distance between the UAV and the indoor receiver, i in .
The building penetration loss L BP depends on the altitude of the indoor receiver z i in ; whilst the free space path loss L FSP is related to the distance between an indoor receiver i in and the UAV d 3D , which depends on the altitude of the indoor receiver z i in . Thus, the altitude of the indoor receiver z i in affects the path loss in terms of L FSP and L BP .

Problem Formulation
Consider a single UAV that is located at (x U AV , y U AV , z U AV ) transmitting data to the outdoor receiver i out located at (x i out , y i out , 0) or to the indoor receiver i in located at (x i in , y i in , z i in ). The data rate for each receiver i that is located either indoors or outdoors is given by: where C is the receiver bandwidth, P t,i is the transmit power of UAV for receiver i, L i is the path loss between UAV and receiver, i, and N p is the noise power. In this work, the interference is implicitly modeled as noise. In this system model, Frequency Division Multiple Access (FDMA) is used as the channel access technique. It is assumed that each UAV allocates equal channel bandwidth to receivers, and in order to avoid interference, each channel is assigned to one receiver.
In this problem, it is considered that the UAV transmits data to N outdoor and M indoor receivers at a desired data rate (r). (M + N) is the total number of receivers inside the coverage subarea, and each receiver has a channel with the bandwidth equal to A/(M + N), where A is the transmission bandwidth of an UAV. The total required transmit power of UAV to satisfy the data rate r for all receivers can be formulated as: The objective is to find an efficient 3D placement of UAV, such that the total required transmit power of UAV is minimized, and at the same time, the desired data rate of all receivers is satisfied. However, the optimal UAV 3D placement that minimizes the total transmit power becomes complicated, as the optimal solutions can be obtained by searching all possible sets of solutions [13], while satisfying the data rate for all receivers. Thus, the efficient solution can be found using the meta-heuristic algorithm.

Efficient UAV 3D Placement Algorithms
Due to the intractability of this problem, two algorithms are proposed to find an efficient 3D placement of a single UAV using two optimization techniques, namely the PSO and KTS algorithms. In PSO, a member of the set of possible solutions is known as a candidate solution, which is referred to as a particle. PSO utilizes the number of particles, and each particle moves around in the search space looking for the best solution. On the other hand, KTS includes two algorithms, namely the K-means algorithm and the Ternary Search (TS) algorithm. K-means is a clustering algorithm that attempts to split a given dataset into a fixed number of k clusters; whilst the ternary search employs the divide and conquer algorithm used to find the minimum or maximum values of a function.

Particle Swarm Optimization Algorithm (PSO)
In 1995, the PSO algorithm was introduced by Kennedy and Eberhart [24]. PSO is an intelligent optimization algorithm that is based on the swarm intelligence paradigm inspired by animals' social behavior such as a flock of birds or a school of fish. In PSO, the optimization problem is solved by iteratively trying to improve a candidate solution based on the local best experience for each particle and the global best experience for all candidates. In PSO, every particle is a candidate solution to the optimization problem.
Algorithm 1 shows the pseudocode of the proposed algorithm using the PSO optimization technique.
Step 1 in Algorithm 1 presents the inputs of the algorithm, namely N_pop defines the population of candidate solutions of the algorithm, W refers to the inertia weight, while r 1 and r 2 denote the two random numbers uniformly distributed randomly in a range between zero and one and c 1 and c 2 are the acceleration coefficients.
In the initialization step, the values of the constriction factor, κ, cognitive parameter, φ 1 , and the social parameter, φ 2 , must be selected, where κ = 1 and φ 1 + φ 2 > 4 [25]. This will help the PSO to find the efficient solution for the optimization problem.
The PSO is initialized with a group of random solutions for all particles' positions, Lo, and the particles' velocity, V, as shown in Steps 7 to 14. Then, in every iteration, for each particle, the best local location, L_BL, and the velocity are updated according to Equations (7) and (8), respectively. Moreover, the global best location, G_BL, is also updated, as in Steps 15 to 25. This will help the swarm of particles to move toward the best solution.

K-Means with Ternary Search Algorithms
K-means is an algorithm that uses an iterative refinement technique to solve the clustering problem. This is accomplished by partitioning a given dataset into k clusters. In this algorithm, each cluster is represented by its centroid. Thus, there are k centroids for k clusters. Each point in the dataset is assigned to a cluster that has the nearest centroid. More specifically, the K-means algorithm can be performed in the following stages:

1.
Initially, random guesses for cluster centroids are made, as shown in Step 3 of Algorithm 2.

2.
The nearest centroid is determined for each data point, by calculating the Euclidean distance between each point and the centroid of the cluster, as shown in Step 5 of Algorithm 2.

3.
In each cluster, the centroid is replaced by a new value. This new value is the means of the points belonging to the cluster, as in Step 6 of Algorithm 2.

4.
Repeat the process in Items 2 and 3 above, until the solution converges. The convergence happens when the centroids and their locations are no longer changed; more specifically, when the cluster mean is not changed as in Steps 4 to 6 of Algorithm 2.
Vmax: Upper bound decision variable.
(κ, φ 1 , φ 2 ): Construction coefficients. 2: Initialization: Moreover, the TS algorithm is used in tandem with the K-means algorithm, to determine the position of a specific value in a dataset. The sorted dataset is divided into three parts, and after that, the ternary search determines in which part the element exists. In Algorithm 2, Step 7 initiates the min and max variables for the left and right edge of the intervals, respectively. Then, in Step 10, the TS algorithm will repeat Steps 11 to 16 while max − min < . In this work, K-means is used to find the 2D placement (x U AV ,y U AV ) of UAV, since the outdoor and indoor receivers in the coverage area are assumed to be non-uniformly distributed. Then, the ternary search algorithm is used to find the UAV altitude (z U AV ). Algorithm 2 shows the pseudocode of the KTS algorithm.
Algorithm 2 K-means with ternary search algorithm.

Providing Wireless Coverage Using Multiple UAVs Equipped with Directional Antennas
In this section, an efficient deployment of multiple UAVs equipped with directional antennas is proposed. More specifically, the UAV deployment strategy to provide wireless coverage for receivers distributed over a large area of square, rectangle and circular 2D regions is investigated. Each UAV employs a directional antenna, which forms a subarea having a circular coverage pattern, referred to as the circle cell. Therefore, the CPT [26] is utilized to divide the total coverage area into non-overlapped circle cells, such that the coverage area is maximized. In order to avoid the interference between contiguous cells, the non-overlapped constraint is considered. Then, for each circle cell, an efficient 3D UAV placement that minimizes the total UAV transmit power required to cover all receivers within the cell is found. Furthermore, in this work, it is assumed that the distribution of the users is not given. Thus, the CPT is more suitable to be utilized in finding the UAV 2D placement, when compared with the K-means algorithm. Moreover, in this work, an efficient deployment of multiple UAVs is determined, such that the coverage area and density are maximized. However, the proposed algorithm presented in Section 2 does not consider this objective.

Case of a Square Region
Gatherings of people in some crowded events can form a square shape, such as the gathering of pilgrims in Arafat. The optimal packing of n identical and non-overlapped circles in a unit square is a widely-explored problem in the literature [27][28][29]. Many researchers have applied the CPT to find the maximum radius of n equal circles that can be packed in a unit square without overlapping.
In this scenario, it is assumed that the ground receivers are located in a geographical area that has a square shape and n UAVs can be deployed to maximize the total wireless coverage and density. More specifically, CPT is utilized to divide the total coverage area into non-overlapped circle cells. Then, for each cell, one UAV can be deployed to provide wireless coverage for receivers within the cell.

Problem Formulation
The packing of equal circles in a square can be formulated in the following way: • P1 : Place n identical non-overlapping circles in a unit square, with the objective function to maximize the radius of the circles r, such that the coverage area and coverage density are maximized.
The problem P1 is a continuous, nonlinear, inequality-constrained global optimization problem [30]. The problem P1 can be formulated as: where (x i , y i ) denotes the center coordinates of the circle i, (x i − x j ) 2 + (y i − y j ) 2 is the Euclidean distance between the centers of circles i and j, 1 ≤ i < j ≤ n and ris the coverage radius of each cell. The objective function is to maximize the radius of the coverage cell. The constraint equations of (9.a) and (9.b) ensure that all packed circles lie inside the square. The constraint equation of (9.c) guarantees no overlapping between circles. There has been a number of optimum solutions proposed for the problem of placing n identical non-overlapping circles in a unit square, such that the circle radius is maximized [27][28][29]. More specifically, in [27][28][29], the optimum arrangement of n identical non-overlapping circles in a unit square was solved by maximizing the distance m between any pairs of circles.
For a unit square, the density of packing n identical non-overlapping circles d n can be defined as the ratio of area occupied by the packed circles to a unit square area. Thus, the density of the packed circles in a unit square is related to the maximum radius of the packed circle r n by the following equation [31]: Table 1 shows the optimal results of the maximum radius of the packed circle r n and the corresponding maximum density d n of packing equal circles in a unit square for 2 ≤ n ≤ 22. In this work, the coverage density of packing n identical non-overlapping circles d n is defined as the ratio of the area occupied by the packed circles to the area of the coverage region.

Case of a Rectangle Region
There are many crowded events where the users can form a rectangular shape, for example the movement of a huge number of attendees through streets or rectangular paths, such as the movement of pilgrims in Mina. The packing of n identical circles into a rectangle can be efficiently used to divide the rectangle region into equal circle cells, and for each cell, one UAV can be deployed to provide wireless coverage for users within the cell. This problem is known as an NP-hard problem, so finding the optimal packing solution is difficult and makes the problem very complicated [26,35,36]. Therefore, meta-heuristics approaches or local exhausted search methods can be used to solve this problem and find an efficient packing solution.
In [37], the authors proposed a heuristic algorithm to find the maximum radius of a specified number of non-overlapped equal circles inside a fixed size rectangle. They used the Formulation Space Search (FSS) approach to solve the circle packing problem.

Problem Formulation
The problem of packing equal circles in a rectangle region can be formulated as follows: • P2: Place n identical non-overlapping circles in a rectangle region L × W, with the Cartesian origin (0, 0) as the rectangle center. The objective function is to maximize the radius of the circles r such that the coverage area and density are maximized.
In problem P2, R overlap is the upper bound radius of the packed circles, and it can be defined by nπR overlap = LW, then R overlap = √ LW/nπ. A mixed Cartesian and polar formulation has been used. Hence, problem P2 is a mixed nonlinear formulation and can be formulated as [37]: where (x i , y i ) is the circle i center. The first four constraints ensure that all packed circles are located inside the rectangle. (x i − x j ) 2 + (y i − y j ) 2 is the Euclidean distance between the centers of circles i and j, ∀ 1 ≤ i < j ≤ n, and this distance must be at least 2r to ensure that there is no overlap between circles. The objective of this problem is to maximize the radius of n identical non-overlapping circle cells. The constraints of Equations (11.a), (11.b), (11.c) and (11.d) ensure that all packed circles are located inside the coverage area. Moreover, the constraints of Equations (11.e) and (11. f ) guarantee no overlapping between cells. Similarly, the coverage density of packing n identical non-overlapping circles d n is defined as the ratio of area occupied by the packed circles to the area of the coverage region, as discussed in Section 3.1.

Algorithms for Packing Circles in a Rectangle Region
Many algorithms were developed for solving the problem of packing equal circles in a rectangle region, using meta-heuristic approaches [36][37][38][39][40]. In [36], two algorithms were developed using the meta-heuristic simulated annealing optimization technique to solve the packing problem of circular-based items into a rectangular pallet, such that the number of circular-based items being packed in the rectangular pallet was maximized without overlaps. These algorithms were referred to as the Cylinder Packing Algorithm using Minimizing Overlap with Simulated Annealing (CPA-MinOSA) and the Cylinder Packing Algorithm using Maximizing the Number of circles with Simulated Annealing (CPA-MaxNSA). The pseudocode of the CPA-MinOSA algorithm is shown in Algorithm 3.
In [37], a heuristic algorithm was proposed to solve the problem of packing n identical non-overlapping circles into a 2D container of fixed size, such that the circle radius was maximized. This algorithm was developed based on the Formulation Space Search (FSS) method. Algorithm 4 illustrates the FSS pseudocode.
Moreover, the Packomania website [40] used heuristic algorithms to find the best known solutions for packing identical circles in fixed size shapes, such that the circle radius and packing density were maximized. give (x, y), x 0 (X − δ,

Case of a Circular Region
In this section, a circular coverage region is considered. The CPT of packing equal circles into a circular region is utilized, to provide wireless coverage using UAVs for receivers within the circular coverage region. Here, the problem to find the maximum radius of equal circles that can be packed into a unit circle with radius of r = 1 is discussed. Moreover, the mathematical formulation and the main algorithm used to solve this problem is presented [37,41,42] in this section.

Problem Formulation
The problem of packing equal circles into a unit circular region can be formulated as: • P3: Place n identical non-overlapping circles into a unit circle with radius of r = 1. The objective function is to maximize the radius of the packed circles, such that the coverage area and density are maximized.
The problem P3 can be formulated as: where (x i , y i ) is the center coordinates of a circle i, and the constraint equation of (12.a) is the Euclidean distance between the center of circles i and j, which is defined as (x i − x j ) 2 + (y i − y j ) 2 and must be ≥2r to ensure no overlapping between circles. The second constraint equation of (12.b) guarantees that the packed circles must fully lie inside the circular region.
Here, the coverage density of packing n identical non-overlapping circles d n is defined as the ratio of area occupied by the packed circles to the area of the coverage region, as discussed in Section 3.1.

Algorithms for Packing Circles into a Circular Region
The authors in [37] presented the formulation of the circle packing problem inside a fixed size container with the objective to maximize the radius of the packed circles. To solve this problem, they used the FSS algorithm with the Sparse Nonlinear OPTimizer solver (SNOPT), as discussed in Section 3.2.2.
In [41], the authors found the densest packing of identical circles inside a unit circle using the general Reformulation Descent (RD) heuristic to maximize the radius of the packed circles into a fixed size container. They proposed nonlinear reformulation to transform one coordinate system to the other, such as switching from Cartesian to polar coordinates and vice versa. This solution was invoked in order to solve local search using the gradient method when it reached a stationary point in the Nonlinear Programming problem (NLP). It was observed that this technique allowed the NLP-solver to find a better solution. Moreover, the Packomania website [40] shows the best known solutions for packing n identical circles inside a unit circle with the maximum common circle radius and maximum packing density.
In addition to the previous convex shapes, CPT can be utilized in the problem of providing wireless coverage for crowded events that have other convex shapes such as an equilateral triangle [43] and a polygon [44]. Moreover, CPT can be also utilized for the problem of providing wireless coverage for crowded events that have non-convex 2D shapes [45].

Providing Wireless Coverage Using a Single UAV
This section presents the simulation results of the proposed algorithms to find an efficient 3D placement of a single UAV for providing wireless coverage for outdoor and indoor users, such that the transmit power is minimized.
In this simulation, an area in Mena city is considered, which is divided into several subareas. The PSO and KTS algorithms are invoked to find the efficient placement of a single UAV for providing wireless coverage for outdoor and indoor users in a rectangle subarea, with the dimensions of 300 m × 150 m. More specifically, for each subarea, it is assumed that 35% of the total number of receivers inside the subarea is active users.
The simulations and analysis of the simulation results are performed using MATLAB R . Moreover, all simulation parameters are listed in Table 2. In this scenario, it is considered that there are two buildings in the subarea, and each building consists of 12 floors. The number of outdoor receivers is 2750, where it is assumed that 35% of the total number of receivers is active users. On the other hand, three different numbers of indoor receivers, namely 600, 750 and 960, are considered. For the case of having 960 indoor users, it is assumed that there are 40 active receivers on each floor. Both outdoor and indoor receivers are non-uniformly distributed based on the beta random distribution function, f (x, y). The total number of active outdoor and indoor receivers is presented in Table 3. For each receiver, the data rate r is 0.5 Mbps.  Figure 2a illustrates the distribution of the outdoor receivers and the coordinates of the two buildings inside the subarea; whilst Figure 2b illustrates the 3D view of the distribution of the indoor and outdoor receivers.
The simulation results of the proposed algorithms, when varying the number of indoor active users, are shown in Table 3. More specifically, this table compares the minimum required total UAV transmit power that satisfies the data rate of all receivers, when invoking the PSO and KTS algorithms. The results show that the proposed algorithm using PSO requires less transmit power when the UAV is at the 3D efficient location; more specifically, for the case of having the greatest number of indoor receivers, 960. The 3D location of UAV using the PSO algorithm is (161.66, 62.79, 62.12), as shown in Figure 2a. This requires a total UAV transmit power of 3.255 watts, as illustrated in Figure 3. On the other hand, the 3D location of UAV using the KTS algorithm is (162.64, 63. 23, 41.11), and the required UAV transmit power is 16.594 watts, as we can see in Figure 2b. Thus, the proposed algorithm using PSO exhibits about five-times improvement in terms of minimum required transmit power, when compared with that using the KTS algorithm. This is because, in PSO, the efficient 3D UAV placement is determined without dividing the problem into subproblems. On the other hand, in KTS, the solution to find an efficient 3D placement of UAV is divided into two stages of solving two subproblems. In the first subproblem, K-means is used to find the 2D placement of UAV (x U AV , y U AV ). Then, in the second subproblem, the ternary search algorithm is used to find the altitude (z U AV ) of UAV.

Providing Wireless Coverage Using Multiple UAVs
This section presents the simulation results of the proposed algorithms for the efficient deployment strategy of multiple UAVs equipped with directional antennas, to provide wireless coverage for receivers distributed over a large area of a square, rectangle and circle. The simulations and analysis of the simulation results are performed using MATLAB R .

Case of a Square Region
In this scenario, the case of providing wireless coverage for gatherings of pilgrims in Arafat city that forms a square shape is considered, as shown in Figure 4. The square region has dimensions of 2 km × 2 km.
In this simulation, the CPT and the results in Table 1 are utilized to find the maximum radius of n identical circles that can be packed inside the targeted coverage area of a square shape. Each circle is referred to as a subarea. Then, the 3D placement for UAVs that minimizes the total UAV transmit power required to provide wireless coverage for all outdoor receivers within the circle subarea is determined. For each circle cell, it is assumed that 35% of the total number of receivers inside the coverage circle is active. All simulation parameters are listed in Table 2. Table 4 presents the simulation results of the proposed algorithms to pack n equal circles inside a 2 km × 2 km square region, for n = 8 to n = 22. More specifically, this table presents the maximum radius of the circle, r, such that the coverage area and density are maximized for each case of packing n identical circles in the square region. Consequently, the corresponding optimal 3D placement of UAVs that minimizes the total transmit power required to provide wireless coverage for all receivers inside the circle is presented in Table 4. The center coordinate of the circle cell is denoted as (x i , y i ), which refers to UAV 2D location, and z i is the optimal UAV altitude. The maximum packing density of the coverage region is also evaluated. The directional antenna half beamwidth θ i /2 is defined as tan −1 (r i /h i ), where r i is the radius of circle i and h i is the optimal UAV altitude. Table 4. Simulation results for packing n circles in a 2 km × 2 km square region (x i ,y i ) I = 1, . . ., n. It can be seen from Table 4 that the highest achievable density of the packed circles in a square occurs when n = 9 and n = 16, when the proposed algorithm is invoked. However, the minimum UAV transmit power when packing n = 9 circles is very high.
For the case of packing n = 16 circles, the maximum circle radius r is 250 m. For this case, the distribution of 2745 active receivers inside the subarea is illustrated in Figure 5. It can be seen from Figure 6a that the corresponding optimal UAV altitude is at 156 m, when the minimum transmit power is at the minimum.  Thus, the distribution of the outdoor receivers and the corresponding optimal 3D placement of 16 UAVs can be illustrated as shown in Figure 6b. Each UAV provides wireless coverage for one circle cell (subarea) with a radius of 250 m. The center of each subarea (x i ,y i ) i = 1, . . . , n is the (x U AV , y U AV ) location, whilst the optimal altitude of all UAVs is at z U AV = 156 m. At this location, the UAV required transmit power for each UAV i 1, . . . , 16 is equal to 20 mwatt, and the maximum coverage density d n is 78.5%. Specifically

Case of a Rectangle Region
In this scenario, the case of providing wireless coverage for gatherings of pilgrims that forms a rectangle shape is considered, which is located in the area between Arafat city and Mena city, as shown in Figure 8. The rectangle region has the dimensions of 6 km × 1.8 km. More specifically, this simulation finds the maximum radius of n identical non-overlapping circles that can be packed inside the targeted coverage area of a rectangle shape. Similar to the case of packing circles into a square region, each circle is referred to as a subarea. Then, for each subarea, the optimal 3D UAV placement is determined. The objective is to minimize the total transmit power required to provide wireless coverage for all receivers inside the circle cell, such that the packing density is maximized. The coverage area contains roads and paths for users' movement, and there are areas where nobody crosses. Therefore, for each circle cell, it is assumed that 25% of the total number of users inside the coverage circle is active.  Table 5 presents the simulation results for packing n identical circles inside a 6 km × 1.8 km rectangular region, for n = 10 to n = 33. The 3D placement of UAV that minimizes the total transmit power required to cover all active users in each circle cell is determined. The center coordinate of the circle cell is denoted as (x i , y i ), which refers to UAV 2D location, and z i is the optimal UAV altitude. The maximum packing density of the coverage region and the antenna half beamwidth for each cell is also evaluated. It can be seen from Table 5 that the highest achievable density of the packed circles in a rectangle occurs when n = 26 and the proposed algorithm is invoked.  Figure 9a illustrates the distribution of the outdoor receivers and the corresponding optimal 3D placement of 26 UAVs. Each UAV provides wireless coverage for one circle cell (subarea) with a radius of 330 m. The center of each subarea (x i ,y i ) i = 1, . . ., n is the (x U AV , y U AV ) location, whilst the optimal altitude of all UAVs is at z U AV = 206 m. At this location, the UAV required transmit power for each UAV i 1, . . . , 26 is equal to 0.246 watts, and the maximum coverage density d n is 82.5%. More specifically, the efficient 3D placements

Case of a Circular Region
In this scenario, the case of providing wireless coverage for gatherings of pilgrims in Masjid Al-Haram and its outer courtyard is considered, which forms a circular shape, with a radius of r = 1.125 km, as shown in Figure 10. This region is divided into identical circles, and each circle (subarea) is served by one UAV. This simulation is performed to find the maximum radius of n identical non-overlapping circles that can be packed inside the targeted area of the circular region. Then, the 3D UAV placement that minimizes the total transmit power required to provide wireless coverage for active receivers inside the circle cell for each circle is determined. For each circle cell, it is assumed that 35% of the total number of receivers inside the circular region is active. The simulation results for packing n identical circles inside a circular region with r = 1.125 km are presented in Table 6. Then, the efficient 3D placement of the UAV that minimizes the total transmit power required to cover all active users inside the circle cell (subarea) for each circle cell is found. The center coordinate of the circle cell is denoted as (x i , y i ), which refers to UAV 2D location, and z i is the optimal UAV altitude. The maximum packing density of the coverage region and the antenna half beamwidth for each case of packing n identical circles is also evaluated. It can be seen from Figure 11 that the best coverage density is 80.3% when n = 19. For the case of packing n = 19 circles in the circular region, it can be observed from Table 6 that the required transmit power is 13 mwatt for each UAV to provide coverage for all active receivers inside the circle cell with a radius of r = 231 m; whilst the optimal 3D placement of 19 UAVs is at 144 m.

Discussion
In this section, the main observations obtained from the simulation results presented in Sections 4.1 and 4.2 are discussed.
It is observed that the UAV required transmit power depends on the number of active receivers within the coverage subarea and, hence, affects the altitude of the UAV. More specifically, in the case of finding the efficient 3D placement of a single UAV in providing wireless coverage for outdoor and indoor users, it is observed that as the number of indoor active users increases, the required transmit power increases, which causes the altitude of the UAV to increase, as presented in Table 3. Similar performance is observed for the case of providing wireless coverage using multiple UAVs. More specifically, for this case, as the number of circle cells (subareas) increases, the number of users inside each subarea decreases. Hence, the UAV required transmit power decreases, and consequently, the altitude of the UAV decreases, as observed in Tables 4-6.
With the specific observation from the case of providing wireless coverage using a single UAV, the problem is formulated to find the efficient 3D placement with the objective to minimize the UAV transmit power. In this case, the choice of the algorithm used to find the efficient 3D placement affects the UAV required transmit power. More specifically, it was observed that the PSO algorithm performed better than the KTS algorithm to find the 3D efficient placement of the UAV, where the 3D location found using the PSO algorithm requires less transmit power, by a factor of 1/5 compared to that when using the KTS algorithm, as discussed in Section 4.1.
Providing wireless coverage for outdoor users within a larger coverage area requires more than one UAV. Thus, it is proposed to utilize CPT in dividing the coverage area into subareas, such that the coverage area and density are maximized. Then, the optimal UAV altitude is determined with the objective of minimizing the UAV required transmit power. It is observed that the large number of UAVs does not necessarily results in the maximum coverage density.
As observed in Table 4, the maximum coverage density is obtained when the number of identical cells n = 4, 9 and 16 for the case of providing wireless coverage for users in a square region with the dimensions of 2 km × 2 km. However, the UAV required transmit power is very high when n = 4, 9. Thus, for this case, the deployment of 16 UAVs is capable of providing the maximum coverage density of 78.5%, with the optimal altitude of all UAVs at z U AV = 156 m. At this location, the required transmit power for each UAV is equal to 20 mwatt, which is not the minimum. The minimum transmit power of 0.124 mwatt is observed when n = 22, which is when the number of active users within the subarea is the smallest.
A similar observation is obtained from Tables 5 and 6, where there is a trade-off between the maximum coverage density and the UAV required transmit power. More specifically, having the maximum coverage density does not guarantee having the smallest value of the UAV required transmit power. However, the required transmit power is reasonably small.

Conclusions
In this paper, a UAV deployment strategy is proposed, for providing wireless coverage for users in small and large coverage areas of massively crowded events. More specifically, in a small coverage area, the problem was formulated to provide wireless coverage for outdoor and indoor users, with the objective to minimize the UAV transmit power. In this case, the PSO and KTS algorithms were used to find an efficient 3D UAV placement that minimizes the total required transmit power and satisfies the data rate for users. It was observed that the PSO algorithm performed better than the KTS algorithm to find the 3D efficient placement of UAV. More specifically, for a single UAV at the 3D location found using the PSO algorithm, the UAV requires less transmit power, by a factor of 1/5 compared to that when using the KTS algorithm.
For the deployment of UAV in providing wireless coverage for users in a large area of massively crowded events, CPT is utilized. More specifically, in a large coverage area, CPT is utilized to find the efficient 3D placements of multiple UAVs by packing identical, non-overlapping and interference free circle cells inside three different 2D shapes' coverage area, namely square, rectangle and circular regions. The problems were formulated with the objective to maximize the total coverage area and coverage density. Then, the efficient altitude placement of each UAV was found using the formulation that minimizes the transmit power in each circle cell.
For the case of a square region, it was found that the deployment of 16 UAVs at the optimum altitude of 156 m is capable of providing 78.5% coverage density of a square region with the dimensions of 2 km × 2 km, where each UAV required 20 mwatt. For the case of a rectangle region, the deployment of 26 UAVs at the optimum altitude of 206 m is capable of providing 82.5% coverage density of a rectangle region with the dimensions of 6 km × 1.8 km, where each UAV requires 0.246 watt. For the case of a circular region, the deployment of 19 UAVs at the optimum altitude of 144 m is capable of providing 80.3% coverage density of the circular region with the radius of 1.125 km, where each UAV requires 13 mwatt. However, it was also observed that an increase of the number of UAVs did not necessarily result in the maximum coverage density.
As future work, it is proposed to generalize the solution techniques to provide wireless coverage using UAVs by utilizing the packing of identical circles inside convex and non-convex regions. Moreover, it is also proposed to study the presence of interference between congruent cells in providing wireless coverage.
Author Contributions: A.S. and H.S. conceived of the idea of this work. A.S. formulated the problems, developed the system model and executed the algorithms. N.S.O. and H.S. analyzed the results. All authors contributed equally of the preparation and writing of the paper.

Funding:
The second author acknowledges the support of the Universiti Tenaga Nasional under UNIIG 2017 grant.

Conflicts of Interest:
The authors declare no conflict of interest.