Zhangyu Guan
email: zguan@ece.neu.edu

Nan Cen
email: ncen@ece.neu.edu

Tommaso Melodia
email: melodia@ece.neu.edu

Scott Pudlewski
email: scott.pudlewski.1@us.af.mil

Self-Organizing Flying Drones with Massive MIMO Networking

This article studies distributed algorithms to control self-organizing swarm drone hotspots with massive MIMO networking capabilities -a network scenario referred to as OrgSwarm. We attempt to answer the following fundamental question: what is the optimal way to provide spectrally-efficient wireless access to a multitude of ground nodes with mobile base stations/aerial relays mounted on a swarm of drones and endowed with a large number of antennas; when we can control the position of many-antenna-enabled drones, access association of ground nodes to drones, and the transmit power of ground nodes?

The article first derives a mathematical formulation of the problem of spectral efficiency maximization through joint control of the movement of many-antenna-enabled aerial drones, access association of single-antenna ground nodes to many-antenna drones, and transmit power of ground nodes. It is shown that the resulting network control problem is a mixed integer nonlinear nonconvex programming problem (MINLP). We then first design a distributed solution algorithm with polynomial computational complexity. Then, a centralized but globally optimal solution algorithm is designed based on a combination of the branch and bound framework and convex relaxation techniques to provide a performance benchmark for the distributed algorithm. Results indicate that the distributed algorithm achieves a network spectral efficiency very close (over 95% on average) to the global optimum.

I. INTRODUCTION

The proliferation of innovative Internet of Things (IoT) applications, including intelligent transportation, wireless video surveillance, augmented reality, among others, has resulted in increasing demand for faster wireless communication links with higher spectral efficiency; as well as techniques to reduce interference between co-located links operating on the same spectrum bands and hence to increase the spectral efficiency [START_REF] Noura | A Survey on Interference Management for Device-to-Device (D2D) Communication and Its Challenges in 5G Networks[END_REF]- [START_REF] Chen | Relay Location Optimization for Differential Amplify-and-Forward Cooperative Relaying[END_REF]. This paper studies OrgSwarm, a new approach to provide ground connectivity by exploring infrastructure-less self-organizing swarm drone hotspots with massive MIMO networking.

With OrgSwarm, a swarm of many-antenna-enabled aerial drones collaboratively provide data collection/forwarding services to a group of single-antenna ground ad hoc nodes, and send to ground nodes control commands that are generated either locally at the drones or in a remote fusion center. This work is based upon material supported in part by Air Force Research Laboratory under Contract FA8750-14-1-0074 and Grant ARMY W911NF-17-1-0034.

We attempt to provide an answer to the following question: What is the best way to provide spectrally-efficient wireless access/networking to a group of ground nodes with mobile base stations/aerial relays mounted on drones and endowed with a large number of antennas; when you can control the movement of drones, access association, and transmit power of ground nodes? We attempt to answer this compelling question by maximizing the spectral efficiency of OrgSwarm by jointly considering three tightly coupled network control strategies:

• Movement control of drones. Compared to infrastructurebased cellular networks with static massive MIMO base stations as in [START_REF] Björnson | Massive MIMO: Ten Myths and One Critical Question[END_REF], [START_REF] Wu | Secure Massive MIMO Transmission With an Active Eavesdropper[END_REF], a peculiar feature of OrgSwarm is that many-antenna-enabled aerial drones need to move dynamically to adapt to changes in the location or traffic demands of ground nodes; nodes leaving or joining the network; as well as time and spatially-varying interference level, among others, to provide coverage with higher spectral efficiency. • Access association of ground nodes to drones. This is jointly determined by several factors, including network topology (i.e., the relative locations of aerial drones and ground nodes), wireless channel quality currently experienced by ground nodes and transmit power of ground nodes, among others. Moreover, in a massive MIMO setting, a transmission is typically conducted in two phases: pilot-based channel estimation and data transmission. The maximum number of ground nodes associated to an aerial drone is constrained by the length (in symbols) of the pilot sequences used in the channel estimation phase [START_REF] Larsson | Massive MIMO For Next Generation Wireless Systems[END_REF], [START_REF] Marzetta | Massive MIMO: An Introduction[END_REF]. • Massive-MIMO power control. In OrgSwarm, all ground nodes are allowed to operate over the entire available spectrum band and at any transmission time. Therefore, it is imperative to regulate the transmit power of ground nodes to eliminate mutual interference in both channel estimation and data transmission phases caused by imperfect channel orthogonalization in the case of limited number antennas at each aerial drone, and hence to maximize the achievable SINR and spectral efficiency. Clearly, the above three network control problems are tightly coupled and should therefore be jointly considered to obtain the optimal network operating point. Because of the binary ground-drone association variables (which will be clearer ISBN 978-3-903176-05-8 c 2018 IFIP in Section III), the resulting OrgSwarm control problem is a mixed integer nonlinear nonconvex programming problem (MINLP), which is generally NP-hard and there is no existing solution algorithm that can be used to obtain the globally optimal solution with polynomial computational complexity. This paper makes the following main contributions:

• OrgSwarm control formulation. We formulate mathematically the OrgSwarm control problem with control objective of maximizing the spectral efficiency of OrgSwarm by jointly controlling the movement of aerial drones, ground-drone association and transmit power of ground nodes.

• Distributed solution algorithms. Distributed solution algorithms are designed to solve the formulated OrgSwarm control problem based on primal decomposition, which decomposes the resulting MINLP problem into three distributed sub-problems: single-antenna-to-many-antennas access association, massive MIMO power control, and drone movement. • Globally optimal solution algorithm. To provide a performance benchmark for the distributed solution algorithm, we design a centralized but globally optimal solution algorithm based on a combination of the branch and bound framework and of convex relaxation techniques that result in an ε-optimal solution of the original MINLP problem, where ε is a predefined level of optimality precision that can be set arbitrarily close to 1. • Performance evaluation. The performance of the proposed distributed solution algorithm is evaluated in terms of network spectral efficiency by comparing it to the global optimum through extensive simulation experiments. The reminder of the paper is organized as follows. We review related work in Section II, and describe the system model and problem formulation in Section III. In Sections IV and V, we present the distributed solution algorithm and the globally optimal solution algorithm, respectively. Performance evaluation is presented in Section VI, and finally we draw conclusions in Section VII.

II. RELATED WORK

There is a large and growing body of literature on unmanned aerial vehicular networking, focusing on UAV-assisted guidance [START_REF] Sharma | A Cooperative Network Framework for Multi-UAV Guided Ground Ad Hoc Networks[END_REF], UAV-based data collection [START_REF] Jawhar | A Framework for Using Unmanned Aerial Vehicles for Data Collection in Linear Wireless Sensor Networks[END_REF]- [START_REF] Daniel | Cognitive Agent Mobility for Aerial Sensor Networks[END_REF] and relaying [START_REF] Ono | A Wireless Relay Network Based on Unmanned Aircraft System With Rate Optimization[END_REF]- [START_REF] Ferranti | Drone Cellular Networks: Enhancing the Quality of Experience of Video Streaming Applications[END_REF], ground-aerial channel measurements [START_REF] Willink | Measurement and Characterization of Low-Altitude Air-to-Ground MIMO Channels[END_REF] as well as tracking and control of UAV networks [START_REF] Morgenthaler | UAVNet: A Mobile Wireless Mesh Network Using Unmanned Aerial Vehicles[END_REF]- [START_REF] Chandhar | Massive MIMO as Enabler for Communications with Drone Swarms[END_REF]. Readers are referred to [START_REF] Gupta | Survey of Important Issues in UAV Communication Networks[END_REF]- [START_REF] Zeng | Wireless Communications with Unmanned Aerial Vehicles: Opportunities and Challenges[END_REF] and references therein for an extensive survey of this research area. Most of these works focus on single-antenna aerial vehicles and conventional MIMO, with very few recent efforts considering massive MIMO [START_REF] Chandhar | Massive MIMO as Enabler for Communications with Drone Swarms[END_REF]. Different from [START_REF] Chandhar | Massive MIMO as Enabler for Communications with Drone Swarms[END_REF], where Chandhar et al. derived the achievable uplink capacity from a many-antennas ground base stations to a set of single-antenna aerial drones, in this paper we maximize the aggregate throughput of single-antenna ground nodes to a set of aerial drones each endowed with a large number of antennas.

Compared to conventional multiuser MIMO, massive MIMO can attain much higher spectral efficiency by using a large number of antennas with low-complexity linear precoding technologies [START_REF] Larsson | Massive MIMO For Next Generation Wireless Systems[END_REF], [START_REF] Marzetta | Noncooperative Cellular Wireless with Unlimited Numbers of Base Station Antennas[END_REF]- [START_REF] Zhang | Spectral and Energy Efficiency of Multipair Two-Way Full-Duplex Relay Systems With Massive MIMO[END_REF]. For example, in [START_REF] Ngo | Multipair Full-Duplex Relaying With Massive Arrays and Linear Processing[END_REF], the authors derived an exact achievable rate expression in closedform for maximum-ratio combining/maximum-ratio transmission (MRC/MRT) processing and an analytical approximation of the achievable rate for zero-forcing (ZF) processing for multi-pair full-duplex massive MIMO relay system. In [START_REF] Jin | Ergodic Rate Analysis for Multipair Massive MIMO Two-Way Relay Networks[END_REF], Jin et al. derived ergodic rates for the case with a finite number of antennas and concluded that the ergodic sum-rate can be maintained while the relay power is scaled down by a factor of the number of the antennas at the relay over the number of users. Please refer to [START_REF] Björnson | Massive MIMO: Ten Myths and One Critical Question[END_REF], [START_REF] Marzetta | Massive MIMO: An Introduction[END_REF], [START_REF] Araújo | Massive MIMO: Survey and Future Research Topics[END_REF] for a good survey and tutorial on massive MIMO networking. These papers are focused on infrastructure-based cellular networks with static many-antenna-enabled base stations, and focus on asymptotic performance analysis with respect to a single network parameter (e.g., power). Our paper, instead, considers aerial drone hotspots with massive MIMO capabilities in an infrastructure-less network scenario.

III. SYSTEM MODEL AND PROBLEM STATEMENT

We consider infrastructure-less wireless ad hoc networks with many-antenna-enabled drones, where a set A of manyantenna-enabled aerial drones are used to collect field information from a set G of single-antenna ground nodes. The goal of the aerial drones is to collect field information, make action decisions either locally at each drone or by sending the information fusion results to a remote control center, and finally send the action commands back to ground nodes. We focus on information uploading over ground-drone links since it causes only low-level traffic load to transmit information fusion results and control commands back to ground nodes. Our objective is to maximize the network-wide spectral efficiency by jointly controlling the movement of aerial drones, associating single-antenna ground nodes to many-antenna drones, and controlling the transmit power of ground nodes. Ground-Drone Association. We consider single-home accommodation for ground nodes in favor of tractable complexity in modeling and theoretical analysis, i.e., each ground node is associated with at most one drone. Denote the association vector as α (α g ) g∈G , where α g (α ga ) a∈A is the association vector of ground node g ∈ G, with α ga = 1 if ground node g is associated with drone a and α ga = 0 otherwise. Denote A g ⊂ A as the set of aerial drones ground node g is associated with, G a ⊂ G as the set of ground nodes associated with drone a, and A g and G a as the number of aerial drones in A g and the number of ground nodes in G a , respectively. Then, we have

A g {a|a ∈ A, α ga = 1}, ∀g ∈ G, ( 1
)
A g |A g | ≤ 1, ∀g ∈ G, ( 2
)
G a {g|g ∈ G, α ga = 1}, ∀a ∈ A, ( 3
)
G a |G a | ≤ G max , ∀a ∈ A, ( 4
)
where G max is the maximum number of ground nodes that can be served by each aerial drone a ∈ A. 1Ground-Drone Link Capacity. In massive MIMO networks, a transmission is typically accomplished in two phases, i.e., channel estimation and data transmission [START_REF] Larsson | Massive MIMO For Next Generation Wireless Systems[END_REF], [START_REF] Marzetta | Massive MIMO: An Introduction[END_REF], [START_REF] Lu | An Overview of Massive MIMO: Benefits and Challenges[END_REF]. In the first phase, ground nodes send a set of orthogonal pilot sequences to their service aerial drones for channel estimation. Denote τ as the length (in symbols) of each pilot sequence. For ground node g ∈ G, let I g represent the set of ground nodes using the same pilot sequences as node g, and a(g) denote the service aerial drone of ground node g. Let p = (p g ) g∈G represent the transmit power vector of ground nodes, with p g being the power of gth ground node in G. Then, the achievable signal-to-interference-plus-noise ratio (SINR) of ground node g ∈ G, denoted as γ g , can be written as [START_REF] Yang | Capacity Performance of Multicell Largescale Antenna Systems[END_REF]

γ g = M -G a(g) τ ρβ 2 gg p g /(1 + τ ρξ gg ) 1 + g ∈G μ g g p g + τ ρ(M -G a(g) ) g ∈Ig \g β 2 g g p g
1+τ ρξgg [START_REF] Chen | Joint Power Optimization for Multi-source Multi-destination Relay Networks[END_REF] where M represents the number of antennas available at each aerial drone, β g g is the slow fading channel gain between ground node g ∈ G to the service aerial drone of ground node g, i.e., aerial drone a(g), ρ is transmit signal-to-noise ratio (SNR), ξ g g l∈I g

β lg , and μ g g is defined as

μ g g = ⎧ ⎨ ⎩ β g g , if a(g ) / ∈ A(I g ) β g g 1 - τ ρβ g g 1+τ ρξ g g , otherwise, (6)
with A(I g ) {a(g )|g ∈ I g } representing the set of service aerial drones of interfering ground nodes in I g . Let x = (x a ) a∈A , y = (y a ) a∈A , x = ( x g ) g∈G and y = ( y g ) g∈G respectively represent the coordinate vectors of aerial drones and ground nodes, with x a , y a being x-and y-axis coordinates of aerial drone a, and x g and y g being x-and y-axis coordinates of ground node g, respectively. Then, the slow fading channel gain β g g in ( 6) can be expressed as

β g g = H g g (x, y, x, y)ζ g g , (7)
where ζ g g represents the log-normal slow fading between ground node g and the service aerial drone of ground node g, i.e., aerial drone a(g), H g g d -χ g g (x, y) is locationdependent path loss with χ being path loss factor and d -χ g g being the distance (in meter) between ground node g and aerial drone a(g). Then, the capacity achievable by ground node g ∈ G, denoted as R g (α, p, x, y, x, y), can be expressed as R g (α, p, x, y, x, y

) = B log 2 (1 + γ g ), (8)
where B is bandwidth (in Hz) and γ g = γ g (α, p, x, y, x, y) is the achievable SINR defined in [START_REF] Chen | Joint Power Optimization for Multi-source Multi-destination Relay Networks[END_REF].

Problem Statement. The network control objective can then be stated as maximizing the sum capacity of all ground nodes in G and hence the network spectral efficiency by jointly determining transmit power vector p for ground nodes, coordinate vectors x and y for aerial drones, and grounddrone association vector α. The problem is a mixed integer nonlinear nonconvex programming problem (MINLP) because of the nonlinear nonconcave mathematical expression of γ g with respect to α, p, x, y in ( 5)-( 7) and binary association variables α. Given an arbitrary such problem, how to design distributed solution algorithms to achieve the global optimum is still an open problem. In the following sections, we first present a distributed solution algorithm in Section IV and then design a centralized but globally optimal solution algorithm in Section V to provide a performance benchmark for the distributed solution algorithm.

IV. DISTRIBUTED SOLUTION ALGORITHM

A key step of the distributed solution algorithm design is to decompose the original network control problem into a series of subproblems, by solving which in a distributed manner the original problem can be solved. In our case, network control variables α, p, x, y are tightly coupled with each other in the SINR γ g expression in (5) and, therefore, in the achievable capacity in [START_REF] Wu | Secure Massive MIMO Transmission With an Active Eavesdropper[END_REF]. This results in the fact that the resulting network control problem is architecturally undecomposable. 2 Therefore, we design our distributed solution algorithm by decomposing the network control problem following a primal decomposition approach. The core idea is to divide the feasible set of the original problem into multiple parts by fixing a subset of variables at each time; and solve the resulting subproblems iteratively [START_REF] Palomar | A Tutorial on Decomposition Methods for Network Utility Maximization[END_REF]. In our case, primal decomposition can be accomplished by fixing the aerial drone location variables x, y, association variables α, and transmit power vector p, resulting in three subproblems: ground-drone association, ground node power control and drone movement. In the following, we design an iterative distributed solution algorithm, where at each iteration the three subproblems are solved sequentially. Massive MIMO Association. Let x (ν) and y (ν) represent the coordinates of aerial drones and p (ν) be the transmit power of ground nodes in current iteration ν. Then in iteration ν + 1 the objective of the ground-drone association subproblem is to maximize the aggregate capacity of all ground nodes in A by determining the association vector α (ν+1) subject to association constraint [START_REF] Ali | Resource Allocation, Interference Management, and Mode Selection in Device-to-Device Communication: A Survey[END_REF], which is a nonconvex combinatorial problem. Next, inspired by the behaviors of interacting agents in real-time auctions [START_REF] Mcafee | Auctions and Bidding[END_REF], [START_REF] Gerkey | Sold!: Auction Methods for Multirobot Coordination[END_REF], we present an effective Distributed Auction-based Association Strategy (DAAS) by viewing the drones as sellers while the ground nodes as buyers.

The core idea of DAAS is to let ground nodes iteratively compete for association by reporting their own association preferences to nearby aerial drones, which then make association decisions based on the received preference information. Denote N (ν+1) g ⊂ A as the set of nearby aerial drones of ground node g ∈ G (i.e., aerial drones in the communication range of the ground node). Then, the association preference of ground node g with respect to aerial drone a ∈ N (ν+1) g , denoted as λ ga , can be computed as

λ ga = log(1 + γ ga ) a ∈N (ν+1) g log(1 + γ ga ) , ( 9
)
where γ ga represents the interference-free single-input-singleoutput (SISO) SINR, i.e., the SINR achievable with single antenna and without interference from other ground nodes [START_REF] Ngo | Energy and Spectral Efficiency of Very Large Multiuser MIMO Systems[END_REF], achievable by ground node g if associated to aerial drone a . This can be defined as γ ga ρβ gg | a(g)=a with ρ being the transmit SNR and β gg | a(g)=a being the fading coefficient. Denote λ a = (λ ga ) g∈ Na as the preference vector aerial drone a receives from its nearby ground nodes in N a {g|g ∈ G, a ∈ N g }. Then, aerial drone a sorts λ a in descending order and sends association offers to a preferred set of maximum G max of ground nodes, as follows,

λ dsc a = (λ g1a , λ g2a , • • • , λ g Gmax a P ref erred Ground N odes , • • • ). (10)
Then, the association offer corresponding to the highest SISO SINR is accepted by a ground node if it receives multiple offers. The above procedures are repeated until no ground node receives multiple association offers. In [START_REF] Larsson | Massive MIMO For Next Generation Wireless Systems[END_REF], the rationale of computing association preference based on interference-free SISO SINR is as follows. In massive MIMO settings, particularly when the number of antennas M is large, the received SINR is dominated by path loss and large scale shadowing and tends to the interference-free SISO SINR as M → ∞. Therefore, the capacity with interferencefree SISO SINR can serve as a good indication of association preference and can be computed with low computational complexity.

Massive MIMO Power Control. Given the coordinates of aerial drones x (ν) , y (ν) and association vectors α (ν+1) , the ground node power control subproblem can be expressed as Given :

x (ν) , y (ν) , α (ν+1) M aximize

p : g∈G R g (p)
Subject to : 0 ≤ p g ≤ p 0 , ∀g ∈ G [START_REF] Sharma | A Cooperative Network Framework for Multi-UAV Guided Ground Ad Hoc Networks[END_REF] where R g (p) is defined in ( 5)-( 8) and p 0 is the maximum transmit power of ground nodes. As stated in the following theorem, this problem can be approximated as a convex programming problem in the setting of massive MIMO networking, and hence the optimal power control can be achieved with poly-nominal computational complexity. [START_REF] Sharma | A Cooperative Network Framework for Multi-UAV Guided Ground Ad Hoc Networks[END_REF] can be approximated as a convex programming problem.

Proposition 1. In a massive MIMO setting, for given coordinates of aerial drones and association strategies, ground node transmit power control problem in

Proof. According to convexity preserving properties, the logarithm of a concave function is a concave function [START_REF] Boyd | Convex Optimization[END_REF]. Therefore, to prove the theorem we only need to show that the achievable capacity R g in ( 8) is a concave function. In a massive MIMO setting, it is reasonable to assume that the achievable SINR γ g in ( 5) satisfies γ g 1 since the negative effects of interference and small scale fading can be largely smoothed out [START_REF] Ngo | Energy and Spectral Efficiency of Very Large Multiuser MIMO Systems[END_REF]. Then, it can be shown that the expression of capacity R g in ( 8) has the form of minus-log-sum-exp with respect to logarithm transformation of transmit power variables p g , which is a concave function [START_REF] Boyd | Convex Optimization[END_REF].

Drone Movement. In this subproblem, each aerial drone determines its own best coordinates to adapt to the changes in association strategies and transmit power of ground nodes caused by solving the previous two subproblems. With newly obtained association vector α (ν+1) and transmit power vector p (ν+1) , the subproblem of aerial drone movement can be written, for each aerial drone a ∈ A, as

Given : α (ν+1) , p (ν+1) M aximize x,y : g∈Ga R g (x, y), Subject to : x min ≤ x a ≤ x max , y min ≤ y a ≤ y max , (12)
where x min , y min , x max and y max represent the minimum and maximum x-and y-coordinates, respectively, and G a represents the set of ground nodes associated with aerial drone a with given association vector α (ν+1) . In this subproblem, the mathematical expression of utility function R g (x, y) in ( 12) has a log-convex form, which is in general a nonconcave function with respect to coordinate variables x and y. In this paper we solve subproblem (12) using an interior point method [START_REF] Boyd | Convex Optimization[END_REF] to search for locally optimal coordinates for each aerial drone in favor of a low-complexity distributed solution.

In Section V, we will design a centralized but globally optimal solution algorithm at the cost of extra computational complexity.

Computational Complexity Analysis. In the distributed solution algorithm, the above three subproblems are solved iteratively and sequentially at each iteration. In the ground-drone association subproblem, the association strategy is determined iteratively. In each iteration, the association variable α ga can be determined for at least one ground node, and therefore the maximum number of associations is |G|, and the overall computational complexity of the association is O(|G|). Both the power control and the aerial drone movement subproblems can be solved in polynomial computational complexity, i.e., O(|G||A|). Therefore, the complexity of the overall distributed solution algorithm is O(|G|(|A| + 1)) for each iteration. The convergence of the proposed distributed solution algorithm will be studied in Section VI based on simulations, where results indicate that the algorithm can converge very fast, i.e., within tens of iterations. So far, we have presented a distributed solution algorithm to jointly control the movement of aerial drones, ground-drone association and power control of ground nodes in infrastructure-less wireless ad hoc networks with manyantenna-enabled aerial drones. A natural question is: How does the distributed solution algorithm compare to the global optimum in terms of aggregate spectral efficiency? In the remaining of the paper we answer this question by designing a centralized but globally optimal solution algorithm to provide a performance benchmark for the distributed solution algorithm.

V. GLOBALLY OPTIMAL SOLUTION ALGORITHM

As stated in Section III, the "social" objective of the network control problem is to maximize the network spectral efficiency by jointly controlling movement of aerial drones, the transmit power of ground nodes and the ground-drone association strategies. The social OrgSwarm control problem can be mathematically represented as α,p,x,y : ,p,x,y)

Given : A, G, G max M aximize
U g∈G R g (α
Subject to : 0 ≤ p g ≤ p 0 , ∀g ∈ G, x min ≤ x a ≤ x max , ∀a ∈ A, y min ≤ y a ≤ y max , ∀a ∈ A, ( 1
), ( 2), ( 3), ( 4).

In ( 13), the individual throughput R g is a nonconvex/nonconcave function with respect to coordinates variables x, y and transmit power vector p. Moreover, the association variables α take only binary values. Therefore the resulting network control problem is a mixed integer nonlinear nonconvex programming problem (MINLP), for which there is in general no existing solution algorithm that can be used to obtain the global optimum in polynomial computational complexity. In this paper, we design a globally optimal solution algorithm based on a combination of the branch and bound framework and of convex relaxation techniques [START_REF] Lawler | Branch-And-Bound Methods: A Survey[END_REF], [START_REF] Sherali | A Reformulation-Linearization Technique for Solving Discrete and Continuous Nonconvex Problems[END_REF].

A. Overall Algorithm

Denote Γ 0 = {α, p, x, y| constraints in (13)} as the feasible set of initial problem [START_REF] Kandeepan | Aerial-Terrestrial Communications: Terrestrial Cooperation and Energy-Efficient Transmissions to Aerial Base Stations[END_REF] and let U * (Γ 0 ) represent the global optimum of problem (13) over Γ 0 , then the objective of our algorithm is to iteratively search for a U so that

U (Γ 0 ) ≥ εU * (Γ 0 ), (14)
where ε ∈ (0, 1] is predefined optimality precision. To this end, the algorithm maintains a set Γ = {Γ i , i = 0, 1, 2, • • • } of subproblems by iteratively partitioning feasible set Γ 0 into a series of smaller subsets (see Section V-C). 3 The algorithm also maintains a global upper bound U glb (Γ 0 ) and a global lower bound U glb (Γ 0 ) on U * (Γ 0 ) so that

U glb (Γ 0 ) ≤ U * (Γ 0 ) ≤ U glb (Γ 0 ) (15)
to drive the iterations of subproblem partitions, as follows.

• Global upper bound U glb (Γ 0 ): For each subproblem Γ i ∈ Γ, the algorithm computes a local upper bound 3 In this paper we use Γ i to refer to both subproblem i and the corresponding feasible set.

U lcl (Γ i ) on network utility function U via convex relaxation (see Section V-B). Then the global upper bound U glb (Γ 0 ) can be updated as

U glb (Γ 0 ) = max Γi∈Γ {U lcl (Γ i )}. (16)
• Global lower bound U glb (Γ 0 ): Similarly, for each subproblem Γ i ∈ Γ a local lower bound U lcl (Γ i ) is computed based on the solution obtained by solving the relaxed convex network control problem. Then the global lower bound U glb (Γ 0 ) can be updated as

U glb (Γ 0 ) = max Γi∈Γ {U lcl (Γ i )}. (17)
The algorithm terminates if

U glb (Γ 0 ) ≥ εU glb (Γ 0 ) is reached and the global optimum U * (Γ 0 ) is set to U * (Γ 0 ) = U glb (Γ 0 )
as a upper-bound benchmark. Otherwise, the algorithm selects a subproblem from Γ and further partitions its feasible set into two smaller subsets, computes local upper and lower bounds and updates the global bounds U glb (Γ 0 ) and U glb (Γ 0 ) as in ( 16) and ( 17), respectively. In our algorithm, we select the subproblem Γ i ∈ Γ with the highest local upper bound to partition, i.e.,

Γ i = arg max Γi U lcl (Γ i ). (18)
Based on the global bounds update criterion in ( 16) and ( 17), the gap between the two global bounds converges to 0 as the partition progresses. Furthermore, from [START_REF] Ono | A Wireless Relay Network Based on Unmanned Aircraft System With Rate Optimization[END_REF], U glb (Γ 0 ) and U glb (Γ 0 ) converge to the global optimum U * (Γ 0 ).

B. Convex Relaxation

For each subproblem Γ i ⊂ Γ, which is MINLP in our case, a key step is to obtain a relaxed but convex version of Γ i so that it is easy to compute a tight local upper bound U lcl (Γ i ). In this paper the convex relaxation is designed following a two-phase approach as follows. Phase 1: In this phase the relaxation is accomplished by assuming i) there is no mutual interference among ground nodes, i.e., interference items in the denominator of ( 5) are set to zero, and that all ground nodes use different pilot sequences in channel estimation and hence ξ gg = 0 in (5); ii) the maximum number of ground nodes that can be associated to an aerial drone is not limited to G max in (4). Then, the objective of the relaxed network control problem is to maximize the aggregate capacity of ground nodes by determining the optimal coordinate x and y of drones, i.e.,

M aximize

x, y

: U g∈G R g (x, y) Subject to : x min ≤ x a ≤ x max , ∀a ∈ A, y min ≤ y a ≤ y max , ∀a ∈ A, (19)
where R g (x, y) = B log 2 (1 + γ g (x, y)) with γ g (x, y) defined in [START_REF] Chen | Joint Power Optimization for Multi-source Multi-destination Relay Networks[END_REF]. Since γ g (x, y)

1, R g (x, y) can be approximated as

R g (x, y) ≈ B log 2 (γ g (x, y)) (20) ≤ B log 2 M τ ρp 0 ζ 2 gg H 2 gg (x, y) (21) = B log 2 M τ ρp 0 ζ 2 gg d χ gg (x, y) (22)
= B log 2 (M τ ρp 0 ζ 2 gg ) -χB log 2 (d gg (x, y)), [START_REF] Chandhar | Massive MIMO as Enabler for Communications with Drone Swarms[END_REF] where the inequality in [START_REF] Morgenthaler | UAVNet: A Mobile Wireless Mesh Network Using Unmanned Aerial Vehicles[END_REF] holds since G a(g) ≥ 0 in (5), χ is path loss factor and d gg (x, y) is distance (in meter) from ground node g to its service aerial drone a(g). Since d gg (x, y) in ( 23) is a convex Euclidean norm with respect to x and y [START_REF] Boyd | Convex Optimization[END_REF], log 2 (d gg (x, y) cannot be theoretically guaranteed to be concave. In this phase, we obtain a convex relaxation of [START_REF] Chandhar | Massive MIMO as Enabler for Communications with Drone Swarms[END_REF] based on linear approximation of logarithmic function. To this end, we first replace d gg (x, y) in ( 23) with t, then log 2 (d gg (x, y) in ( 23) can be represented as log 2 (t) subject to t ≥ d gg (x, y). Then, log 2 (t) can be further relaxed using a set of linear functions using a segment and three tangent lines.

Phase 2: Phase 2 of relaxation is invoked if the algorithm is done with partitioning coordinate variables x and y, i.e., for each aerial drone a ∈ A,

x max,ax min,a ≤ Δx, (24) y max,ay min,a ≤ Δy, [START_REF] Samira Hayat | Survey on Unmanned Aerial Vehicle Networks for Civil Applications: A Communications Viewpoint[END_REF] where x max,a and x min,a (y max,a and y min,a ) are upper and lower bounds of x-axis coordinate x a (y-axis coordinate y a ), and Δx and Δy are predefined movement step size of aerial drones in x-and y-axis, respectively. The objective in this phase is to determine the optimal association vector α with given aerial drones coordinates vectors x * and y * and without considering mutual interference among ground nodes as in Phase 1 relaxation. Let R ga represent the capacity achievable by ground node g ∈ G if g is associated to aerial drone a ∈ A, then the optimal association can be obtained by solving the following linear optimization problem:

Given : x * , y * M aximize α : a∈A g∈G α ga R ga (α, x * , y * ) Subject to : 0 ≤ α ga ≤ 1, ∀a ∈ A, g ∈ G, g∈G α ga ≤ G max , ∀a ∈ A, a∈A α ga ≤ 1, ∀g ∈ G. (26)
As variable partition progresses, the association variable α ga becomes fixed either to 0 or 1 in all subproblems, for which the optimal transmit power p can be obtained by solving a convex optimization problem as in Section IV.

C. Variable Partition

Variable partition can be conducted by partitioning association variable α and movement variables x and y. For example, given a subproblem Γ i ∈ Γ, by fixing association variable α ga subproblem Γ i can be partitioned into two subproblems with feasible sets Γ i,1 = {(α, p, x, y) ∈ Γ i |α ga = 0} and Γ i,2 = {(α, p, x, y) ∈ Γ i |α ga = 1}, respectively. For movement variables, say x a ∈ [x min,a x max,a ] for aerial drone a ∈ A, the partition can be conducted by splitting x a from the half, resulting in two subproblems with feasible sets

Γ i,1 = {(α, p, x, y) ∈ Γ i |x a ∈ [x min,a x mid,a ]}, (27)
Γ i,2 = {(α, p, x, y) ∈ Γ i |x a ∈ [x mid,a x max,a ]}, (28)
where x mid,a xmin,a+xmax,a 2

.

VI. PERFORMANCE EVALUATION

We evaluate the performance of the proposed solution algorithms by considering a network area of 1000 × 1000 m 2 , with N g = {4, • • • , 15} ground nodes served by N a = {2, 3} aerial drones. The altitude of the drones are set to 100 meters. The number of antennas of each aerial drone is set to N ant = {20, 40, • • • , 160}, the maximum transmit power of each ground node p 0 is set to {100, 200, • • • , 600} mW, path loss factor is set to χ = 2, and the average noise power is set to 10 -8 mW. Results are obtained by averaging over 20 independent simulation instances with network topology randomly generated.

Performance of the proposed distributed solution algorithm is evaluated by comparing it with the global optimum obtained by the centralized solution algorithm designed in Section V. Figure 1 reports the network spectral efficiency achievable by the distributed solution algorithm and the centralized global optimum. It can be seen that the distributed solution algorithm achieves a performance very close to the global optimum in all the tested 20 instances, where 12 ground nodes and 2 aerial drones each with M = 100 antennas are considered. The average performance of the achievable network spectral efficiency is reported in Fig. 2 with different number of ground nodes. Results indicate that around 95% percent of the global optimum can be achieved by distributed solution algorithm, e.g., 94.9% and 95.4% with 6 and 10 ground nodes, respectively.

The network spectral efficiency is reported in Fig. 3 with the number of antennas at each aerial drone varying from 20 to 160 at step of 20. It can be seen that over 95% of the global optimum can be achieved by the distributed solution algorithm. Results also indicate that the achievable network spectral efficiency monotonically increases with the number of antennas but at a decreasing speed. With more than 100 antennas, the effect of increasing the number antennas is marginal.

In Fig. 4 we plotted the achievable network spectral efficiency against the maximum transmit power of ground nodes. Over 95% of the global optimum can be achieved by the distributed solution algorithms. We noticed that the network spectral efficiency rises only around 10% by increasing the maximum transmit power by 6 times from 100 to 600 mW. This is because in massive MIMO setting the network basically operates at high SINR regime, i.e., in bandwidth-limited regime. We further studied the effects of different network control strategies on the achievable network spectral efficiency in Fig. 5, where the joint network control strategy is compared to the other three strategies, in which aerial drone movement is randomly generated in "w/o Aerial Drone movement", power control and association are executed only once in "w/o" power control and "w/o" association, respectively. It can be seen that the joint network control achieves the highest spectral efficiency in almost all tested instances. Results also indicate that a spectral efficiency close to that of joint control can be achieved by "w/o power control" in most network instances since the network operates in high SINR regime in massive MIMO setting.

VII. CONCLUSIONS

We studied infrastructure-less wireless ad hoc networking with many-antenna-enabled aerial drone hotspots. The network control objective is to maximize network-wide spectral efficiency by jointly controlling movement of drones, associating single-antenna ground nodes to many-antenna drones, and adapting transmit power of ground nodes. The network control problem was formulated as an integer nonlinear nonconvex programming problem (MINLP). Both distributed and globally optimal solution algorithms have been designed and evaluated with extensive simulation results. Results indicated that the distributed solution algorithm can achieve around 95% of the global optimum.

Fig. 1 :

 1 Fig. 1: Aggregate network spectral efficiency achievable by distributed solution algorithm and global optimum.

Fig. 2 :

 2 Fig. 2: Network spectral efficiency in the case of different number of ground nodes.

Fig. 3 :

 3 Fig. 3: Network spectral efficiency in the case of different number of aerial drone antennas.

Fig. 5 :

 5 Fig. 5: Achievable network spectral efficiency with different network control strategies.

The maximum number of served ground nodes cannot exceed the number of antennas available to each aerial drone and the length of pilot sequences used in channel estimation[START_REF] Ngo | Energy and Spectral Efficiency of Very Large Multiuser MIMO Systems[END_REF].

A problem is architecturally decomposable if its dual problem obtained by introducing Lagrange multipliers can be rewritten into a set of subproblems, each of which can be solved locally in a single protocol layer and network device[START_REF] Chiang | Layering as Optimization Decomposition: A Mathematical Theory of Network Architectures[END_REF].


Foroogh Mohammadnia

Marco Fiore

Marco Ajmone Marsan

Adaptive densification of mobile networks: exploring correlations in vehicular and telecom traffic

Keywords: Mobile network traffic, Vehicular traffic, Dense radio access network, Moving base station, Small cells

In this paper we use real data to investigate the correlation in time and space between vehicular traffic and telecom traffic in cellular networks. The analysis of such a correlation is critical to assess the feasibility of cellular network architectures where densification is achieved with small-cell mobile base stations carried by vehicles. This is an attractive possibility for the provisioning of on-demand capacity through temporary dense small cell deployments where and when needed.

Our results indicate that vehicular traffic at penetration rates expected for small-cell-carrying vehicles is much more bursty than telecom traffic. Yet, some correlations between the two emerge, even when considering an entire large metropolitan area. More importantly, correlations tend to be stronger in densely urbanized areas and during high-demand time periods, i.e., where and when radio access densification is most needed. Overall, our results indicate small-cell base stations carried by vehicles as a promising and cost-effective approach to dynamic provisioning of capacity in future-generation cellular networks.

I. INTRODUCTION

The use of large numbers of Small Cells (SCs) for the densification of Radio Access Networks (RANs) is regarded as a most promising approach for the provision of broadband services to large numbers of mobile network end-user terminals. However, the deployment of the SC base stations (BSs) necessary to make dense RANs a reality implies huge investments, which Mobile Network Operators (MNOs) are reluctant to approve in a period of declining profits.

A major reason for the high SC deployment cost is that peaks of bandwidth demand vary significantly in time and space [START_REF] Furno | Mobile Demand Profiling for Cellular Cognitive Networking[END_REF], [START_REF] Furno | A comparative evaluation of urban fabric detection techniques based on mobile traffic data[END_REF], thus making the SC utilization highly variable, and low on average. For instance, users tend to commute from home to work in the morning of working days, making business districts crowded during working hours. In this period, MNOs need high-capacity RANs, hence a dense presence of SCs, to be located in business areas. However, after work, users move out of their offices, so that the demand for RAN capacity becomes much lower in business districts, and many of the installed SCs become redundant. The opposite is true for residential areas, where capacity is mostly necessary in the evening, so that a dense RAN layout becomes useful then, not during working hours. In addition, during commuting and in the event of traffic jams or accidents, dense RAN layouts are needed where road traffic is congested. The cost of a RAN which is dense everywhere can be extremely high, and the associated long periods of low utilization imply a poor return on investment for the MNO.

A possible approach to make high RAN capacity available where and when needed is to use mobile SC BSs carried by vehicles. This approach leads to a network where not only end users, but also (part of) the network infrastructure is mobile. We name such an architecture a mobile mobile network, or MOMONET. The idea of mobile BSs is not new. Operators have exploited truck-mounted BSs for the quick provision of coverage in areas where service was not otherwise available, or where additional capacity was needed [START_REF] Dael | [END_REF]. Some truckmounted BSs were periodically moved from beach areas in summer to skiing resorts in winter. However, in this paper we discuss a new, much broader use for mobile BSs, looking at the possible impact of vehicle-mounted SC BSs on the implementation of dense RANs based on SCs.

The feasibility of RAN densification with vehicle-mounted SC BSs raises a number of research questions that will require accurate investigation. For example: the selection of algorithms for the association of end users to BSs and for handover mitigation; the management of interference among transmissions from traditional BSs, from small cell BSs, and from end users; the selection of technologies for the implementation of small cell BSs on vehicles, and for the implementation of backhaul wireless connections; the reduction of the impact of the energy consumption of small cell BSs on vehicles power systems; the impact of self-driving cars.

Vehicles have the desirable property of moving together with end users 1 , so that in a business district during working hours we normally have both many end users and many vehicles. If a significant fraction of those vehicles carry a SC BS, a temporary dense SC deployment is created. Quite nicely, this temporary dense SC deployment will be recreated in residential districts when drivers return home with their cars, as well as in traffic jams.

Note that the correlation between the number of vehicles and the number of MNO end user terminals, hence between telecommunication traffic and vehicular traffic, was assumed in several works (see, e.g., [START_REF] Janecek | The Cellular Network as a Sensor: From Mobile Phone Data to Real-time Road Traffic Monitoring[END_REF] and [START_REF] Caceres | Traffic flow estimation models using cellular phone data[END_REF]), starting with the Real Time Rome project of the MIT Senseable City Lab [7], [8], but has never been carefully measured in real-world settings. In this paper, we aim at filling this gap. We use recent measurement data on two months of telecom traffic and vehicular traffic in the city of Milan, Italy, and evaluate their spatio-temporal correlation.

The rest of the paper is organized as follows. In Section II we briefly overview previous works on mobile BSs. In Section III we present our reference dataset. In Section IV we provide a preliminary exploration of the data. In Section V we illustrate and comment correlation results. In Section VI we investigate relationships further, by assessing the average distance between end user terminals and SC BSs carried by vehicles. Finally, Section VII concludes the paper.

II. RELATED WORK

The idea of using cars as active elements of a telecommunication network is not new. Several studies in the domain of vehicular networks considered the possibility of integrating a Wi-Fi Acces Point (AP) in cars, and possibly also a cellular interface, so that a vehicle can provide connectivity to the sourrounding cars. For example, [START_REF] Sadek | A Robust Multi-RAT VANET/LTE for Mixed Control & Entertainment Traffic[END_REF] proposed the use of mobile vehicular gateways that exploit Wi-Fi for vehicleto-vehicle (V2V) communications and LTE for vehicle-toinfrastructure (V2I) communications. The study in [START_REF] Frangiadakis | Virtual Access Points for Vehicular Networks[END_REF] introduced the concept of Virtual APs, which allow extending the reach of roadside access points: vehicles that receive a message, store it, and rebroadcast it into non-covered areas. In [START_REF] Dressler | Towards a Vehicular Cloud Using Parked Vehicles as a Temporary Network and Storage Infrastructure[END_REF]- [START_REF] Malandrino | The Role of Parked Cars in Content Downloading for Vehicular Networks[END_REF], parked vehicles were exploited, in addition to roadside units, to improve the performance of video downloads and other services toward travelling cars. In [START_REF] Borsetti | Virtual data mules for data collection in road-side sensor networks[END_REF], mobile nodes were used for the collection of data from sensors.

These approaches are different from a MOMONET, which leverages vehicles as a support for SC BSs of a dense RAN. MOMONETs allow the seamless integration of the vehiclemounted SC BSs with the "normal" macro-and micro-BSs at fixed locations; specifically, the moving SC BSs provide additional capacity, when and where needed, to the end users of an otherwise traditional RAN. End-user terminals can thus freely and transparently transfer their services between a normal and a vehicular BS. MOMONETs generalize the concept of early deployments of SC BSs onboard public transport vehicles to serve passengers in ultra-dense urban environments [START_REF] Sui | On the deployment of moving networks in ultra-dense urban scenarios[END_REF]- [START_REF] Yasuda | A Study on Moving Cell in 5G Cellular System[END_REF], and that of using vehicle-carried SC BSs in public safety networks [START_REF] Shin | Moving small cells in public safety networks[END_REF]. A concept similar to MOMONET has also been briefly hinted at in recent works [START_REF] Andreev | Future of Ultra-Dense Networks Beyond 5G: Harnessing Heterogeneous Moving Cells[END_REF], [START_REF] Pous-Fenollar | Mobile Crowdcell[END_REF], whose authors mention the possible use of 4G or 5G small cells within moving and parked cars, or on drones, so as to serve users both indoor and outdoor. can be seen on the map. The time granularity of the dataset is 15 minutes (96 intervals per day). For each time interval, the dataset reports the total amount of mobile network data traffic (voice and data, no signaling) recorded in every square by Telecom Italia Mobile (TIM), the incumbent Italian MNO, as well as the GPS trajectories of active (i.e., not parked) vehicles tracked by InfoBlu (about 4 million, on a national scale), a large Italian infomobility services provider. The total number of time intervals is 61 × 96 = 5, 856. Considering the number of area squares for which data are reported separately, we have a total of 3,373,056 data points.

It is important to remark that the penetration rates of the two measurement technologies are different. While the telecom traffic in the dataset corresponds to more than one third of the overall mobile network demand, the fraction of the monitored vehicles is less than 1%, as it mainly encompasses probe vehicles and commercial fleets monitored by InfoBlu. We stress that these settings are consistent with the goal of our investigation: we expect a small percentage of vehicles to act as SC BS carriers, and would like them to cover the largest fraction possible of the total telecom traffic.

IV. DATA EXPLORATION IN REPRESENTATIVE CASES

We start by showing the average telecom and vehicular daily traffic patterns in a representative selection of squares. We consider squares of different types, based on several major land uses in the area: a) residential, b) business, c) train station, d) airport, and e) soccer stadium. Figure 2 shows the average telecom and vehicular traffic patterns for all days of the week for the five selected squares (telecom traffic in the left plots and vehicular traffic in the right plots). We observe a number of interesting features in both telecom and vehicular daily traffic patterns. First of all, while in telecom traffic clear behaviors exist, the vehicular traffic is highly irregular. In a residential area, plots (a) and (b) in Figure 2, we can clearly see that the telecom traffic is very low at night (roughly from 2 am to 7 am), it increases (sharply early in the morning, then more mildly) until 10 pm, and finally starts declining again. No big differences exist among days, except for the fact that the traffic on Sundays is higher during the night and around midday.

(c) (d) (e) (f) (g) (h) (i) (j)
Fig. 2: Average telecom and vehicular traffic (left and right, respectively) for all days of the week in selected square areas. Plots (a) and (b) refer to a residential area, plots (c) and (d) to a business area, plots (e) and (f) to a train station, plots (g) and (h) to an airport, and plots (i) and (j) to a soccer stadium.

Vehicular traffic exhibits instead many ups and downs, with the highest peaks between 7 am and 10 am; also, vehicular traffic on Sundays is lower than in other days. In a business area, plots (c) and (d) in Figure 2, telecom traffic is very low at night (roughly from 2 am to 7 am), it increases sharply around 8 am in working days, and then gradually decreases after 2 pm. The peak of telecom traffic in weekends is about half the one in working days. Vehicular traffic again exhibits high variability, with peaks at around 8 am and between 3 pm and 5 pm, and is lower on Sundays.

In a train station area, plots (e) and (f), telecom traffic is once more very low at night, then sharply increases around 8 am in working days, and remains close to its peak value from 9 am to 9 pm, then sharply decreases. The peak of telecom traffic in weekends is about one third of that in working days. Vehicular traffic is very bursty, as usual, peaking in the afternoon and evening.

In an airport area, plots (g) and (h), telecom traffic is very low at night, sharply increases around 5 am, and remains close to its peak value from 8 am to 10 pm, then sharply decreases. The peak of telecom traffic in weekends is about two-thirds of the one in working days. Vehicular traffic as usual exhibits many ups and down along most of the day.

Finally, in a soccer stadium area, plots (i) and (j), telecom traffic exhibits very high peaks around 9 pm in Saturdays and Sundays, when soccer matches are played (and concerts take place). Otherwise, traffic is close to zero. Vehicular traffic is once more irregular, peaking at around 7 am and early in the afternoon.

The inspection of the raw data about telecom and vehicular traffic that we presented so far does not provide many clues to justify correlation; rather, it reveals that telecom traffic dynamics vary widely in space and time, and our reference vehicular traffic is generally much more bursty than telecom traffic. This latter effect is mainly due to the reduced fraction of vehicles carrying SC BSs that will characterize MOMONET deployments.

We provide a more global overview of the evolution of vehicular and telecom traffic in Figure 3. There, the variation of traffic in space is visualized through heat maps, where for each square the color indicates the (telecom or vehicular) traffic level (dark blue is lowest and dark red is highest -the color scale is logarithmic). The different plots portray heat maps for a same day (March 2, 2015) at five different times: 5 am, 10 am, 3 pm, 8 pm and midnight. We can see from these plots (and for similar ones at all time intervals in the day) that telecom traffic, in comparison to vehicular traffic, rises earlier in the morning, then it remains high later in the evening, and is more pervasive over the entire Milan area. Also the simple visual inspection of the geographical distribution of vehicular and telecom traffic above does not seem to highlight evident correlations, except for the fact that traffic peaks coincide with the city center in both cases.

V. CORRELATION ANALYSIS

As the qualitative analysis in Section IV did not reveal evident relationships, we carry out a comprehensive quantitative study of the correlations between telecom and vehicular traffic. Specifically, we separately explore the spatial, temporal, and spatiotemporal dimensions of the problem.

A. Spatial correlation

The results about correlations from a spatial perspective are summarized by the scatterplots in Figure 4. There, the abscissa refers to the telecom traffic value, and the ordinate reports the vehicular traffic value for each dot corresponding to one square in the Milano area (thus 576 dots), and traffic values are averaged over all time intervals of working days (Monday through Friday; 44 × 96 instances -top plot) and weekends (Saturday and Sunday; 17 × 96 instances -bottom plot). Therefore, Figure 4 outlines correlations emerging across the geographical areas in the target metropolitan region. Plots also show the least square linear regression line, whose slope is 0.018 for working days and 0.016 for weekends.

The positive slope values indicate that a positive correlation exists between telecom and vehicular traffic across space, despite the dispersion of points around the regression line.

B. Temporal correlation

Correlations in the time dimension are considered in Figure 5, where each dot of the scatterplots represents one daily time interval (for a total of 96 dots). Results are aggregated over the whole metropolitan area: the coordinates of each dot thus represent the average telecom and vehicular traffic in the time slot, averaged over 576 squares and either 44 working days in plot (a) or 17 weekend days in plot (b). The shape of these plots is quite interesting, if we look at them while considering the sequence of time slots. For this, plots (c) and (d) report the same data as plots (a) and (b), respectively, in a 3D plot, where the z-axis represents time of the day. We see that a significant amount of telecom traffic is present at 12:15 am, while vehicular traffic is close to zero. Telecom traffic then declines, and both types of traffic are at their minimum around 3 am. Vehicular traffic then starts to increase, while telecom traffic remains low (time runs clockwise along the dot sequences in Figure 5, plots (a) and (b)). Then, telecom traffic increases, and remains high until late in the evening, when vehicular traffic has become low. Finally, also telecom traffic dies out. Note that the effects are more pronounced during working days than in weekends.

Linear regression lines highlight a positive correlation also in this case, although points are typically far from the behavior predicted by a simple linear model. The limited goodness of fit can be imputed to the temporal offset between the high activity periods of vehicular and telecom traffic explained above.

C. Spatio-temporal correlation

We delve deeper in the analysis by disaggregating traffic over both space and time. Figure 6 shows plots where each dot refers to one square in one 15-minute interval of the day (hence 576 × 96 = 55, 296 dots in each plot). The coordinates of each dot are the average telecom and vehicular traffic recorded during either 44 working days or 17 weekend days. Also in this case, the least square linear regression lines have positive slopes, equal to 0.014 in the case of working days, plot (a), and 0.013 in the case of weekends, plot (b), which entail positive correlation for both. However, dots are now even more scattered over the quadrant, indicating that the linear regression is not representative of many combinations of squares and daytimes.

D. When and where is correlation higher?

The results above yield the following conclusions: i) there exist a general positive correlation between vehicular and telecom traffic, but ii) the correlation is low, due to a high variability of relations over time and space. A sensible question is then whether some hours and geographical areas show especially strong correlations, and are thus suitable periods and locations for MNOs to take advantage of the MOMONET paradigm. To answer this question, we leverage three different correlation coefficients:

• The Pearson product-moment correlation coefficient r, defined as:

r = n n i=1 x i y 1 -   n i=1 x i n j=1 y j     n n i=1 x 2 i - n i=1 x i 2   n n i=1 y 2 i -( n i=1 y i ) 2 ,
(1) where n is the number of observations for both telecom and vehicular traffic, and x i , y i are the telecom and vehicular traffic observations, respectively.

• The Spearman rank correlation coefficient s, defined as:

s = 1 - 6 n i=1 y 2 i d 2 i n(n 2 -1) (2)
where d i is the difference between the ranks (ordered positions) of corresponding telecom and vehicular traffic samples.

• The Kendall rank correlation coefficient q, defined as:

q = n c -n d n(n -1)/2 (3)
where n c is the number of corresponding telecom and vehicular traffic samples with the same rank and n d is the number of corresponding telecom and vehicular traffic samples with different rank. We first investigate which are the periods during the day when the correlation is higher. Figure 7 (a) shows the curves of the Pearson correlation coefficient for the average working day and for the average weekend day. The curves for the Spearman and Kendall correlation coefficients are not reported, but are very similar. We observe that correlation is low overnight in working days and in the early morning of weekends. On the contrary, the correlation coefficient values are higher -in the range from 0.3 to 0.5 -during high-demand periods, i.e., in working and evening hours. These are good news, since the correlation is low in the night periods, when the demand for mobile services is typically low, and network densification is not necessary. Instead, MOMONETs are the most effective exactly when densification is needed, that is, in presence of mobile network traffic demand peaks.

Concerning the geographical dispersion of correlation, Figure 7 (b) portrays the cumulative distribution functions (CDFs) of the values of the three correlation coefficients for all spatial squares in the considered metropolitan area. We can see that negative correlations are limited to a few cases, and correlation coefficient values are positive in the large majority (90-95%) of areas. Differences across coefficient types are also small, as the Pearson and Spearman correlation coefficients yield very similar CDFs, while the CDF of the Kendall correlation coefficient values indicates some prevalence of lower values of correlation coefficients. However, the important observation is that the three CDFs confirm the existence of spatial heterogeneity of correlation: depending on the considered square, the correlations vary significantly, and can reach values as high as 0.84.

In order to understand where such high-correlation squares lie, Figure 7 (c) shows a heatmap where the color of each square corresponds to the Pearson correlation coefficient computed over an average day. We can clearly see that correlation is stronger in the city center. On the one hand, this is expected, since both vehicular and telecom traffic are higher in those squares, hence they follow more predictable patterns. On the other hand, and more importantly, those are also the squares where the demand for radio access network capacity is higher, hence densification is needed the most. We conclude that, in our reference scenario, MOMONETs prove to be mainly exploitable not only when, but also where network densification is needed the most.

VI. AVERAGE DISTANCE

As a final complement to our analysis, we comment on the average distance between SC BS-equipped vehicles and the location of end user terminals that can use the SC BSs carried by vehicles to transmit and receive mobile network traffic. Since our dataset reports exact vehicle positions within each square, but only the origin square (and not the exact end user terminal location) for telecom connections, we randomly place end user terminals within the area covered by each square. Of course, this is an approximation, but we have no information to drive a more realistic placement of end users. The average distance between the end user terminal and the closest SC BS (the one to which the end user terminal normally associates) is a metric of great interest for MOMONETs, since it is a proxy for the connection data rate.

Results in Figure 8 (a) show the evolution of such average distance over time, averaged over all spatial squares and averaged over all 61 days. The average distance is high (it can grow over 190 m) in periods of low telecom and vehicular traffic, and it stays at around 120 m during periods characterized by high network traffic, between 7 am and 6 pm roughly. These distances risk to make connections to SCs frequently unviable, due to poor signal strengths (consider that the power emitted by a SC is typically less than 1W, and reach depends on propagation conditions, but normally is less than 100 m).

The same observation holds when considering the geographical distribution of distances, averaged over the whole day. In However, when focusing on the city center, the outcome is different. Plot (c) in Figure 8 depicts the average distance between the position of end user terminals and the closest SC BS carried by a vehicle, in the 158 squares in the center of Milan averaged over all 61 days. Between 9 am and 8 pm average distances are around 50 m, so that in this portion of the city we can expect a large fraction of terminals to be within reach of a SC BS carried by a vehicle. These values nicely complement the results in Section V, and confirm that MOMONETs represent a suitable network paradigm to enable the adaptive densification of radio access envisioned for nextgeneration 5G systems, at least in the highly urbanized areas where mobile network capacity will be needed the most.

VII. CONCLUSIONS

We used real traces of vehicular traffic and telecom traffic to investigate the correlation in time and space between the two types of traffic. The motivation for the study lies in the assessment of the feasibility of cellular network architectures where densification is achieved with small-cell mobile base stations carried by vehicles, according to a novel mobile mobile network, or MOMONET, architecture. Our results prove that correlation between vehicular and telecom traffic is low when we consider an entire metropolitan area, but becomes stronger when we look at the areas and the time periods where and when densification is necessary to serve peaks of telecom traffic. Similar considerations hold when looking at the distance between mobile SC BSs and the users they could serve. Ultimately, these observations support small-cell base stations carried by vehicles, and the associated MOMONET paradigm, as a promising and cost-effective approach to the adaptive densification of cellular networks.

Fig. 1 :

 1 Fig. 1: The subdivision of the Milan area into 576 squares.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Heatmaps of Telecom (first row) and vehicular (second row) traffic on March 2, 2015. The first column shows data at 5 am; the second column at 10 am, the third at 3 pm, the fourth at 8 pm, the last at midnight. The plotted area is about 18 × 21 km.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Correlation between average telecom and average vehicular traffic in each timeslot over 44 working days in plot (a) and over 17 weekend days in plot (b); plots (c) and (d) report the same data as (a) and (b), with time increasing along the z-axis.

Fig. 7 :

 7 Fig. 7: Pearson correlation coefficient between average telecom and vehicular traffic over all squares and over 44 working days and 17 weekend days in plot (a); cumulative distribution function of the correlation coefficients between average (over 61 days and over timeslots) telecom and vehicular traffic in each square in plot (b); Pearson correlation coefficient heat map between telecom and vehicular traffic for each Milano square, averaged over 61 days in plot (c).

Fig. 8 :

 8 Fig. 8: Average distance between vehicles and origins of telecom traffic versus time of day averaged over 61 days and over all Milano squares in plot (a); average distance between vehicles and origins of telecom traffic for all square areas over one day (March 2, 2015) in plot (b); and average distance between vehicles and origins of telecom traffic versus time of day for 158 squares in the center of Milan averaged over 61 days in plot (c).

Here we mostly refer to humans, but vehicles can also be host to several connected devices; for instance, a recent study by HIS[START_REF] Ihs | IHS Connected Car consumer survey finds strong appetite for ADAS features and autonomous cars[END_REF] predicts that by

the fraction of connected cars in Europe will reach 50%. ISBN 978-3-903176-05-8 © 2018 IFIP

ACKNOWLEDGMENTS

This work was supported by the Siebel Energy Institute through the Seed Grant Award 59/2016. This research received funding from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme under REA grant agreement n.630211, ReFleX.