Christelle Caillouet

Tahiry Razafindralambo

Dimitrios Zorbas

Optimal placement of drones for fast sensor energy replenishment using wireless power transfer

Lifetime is the main issue of wireless sensors networks. Since the nodes are often placed in inaccessible places, the replacement of their battery is not an easy task. Moreover, the node maintenance is a costly and time consuming operation when the nodes are high in numbers. Energy harvesting technologies have recently been developed to replenish part or all of the required energy that allows a node to function. In this paper, we use dedicated chargers carried by drones that can fly over the network and transmit energy to the nodes using radio-frequency (RF) signals. We formulate and optimally solve the Optimal Drone Placement and Planning Problem (OD3P) by using a given number of flying drones, in order to efficiently recharge wireless sensor nodes. Unlike other works in the literature, we assume that the drones can trade altitude with coverage and recharge power, while each drone can move across different positions in the network to extend coverage. We present a linear program as well as a fast heuristic algorithm to meet the minimum energy demands of the nodes in the shortest possible amount of time. Our simulation results show the effectiveness of our approaches for network scenarios with up to 50 sensors and a 50 × 50m terrain size.

I. INTRODUCTION

A. Context

Wireless sensor networks have brought a tremendous amount of civil, military and industrial applications in the context of the Internet of Things. They are capable of periodically monitoring their vicinity and therefore generating qualitative and quantitative measurements about their environment. Their main advantage is their low size and the ability to take measurements in real-time and forward them for further analysis. However, one of their main drawbacks is their limited battery lifetime.

Optimizing the battery consumption of wireless sensors has been the focus of many works in the literature. Since wireless sensors are autonomous objects powered by batteries, the performance of the underlying application strongly depends on their energy capacity. The data generation and transmission frequency of the sensors are always fundamental for the application. These two services consume the majority of the sensors' energy, thus, limiting the application's lifetime. Optimizing battery consumption from data generation and transmission frequency is important but not at the expense of the application.

Even though applications are already optimized, prolonging the lifetime of a wireless sensor network is still a fundamental issue. Since the deployment of the wireless sensors may happen in harsh environments where the sensors are not accessible, it is not always possible to replace the batteries with a low cost. To tackle or to alleviate this issue energy harvesting solutions have been proposed.

There are multiple energy harvesting solutions to increase the lifetime of sensor networks. On one hand, it is possible to harvest energy using, for example, solar panels or other types of harvesting. However, even if this type of harvesting is very effective and can greatly extend the sensors' lifetime, the amount of power is not constant and most of the time requires a big in size extra equipment. On the other hand, using a technology such as the RF-power harvesting [START_REF] Lu | Wireless networks with rf energy harvesting: A contemporary survey[END_REF], we can transmit a(n) (on-demand) constant amount of energy to the nodes, and recharge multiple nodes at the same time.

In this paper, we use the energy harvested by transmitted RF signals from flying drones to power up energy-constrained sensor nodes. We examine the scenario where the nodes are not always on, but we periodically visit them in order to (partially) recharge their battery and collect the data. We study the problem of optimally deploying flying drones over the sensor area. These flying drones carry dedicated chargers facing to the ground and are capable of emitting RF signals towards the sensors to recharge their battery as it is illustrated in the example of Fig. 1. We develop a Mixed Integer Linear Program (MILP) and an heuristic to optimize drones utilization.

B. Contributions

Unlike previous works in the literature which do not consider realistic harvesting values neither drone mobility, (a) we take into account realistic harvesting values (for both transmitters and receivers) obtained from a RF-power harvesting module manufacturer, (b) we place the drones in 3D-space while taking into account the effect of the altitude adjustment to the harvesting power, and (c) we consider that the drones can change positions in the 3D-space to improve coverage and, thus, minimize the total energy replenishment time. To this end, we formulate and optimally solve the optimal drone placement and planning problem (OD3P) for the nodes energy replenishment using a mixed integer linear programming formulation. A sub-optimal greedy heuristic algorithm is also proposed as a time-efficient solution for large scale networks. We present extensive simulation results for network instances with up to 50 nodes and 50 × 50m terrain sizes and we study the effectiveness and scalability of the approaches.

C. Outline

The rest of this paper is organized as follows. Section II briefly reports the most recent related work. In Section III we describe in detail the OD3P problem. Then, in Section IV, we present the optimization model followed by an heuristic to solve OD3P in Section V. We then evaluate the performance of OD3P and we discuss its scalability in Section VI before concluding the paper in Section VII.

II. RELATED WORK

Recent works have already considered using vehicles as chargers to enhance wireless sensor network lifetime. However, they often limit the problem in a two-dimensional space [START_REF] Zhang | P3: Joint optimization of charger placement and power allocation for wireless power transfer[END_REF] which decreases the accessibility and the recharge opportunities. Besides, some works focus only on a recharge model where the nodes are recharged one after the other [START_REF] Najeeb | Extending wireless rechargeable sensor network life without full knowledge[END_REF]. Given the vehicle's lifetime restrictions, this solution may not be scalable since a considerable amount of time is spent to travel from one node to another. In our work, we can cover multiple nodes by a single drone depending on its position and the density of the network.

Optimal dedicated charger positioning problems are also examined in [START_REF] Zorbas | On optimal charger positioning in clustered rf-power harvesting wireless sensor networks[END_REF] and [START_REF] Pang | Charging coverage for energy replenishment in wireless sensor networks[END_REF]. The first work presents a localized solution as well as centralized approximation algorithm to place mobile chargers among clustered wireless sensor nodes. The objective of the problem is to extend the lifespan of as many nodes as possible. An approximation algorithm for the same problem is presented in [START_REF] Pang | Charging coverage for energy replenishment in wireless sensor networks[END_REF]. In [START_REF] Zorbas | On the optimal number of chargers in battery-less wirelessly powered sensor networks[END_REF], the potential of using battery-free nodes is examined. Dedicated chargers with fixed positions and omni-directional antennae are placed across the network to continuously transmit power to the nodes. A number of solutions for the charger placement problem while taking into account the electromagnetic radiation activity is, also, proposed in the literature [START_REF] Dai | Radiation constrained wireless charger placement[END_REF], [START_REF] Wang | Energyefficient secure transmission for wireless powered internet of things with multiple power beacons[END_REF], [START_REF] Rault | Avoiding radiation of on-demand multi-node energy charging with multiple mobile chargers[END_REF]. However, all the previous works consider only 2-dimensional charger deployments and a non-realistic high-power omni-directional energy emission.

The financial cost of the deployment is an important parameter and [START_REF] Zorbas | Computing optimal drone positions to wirelessly recharge iot devices[END_REF] tries to minimize the number of drones used to maintain an operating wireless sensor network for a given amount of time. The authors provide a integer linear model related to the set cover problem and heuristics, to minimize the number of drones needed to wirelessly recharge the network. However, this work does not consider drone mobility and, thus, its practicality is limited only to very dense networks. In sparse networks, a number of drones that equals the number of nodes would be required.

III. PRELIMINARIES

Given a set of nodes I with known and fixed 2D positions, the goal is to optimally deploy and plan the positions of a set of drones N in the 3D-space in order to provide an amount of energy at least equal to E to all the nodes in the network in the shortest possible amount of time. This amount of energy is an empirical value that has to be sufficient to power-up a node, to allow it to take a measurement, and transmit the data to the drone (or to a base station). We assume that the transmission time is negligible compared to the recharge time of the nodes.

Each node is equipped with an RF-power harvesting module capable of converting the RF-power from transmitted signals to DC power. A charger with a directional antenna facing the ground is used to continuously emit power to the network and recharge the nodes. A drone can adjust its altitude to recharge multiple nodes at the same time.

The amount of power each node receives is affected by its distance to the transmission source (charger on the drone) and the environmental conditions. Formula (1) describes the power harvested by a node i when it is located within the line of sight of the charger [START_REF] Zorbas | Assessing the cost of rf-power harvesting nodes in wireless sensor networks[END_REF]:

P hi = P dij rx f dij , (1)
where P dij rx is the received power and f dij is the efficiency of the harvesting antenna at distance d ij . The received power at distance d is given by the following propagation model [START_REF] Win | A mathematical theory of network interference and its applications[END_REF]:

P d rx = P 0 e 2σG d 2b , (2)
where e 2σG d 2b has a log-normal distribution with a shadowing coefficient σ (G ∼ N (0, 1)) and b is the amplitude loss exponent. P 0 is the received power at reference distance.

The minimum harvested power received by a node to recharge its battery depends on the efficiency of the converter for the corresponding received power at this particular distance. We denote by Γ this minimum harvesting threshold, which is a hardware-depended constant.

Formula (3) defines the harvested energy of a node i for a given time period t:

H t i = t 0 P hi tdt. ( 3
)
This energy is stored directly in a super-capacitor with some leakage properties expressed as ηH t i , where η ∈ (0, 1).

IV. THE OD3P OPTIMIZATION MODEL

We set J the set of possible drone 3D-positions and denote with E i the energy needs for each node i in I. We define as β n j the binary variable stating whether a drone n ∈ N is deployed at position j ∈ J . Let α n be a binary variable stating if drone n ∈ N is used or not. We also add into the model the continuous variables t n j representing the associated time a drone n ∈ N is located at position j ∈ J . We finally set τ the maximum total recharge time of a drone. τ practically represents the drone's maximum flying time.

The goal is to minimize the time needed to harvest and store the required quantity of energy from the deployed drones to the sensors, while not exceeding the budget (in terms of number of drones) B ≤ |N |. We specify that at most one drone can be placed at each 3D-position, determining one possible altitude. The Mixed Integer Linear Program (MILP) is given as follows:

min λ (4) j∈J t n j ≤ λ, ∀n ∈ N (5) n∈N j∈J (1 -η)H t n j i ≥ E i , ∀i ∈ I (6) n∈N α n ≤ B ( 7
) j∈J t n j ≤ τ α n , ∀n ∈ N (8) n∈N β n j ≤ 1, ∀j ∈ J (9)
t n j ≤ τ β n j , ∀n ∈ N , j ∈ J (10) n∈N t n j + j ∈J ,(xj ,yj )=(x j ,y j ) ( n ∈N t n j ) ≤ λ, ∀j ∈ J (11) α n ∈ {0, 1}, β n j ∈ {0, 1}, t n j ∈ R, λ ∈ R (12)
The objective function seeks to minimize the longest drone flying time min max n∈N j∈J t n j , which corresponds to the minimum recharge time that ensures the harvested energy needs per node. It has been decomposed into one objective (4) minimizing a continuous variable λ. This variable is, thus, the upper bound of the flying time of each drone (Constraints (5)).

Constraints [START_REF] Zorbas | On the optimal number of chargers in battery-less wirelessly powered sensor networks[END_REF] ensure that the harvesting energy is greater than the energy requirements of each sensor. To do so, we compute for each position j ∈ J the harvesting power for node i ∈ I. We then verify that the sum over the positions harvesting power greater than Γ to i is enough to fulfill the node's requirement E i . Constraint [START_REF] Dai | Radiation constrained wireless charger placement[END_REF] ensures that the maximum budget is not exceeded. We can place each drone at different locations during the operation (Constraints (8)), but the total length of its flying time cannot exceed the imposed time limit (i.e., τ ). Only one drone can be placed at a given position (Constraints 9) in order to limit the movement of the drones through time. Constraints [START_REF] Zorbas | Assessing the cost of rf-power harvesting nodes in wireless sensor networks[END_REF] verify that if a 2Dlocation with a corresponding altitude is chosen for a drone, then no other drone can be placed above or under it. In this way we avoid multiple drones to be aligned at different altitudes and, thus, limit physical collisions between the drones as well as interference.

OD3P problem is N P -hard since it is related to the set covering problem [START_REF] Zorbas | Computing optimal drone positions to wirelessly recharge iot devices[END_REF]. There is no algorithm that solves exactly OD3P in polynomial time. Due to the limited size of this paper, we omit a formal proof of hardness.

Solving the proposed model is challenging since there is no straightforward solution to determine the optimal recharge time while minimizing the drone's altitude. Indeed, placing the drones at low altitude the total recharge time can be improved since the harvesting power is increased. When we have enough available drones to recharge the sensors (i.e., Constraint [START_REF] Dai | Radiation constrained wireless charger placement[END_REF] is not tight), the model will use them to reduce the average altitude of the selected positions. However, given the number of available drones, a trade-off exists here; by reducing the altitude of the drones, the coverage is limited and, thus, more successive drone positions are needed. On the contrary, due to lifetime limitations of the drones, the number of visiting positions may be also limited. This trade-off is captured by the model when computing the time associated to the position (i.e., the value of variable t n j ). As far as more drones are utilized, the model will find the best trade-off between the number of positions visited by a drone (minimizing the sum of t n j for all visited j), and the harvesting efficiency of each selected position (maximizing the sum of H t n j i for all j by selecting the appropriate altitude).

V. OD3P GREEDY HEURISTIC ALGORITHM

The complexity of MILP solution requires the use of a faster algorithm when dealing with topologies with a high number of nodes or drone positions. Hence, we propose a greedy heuristic algorithm to overcome this issue and approach the optimal solution.

The greedy algorithm (see Algorithm 1) computes a set O containing the set of sensors in decreasing order according to their energy needs E i . The algorithm, then, selects the first encountered sensor in O that has the smallest set of possible harvesting 3D positions. This set is computed by comparing the received harvested power at node i to threshold Γ. Every time a sensor is chosen, the 3D position with the highest harvesting power is selected. Given this position, the algorithm computes the maximum time length that the drone can spend at this position. We, thus, obtain a triplet (i * , j * , t * ) composed of a sensor, a position, and an associated harvesting time.

If the position can be associated to a drone n * , then the energy costs of the reached sensors from j * are updated and the position is marked as unavailable so that it cannot be further chosen. The sensors that are fully recharged are removed from the ordered set O. T j is the available time for each position j ∈ J . It depends on the time associated with same 2D-positions in order to ensure Constraints [START_REF] Zorbas | Assessing the cost of rf-power harvesting nodes in wireless sensor networks[END_REF]. If the cumulative flying time among same 2D positions reaches the time limit, then all these positions are removed so that they cannot be chosen. If no drone can go to position j * with associated time t * (i.e. if there is not enough time left), then position j * is simply discarded. T n is the available time for drone n ∈ N . It is bounded by the maximum time period τ , and seeks to be minimized. At the end of the first step of the algorithm, the flying time of the drone n * is updated. At the end of the algorithm, we compute the minimum value λ of the recharging time phase.

Algorithm 1 Greedy heuristic algorithm for OD3P

1: Input: sensors I with associated energy needs E i , i ∈ I, drones N , positions J , budget B 2: Output: a subset S n ⊂ J of positions with associated time duration t n j for each drone n ∈ N minimizing max n∈N j∈J t n j 3: O : set of sensors in decreasing order of E i 4: J i = {j ∈ J : P hi ≥ Γ} in decreasing order of P hi 5: T j = 0 for all j ∈ J 6: T n = 0 for all n ∈ N 7: while O = ∅ do

8: i * = { first sensor in O with minimum |J i |} 9: j * = first position in J i 10: t * = min{τ -T j , E i * /P h i * } 11: n * = min{T n for all n with T n + t * ≤ τ } 12:
if n * has been found then 13:

T n * = T n + t * 14:
Add j * and t * to S n 15:

E i = max{0, E i -P hi * t * } for all i ∈ O 16: O = O \ {i ∈ O, E i = 0} 17:
Re-order O in decreasing order of the updates E i 18:

for all j ∈ J , (x j , y j ) = (x j * , y j * ) do 19:

T j = T j + t * 20: if T j = τ then 21:
Remove j from J i for all i ∈ I 22: Remove j * from J i for all i ∈ I 23: return S n for all n ∈ N and λ = max n∈N j∈Sn t n j The computational complexity of the greedy algorithm depends on the number of nodes and the complexity of the implementation of the sorting algorithm (O(n log n) in best case). Since the algorithm seeks one node at each iteration (see line 7), and every iteration is computed in O(m log m) run time, where m the number of nodes remaining in set O, the total running time will be |I| i=1 i log i.

VI. PERFORMANCE EVALUATION & DISCUSSION OF THE

RESULTS

A. Simulation setup & evaluation metrics

We assume that the nodes are randomly placed on a 2dimensional grid of size 50m × 50m. We vary the number of nodes between 5 and 50, and the budget (in terms of number of drones) from 3 to 10. The number of possible drone positions is chosen between 125 and 500 (i.e., 5 × 5 to 10 × 10 grid points × 5 altitudes). The energy requirement of each sensor is set to 150mJ. For instance, this is the minimum amount of energy a Zolertia Z1 node 1 needs in order to bootup, take a measurement, and transmit it to the gateway. The 1 https://github.com/Zolertia/Resources/wiki/The-Z1-mote drones can cover the nodes' energy demands by transmitting 3W EIRP power (915MHz) to the nodes, which according to our adopted path-loss model and the manufacturer's data results in approximately Γ = 6.3 × 10 -5 , 0.01 W of harvested power depending on the distance. The power harvesting values are computed given the efficiency of the Powercast P2110B modules2 operating at 915MHz. The efficiency varies from 0 to about 65% depending on the received power. The chosen parameters for the harvesting model are P 0 = 10mW , b = 1.05, σ = 1, and η = 0.3. The maximum recharging time of a drone is set to 1 hour.

For each size of I, J and budget B ≤ |N |, we compute 10 different random topologies. We summarize the results with the mean value for each topology size. We did not plot the 95% confidence interval to keep the plot readable since they are very close to the mean value. The model is implemented in Java language and solved using the IBM CPLEX solver 12.8 on an Intel Core i7-5500U CPU, 2.40 GHz, 16 Gb RAM computer under Microsoft Windows 10 operating system.

B. Study of the optimal energy replenishment time

In this study we compute the optimal time needed to replenish the energy of the sensors for different numbers of used drones, sensors, and 3D-positions (grid positions). The results are computed by solving the OD3P optimization problem and are depicted in Fig. 2. They reveal the model's ability to compute optimal values even for a high number of sensors and drone positions. On one hand, as it was expected, the recharge time decreases as more drones are added in the field or more drone positions are taken into account. More drones are also needed for higher node populations. Moreover, the optimal number of drones varies between 7 and 9 drones. On the other hand, the recharge time tends to remain constant or increase as more drones are added above the optimal point. This happens because once the optimal solution has been found, the model does not allow placing drones at different altitudes while covering the same number of sensors. As a consequence, some drones are slightly moved from the optimal position so that the rest of the drones can fit in the 3D field. This results in a constant or sometimes higher recharge time.

C. Study of the drone positions

In the next study, we compute the number of successive positions of the drones throughout the process as well as the mean recharge time per position. We remind that each drone is allowed to change positions in order to recharge as many sensors as possible within the imposed time limit. As we can observe from Fig. 3a, a drone moves on average to several positions in order to fulfill the energy requirements of the nodes (Constraints (6)). More positions are selected as more nodes are added in the field. Moreover, the number of positions decreases when increasing the drone population. Despite that, it is interesting to observe that the number of positions converges to 2 or 3. This number is related to the minimum time required to recharge a single node at the lowest altitude as well as to the imposed energy replenishment time. This can be also explained by dividing the results of Fig. 2a by those of 3b, which will give us the mean number of drone positions.

Placing the drones at the lowest possible altitude, the harvesting energy is maximized due to the increased power and the high efficiency of the harvesting module. However, this action decreases coverage, thus, more drones are needed in the field. In other words, as explained in Section IV there is a trade-off between coverage and recharging time that the computational model has to solve. Fig. 3c reflects the behavior of the model while dealing with this trade-off. The figure reports the mean and the maximum altitude of the drones for the examined scenarios. The results witness that even if there are selected positions at 3m (for |J | = 125 or 250), or 5m (for bigger |J |), the average altitude of all the selected positions barely exceeds 1m. This shows that low altitude positions are promoted by the model.

D. Scalability study & comparison

In this study we examine the computational model's scalability and we compare MILP's performance to that of the greedy algorithm.

Fig. 4b presents the execution time of MILP for different drone positions and populations. The MILP optimally solves all the examined instances. However, the resolution time presents an exponential increase as the number of drone positions and the number of drones increases. However, as depicted in Fig. 4a and4b, it needs only a few seconds for scenarios with up to 20 nodes and up to 245 drone positions, respectively. On the contrary, the greedy algorithm solves quasi instantly all the topologies. Indeed, the maximum captured running time of the algorithm is 72ms. We must note that, in terms of total recharge time, the greedy approach achieves the optimal value for only 8.2% of the instances, while about 17% of the solutions are up to 10% longer than the optimal. The optimal value is usually reached for small number of drone positions (i.e., 125, 180, and 245). When the problem complexity increases greedy is on average 26.07% far from the optimal (Fig. 5a). Moreover, as we can observe from Fig. 5b and5c, the solutions computed by the heuristic impose on average a lower drone altitude as well as a lower recharge time per position. This result reveals the ability of the model to provide a better trade-off between the altitude and the optimal recharge time compared to the heuristic algorithm.

VII. CONCLUSION AND FUTURE WORK

In this paper, we presented an optimal linear formulation of the drone placement and planning problem in order to minimize the energy replenishment time of ground sensors. We developed and optimally solved a Mixed Integer Linear Program (MILP) to tackle the problem. We also proposed a heuristic algorithm capable of solving large instances of the problem in a very short amount of time. The originality of our work lies on the focus of our approach which is the (a) Variable number of sensors. minimization of the recharging time, given a number of mobile drones. The simulation results showed the effectiveness of our strategy both in terms of deployment and computation cost. They also showed that the drones altitude is minimized to ensure good harvesting capabilities which gives good insights to the analysis of practical distributed algorithms and industrial scenarios. In the next steps of this work, we would like to investigate a fully localized algorithm for the drone placement and planning problem. We will also focus on integrating the drones' trajectories inside the MILP and the heuristic formulation.

Fig. 1 :

 1 Fig. 1: An example employing 7 nodes and 2 drone positions.

  Variable number of drones & sensors.

  Variable drone positions & populations.

Fig. 2 :

 2 Fig. 2: Optimal energy replenishment time for (a) variable drone and sensor populations, and (b) variable drone positions and drone populations.

  Number of successive positions.

  Mean recharge time per position.

  Maximum and mean altitudes.

Fig. 3 :

 3 Fig. 3: Analysis of the computed positions per drone & time per position.

  Variable number of drone positions.

Fig. 4 :

 4 Fig. 4: Resolution time for variable number of nodes and drone positions.

  Mean number of successive drone positions.

Fig. 5 :

 5 Fig. 5: Comparison with the greedy algorithm.

http://www.powercastco.com/wp-content/uploads/2016/12/P2110B-Datasheet-Rev-3.pdf

ACKNOWLEDGMENTS

This publication has emanated from research conducted with the financial support of Science Foundation Ireland (SFI) and is co-funded under the European Regional Development Fund under Grant Number 13/RC/2077.

It has also received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skodowska-Curie grant agreement No. 713567, and by the ANR program Investments for the Future under reference ANR-11-LABX-0031-01.


Christelle Caillouet
email: christelle.caillouet@unice.fr

Frédéric Giroire
email: frederic.giroire@cnrs.fr

Tahiry Razafindralambo
email: tahiry.razafindralambo@univ-reunion.fr

Efficient Data Collection and Tracking with Flying Drones

Keywords: UAV, Coverage, Wireless Sensor Network

HAL is

Introduction

Wireless sensor networks (WSNs) have tremendous applications in environmental monitoring or surveillance. Within the context of the Internet of mobile things, this paper focuses on observing a set of mobile targets, i.e., sensors, using a fleet of Unmanned Aerial Vehicles (UAVs), i.e., flying drones. Specific applications such as wildlife monitoring [START_REF] Dyo | Wildsensing: Design and deployment of a sustainable sensor network for wildlife monitoring[END_REF], vehicle observation and tracking [START_REF] Guido | Evaluating the accuracy of vehicle tracking data obtained from unmanned aerial vehicles[END_REF] in smart cities and monitoring sporting events [START_REF] Hayat | Survey on unmanned aerial vehicle networks for civil applications: A communications viewpoint[END_REF][START_REF] Cambra | Ad hoc network for emergency rescue system based on unmanned aerial vehicles[END_REF] can benefit from the results presented in this paper. In our context, the flying drones act as a wireless mobile backbone, covering the targets on the ground and collecting information from them. Gathered data is continuously sent to a base station for further processing using a multi-hop communication paradigm through the wireless backbone. Our problem can be described as an optimization of spacial and temporal coverage with mobile flying drones and is named the Aerial Data Collection Problem (ADCP).

The use of mobile flying devices to cover mobile ground targets has become an important topic in the past few years. The problem tackled in this paper is crucial since it allows the monitoring system to optimally adapt to the evolving space and time placement of the ground targets, but also to the changing requirements of the application itself. Therefore, deployment efficiency and cost can be jointly optimized for a wide adoption of wireless mobile monitoring systems. Figure 1 shows the studied multi-tier network architecture composed of mobile sensors, flying drones, and a fixed base station. The sensors are located on the ground and produce data. Their mobility pattern is usually unknown. Data is then gathered by a fleet of UAVs whose mobility is fully controlled in order to track the sensors. Each ground sensor must be covered at any time by at least one UAV for ground-to-air communications. Each UAV moves in a 3D-space and their altitude must be managed for the coverage and data collection. When a UAV is at a high altitude, the covered area is wider on the ground. However, altitude increases the distance between UAVs and ground sensors and thus reduces the ground-to-air communication transmission quality [START_REF] Al-Hourani | Modeling air-to-ground path loss for low altitude platforms in urban environments[END_REF]. Finally, data gathered by UAVs are sent to a base station in a multi-hop fashion. A connected backbone of UAVs and the base station must be formed and maintained at anytime to transmit the data.

Deployment of UAVs is studied both from theoretical and practical considerations. From a practical point of view, research focuses on distributed UAV deployment, especially the optimal way of moving one UAV from one point to another. The robotics community is strongly involved in this research topic. From theoretical and practical points of view, the networking community and especially the wireless sensor network community focus on designing the appropriate communication protocol for UAV mobile networks. From a theoretical point of view, research focuses on optimal or approximated computations of the positions of each UAV. Multiple works propose exact mathematical formulations or heuristic algorithms for the coverage and positioning problems using UAVs with several objective functions. The work presented in this paper lays in this last category. However, in this paper we try to encompass positioning aspects for coverage, communication aspects, especially connectivity, and constraints related to UAV movement. This paper tries to fill the gap between the networking approach and the theoretical approach by providing and describing a tool and framework for the evaluation and comparison of data collection and tracking strategies algorithms for wireless flying drones. Indeed, the results of this paper can be used and modified to serve as a basis for the networking community to evaluate the performance and enhancement an algorithm can bring.

This paper is organized as follows. In Section 2, we present an overview of related works and compare our paper to specific results in the literature. In Section 3, we describe the studied ADCP with its assumptions and the considered scenario. In Section 4, we present a mixed-integer linear program optimally solving ADCP. However, the model does not allow us to solve large instances of the problem. We thus propose a scalable decomposition model in Section 5. This new formulation can be solved efficiently using column generation. We present the results in Section 5.4 validating our model and studying the evolution of the UAV deployment over time. A heuristic enhancing the resolution time of the column generation process is presented and validated in Section 6. Section 7 discusses the model assumptions and possible applications, together with potential future work. We finally conclude this paper in Section 8.

State of the art

This section is dedicated to the position of our paper regarding different aspects of the literature: global fundamental approaches (collection or dissemination), types of mobility patterns (random, constrained, fully controllable), community specificities (robotics, networking, operational research), and finally specific assumptions (3D-positioning, coverage, data gathering, and objectives).

Global positioning of our contribution

Mobile nodes such as UAVs (Unmanned Aerial Vehicles) have been introduced into WSNs (Wireless Sensor Networks) for multiple purposes. There are two main approaches: 1) The first one is the up-link process in which the set of mobile nodes help the WSN in the data collection/gathering process from sensors to a central entity [START_REF] Macharet | Bi-objective data gathering path planning for vehicles with bounded curvature[END_REF]. 2) The second approach is to help the down-link process in which the set of mobile nodes help a central entity, for example the base station, disseminate/deliver information, update software or recharge sensors' batteries [START_REF] Chen | Charge me if you can: Charging path optimization and scheduling in mobile networks[END_REF] which implies visiting all the sensors for a given period of time or at a given time. In this paper, we focus on the data collection/gathering process and cover this aspect from different perspectives. One important aspect is the type of mobility.

The literature is very rich regarding mobility management. There are three main mobility categories including random mobility, partially controllable mobility and fully controllable mobility of mobile nodes. In this paragraph, we provide some examples of each category.

• Random mobility (not controllable). In this category, results cannot provide guarantees on data collection/gathering and solutions are often subject to long delays. Early papers such as the one by Shah et al. [START_REF] Shah | Data mules: modeling and analysis of a three-tier architecture for sparse sensor networks[END_REF] and Hamida et al. [START_REF] Hamida | A line-based data dissemination protocol for wireless sensor networks with mobile sink[END_REF] lie in this category and provide a first approach of data collection/gathering in wireless sensor networks using mobile nodes. Unlike the results from [START_REF] Shah | Data mules: modeling and analysis of a three-tier architecture for sparse sensor networks[END_REF] and [START_REF] Hamida | A line-based data dissemination protocol for wireless sensor networks with mobile sink[END_REF], we are interested in data gathering guarantees and optimality.

• Partially controllable mobility. Trajectories and/or stop positions are constrained. This case is often environment specific, for example bus trajectories in a city. In this case, data gathering/collection can suffer from high delays. The paper from Gandham et al. [START_REF] Gandham | Energy efficient schemes for wireless sensor networks with multiple mobile base stations[END_REF], from Basagni et al. [START_REF] Basagni | Moving multiple sinks through wireless sensor networks for lifetime maximization[END_REF] and from Luo et al. [START_REF] Luo | Joint sink mobility and routing to maximize the lifetime of wireless sensor networks: The case of constrained mobility[END_REF], were among the first to provide some solutions in this category. Unlike [START_REF] Gandham | Energy efficient schemes for wireless sensor networks with multiple mobile base stations[END_REF], [START_REF] Basagni | Moving multiple sinks through wireless sensor networks for lifetime maximization[END_REF] and [START_REF] Luo | Joint sink mobility and routing to maximize the lifetime of wireless sensor networks: The case of constrained mobility[END_REF], we want a data gathering process with minimal delay. Each sensor has to be permanently in contact with a mobile node in order to transmit the data directly after their production.

• Fully controllable mobility. It widens the application usage of wireless sensor networks. It is interesting to focus on fully controllable mobility since it is possible to compute the node's trajectory to optimize the collection process [START_REF] Wichmann | Smooth path construction for data mule tours in wireless sensor networks[END_REF]. The contributions of this paper lie under this category. Unlike the paper from Wichmann et al. [START_REF] Wichmann | Smooth path construction for data mule tours in wireless sensor networks[END_REF], we use multiple mobile nodes for data gathering. Unlike other works such as the one presented by Huang et al. [START_REF] Huang | Viable path planning for data collection robots in a sensing field with obstacles[END_REF], the data gathering process in this paper is time and space dependent since the data producers (sensors) are mobile and should therefore, be tracked and covered during the deployment procedure.

Community positioning of our contribution

The data collection process with mobile nodes or UAVs with fully controllable mobility have already been studied by different communities such as the wireless sensor network community, the robotics community and the operational research community. Each community has its own view and way to tackle the data gathering process, focusing on different issues,

h h h h h h h h h h h h h h h

Constraints/Objectives Papers [START_REF] Mozaffari | Efficient deployment of multiple unmanned aerial vehicles for optimal wireless coverage[END_REF] and thus, putting more importance on some constraints, objectives or characteristics. i) The wireless sensor network community focuses on the communication protocol quality which may include connectivity of the network [START_REF] Razafindralambo | Connectivity preservation and coverage schemes for wireless sensor networks[END_REF]. For example, given a location for each device, the objective is to implement protocols above the deployed network [START_REF] Becvar | Performance of mobile networks with uavs: Can flying base stations substitute ultra-dense small cells?[END_REF][START_REF] Razafindralambo | Spread and shrink: Point of interest discovery and coverage with mobile wireless sensors[END_REF]. In these papers, the main focus is to implement algorithms for data gathering. The deployment and its optimality are not considered. An interesting point from [START_REF] Razafindralambo | Connectivity preservation and coverage schemes for wireless sensor networks[END_REF][START_REF] Becvar | Performance of mobile networks with uavs: Can flying base stations substitute ultra-dense small cells?[END_REF][START_REF] Razafindralambo | Spread and shrink: Point of interest discovery and coverage with mobile wireless sensors[END_REF] is the fact that permanent connectivity in the UAV or mobile robots network is important and should be guaranteed. In our work, we also consider this point and add another constraint related to optimal coverage and position cost. ii) The robotics community mostly focuses on the robots' physical and mobility constraints such as turn curvature or obstacle avoidance [START_REF] Ny | On the dubins traveling salesman problem[END_REF], and collaboration among mobile devices for a specific task [START_REF] Zavlanos | Graph-theoretic connectivity control of mobile robot networks[END_REF]. Moreover, it focuses on the deployment itself or how to drive each device to its specific location while maintaining connectivity among them and ensuring a good coverage of the area to monitor [START_REF] Yanmaz | Connectivity versus area coverage in unmanned aerial vehicle networks[END_REF]. In [START_REF] Ny | On the dubins traveling salesman problem[END_REF][START_REF] Zavlanos | Graph-theoretic connectivity control of mobile robot networks[END_REF][START_REF] Yanmaz | Connectivity versus area coverage in unmanned aerial vehicle networks[END_REF], authors put a focus on velocity, obstacle avoidance, and UAVs or robot movement. Unlike our work, these papers do not focus on the application itself and how to optimally perform the data gathering process. However, it would be interesting to consider highlevel constraints on movement in future works. iii) The operational research community is mostly focused on high-level path planning, looking for optimal trajectory [START_REF] Dixon | Optimizing cascaded chains of unmanned aircraft acting as communication relays[END_REF][START_REF] He | A progressive approach to reducing data collection latency in wireless sensor networks with mobile elements[END_REF], optimal topology construction [START_REF] Alduraibi | Coverage-based node placement optimization in wireless sensor network with linear topology[END_REF] or to maximize end-to-end throughput in UAV networks [START_REF] Chen | Optimal positioning of flying relays for wireless networks: A los map approach[END_REF]. This community mostly proposes mathematical formulations for the coverage (which can be coupled with data gathering) and positioning problems using UAVs to find optimal or approximated solutions. Most contributions from this community simplify the network and robotics aspects and do not consider the network's connectivity during the data gathering process nor the 3D-placement of mobile robots/UAVs ( [START_REF] Dixon | Optimizing cascaded chains of unmanned aircraft acting as communication relays[END_REF][START_REF] He | A progressive approach to reducing data collection latency in wireless sensor networks with mobile elements[END_REF][START_REF] Alduraibi | Coverage-based node placement optimization in wireless sensor network with linear topology[END_REF][START_REF] Chen | Optimal positioning of flying relays for wireless networks: A los map approach[END_REF]).

Some contributions from the literature are very interesting since they try to mix the constraints, objectives, characteristics, and needs from different communities. For example, [START_REF] Huang | Viable path planning for data collection robots in a sensing field with obstacles[END_REF] mixes approaches from the operational research community and the robotics community by providing an optimal path planning solution for data gathering with multiple robots while taking into account the physical constraints such as velocity or turning curvature. In the same way, [START_REF] Razafindralambo | Connectivity preservation and coverage schemes for wireless sensor networks[END_REF] takes into account the physical constraints of a robot while guaranteeing wireless mobile sensor network connectivity. In [START_REF] Quintero | Optimal uav coordination for target tracking using dynamic programming[END_REF], dynamic programming is used to find an optimal routing for two camera-equipped UAVs, cooperatively tracking a single target moving on the ground. In [START_REF] Kar | Node placement for connected coverage in sensor networks[END_REF], the authors provide a study on node placement including connectivity constraints, but the model is not optimal and is based on approximation. Moreover, the algorithms are designed for 2D-space with fixed coverage radius. In [START_REF] Mozaffari | Efficient deployment of multiple unmanned aerial vehicles for optimal wireless coverage[END_REF], the authors consider maximizing the total coverage area of the UAVs and their lifetime. But, in the models, the UAVs are not connected and assumed to be placed at the same altitude and thus, do not consider real 3D-positioning. [START_REF] Pugliese | Modelling the mobile target covering problem using flying drones[END_REF] uses the same assumptions as we do in this paper. However, they to not consider connectivity between UAVs. In our paper, following the previous work from [START_REF] Caillouet | Optimization of mobile sensor coverage with uavs[END_REF], we provide a solution that considers constraints from the wireless sensor community, namely connectivity and communication quality, while computing the optimal placement of each UAV. Moreover, we show that the classical Aerial Data Collection Problem (ADCP) formulation is outperformed by the connected set formulation for bigger instances of ADCP, and we propose a heuristic algorithm to enhance resolution time while relaxing the optimality condition.

Specific positioning of our contribution

In this section, we position our paper regarding relevant specific objectives and constraints of ACDP, and with the related literature. We briefly expose the originality of our paper compared to the state of the art. The following points highlight the specificity of our work and Table 1 shows the positioning compared to the literature. This table highlights that, compared to the literature, the models we described here include different assumptions and constraints. Some of the assumptions and constraints were individually tackled in the literature but our model considers them together in a complete framework.

• application: The application we consider is the coverage of ground sensors by UAVs. Therefore, we consider the 3D-positioning of each UAV which affects the coverage range on the ground since ground coverage depends on the UAVs' altitudes.

• sensors' characteristics: We assume that ground sensors are mobile (random unknown mobility). We want to ensure that all sensors are covered by a UAV at any time.

• data gathering: Data gathered by the UAVs should be transmitted by the UAVs to a base station. Since the distance between the UAVs and the base station can exceed the communication range of a UAV, the network of UAVs should be connected in a multi-hop fashion to the base station at any time.

• objectives: We want to minimize the number of UAVs for global cost minimization of ADCP, and for better management and coordination of the UAVs. We also consider the minimization of their respective altitudes and positions to keep them as close as possible to the base station in order to sustain the traffic. Minimizing the number of UAVs in real deployment is important, since it reduces the deployment cost. From a networking and communication point of view, reducing the number of UAVs also limits the possible communication interferences. It is important to notice that reducing the number of UAVs may affect robustness, but this point is left for future work.

Problem definition

From an optimization point of view, data collection from a set of nodes to a central entity has generally been modeled by two different approaches [START_REF] Somasundara | Mobile element scheduling for efficient data collection in wireless sensor networks with dynamic deadlines[END_REF][START_REF] Hung | Energy-aware set-covering approaches for 13 approximate data collection in wireless sensor networks[END_REF]. First, the problem can be viewed as a vehicle routing problem in which the base station is the collecting point, and the goal is to compute the route of each vehicle, visiting all the customers once [START_REF] Bullo | Dynamic vehicle routing for robotic systems[END_REF]. Second, the problem can be viewed as a set covering problem in which a subset of nodes is linked to all the nodes in a graph [START_REF] Popescu | Multirobot patrolling in wireless sensor network using bounded cycle coverage[END_REF]. In our case, we choose to model our aerial data collection problem as a set covering problem since the UAVs evolve in the air with an associated altitude. Moreover, since the ground targets are mobile, we consider its dynamic version [START_REF] Chrissis | The dynamic set covering prblem[END_REF], minimizing deployment cost and data collection delay. We also use equations related to flow problems in order to ensure the connectivity and coverage for ADCP, as seen in the next sections.

ADCP assumptions

Each UAV u is located in the three-dimensional space and P is the set of possible 3D locations. Let p = (x u , y u , h u ) ∈ P be the position (x u , y u ) of UAV u in the 2D-plane and h u its altitude. Time is discretized so that the target positions at each time are estimated and given by the successive sets N t , for t ∈ [0, T ], where T is the monitoring time period length. Each target n ∈ N t is located on the ground and associated with twodimensional coordinates (x t n , y t n ). We compute the observation radius r h u of UAV u depending on its altitude h u and its directional antenna half beam-width, or visibility angle θ [START_REF] Al-Hourani | Modeling air-to-ground path loss for low altitude platforms in urban environments[END_REF]. The coverage area of UAV u on the 2D-plane is represented by a disk of radius bounded by:

r h u ≤ h u • tan θ 2 .
We assume that a UAV u deployed at location p ∈ P covers target n ∈ N t if the euclidean distance between its projection on the 2D-plane and the target d(u, n) is below the observation radius of the UAV:

d(u, n) = (x u -x t n ) 2 + (y u -y t n ) 2 ≤ r h u .
Similarly, for air-to-air communications, a UAV u can communicate with another UAV v if their distance on the 3D plane

D uv = (x u -x v ) 2 + (y u -y v ) 2 + (h u -h v ) 2 ≤ R u ,
where R u is the communication range of UAV u. For efficient data collection, we enforce the deployed UAVs to be connected with each other (through possible multi-hop communications) and with the base station b located on a corner of the area.

Dynamic graph model

Given the target sets N t , the coverage radius, and the communication range of the UAVs, one can derive a dynamic representation of the topology of our problem. Due to their dynamic and evolving characteristics, these networks are best represented by dynamic graph models [START_REF] Zaki | Comprehensive survey on dynamic graph models[END_REF]. For ADCP, we define a set of graphs G t = (V t , E t ) in which the topology changes at each time t ∈ [0, T ]. Definition 1. At time t ∈ [0, T ] of the observation time period, the current topology of ADCP is modeled by the directed graph G t = (V t , E t ) where:

• V t = {b} ∪ P ∪ N t ,

and

• E t = {(b, p), p ∈ P and (x p , y p ) = min q∈P d(b, q)} ∪ {(p, q), p, q ∈ P and R p ≥ D pq } ∪ {(p, n), p ∈ P, n ∈ N t and r h p ≥ d(p, n)}.

The set V t contains the base station, all the positions for deploying drones, and the location of the targets at time t. Edges in E t are defined to ensure (i) the connectivity of the backbone induced by the UAVs and the base station, and (ii) the coverage of all the targets at their current location. We assume that the base station has an edge to the positions p for which the projection on the 2D plane is the closest to the base station.

ADCP objective

In this section, we define a particular cost representation used in the objective function of ADCP involving considerations about cost, wireless communications, and delay. First, the management and coordination of the UAVs represent a specific difficult task to optimize. We thus seek to minimize the global cost of using a fleet of UAVs by minimizing the number of deployed UAVs over time. Moreover, wireless communication quality depends on various parameters such as environmental characteristics, wireless interferences, distance between devices, etc. Deploying the UAVs at low altitude maximizes the communication quality between the targets and the UAVs. However, it leads to more deployed UAVs, since the coverage area is smaller. The trade-off between these two antagonistic objectives has been studied in [START_REF] Caillouet | Efficient Deployment of Connected Unmanned Aerial Vehicles for Optimal Target Coverage[END_REF] for static targets.

In this paper, we define a cost parameter combining the different targeted optimization objectives: Definition 2. To each position p ∈ P, we associate a cost depending on the 3D-location:

C p = D bp = (x u -x b ) 2 + (y u -y b ) 2 + (h u ) 2 ,
which is the distance in the 3D-plane between the position p and the base station b.

Using this cost in the objective function of our optimization model jointly optimizes the number of UAVs, their distances from the base station, and their altitudes (Figure 2). This allows for a better management of the fleet of flying drones, limiting their energy consumption over time and the delay of data collection. With the objective of determining a UAV positioning at each time step and minimizing the sum of their associated cost C p , our optimization model of ACDP ensures that, (i) the UAVs stay close to the base station for efficient data collection, (ii) the UAVs are at the lowest altitude possible to ensure good connectivity with the sensors, and (iii) we use the minimum number of UAVs.

Scenario and tools description

Models developed in this paper have been implemented in the Java language using the JGraphT library 1 , solved using IBM Cplex solver 12.7.1 on an Intel(R) Core(TM) i7-5500U CPU, 2.40 GHz, 16 Gb RAM machine, running Microsoft 8.1 Professional operating system.

Instances are deployed in a square area of size 100 m × 100 m. We consider 10 time slots for target mobility. The 2Dcoordinates of the targets are initially chosen randomly in the monitored area. Their final destination and their velocity are also randomly chosen. At each time step, the next position of each target is computed based on their previous location, their velocity, and their direction. If a target reaches its destination coordinates before the last slot, it stays at its final location for the remaining slots.

We divide the monitored area into equal squares. One possible 2D-position of a UAV is located in the center of each square. In this way, the candidate sites for placing a UAV form a regular square grid. For each 2D-coordinate of point p, (x p , y p ), we set the allowed altitudes to {10 m, 25 m, 45 m}. The base station is placed on the bottom left corner at coordinates (0.0, 0.0) with an altitude of 0 m.

Our instances contain between 5 and 50 targets and between 75 and 300 possible UAV locations (i.e., from 5 × 5 2Dcoordinates ×3 altitudes to 10×10 2D-coordinates ×3 altitudes). The UAV visibility angle θ is set to 60 degrees, and their communication range to 30 m. Table 2 summarizes the parameters used in our scenario.

Classical ADCP formulation

Mixed integer linear programming formulation (MILP)

Let z t p be a binary variable stating if position p ∈ P is selected at time t to deploy a UAV. Classical ADCP can be modeled as a flow problem, given the temporal graphs of Definition 1. The existence of a flow in the graph G t from a source node, the base station b, to destination nodes, the targets of N t , ensures that we obtain a connected set of flying drones with the base station, and such that the selected locations for the UAVs cover all the targets at time t. Indeed, there exists an edge between a UAV position and a target if and only if the target is located within the coverage area of the UAV. We thus define f t uv ∈ R as the amount of flow going through link (u, v)

∈ E t at time t ∈ [0, T ].
The linear formulation of ADCP, called ADCP MILP (for ADCP Mixed Integer Linear Program), is the following :

min         T t=0 p∈P C p z t p + max t∈[1,T ] | p∈P C p z t p - q∈P C q z t-1 q |         (1)
q∈V t ,q p f t pq -

q∈V t ,q p f t qp =          |N t | if p = b 0 if p ∈ P -1 if p ∈ N t , ∀t ∈ [0, T ], p ∈ V t (2) f t pq ≤ z t p • |N t |, ∀t ∈ [0, T ], (p, q) ∈ E t with p ∈ P (3) f t pq ∈ R, z t p ∈ {0, 1} (4)
Objective ( 1) not only minimizes the global deployment cost over time, but also the changes between two consecutive time slots. The metric used here is the one introduced in Definition 2. Recalling Section 3.3, the goal is to maintain a minimum number of UAVs to monitor the mobile targets, while placing them at a minimum distance from the base station. The first part of the objective function (1) ensures that we deploy the minimum number of UAVs as close to the base station as possible, ensuring the lowest altitude to cover the sensors and the fastest delay to gather the monitored information. The second part seeks to minimize the changes of selected positions. We want to minimize the difference in terms of cost (Def. 2) for two consecutive time slots, in order to minimize the delay of data collection between the selected positions over time. If positions need to change due to target mobility, then we ensure that the UAVs remain the closest possible to the base station, corresponding to the departure location of the drones. Then, minimizing their deployment cost over time gives a solution with limited UAV moves over time.

To get rid of the absolute value in the objective function, we introduce a new variable λ ≥ 0. Objective (1) becomes : min

        T t=0 p∈P C p z t p + λ         , (5)
and we add the following constraints to the program:

λ ≥ p∈P C p z t p - q∈P C q z t-1 q , ∀t ∈ [1, T ] (6)
λ ≥ p∈P C p z t-1 p - q∈P C q z t q , ∀t ∈ [1, T ] (7)
Minimizing λ thus gives the minimum absolute value of the cost evolution between two consecutive time slots. Constraints (2) model the flow conservation. We want to ensure the existence of a flow of size |N t | sent by the base station b to the targets of N t . The flow is forwarded by the UAVs until the targets receive one unit each. Since, a priori, we do not know the locations of the UAVs, we define constraints (3) ensuring that no flow can go through a link incident to an unchosen position. Consequently, if a location p ∈ P is not chosen, i.e., z t p = 0, then no flow exists on its incident links.

Results and discussions

In Table 3, we can see that even for small topologies with 5 targets and 75 possible locations for the UAVs (i.e., 5 × 5 coordinates ×3 altitudes), ADCP MILP takes several thousands of seconds (∼ 4500s) to compute the optimal solution. However, when the number of targets is doubled from 5 to 10, then resolution time is multiplied by 3 for |P| = 75. When the number of possible 3D-locations is increased from 75 to 108, the resolution time is multiplied by 30 for |N t | = 5. Consequently, for larger numbers of targets or possible UAV positions, ADCP MILP does not scale and is unable to give a solution, as the CPLEX solver runs out of memory.

Indeed, this mixed-integer formulation suffers from the number of binary variables of the deployment constraints and time period length. Even for small networks and considering few 3D positions and time slots, the program generates MILPs with a huge number of constraints and integer variables. This combinatorial hardness makes the problem intractable for topologies with more than one hundred 3D-positions. In order to tackle larger topologies, we introduce an approach based on a decomposition model in the following sections.

P

Connected set (CS) formulation

We now present a new formulation of ADCP, called ADCP CG OPT (for ADCP Column Generation, with Optimal pricing policy), involving sets of selected locations for the UAVs. This formulation allows us to separate our problem into two parts: (i) a master problem in which we deal with time and ensure a complete coverage of the targets at each time step, and (ii) a pricing problem generating connected subsets of UAVs with the base station, covering the targets at a given time.

Variables of the master program do not depend on one specific UAV location, but on more sophisticated structures, i.e., subsets of deployed UAVs with specific characteristics. Therefore, we do not deal with possible positions for each UAV, but only seek to select one subset of connected UAVs at each time in order to cover the sensors. We give a formal definition of the connected set considered in this new model. Definition 3 (Connected set (CS)). Let S ⊆ P be a subset of 3D-positions with the following properties :

• One deployed UAV is associated with each 3D-position (x p , y p , h p ) ∈ S ;

• Nodes of S form a connected graph with the base station b, i.e., there exists a path between b and every u ∈ S .

We thus adapt the cost associated with each connected set, using Definition 2, in the following: Definition 4 (CS cost). For each subset of UAVs S , we compute an associated cost depending on the 3D-location of the UAVs:

C S = u∈S D bu ,
where D bu is the distance on the 3D-plane between UAV u ∈ S and the base station b introduced in Definition 2.

Master program

Let S denote the set of all possible CS S . S is exponential in size. The goal here is to prevent the enumeration of all possible subsets of UAVs using a known technique of optimization called column generation. Let z t S be a binary variable indicating if CS S is selected at time t. We introduce another set of binary variables χ t n to model the coverage of target n at time t. The linear master program modeling ADCP with S is the following:

min         T t=0 S ∈S (z t S • C S ) + λ         (8)
S ∈S

z t S = 1, ∀t ∈ [0, T ] (9) χ t n ≥ 1, ∀t ∈ [0, T ] , n ∈ N t (10) χ t n ≤ S ∈S         z t S • p∈S |{(p, n)} ∩ E t |         , ∀t ∈ [0, T ] , n ∈ N t (11) λ ≥ S ∈S (z t S • C S ) - S ∈S (z t-1 S • C S ), ∀t ∈ [1, T ] (12) λ ≥ S ∈S (z t-1 S • C S ) - S ∈S (z t S • C S ), ∀t ∈ [1, T ] (13)
The objective function ( 8) is an adaptation of the previous equation ( 5) to deal with S. Constraints ( 9) state that we must select exactly one subset S of UAVs at each time t. Constraints ( 10) and ( 11) ensure that all the targets are covered at each time t if they are within the coverage area of at least one UAV of set S . For each position p ∈ S , we verify if the link (p, n) exists in the set of edges E t of graph G t . Given Definition 1, it means that the target n ∈ N t is within the coverage area of a UAV located at position p ∈ P. Therefore, target n ∈ N t is covered at time t if set S is selected at time t and at least one point p of S covers n. Constraints ( 12) and ( 13) concern the absolute value in the original objective function as constraints ( 6) and ( 7) do in the classical formulation of Section 4.

Pricing program

The sub-problem, or pricing problem, aims at generating new CS fulfilling the property of forming a connected backbone of UAVs with the base station. It is based on the dual formulation of the master problem.

Given the master problem ( 8)-( 13), let β (i) be the dual variables associated with constraints (i). For each S ∈ S, t ∈ [0, T ], the associated dual constraint is of the form:

C S • (1 -γ) - n∈N t
β (11) nt p∈S

|{(p, n)} ∩ E t | ≥ β (9)
t where γ is a coefficient involving dual variables β (12) and β (13) depending on the value of t (i.e., γ = f (β (12) , β (13) )). Indeed, constraints ( 12) and ( 13) are defined for t ≥ 1 and involve variables z t S and z t-1 S . We get the following definition of γ:

γ =           
β (12) t+1 -β (13) t+1 if t = 0 -β (12) t + β (13) t if t = T β (12) t+1 -β (13) t+1 -β (12) t + β (13)

t if 0 < t < T
The pricing problem of ADCP thus seeks to compute new CS violating the dual constraints, to add to the set of columns of the master problem. In other words: Definition 5 (Minimum Weighted Connected Subset (MWCS)). Given weight functions β (i) and a time slot t ∈ [0, T ], the Minimum Weighted Connected Subset (MWCS) consists in finding a subset S ∈ S of UAVs for which

C S • (1 -γ) - n∈N t
β (11) nt p∈S

|{(p, n)} ∩ E t | is minimum.
Thus, a minimum weighted connected subset generation either gives a good candidate to add to the set of variables of the master problem, or proves that no such column exists. If the cost computed by the pricing problem is smaller than β (9) , the generated CS is added to the set of considered connected sets S, and the corresponding variables z t S are added for all t. Constraints of this sub-problem define the structure of the CS of UAVs, fulfilling Definition 3.

Given the set P of possible 3D-locations for the UAVs, the goal is to select a subset of locations to deploy UAVs, such that they form a connected backbone with the base station efficiently collecting data from the mobile sensors on the ground.

In order to improve the efficiency of the computed subset, we define a pricing problem for every t ∈ [0, T ], and we specify that the connected subset must cover the sensors of N t .

An optimal linear formulation of MWCS is similar to the classical formulation presented in Section 4, except that we remove the subscript t on each variable and only look for a single flow between the base station and the current target set N t . To do so, we introduce variables y p ∈ {0, 1} determining if location p ∈ P is selected to be included in the CS. The objective of the pricing program is to minimize the weighted connected subset cost depending on the dual constraint presented above. We ensure the existence of a flow between the base station and the targets n ∈ N t (Constraints (15)), and ensure that no flow can go through incident links of unchosen locations (Constraints ( 16)).

min p∈P y p         C p (1 -γ) - n∈N t
β (11) nt p∈P

|{(p, n)} ∩ E t |         (14) q∈P,q p f pq - q∈P,q p f qp =          |N t | if p = b 0 if p ∈ P -1 if p ∈ N t , ∀p ∈ V t (15) f pq ≤ y p • |N t |, ∀(p, q) ∈ E t with p ∈ P ( 16
)
f pq ∈ R, y p ∈ {0, 1} (17)

Resolution method

Column generation is a decomposition method that combines the resolution of a restricted master problem, i.e., the master problem with a limited number of variables/columns allowing the existence of at least one feasible solution of the linear relaxation of the problem, and a pricing problem or subproblem, generating new columns to add to the master problem in order to improve its objective value. When the pricing problem becomes infeasible, then by the separation/optimization theorem we know that the optimal value of the relaxed master problem has been reached.

We use column generation to optimally solve the linear relaxation of the master problem with a restricted set of initial variables. Indeed, we know that there exists an exponential number of CS as indicated in Section 5. At the beginning of our resolution process, we only generate one CS S 0 containing all the possible 2D-locations for UAVs, and placing them at the highest altitude available at these locations. This set S 0 is actually the most covering set, ensuring an initial coverage of the sensors at each time. S 0 is also guaranteed to be connected since there is one UAV placed at each possible position. However, this solution is not efficient in terms of objective value corresponding to the deployment cost. Placing a UAV at every possible location is very costly, and assigning them the highest possible altitude degrades the air-to-ground communication quality. Consequently, the column generation process iteratively solves the master problem and the pricing problem for every t ∈ [0, T ], generating new CSs S of UAVs with lower associated cost C S , optimizing the MWCS objective involving the dual values obtained from the master problem resolution.

After each resolution of the master program, we get the dual costs associated with columns z t S and solve the pricing program t times to find new improving CSs. When a new CS S has been found by the pricing program for a specific time, we choose to generate t new variables z t S to add to the master program in order to be able to reuse the CS for other time slots than the one for which it has been found. This method allows us to accelerate the resolution of ADCP CG OPT as presented in the next section.

At the end of the column generation process, we thus obtain the optimal solution of the linear relaxation of the master program, giving a lower bound for our problem. We finally run the master program again, with the set of CSs found by the column generation process, but with binary variables z and χ. This gives us an effective integer solution of ADCP that we evaluate in the following.

Results for ADCP CG OPT

We use here the same scenario as presented in Section 3.4.

Performance of the model

Table 4 summarizes the performance of ADCP CG OPT. We present results for the tested topologies with different numbers of mobile sensors |N t | and different number of 3D-locations for the UAVs |P|, listed in columns and . Note here that, |P| = 147 means (7 × 7) 2D-locations ×3 altitudes.

Column presents the number of generated connected subsets of positions during the column generation process. We remark that this number is very low and always below 10 CSs. This means that the number of generated columns needed to reach the optimal value of the relaxed master problem is small compared to the exponential size of the set of possible CSs S. In our model, when we run the pricing program for one particular value of t, the new subset found involves adding t columns in the master program, one for each time slot t. Consequently, a new subset of UAVs found for a specific time t can be reused for other time slots in the master program. This greatly improves the column generation process. However, the pricing program is ran for each time slot t in order to determine if there still exists a new subset to consider in the master program. This condition degrades the total resolution time as we will see in the next section.

Column gives the optimal value of the linear relaxation of the master program (z * LP ) obtained at the end of the column generation process. Column gives the value of the integer master program with the set of variables of the last iteration of the column generation (z ILP ). Column provides the accuracy of the integer solution which is computed as = 1-(z * LP /z ILP ). can be seen as the worst-case distance between z ILP and the optimal integer value of ADCP. A small value of thus means that z ILP is very close to the optimal integer value. When = 0, then the found integer solution is optimal for our problem. This optimum is reached for different numbers of sensors (5 to 20, and 35) and for different numbers of possible positions (75 to 300). We can see from column of Table 4 that our solution is always at most at 2% away from the optimal value, validating our approach.

Quality of the solution computed with ADCP CG OPT

We investigate the number of deployed UAVs for different numbers of mobile sensors. It is important to note here that the number of deployed UAVs depends on the evolving position of the sensors on the ground and thus, depends on time. We can see from Figure 3 that the number of deployed UAVs increases with the number of sensors to cover, for all numbers of 3Dpositions |P|. We present the mean number of deployed UAVs over time by a bar inside each boxplot in Figure 3. We can remark that the number of UAVs is less than or equal to the number of sensors. Indeed, a UAV can cover more than one sensor, optimizing the deployment cost in an effective way.

Figure 3 also depicts that the number of UAVs increases or is reduced on average by 1 unit over time (except for |N t | = 35 and |P| = 108). This information is given by the size of the boxes in the boxplot. For some scenarios, the number of UAVs remains constant over time, especially for |P| = 147, or when the number of sensors increases (|N t | = 45 or 50 for |P| = 75). This behavior is interesting for real deployment because it limits the use of new UAVs from the base station and enforces re-deployment of already placed UAVs.

In Figure 4, we depict the evolution of the cost of the selected sets over time. We seek to keep the cost as constant as possible due to the second part of the objective function of our optimization problem. However, the sensors' mobility enforces the UAVs to often change positions. In particular, when the number of sensors increases, then a configuration when several sensors are located close to each other can drastically change the cost of the subsets. For instance, the case for |N t | = 40 and |P| = 108 has a drop of the cost at time t = 1. This is due to the removal of a UAV located far away from the base station. The UAV is then needed again, making the cost become higher again. Moreover, the evolution of cost is optimized only for consecutive time slots. Due to the mobility pattern of the sensors, small changes are needed at each time, either by adjusting the altitude of the UAVs, or by changing the location of the most remote UAV.

Resolution time

Table 4 provides raw value of ADCP CG OPT computation times in column . We report in Figure 5 the evolution of the computation time depending on both the number of sensors and the number of possible 3D-positions. We remark that it takes at most 1303.261s (∼ 20 minutes) to compute solutions for instances with less than 147 possible 3D-positions for the UAVs, combined with any number of sensors. From Figure 5, we can see that the resolution time increases exponentially with the number of positions beyond 147 possible 3D-positions. When the number of sensors increases, computational time also increases to several hours of resolution, for example with 15 sensors and 192 possible 3D-positions for which the resolution time is 52511s (∼ 15h). This is due to the complexity of the pricing program which is NP-hard since it is related to a Steiner tree problem. Moreover, we have here the additional constraint that we do not know a priori the set of nodes to include in the tree. It is worth noting that when the number of time slots is increased, the time for the column generation process also increases since we need more pricing resolutions.

Compared to ADCP MILP, ADCP CG OPT shows better computational performance and provides solution quality close to the optimal. However, when both the number of 3D-positions and sensors increases, the computation time becomes unrealistic for a reactive system as the one we describe. It is important to find a trade-off between the computation time and the solution quality. This trade-off is studied in the following sections.

Heuristic pricing policy

In order to improve the pricing program resolution and thus, reduce the computation time, we change our pricing program policy and propose a fast heuristic algorithm instead of the optimal pricing computation presented in Section 5.2 for MWCS optimization. The algorithm is presented in the next section, and we analyze the results found on our tested scenarios in Section 6.2.

Heuristic MWCS algorithm

Algorithm 1 Heuristic MWCS We use Dijkstra shortest paths algorithm (ADCP CG DIJ) to find a new CS (see Algorithm 1) as the pricing program in our column generation process. We iteratively compute the weighted shortest paths from the base station to the targets in the graph G t for time t, in which the weight of link (p, q) ∈ E t corresponds to the dual cost associated with Definition 5. We run |N t | times the Dijsktra algorithm. At each step, the longest path found is selected and we record the positions of the UAVs 4 and describe here the results obtained for the column generation model with the heuristic pricing policy (ADCP CG DIJ). First, we can remark that the number of generated subsets during the column generation process (Column z of Table 5) is very low and always below 10 as has already been noticed for ADCP CG OPT. In columns {, |, and }, as for ADCP CG OPT, we recall the optimal value of the linear relaxation of the master program (z * LP ) and the value of the integer master program with the set of variables of the last iteration of the column generation (z ILP ) together with the accuracy . Solutions for ADCP CG DIJ are less accurate than ADCP CG OPT since their accuracy has values between 8 and 26% away from optimal (see column }), however it generates less subsets for small topologies to obtain very efficient solutions in a fast resolution scheme.

Require: Time t, G t = (V t , E t ), dual costs β (i) G = G t , N = N t , S = ∅ for all e = (p, q) ∈ E(G) do if q ∈ N t then cost e = 0 else cost e = D bq • (1 -γ) - n∈N t β (11) nt |{(q, n)} ∩ E(G)

Resolution time

The time needed to compute the solution for ADCP CG DIJ is below ∼ 1.6s for all instances with any combination of different numbers of 3D-positions |P| and sensors |N t |. The raw data for the computation time of ADCP CG DIJ are given in Table 5 column ~. This validates our policy to improve the computational time of ADCP.

Quality of the solution computed with ADCP CG DIJ

In Figure 6, we investigate the number of deployed UAVs over time for different numbers of mobile targets. We present the mean value over time by a bar inside each boxplot. We can see that the mean size of selected subsets is similar to ADCP CG OPT (see Figure 3a), and this mean value can also be lower than the one obtained with the optimal pricing policy, for example when |P| = 75 and |N t | = 20. However, for a given topology, the evolution of the size during the time period changes more than ADCP CG OPT. For each topology, the variation of the CS size stays usually bounded by 3. And for a large number of sensors, we see that ADCP CG DIJ generates bigger CSs than ADCP CG OPT (12 UAVs needed for the heuristic while the optimal pricing uses only 10 UAVs for |N t | = 50).

We investigate more specifically the solutions of ADCP CG DIJ in Figure 7. They globally use more UAVs than ADCP CG OPT during the time period (see Figure 7a), and the placed UAVs are located at higher altitudes due to Algorithm 1. However, the selected sets of ADCP CG DIJ maintain the positions more efficiently than ADCP CG OPT (see Figure 7b). This is an advantage for maintaining data collection and maximizing the aerial network lifetime.

Discussion and future works

In this section, we discuss the limitations and possible extensions of our work. Especially, we discuss the assumptions made in our model and the possible usage of the results.

• Application: This work will be used to solve, design and benchmark a real application that will be developed for wild animal tracking designed by biologists for marine turtle tracking. The turtles are tracked by a fleet of flying drones when they are in a bay area close to the base station. Related to the assumptions given by the biologists, the turtles are equipped with sensors (temperature, cardiac frequency, etc.), a GPS, and two communication interfaces. One interface (satellite) is used to send the position of the turtle to the base station when they move to the surface in the bay. This position is used as a trigger to deploy the fleet of drones. The second interface (not decided so far) will be used to transmit data from the turtles to the drones.

• Computation time: Figure 8 compares the resolution time of the three models. In Figure 8a the number of 3Dpositions is fixed with |P| = 75. In Figure 8b the number of sensors is fixed with |N t | = 5. These figures show the scalability of ADCP CG DIJ compared to the other methods. However, it is important to notice that choosing a method such as ADCP CG DIJ reduces the accuracy of the solution compared to the optimal pricing policy.

• Uplink and Downlink transmissions: In this paper, we focus on uplink transmissions from sensors to sink. However, the model can be used for data collection but also software updates. The permanent connectivity provided by the deployment of the UAV network allows such a behavior.

• Connectivity: In this paper, we focus on permanent connectivity of the network composed by UAVs. This constraint is strong and could be relaxed for applications where data from sensors are not required to be sent directly to the base station. An intermittent connectivity could be sufficient. For example, network connectivity should be established at regular periods and for a given duration. This can be solved by modifying our model, by relaxing the connected set problem at specific time slots. However, in our case, for the purpose of observation, data have to be sent in real time since different data are produced when the turtles are in the bay.

• Scalability: Depending on the size of the problem, the number of sensors, the size of the area to be observed, and the number of drones, the solution computed by the model may not be fast enough to be directly used to drive the drones. We are already working on some heuristics and some distributed models that could scale. In this case, the model described here will be used as a reference. However, based on the application description, the model could be used as is.

• Path planning: In this paper, we do not optimize the path of each drone but only provide successive positions that should be occupied. In future work, we will include an optimal path computation for each drone. Moreover, we will also include drone properties such as angle of curvature, etc.

• Communication protocol: The next step of this paper is to design or adapt communication protocols above the UAV network or to use existing protocols and modify the topology of the UAV network to obtain optimal performance from the network protocols. This second case is interesting and is left to future work since protocol properties can be injected into the model.

• Practical considerations: It is important to link the results from this paper to practical implementation. In this work, we try to be generic enough to be able to easily modify the model to fit any technology. We think that the model presented in this paper could be associated with communication technology such as WiFi for the air-to-air communication, LoRa for the ground-to-air communication and that off-the-shelves drone technologies fit the as-sumptions of our model. However, further investigations should be made regarding this subject.

Conclusion

In this paper, we address the aerial data collection problem from a set of mobile targets using a fleet of UAVs. We present linear programming formulations: a MILP and a column generation approach with two different pricing policies. We show that ADCP MILP does not scale for large instances, while ADCP CG OPT solves ADCP almost exactly with an NP-hard optimal pricing problem that becomes hard to solve with larger possible 3D-positions. We thus propose a heuristic pricing program generating larger subsets of UAVs but solving ADCP in less than 2 seconds with a small number of columns. Taking into account the work of the literature, this is the first model that scales with an increasing number of 3D-positions and targets, while considering mobility and connectivity.

Bibliography

Figure 1 :

 1 Figure 1: Multi-tier network architecture for the aerial data collection problem (ADCP).

Figure 2 :

 2 Figure 2: ADCP objectives.

  1 http://jgrapht.org/ Properties Number of Sensors |N t | [5;50] Sensors' mobilities random Terrain size 100 m×100 m 2D-positions (grid) [5×5] -[10×10] Possible altitudes {10 m, 25 m, 45 m} Base station coordinates (0,0) -alt. 0m Observation angle θ 60 degree Air-to-air communication range R u 30 m Number of time slots 10Table 2: Summary of instance generation.

  | s = b : source node (at the beginning it is the base station) while N ∅ do L = {L = (s, l 1 , ..., l k , n) shortest paths from s to n ∈ N} Select L * s.t. l * 1 = max L∈L h l 1 and cost L * = max L l * 1 e∈L cost e if length(L * ) = 1 then N = N \ {n} : path to n has been found else if e∈L * cost e = 0 then S = S ∪ {l * 1 , ..., l * k } and N = N \ {n} else S = S ∪ {l * 1 , ..., l * k } and N = N \ {n * = l * k } for e ∈ L * do cost e = 0 : encourage using same intermediate nodes s = l * 1 : new source for shortest paths return new CS S ⊂ P covering N t and connected with b

  (a) |P| = 75. (b) |P| = 108. (c) |P| = 147.

Figure 3 :

 3 Figure 3: Distribution of the size of selected subsets of UAV positions over time.

Figure 4 :

 4 Figure 4: Evolution of deployment cost over time.

Figure 5 :

 5 Figure 5: Resolution time in function of the number of 3D-positions for ADCP CG OPT. (y axis is in log scale)

6. 2 .

 2 Results for ADCP CG DIJ 6.2.1. Performance of the model

Figure 6 :

 6 Figure 6: Distribution of the size of selected CSs over t for |P| = 75 with ADCP CG DIJ.

  (a) Sum of |S |. (b) Common positions.

Figure 7 :

 7 Figure 7: Evaluation of Dijkstra pricing for |P| = 75.

5 Figure 8

 58 Figure 8

Table 1 :

 1 Contributions of our paper and related SOA.

Table 3 :

 3 Resolution time in seconds of ADCP MILP.

	P |N t |	P	P P P |P|	75	108
	5			4575.593	129764.829
	10			12543.451

Table 4 :

 4 Computational results of ADCP CG OPT.

	|N t |	|P|	# generated	z * LP	z ILP		Time (s)
			CSs
	5	75	5	3424.389 3424.389	0	25.668
	10	-	6	3905.826 3927.719 5.6 × 10 -3	31.473
	15	-	8	5657.524 5707.063 8.7 × 10 -3	136.312
	20	-	9	6096.653 6146.192 8.1 × 10 -3	249.439
	25	-	10	6468.669 6484.211 2.4 × 10 -3	233.053
	30	-	10	6574.307 6598.206 3.6 × 10 -3	276.054
	35	-	9	7081.05	7081.05	0	537.511
	40	-	8	7184.666 7203.332 2.6 × 10 -3	371.639
	45	-	10	7501.873 7536.078 4.6 × 10 -3	517.278
	50	-	10	7705.093	7718.27	1.7 × 10 -3	650.059
	5	108	6	3227.604 3235.296 2.4 × 10 -3	144.652
	10	-	9	3522.636 3545.914 6.6 × 10 -3	120.476
	15	-	8	4958.595 5046.899 1.8 × 10 -2	935.043
	20	-	7	5183.688 5183.688	0	888.504
	25	-	8	5312.243 5330.563 3.4 × 10 -3	745.25
	30	-	9	5345.651 5403.246 1.1 × 10 -2	769.817
	35	-	8	5885.56	5917.416 5.4 × 10 -3	1315.132
	40	-	10	6106.151 6188.343 1.3 × 10 -2	1150.683
	45	-	10	6360.619 6376.655 2.5 × 10 -3	911.651
	50	-	10	6559.157 6588.797 4.5 × 10 -3	1303.261
	5	147	6	2889.468 2889.468	0	365.313
	10	-	7	3128.349 3128.349	0	351.032
	15	-	8	4606.988 4606.988	0	3553.59
	20	-	8	4769.846 4774.719 1.0 × 10 -3	6670.563
	25	-	9	5031.934 5047.334 3.1 × 10 -3	6415.417
	30	-	9	5088.588	5092.94	8.6 × 10 -4	6440.275
	35	-	8	5472.193 5517.763 8.3 × 10 -3	8352.025
	40	-	9	5638.079 5674.783 6.5 × 10 -3	7682.315
	45	-	10	5899.808 5922.879 3.9 × 10 -3 11361.689
	50	-	9	5973.393 5996.662 3.9 × 10 -3	8723.237
	5	192	6	2782.009 2782.009	0	2000.798
	10	-	8	3174.343 3184.558 3.2 × 10 -3	1938.628
	15	-	10	4774.253 4789.456 3.2 × 10 -3 52511.092
	20	-	10	4879.655 4906.987 5.6 × 10 -3 69763.808
	25	-	10	5077.339 5095.287 3.5 × 10 -3 68535.967
	5	243	6	2652.332 2690.736 1.4 × 10 -2 11764.966
	10	-	8	3001.447 3010.627 3.0 × 10 -3	9415.008
	5	300	8	2576.889 2576.889	0	32964.699
	10	-	8	2761.192 2776.259 5.4 × 10 -3 15286.829

Table 5 :

 5 Computational results of ADCP CG DIJ.

	x	y	z	{	|	}	|N
	t | |P|	# generated	z * LP	z ILP		Time (s)
			sets
	5	75	5	3424.389 3909.936 1.2 × 10 -1	0.204
	10		5	3905.826 4510.616 1.3 × 10 -1	0.159
	15		6	5657.524 6339.536 1.0 × 10 -1	0.122
	20		6	6096.653 6415.276 4.9 × 10 -2	0.186
	25		6	6468.669 6763.646 4.3 × 10 -2	0.255
	30		5	6574.307 6866.951 4.2 × 10 -2	0.258
	35		6	7081.05 7434.808 4.7 × 10 -2	0.292
	40		8	7184.666 8000.568 1.0 × 10 -1	0.313
	45		8	7501.873 8413.125 1.0 × 10 -1	0.372
	50		9	7705.093 8536.661 9.7 × 10 -2	0.431
	5	108	5	3227.604 3878.835 1.6 × 10 -1	0.215
	10		8	3522.636 4463.87 2.1 × 10 -1	0.246
	15		10	4958.595 5970.774 1.6 × 10 -1	0.361
	20		10	5183.688 5799.398 1.0 × 10 -1	0.407
	25		8	5312.243 6240.422 1.4 × 10 -1	0.485
	30		7	5345.651 6368.844 1.6 × 10 -1	0.528
	35		9	5885.56 7606.585 2.2 × 10 -1	0.539
	40		10	6106.151 7608.826 1.9 × 10 -1	0.563
	45		9	6360.619 8085.365 2.1 × 10 -1	0.641
	50		10	6559.157 8518.769 2.3 × 10 -1	0.715
	5	147	8	2889.468 3205.161 9.8 × 10 -2	0.280
	10		7	3128.349 3369.58 7.1 × 10 -2	0.450
	15		9	4606.988 5553.485 1.7 × 10 -1	0.510
	20		8	4769.846 5526.248 1.3 × 10 -1	0.588
	25		8	5031.934 5756.901 1.2 × 10 -1	0.919
	30		9	5088.588 6080.205 1.6 × 10 -1	0.868
	35		10	5472.193 6344.368 1.3 × 10 -1	1.135
	40		9	5638.079 6320.016 1.0 × 10 -1	1.144
	45		10	5899.808 6743.387 1.2 × 10 -1	1.147
	50		9	5973.393 6907.604 1.3 × 10 -1	1.218
	5	192	8	2782.009 3359.266 1.7 × 10 -1	0.799
	10		10	3174.343 3979.738 2.0 × 10 -1	0.800
	15		9	4774.253 5719.449 1.6 × 10 -1	1.103
	20		9	4879.655 5815.745 1.6 × 10 -1	1.244
	25		9	5077.339 6416.000 2.0 × 10 -1	1.206
	5	243	9	2652.332 3371

.896 2.1 × 10 -1 1.396 10 10 3001.447 3584.85 1.6 × 10 -1 1.211 5 300 10 2576.889 3059.354 1.5 × 10 -1 1.412 10 10 2761.192 3663.003 2.4 × 10 -1 1.542