Thierry Pécot

Jérôme Boulanger

Charles Kervrann

Patrick Bouthemy

Jean Salamero

Estimation of the flow of particles within a partition of the image domain in fluorescence video-microscopy

Automatic analysis of the dynamic content in fluorescence video-microscopy is crucial for understanding molecular mechanisms involved in cell functions. In this paper, we propose an original approach for analyzing particle trafficking in these sequences. Instead of individually tracking every particle, we estimate the particle flows between predefined regions. This approach allows us to process image sequences with a high number of particles and a low frame rate. We investigate several ways to estimate the particle flow at the cellular level and evaluate their performance in synthetic and real image sequences.

I. INTRODUCTION

The molecular mechanisms that coordinate transmembrane protein trafficking at the cellular level are still unclear. Fluorescence microscopy allows us to localize the distribution of proteins of interest along time. The analysis of the dynamic behavior of these proteins enables to better understand the underlying mechanisms. A common and established approach for extracting this type of information consists in individually tracking particles in fluorescence video-microscopy. The most widespread tracking concept is the correspondence approach [START_REF] Cheezum | Quantitative comparison of algorithms for tracking single fluorescent particles[END_REF], [START_REF] Sbalzarini | Feature point tracking and trajectory analysis for video imaging in cell biology[END_REF] which consists in detecting particles independently in each frame and then associating the detected objects over time. The object association is particularly difficult when the density of particles is high, their appearance is similar and their trajectories interact. Kalman filtering [START_REF] Genovesio | Multiple particle tracking in 3-d+t microscopy: method and application to the tracking of endocytosed quantum dots[END_REF], particle filtering techniques [START_REF] Smal | Multiple object tracking in molecular bioimaging by Rao-Blackwellized marginal particle filtering[END_REF], graph-theory based methods [START_REF] Jaqaman | Robust single-particle tracking in live-cell time-lapse sequences[END_REF] and multiple hypothesis tracking [START_REF] Chenouard | Multiple hypothesis tracking in microscopy images[END_REF] have been developed to improve temporal matching. Practically, these approaches require a high temporal sampling rate to successfully track particles in fluorescence video-microscopy. Unfortunately, the phototoxicity of the acquisition process with this modality limits the number of images that can be acquired. Consequently, a high temporal sampling rate implies a short acquisition time. Another line of work in computer vision considers the object flows between predefined regions in the observed scene to recover the individual object tracks along time [START_REF] Berclaz | Multiple objects tracking using K-shortest paths optimization[END_REF]. Unfortunately, this method cannot be directly applied to fluorescence microscopy images because of the higher number of objects to track and their possible appearances and disappearances along time while only few occlusions happen in natural scene videos.

The particle flow between cell compartments reveals the relationships between the different cell domains, a generic problem in the study of the cell communication (endocytosis, recycling, motility, etc). Another possible application is the determination of the number of objects interacting with a region of interest over a given period allowing the evaluation of the trigger effect. This approach can also be applied to object flow estimation between non connected regions following the so-called network tomography as described in [START_REF] Pécot | Minimal paths and probabilistic models for origin-destination traffic estimation in live cell imaging[END_REF].

In this paper, we describe an original approach to extract global information about trafficking by analyzing the particle flow between connected regions. These regions are arbitrarily defined to form a regular cell tessellation or correspond to a specific cell segmentation such as cell compartments. The flow will be inferred from the variation of the number of particles in each region. In addition to the flow computed between regions in the cell, the number of appearances and disappearances is also estimated by introducing a virtual region. The particle flow estimation proceeds in three steps: i) the cell is partitioned into regions; ii) the particles are detected in the image sequence; iii) the particle flow is estimated at each time step from the particle number variation in each region. Except some trivial cases, flow estimation is an ill-posed problem. Therefore, we investigate two different regularization schemes [START_REF] Lawson | Solving Least Squares Problems[END_REF], [START_REF] Combettes | Fixed-Point Algorithms for Inverse Problems in Science and Engineering[END_REF] and compare them to a simple correspondence approach. The particle flow estimation provides less information than a tracking framework as it only describes the global particle paths in the cell compared to individual tracks. However, it offers an efficient way to extract relevant information from the image sequences while further keeping cells alive by reducing the frame rate, which leads to observe biological processes over a longer period of time.

The remainder of the paper is organized as follows. In Section II, we describe our particle flow estimation approach. In particular, we propose two regularizing methods to solve the flow estimation that is an ill-posed problem. In Section III, we compare the performance between flow estimation methods and tracking methods on synthetic and real fluorescence image sequences. Section IV contains concluding remarks.

II. FLOW ESTIMATION

In this section, we introduce the partitioning of the cell and define the associated graph. Then, we present the particle flow estimation approach and we propose two regularizing methods to solve this ill-posed problem.

A. Image partition

Let us consider a partition of the image domain Ω. An oriented graph G(E, V ) with |E| edges and |V | vertices is associated to the partition. The vertices correspond to the regions of the partition and the edges to the boundaries between the regions. A virtual region is added and connected to all the other regions to handle the particles that appear or disappear over time. The graph is oriented since only positive flows between regions are considered. An example of an image partition and the associated graph are shown in Fig. 1. In Fig. 1(b), the virtual region is represented as a red disk while the other regions appear as blue disks.

Several partition strategies are possible. A very simple partition consists in dividing the image domain into a regular grid of size n × m. In this case, the graph is defined with |V | = nm + 1 vertices and |E| = 6nm -2n -2m edges. If the particle trafficking is known to happen between cell compartments, the segmentation of these compartments supplies an adequate image partition.

B. Flow estimation

Matching particles is a hard task in fluorescence microscopy, especially when the number of particles is high, appearances and disappearances are frequent and when the frame rate is low. In contrast, we propose to infer the particle flow between regions by evaluating the number of particles in each region at each time step. Let us define

Y (t) = {y i (t), i ∈ [1, |V |]} the number of detected particles in each region i ∈ [1, |V |] at each time step t.
The following conservation equation (Kirchhoff law) can be written:

∂Y ∂t -M X = 0. (1)
This equation allows us to associate the particle flow

X(t) = {x j (t), j ∈ [1, |E|]} on each edge j ∈ [1, |E|]
with the temporal variation of the particle number in each region. In this context, the divergence operator M corresponds to a |V | × |E| matrix such that for each vertex i and edge j:

m ij =   
1 if edge j points to vertex i, -1 if edge j originates from vertex i, 0 otherwise.

(2)

As the number of edges is generally higher than the number of vertices, the matrix M has more columns than lines. The estimation of the flow X(t) from the measures Z = ∂Y (t)/∂t is performed through the inversion of M or the minimization of Z -M X 2 . The number of particles that appear or disappear is estimated by introducing the virtual region i * and by computing the total variation of particles between two time steps to be attached to that virtual region:

z i * (t) = -   i∈V \i * y i (t + 1) - i∈V \i * y i (t)   . (3)
As the graph is oriented, the flow on the edges is positive and this constraint can be exploited to constrain the solution. Then, we resort to a least-mean square algorithm with a positivity constraint to estimate the flow if the number of edges is inferior to the number of vertices. As the problem is ill-posed in most cases, we need to regularize the problem. First, we investigate a spectral cut method. More precisely, we estimate the particle flow with the technique proposed by Lawson and Hanson [START_REF] Lawson | Solving Least Squares Problems[END_REF], the so-called NNLS algorithm.

The probability that a particle crosses the boundary between two regions is relatively low. Then, we can assume that the flow estimates must be sparse, which leads to another regularization method. As a consequence, we minimize the following global function:

J(X, Z) = 1 2 Z -M X 2 + κ|X| + 1 2 d 2 C (X), (4
) where d C is the distance to the space C of the positive solutions and κ is a real positive parameter. To estimate the flow, we adopt the PPXA algorithm [START_REF] Combettes | Fixed-Point Algorithms for Inverse Problems in Science and Engineering[END_REF] well suited to non differentiable convex minimization problem. It uses the proximal operators associated to the different energy terms involving the proposed constraints (see Algorithm 1). In particular, the operator associated to the square distance function is given by [10]: 1 2 (x + P C (x)), where P C (x) = x if x > 0 and P C (x) = 0 otherwise.

III. EXPERIMENTAL RESULTS

In this section, we compare on synthetic and real data the performance of the two object flow estimation methods described in the previous section along with a nearestneighbor matching (NNM) algorithm.

A. Synthetic data

Using synthetic data, we aim at evaluating the influence of the number of particles involved, the particle speed, the appearance and disappearance rate, along with the number of regions of the partition. In addition to the two particle flow estimation methods described in Section II, we also consider Data:

M, Z, ε ∈ [0, 1], γ > 0, κ > 0 Result: X N U 1,1 = 0, U 2,1 = 0 for n ∈ [1, N ] do P 1,n = (I + γM T M ) -1 (U 1,n + γM T Z) P 2,n = 1 2 (U 2,n + P C (U 2,n )) P 3,n = max(U 3,n -1 κ , 0) sign(U 3,n ) P n = P 1,n + P 2,n + P 3,n ε ≤ λ n ≤ 2 -ε for i ∈ [1, 3] do U i,n+1 = U i,n + λ n (2P n -X n -P i,n ) end X n+1 = X n + λ n (P n -X n )
end Algorithm 1: PPXA minimization algorithm used to estimate the particle flow with proximal operators associated to the different constraints. a simple Nearest-Neighbor Matching (NNM) algorithm. The latter consists in pairing each particle at time t to the nearest particle at time t + 1. The particles that cross a boundary between two regions in the partition contribute to the flow. The particle flow estimated for each method is compared to the ground truth. For each simulation, the true positive probability (sensibility) P d is defined as the ratio between the number of correct detections and the total number of particles that cross a boundary between two regions over the sequence. The probability of false alarm (specificity) P f a is defined as the ratio between the number of false alarms and the number of boundary crosses.

Four different scenarios with 50 particles are simulated. All the scenarios correspond to a directed motion combined with a random walk. In the second scenario, the particle appear and disappear with a probability equal to 0.1. In the two last scenarios, on top of appearances and disappearances, the particles can stop at a given location with a probability equal to 0.8, and start over with a probability equal to 0.2. In the three first scenarios, the particle speed is equal to 3 pixels per frame while it is equal to 10 pixels per frame in the last scenario. The number of regions in each simulation varies in the range [START_REF] Smal | Multiple object tracking in molecular bioimaging by Rao-Blackwellized marginal particle filtering[END_REF]625]. For each setting, we report the average result over 10 simulations for each of the three methods. To supply a ground truth, we compute the crossing of the Evaluation of the influence of the number of objects for a directed motion combined with a random walk, with particle appearances/disappearances and start/stop. simulated trajectories with the boundaries of the regions. The time computation for the NNM algorithm increases as K 2 . The time computation for the NNLS and the PPXA algorithms does not depend on the number K of particles but on the number of edges in the graph.

As shown in Fig. 2, the NNM algorithm is almost insensitive to the image partition, except for the P f a value that increases slowly with the number of regions. But the P d value decreases when vesicles can appear/disappear and start/stop. This value also decreases a little when the speed increases. Except in the first scenario, the PPXA algorithm shows higher P d values and lower P f a values than the NNM algorithm. The particle appearances/disappearances and start/stop necessitate to consider a higher number of regions to obtain satisfying results. The NNLS algorithm shows a similar behaviour than the PPXA algorithm but the P d values (resp. P f a values) are lower (resp. higher). This demonstrates the importance of the sparsity constraint on the particle flows.

The influence of the number of particles on the results is evaluated and illustrated in Fig. 3 for the scenario shown in Fig. 2 (c). These results demonstrate the higher P d values and lower P f a values obtained with the PPXA algorithm when compared to the NNM and NNLS algorithms whatever the number of particles is.

B. Real TIRF image sequence

We compare the PPXA algorithm and the NNM algorithm on a real image sequence (Fig. 4(a)) of 237 frames of size 337 × 400 acquired in TIRF microscopy and showing the protein clip170 fluorescently labeled. This sequence is particularly challenging because of the high number of vesicles with varying speed. Furthermore, many vesicle appearances and disappearances happen as the acquired volume is thin in TIRF microscopy. The ground truth is defined by visually inspecting the image sequence. The results obtained with the NNM and PPXA algorithms for different sampling rates are reported in table I. We consider a simple partition of 5 regions (see Fig.

IV. CONCLUSION

We have proposed an original approach that determines the flow of particles between regions issued from a partition of the image domain. In contrast to usual approaches, we do not track individual particles, we only need to locally count particles on regions over time and minimize a global energy function. We have specified three methods to determine the particle flow. We have conducted comparative experiments on synthetic and real fluorescence image sequences. We have shown that adding a sparsity constraint on the number of detected events allows us to reduce the number of false alarms. Compared to usual tracking methods, our approach is simpler and the results are very stable with respect to the only two parameters. We have thus demonstrated the efficiency and the contributions of this approach.

Fig. 1 .

 1 Fig. 1. (a) Image partition example; (b) graph associated to the image partition (a).

Fig. 2 .

 2 Fig. 2. Performance evaluation of the NNLS, PPXA and NNM algorithm for four different scenarios. (a-d) Directed motion combined with a random walk, (b-d) object appearances/disappearances,(c-d) object start/stop. The particle speed is equal to 3 pixels per frame in the first three scenarios while it is equal to 10 pixels per frame in the last scenario.

  Fig. 3.Evaluation of the influence of the number of objects for a directed motion combined with a random walk, with particle appearances/disappearances and start/stop.

Fig. 4 .

 4 Fig. 4. (a) Image #1 of the sequence; (b) edge numbers between the regions.

4

 4

  (b)) to have a global idea about vesicle trafficking over the sequence. As the PPXA algorithm needs a high number of regions to obtain satisfying results, a partition of 730 regions is used and the final results are converted and reported in the partition with 5 regions. These results indicate that the PPXA algorithm better estimates the vesicle flows between the regions.

V. ACKNOWLEDGMENTS

The authors would like to thank Hélène de Forges and Franck Perez (Institut Curie, Paris) for providing the data set. They also acknowledge France-BioImaging infrastructure supported by the French National Research Agency (ANR-10-INBS-04-07, "Investments for the future").


Tan Viet

Anh Truong

Petr Benda
email: petr.benda@agents.fel.cvut.cz

Jiri Vokrinek
email: jiri.vokrinek@agents.fel.cvut.cz

The Task Agent Resource Function application in UAV domain

The goal of the paper is to present an application of TARF (Task Agent Resource Function) in UAV domain. This TARF is used to optimize the tasks allocation in case of cooperation between multiple Unmanned Aerial Vehicles. It plays the role of a generic mission planner for cross domain such as UAV, maritime, automotive and manned aerial vehicle (MAV). This paper concentrates on how to integrate TARF in Paparazzi platforman UAV platform based on PPRZ autopilot -in order to perform a swarm of aircrafts doing different tasks. The TARF contributes a new module to this platform for the simulation and real demonstrations.

I. Introduction

The Unmanned Aerial Vehicles (UAVs) have been developed more and more since the last decade. They are used in varied applications such as surveillance missions, gas pipelines monitoring, 1 localization and mapping, 2 meteorology 3 and photography. [START_REF] Hardin | An unmanned aerial vehicle for rangeland photography[END_REF] Several UAV open source autopilot projects developed recently [START_REF] Beard | Autonomous vehicle technologies for small fixed-wing UAVs[END_REF][START_REF] Tu | The Design of Small UAV Autopilot Hardware System Based on DSP[END_REF][START_REF] Lim | Build your own quadrotor: Open-source projects on unmanned aerial vehicles[END_REF] are enable to manage the navigation and monitoring before and during the flight. However, these autopilots are still far from controlling and optimizing the multi-tasks mission for an UAV swarm. This paper proposes a TARF solution for Paparazzi autopilot. [START_REF] Brisset | The paparazzi solution[END_REF][START_REF]Paparazzi, Open-source hardware and software project[END_REF] Paparazzi is an open source platform used to control and simulate offline/online one/multi UAV for different types of UAVs as well as tasks.

The task allocation issue has been studied longtime ago, and recently in UAV domain (see [START_REF] Hu | Dynamic group method based on health management applied to mission tasking for multi UAVs[END_REF] and [START_REF] Palmer | Decentralized cooperative auction for multiple agent task allocation using synchronized random number generators[END_REF] for examples). They focused not only on trajectory planning for agents [START_REF] Gancet | Task planning and control for a multi-UAV system: architecture and algorithms[END_REF] but also on how to distribute tasks for them. [START_REF] Bellingham | Multi-task allocation and path planning for cooperating UAVs[END_REF] The sequence of task is given at the desired time for one agent [START_REF] Chanthery | Planning algorithms for autonomous aerial vehicle[END_REF] and more complicated for multi agents. [START_REF] Gerkey | A formal analysis and taxonomy of task allocation in multi-robot systems[END_REF] Recently, the online task planning [START_REF] Pongpunwattana | Real-time planning for multiple autonomous vehicles in dynamic uncertain environments[END_REF][START_REF] Bertuccelli | Real-time multi-UAV task assignment in dynamic and uncertain environments[END_REF] is also conducted after dealing with the offline one. [START_REF] Grotli | Task assignment for cooperating UAVs under radio propagation path loss constraints[END_REF] The goal of this paper is to use the TARF as an online mission planner for multi-UAVs. This mission planner is then integrated in Paparazzi and contributing a new module to this platform. The Paparazzi community could perform different kinds of the cooperation of multiple UAVs afterward due to this module.

Among different problems of cooperation of multiple UAVs, the mission planning is merged as the main activity. The principal question is how to allocate available tasks to available aircrafts at one moment of time. Many algorithms have been proposed for mathematical and applied study. In this paper, we describe the mission planner TARF developed by ENAC (Ecole National Aviation Civil) and CTU (Czech Technical University). This mission planner is used for planning the cooperation between multiple UAVs. The UAV can be the fixed-wing aircraft or rotor-craft depending on the goal type. All the simulations were done by using the Paparazzi simulation platform. In the first simulation, we dealt with the scenario of three aircrafts doing three surveillance tasks.

The paper is organized as following. We first describe our TARF mission planner in section II. The demonstrator architecture and scenario are presented in section III and IV respectively. Finally, we show the simulation results in section V.

II. The TARF integration in Paparazzi

II.A. TARF description

The TARF algorithm (see in Figure 1) based on the multi-agents solver that minimizes the insertion cost computed by resource agents. [START_REF] Vokrinek | Abstract architecture for task-oriented multi-agent problem solving[END_REF] Figure 1. The TARF algorithm of task allocation, [START_REF] Vokrinek | Abstract architecture for task-oriented multi-agent problem solving[END_REF] where a and t are agent and task, respectively and CNP is contract net protocol. [START_REF] Smith | The contract net protocol: High-level communication and control in a distributed problem solver[END_REF]

II.B. Mission formulation and TARF API (Application Programming Interface)

We define the mission as the following data structures: Mission = {Waypoint, Task, Aicraft} Aircraft = {Aircraft Type, Aircraft Status, Speed} Task = {Task Type, Task Position, Task Status, Duration} Waypoint = {x, y, altitude} TARF connects to other modules via the API as in Figure 2. The input of TARF (TARFin) is all the information about the mission including agents (aircrafts), tasks and other resources. The output of TARF (TARFout) is the list of agents (aircrafts) and the corresponding (assigned) tasks of each aircraft.

II.C. TARF-Paparazzi google protocol buffer

Because of the different programming languages used in TARF (in java language) and Paparazzi at ground station level (in c language), a TARF-Paparazzi protocol buffer is proposed as in Figure 3. This protocol buffer developed by Google uses the server/client mechanism. The TARF client (in Paparazzi platform) sends information about the tasks and aircrafts to TARF server (Task allocator) and receives the allocation update including allocated tasks. The "mission planner" in Paparazzi monitors and sends all the feedbacks during

of 12

flight to TARF client. There is also the low level planning of collision avoidance in Paparazzi platform. When there is a conflict between two UAVs, they will change their altitudes to guarantee their desired trajectories. The Figure 4 shows the demonstration architecture of one UAV in the cooperation between UAVs. This architecture consists of the Paparazzi modules in dark green and the new module TARF in pink.

III. UAV Demonstration architecture

Each aircraft is equipped with an autopilot (AP) which includes the natigation, stabilization during flight as well as a Collision Avoidance module. This collision avoidance use the "TCAS" technique to detect and avoid the collision between UAVs (see on the Paparazzi website [START_REF]Paparazzi, Open-source hardware and software project[END_REF] for more details about this module).

Aircrafts connect to the Ground Control Station (GCS) by the LINK agent. All the telemetry data and uplink (control) data are stored and dispatched by Server via a system of particular messages (in module Messages). All these agents are interconnected by IVY bus, [START_REF]A simple protocol and a set of open-source (LGPL) libraries[END_REF] an open-source software bus that can plug modules from many different languages such as C, C++, java, python, etc.

The mission is defined by the operator in a XML file as an input of the TARF. TARF collects the data of tasks, agents and resources of the mission. Then it allocates these tasks to aircrafts by sending the commands to these aircrafts via IVY. Paparazzi ground segment monitors the flight and send the feedbacks to TARF.

IV. Scenario

The scenario (in Figure 5) shows a typical issue in case of multi-UAVs. When an aircraft in an UAV swarm is lacking of battery, all the tasks will be redistributed to the rest aircrafts. The mission is composed of the surveillance task in a three different zones. For this mission, we have three UAVs (of the same type -fixed wing aircraft for this simulation) and one operator. Basically, the main task is the surveillance of the large section with several UAVs. This task is decomposed to the surveillance tasks in three different sections each for only one UAV. The reconfiguration phase happens when an UAV is lacking of battery, one of the other UAVs will replace this one to finish its on-going task. The background tasks include then the separation phase -allocation of task and the deconflict phase -collision avoidance.

V. Results

The simulation is done by using Paparazzi platform and results are shown by the following video: https://www.dropbox.com/s/y257cezyajioca8/demotarf.mp4?n=141175626 Three surveillance sections are respectively defined by (200 × 350(m)) S1-S2, (220 × 390(m)) S2-S3 and (200 × 370(m)) S4-S5 rectangle. The flight zone in the simulation is taken from the Muret (a small region near Toulouse, France) where our real demonstration will be carried out afterward. The google earth map of this region is loaded and served in all the figures of simulation as the background of GCS. All three aircrafts in the simulation use the Funjet aerodynamic model and the typical fixed-wing aircraft airframe and setting parameters. They are launched from the HOME waypoint and reaching first their standby points before doing their tasks. The other waypoints such as CLIMB, AF, etc are used for the the take-off and landing procedure. The TARF planner is running as following steps:

• Launching three aircrafts: They are launched from the HOME point one by one (see Figure 6).

• Initializing TARF: These three aircrafts are sent to their three standby points STDBY1, STDBY2, STDBY3 after taking-off. Normally, they are circling around these points and waiting for the next commands. In the simulation, all these procedures are reduced and the initiation step begins when the TARF-Paparazzi interface is called (see Figure 7). In this stage, TARF read all the data from mission XML file and parse them to the particular messages of waypoints, tasks, aircrafts.

• Starting TARF: TARF allocates three surveillance tasks to the three aircrafts (see Figure 8). Basically, TARF distributes respectively tasks to available agents. At the beginning, three aircrafts are available and the first task (S1-S2 surveillance task) is allocated to its nearest aircraft -the aircraft AC D3CoS 3 in red. The second task (S2-S3 surveillance task) is then allocated to the aircraft AC D3CoS 2 in yellow and the third task (S4-S5 surveillance task) is allocated to the aircraft AC D3CoS 1 in blue, respectively. During the flight, Collision Avoidance module is always triggered each 500 (ms) in order to detect the conflict between UAVs (see Figure 9 for an example). The trajectories of the surveillance tasks is tracked as in the Figure 10. These trajectories are defined by dividing the zone as the parallel lines. Aircrafts follow all these lines to cover their surveillance zones. The radius of the circle used for returning aircrafts once they goes out of the zone is determined in the flight plan before flying and can not be changed during the flight.

• Updating TARF: Three aircrafts are continuing doing their current tasks. One agent of GCS is detecting any problem during flying including the battery, connection, etc and send to TARF. The "Kill AC" button in the TARF-Paparazzi interface is used to simulate a problem concerning the battery of the aircraft AC D3CoS 2 (see Figure 11). The yellow aircraft is then supposed to land and being considered as dead (Figure 12). The TARF receives the feedback of on-doing tasks, recalculates and sends the allocation update including the list of new tasks corresponding to aircrafts. In this case, the S2-S3 surveillance task is still not done and the two rest aircrafts (AC D3CoS 1 and AC D3CoS 3) are continuing to finish their own tasks (see Figure 13). This task is always waiting for one of two these aircrafts once it accomplishes its assigned task. Actually, the surveillance task is a cycle task. Aircrafts cover the zone as much as possible to obtain the goal which might be detecting targets or measuring the atmosphere inside the zone, etc. Our mission XML defines the duration for each surveillance task so that the aircraft could finish its task and go to another task if any. In our simulation, the durations of task 1 (S1-S2 surveillance task) and task 3 (S4-S5 surveillance task) are 400(s) and 700(s), respectively. As the result, TARF allocates the task 2 (S2-S3 surveillance task) to aircraft AC D3CoS 3's list of tasks. Once aircraft AC D3CoS 3 terminates task 1, it goes to task 2 and continue this task of AC D3CoS 2 at the beginning (see Figure 14). These two aircrafts are doing their tasks until finishing the tasks durations (see Figure 15).

VI. Conclusion

The paper presents the TARF planner module integrated in Paparazzi platform. This planner allows us to allocate tasks to aircrafts in the context of cooperation between multiple UAVs. The simulation of the proposed scenario gave the good results. Aircrafts could cooperate together to perform the desired tasks set. The scenario in the paper deals with three aircrafts with the same type (fixed-wing aircraft). The replacement/substitution of multiple UAVs when there are problems of battery, connection missing, etc could be resolved with this TARF planner. For the next step of this work, we will simulate with other scenarios for the different tasks and the different types of aircrafts as well (that could be the combination between aircrafts and rotorcrafts). We are going to apply this TARF module within the D3CoS (Designing Dynamic Distributed Cooperative Human-Machine Systems) 22 project to make a real flight demonstrator of an UAV swarm. The TARF could also be used in other domains in the project such as automotive, manned aerial vehicle and maritime.

Figure 2 .

 2 Figure 2. The API of TARF.

Figure 3 .

 3 Figure 3. The architecture of TARF-Paparazzi google protocol buffer.

Figure 4 .

 4 Figure 4. The simulation architecture based on Paparazzi platform.

Figure 5 .

 5 Figure 5. The scenario with three aircrafts, three tasks and one operator.

Figure 6 .

 6 Figure 6. Launching stage: three aircrafts in red, blue and yellow are at the HOME point at the beginning. All the waypoints are in the UTM local coordinate relative to the the HOME point. We can see in the upper left corner the GCS interface of Paparazzi platform with strips of button for controlling/monitoring/editing before and during flight. The warning messages system is found on the lower left corner. The blocks of tasks and the altitude monitoring interface are on the lower side of the figure.

Figure 7 .

 7 Figure 7. TARF initialization: three aircrafts are going their standby points STDBY1, STDBY2, STDBY3 and the TARF-Paparazzi interface is called on the left of the screen

Figure 9 .

 9 Figure 9. The conflict happens between the blue and yellow aircraft. As the result of using "TCAS" technique, the yellow aircraft is climbing and the blue one is descending for avoiding the collision between them. This phase of conflict is shown in the altitude interface on the lower side of the figure.

Figure 10 .

 10 Figure 10. The trajectory shapes of three aircrafts doing their surveillance tasks. The grids is pre-calculated in autopilot

Figure 11 .

 11 Figure 11. Kill aircraft simulation: at the beginning of the command, AC D3CoS 2 starts to land

Figure 13 .

 13 Figure 13. Aircrafts are continuing doing their surveillance tasks until finishing the durations

Figure 15 .

 15 Figure 15. Two aircrafts have almost finished their tasks and they are ready to land

of 12

Acknowledgments

This work is financed by the europeen project "D3CoS". We thanks Gautier Hattenberger for his valuable help.


Murat Bronz
email: murat.bronz@enac.fr

Jean-Philippe Condomines
email: jean-philippe.condomines@recherche.enac.fr

Gautier Hattenberger
email: gautier.hattenberger@enac.fr

Development of an 18cm Micro Air Vehicle : QUARK

This paper describes the development of the QUARK micro unmanned air vehicle. The main objective of the study is to show the feasibility of designing an autonomous MAV smaller than 20 cm by using open-source programs and off-the-shelf components. Effort is given to show how to choose the system components and the design variables correctly in order to end up with an optimum design. An open-source multi-disciplinary conceptual aircraft design program called CDSGN is used in order to select the appropriate wing planform and the optimum battery capacity. Propeller and electric motor are selected throughout an optimization procedure, which analyze the performance of each suitable off-the-shelf propeller and motor combination for the specific mission requirements of the new design. A 2 cm × 2 cm custom autopilot board is specifically designed according to the requirements of the QUARK. The modelization of the MEMS sensors (Micro-Electromechanical Systems), that are used on the autopilot, is described and analysed in order to improve the stabilisation and the navigation of the aircraft. Longitdinal dynamics of the vehicle are examined by using an innovative technique, and following that, corresponding flight control gains are tunned prior to the outdoor flight test.

The first prototype has been flown and a smaller 15.9 cm version is designed according to the obtained results. The expected endurance for the final version is over 35 minutes with reaching up to 35 km of range.

Introduction

The birth of Micro Air Vehicles goes back to early 1990s, when the RAND Corporation performed a feasibility study. Following that interest, in 1996 Aerovironment showed the feasibility of flying a six-inch MAV in 1996, and later improved the performance of the existing design and performed a 16 minutes of flight in 1997 by manual control. After this achievement, Aerovironment was awarded a DARPA founded Phase II Small Business Innovation Research (SBIR) contract. As a result of this contract, the Black Widow MAV was designed, which is still one of the smallest fixed-wing MAV. Black Widow MAV had 30 minutes of endurance and was able to transmit real-time video to the ground within 1.8 km range. There was an autopilot for stability augmentation, however there was no GPS on-board for autonomous navigation. Recent progresses in the miniaturization of the inertial instruments i.e. Inertial Navigation System (INS) and Global Positioning System (GPS) hardware have led to the design of small, low-cost integrated navigation systems. Despite all these improvements, in the year of 2013, there is no fixed-wing MAV with autonomous navigation capability, which remains a big shortfall. Figure 1 shows INS relative strap-down system cost "projections" as a function of inertial instrument technology and performance.

The performance of MEMS-based gyroscopes and accelerometers is significantly affected by complex error characteristics that are stochastics in nature. Therefore, the modelisation of the sensors remains a challenge in order to improve the navigation.

Design Methodology

The design method was kept simple. The idea was to gather the most appropriate off-the-shelf components that are required for the flight, control, and communication of the aircraft and decide the minimum core size while taking into account the interference and center of gravity issues. Then by using a multi-disciplinary parametric conceptual design program called CDSGN [START_REF] Bronz | A Contribution to the Design of Long-Endurance Mini Aerial Vehicles[END_REF], the wing planform, battery capacity, cruise speed, motor and propeller combination was optimized while the cost function being the endurance and range performance maximisation. Certain stability derivatives, such as pitch damping (C mq ) and spiral stability (C l b C nr /C lr C n b ), are used as constraints during the optimization phase.

Off-the-Shelf Component Search

One of the main objectives of this study was to use completely off-the-shelf commercially available components for the realization of the aircraft.

Mini-micro Scale Propellers. Among the available mini-micro scale propellers, seven different candidates were found and taken into consideration for the final propulsion system selection. Figure 2 shows the pictures and Table 1 presents the specifications of the candidate propellers. Each of the propellers modelled geometrically by measuring their chord length, and twist angle at several sections throhout the blade span. The corresponding airfoil geometry is extracted by cutting the propeller blade into pieces and then aerodynamic characteristics of each section airfoil is examined by using XFOIL [START_REF] Drela | An analysis and design system for low reynolds number airfoils[END_REF]. Following that, the performance modelisation of the propellers done by using QPROP [START_REF] Drela | QPROP Formulation[END_REF]. 4 shows the main parameters that were variable for the wing planform while optimizing performance of the QUARK MAV.

System Integration and Design

It is crucial to take into account the system integration during the design phase. Finally, the limitation on the battery shape, servo choice and GPS module (because of their availability and ordering issues by a research lab) restricted the minimum dimensions of the parameters and created a solid core system that has to be fitted into the vehicle. According to the center of gravity and the wing thickness issues, the placement of the components have been fixed without too much option. Figure 9 shows the internal placement of the components in the final design. Note that the requirement of placing every component as close as possible to the leading edge comes from the center of gravity alignment. Figure 5 shows the first prototype that is build from EPP foam by using a CNC hot-wire cutting machine and the custom designed smallest Paparazzi autopilot for the QUARK MAV.

Sensor modelling

All the sensors embedded in the Quark autopilot are low-cost and so have imperfection. The major error sources in the navigation system are due to : -the inertial instrument imperfections (the sensor itself); -the numerical integration process in the INS mechanization. Inertial sensors have both deterministic and random errors. Deterministic errors are due to manufacturing (characteristized by the manufactuers and generally over-estimated with respect to effective values [START_REF] De Agostino | Performance Comparison of Different MEMS-based IMUs[END_REF]) and mounting defects. This can be obtained using calibration procedure. On the other hand, the stochastic errors are random errors that occur due to random variations of bias or scale factor over time. The accurate modeling of the MEMS sensor errors is one of the most challenging tasks in the design of low-cost navigation. Indeed, to overcome the disadvantages associated with autonomous operation of the GPS and INS, the two systems are often paired together and the INS/GPS data fusion is commonly permormed using Kalman Filter (KF) [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. There are more recent derivatives of KF known as GMEKF [START_REF] Martin | Generalized multiplicative extended kalman filter for aided attitude and heading reference system[END_REF] and π-IUKF [START_REF] Condomines | Nonlinear state estimation using an invariant unscented kalman filter[END_REF]. KF has widely been used for data fusion but there are, however some inadequacies, including the following :

• the requirement for a priori information of the system and measurement covariance matrices, R and Q, for each sensor;

• the necessity to have a gaussian assumption on probability distribution and a white bandwith noise.

A general model for the sensor output can be described by the following equation [START_REF] Wenger | Systems Concepts and Performances Analysis of Multi-Sensor Navigation System for UAV Applications[END_REF] :

ζ m = (1 + S f )ζ t + b(t) + n(t) (1)
where ζ m is the mesured quality at the sensor output and ζ t is the true quatity. S f represent scale factor error (In this work, we will ignore the scale factor error S f = 0 i.e, calibrated out at the factory), b(t) represents the time varying bias (random drift) and n(t) is observation noise, which is assumed to be band limited white noise with spectral density R. This can be characterized by taking the standard deviation of the sensor's output over a short period of time with no input applied. The bias term b(t) is a random process consisting of the integral of a white noise process ω(t) with spectral density Q. The bias signal is modeled as follows : In order to extract the covariance and static bias of each sensor noise afecting the ouput, we have analyzed 1000 samples of each sensor using a PDF (Probability Density Function). These data

ḃ(t) = ω b , E{ω(t 1 )ω(t 2 )} = Qδ(t 1 -t 2 ) (2) Static Bias σ R = √ σ 2 • dt R (datasheet) Q X-

International Micro Air Vehicle Conference and Flight Competition (IMAV2013)

17-20 September 2013, Toulouse, France coming from hight frequency dataloger system working at 100 Hz (dt=0.01s) for INS sensor and 5 Hz (dt=0.2s) for GPS. Remark that, a simplest Allan variance analysis [START_REF] El-Sheimy | Analysis and modeling of inertial sensors using allan variance[END_REF] with a record of 2 hours of data has been conducted to obtain spectral density Q. The determined models for each of the sensors ( gyrometers and accelerometers) are tabulated in Table 3, and shows that spectral density values R (µg or • /s/ (Hz) ) given by the manufacturer are generally overestimated or underestimated with respect to effective values. For example, a spectral density value of 0.0082 Inertial sensor always contain a certain amount of measurement noise (e.g, vibration). On the sensor level, this noise is usually separated into long-term (low-frequency) and short-term (highfrequency) noise. Wavelet denoising is frequently used to process inertial sensors measurements to improve the signal-to-noise ratio (SNR). However, in the case of QUARK Paparazzi autopilot, it will not be possible to use wavelet denoising because it is computationally costly. Frequency Fourier Transformation (FFT) of both gyrometer and accelerometer, recorded with an on-board datalogging system, for main flight phases (waiting, takeoff, cruise an landing) are represented in Figure 7. FFT which depends on the throttle, shows similar characteristics on both sensors after 20 Hz which can be considered as vibration noise. Thus, a Low Pass Filter (LPF) with a cut frequency at 22 Hz can be used to remove high frequency vibration. This LPF has been implemented to the autopilot. The difficulty of flying such a small aircraft is obvious, so that prior to the outdoor flight tests, it is decided to tune the longitudinal control loop gains with a restricted flight method, that is shown in Figure 8, and tune the longitudinal control loop gains prior to the maiden flight. The experiment is done in a closed place, where a pole is used as the center fixation point. The aircraft is connected from its wing-tip to an encoder which is placed on top of the fixation pole. The encoder was turning freely and measuring the exact rotation rate during the experiment. Additionally, a high frame rate camera is placed on top of the encoder directed to the aircraft which will be recording the vertical movement of the aircraft during the flight with a fixed reference.

Support Pole

The length of the line has been increased gradually. The center of gravity was located at aircraft's aerodynamic center because of component integration issues resulting a highly agile vehicle. Therefore, the experiment started directly with autopilot assisted mode where the pilot controls the attitude of the vehicle through autopilot (fly-by-wire). The longitudinal control loop gains tuned while watching the oscillation of the vehicle and the line length has been increased up to the wall limits.

As a result of the experiment, several important flight parameters are obtained such as, the longitudinal control loop gains, the minimum flight speed, approximate cruise throttle setting, flight angle of attack at different speeds. Additionally, the control surface travels obtained for the pitch control have been copied for the roll control travels as a starting point for the outdoor flight.

Outdoor Flights

The first parameters and gains that are obtained from the pole flight tests are applied for the maiden flight. The flight has been started directly in autopilot assisted mode as explained in the previous section. After a short attempt, a small angle of attack increment was required in order to sustain the level flight. The pitch proportional gains has been increased (10%) in comparison to the values that are obtained from the pole flight tests. However the lateral behavior of the vehicle was satisfactory directly from the beginning (in no-wind condition). Although the maiden flight took less then 5 minutes, the amount of battery energy used for the flight gives an estimated 20 minutes of flight time for the current prototype with small batteries.

The second flight has been made in a windy day with more than 10 m/s wind gusts. The main objective of this flight was to tune the navigation control loop gains in order to fly completely autonomously. Unfortunately, 20 m after the launch the telemetry link started to degrade and as the safety pilot radio command link was transmitted through the telemetry, it was no more possible to take back the flight control of the vehicle. However, the starting navigation control loops were sufficient enough to keep the vehicle in the air following the predefined flight circle on top of experiment area. Therefore, we had the chance to change some of the stabilisation and navigation control gains whenever we catch telemetry connection.We have observed severe rolls because of the strong wind gusts several times, and the counter force applied in order to correct them was seemed to be not sufficient and fast enough. That leads us to increase the lateral gains which should theoretically improve the response for the wind-gusts, unfortunately before being able to change and tune the gains (lack of constant communication), the vehicle crashed because of a strong wind gust. Thanks to its light weight and low inertia, beside a broken propeller there was no other major damage.

The third flight has been made in calm weather conditions with some modifications to the telemetry link with the hopes of improved the communication. Unfortunately, the telemetry behavior was the same as the previous flight limiting the further tuning of the vehicle. As the weather was calm with less wind, the vehicle flew much better and follow the desired trajectory closely compared to the previous flight. After 5 minutes of autonomous flight, the plane started to loose altitude at a point of the trajectory where the communication was lost. It finally crashed in a field where it was not possible to retrieve it. Without any telemetry data, the reason for this crash remain unknown.

Despite having lost the only prototype, we have obtained precious information from the flights. The prototype showed the feasibility of flying smaller than 20 cm vehicle with cheap off-the-shelf

International Micro Air Vehicle Conference and Flight Competition (IMAV2013)

17-20 September 2013, Toulouse, France components. It was also a good opportunity to verify and correct our aerodynamic codes for this scale which gives a big confidence for the future designs.

On-going Work

According to the experience gained from the prototype design and flights, a new version of the QUARK MAV is designed with 15.9 cm wing span. This time the manufacturing of the vehicle is going to be done by using precisely machined molds by CNC. The total take-off weight is 67 g.

Figure 9 shows the CAD model of the new QUARK with the already manufactured aluminium molds. As a continuation and future work of this study, the telemetry integration will be taken into consideration more seriously in order to improve the communication. Additionally in order to make this scale of MAVs useful for real-life applications such as surveillance and information gathering from an hazardous place, a micro camera and its transmitter will be integrated into system. The camera integration has been already included on the final QUARK design.

Conclusion

Present work explains the efforts given for the realisation of an MAV with a maximum size of smaller than 20 cm. The first prototype (18 cm) is designed and manufactured by using off-the-shelf components. Complete autonomous flight including stabilisation and navigation has been achieved with the prototype. The results obtained from the limited number of flights shows a good confidence for the final design (15.9 cm) to fly over 35 minutes at a flight speed of 15 m/s and reach up to 30 km for the straight line range at a speed of 16 m/s.

Figure 1 :

 1 Figure 1: Strapdown INS cost as a function of instrument technology

Figure 2 :

 2 Figure 2: Pictures of candidtae commercially available off-the-shelf mini and micro propellers

Figure 4 :

 4 Figure 4: Main parameters for the wing planform that were variable during the multi-disciplinary optimization of the QUARK MAV by CDSGN program

Figure

Figure 5 :

 5 Figure 5: The prototype QUARK MAV wit the new developed QUARK Paparazzi Autopilot (October 2012)

  rotation rate (degrees/s) Mean : -1.55 deg/s RMS error : 0.082 deg/s Mean : -1.13 deg/s RMS error : 0.167 deg/s Mean : -1.70 deg/s RMS error : 0.221 deg/s

Figure 6 :

 6 Figure 6: Distribution and RMS error of the used gyrometer

Figure 7 :

 7 Figure 7: Frequency Fourier Transformation of roll-axis gyrometer and X-axis accelerometer values with the flight altitude information on the top.

Figure 8 :

 8 Figure 8: The longitudinal control test setup for restricted flight tests

Figure 9 :

 9 Figure 9: The CAD representation of the final QUARK MAV with the internal components and the aluminium CNC molds on the right

Table 1 :

 1 Coefficients of candidate commercially available off-the-shelf mini and micro electric motors

	Propeller Name	Blade number Diameter Pitch Mass
		[-]	[inch]	[inch]	[g]
	Three Blade Prop	3	3.8	-	1.9
	Graupner CAM Slim	2	3	3	0.8
	Graupner CAM Speed	2	4	3	0.95
	GEMFAN mini 3x2	2	3	2	1.0
	GEMFAN mini 4x2.5	2	4	2.5	1.4
	Red propeller	2	3	-	0.7
	HobbyKing nano prop	2	2.2	-	0.3

Table 2 :

 2 Coefficients of candidate commercially available off-the-shelf mini and micro electric motors

	International Micro Air Vehicle Conference and Flight Competition (IMAV2013)
	17-20 September 2013, Toulouse, France
	Centerline
	Generic Fuselage

Table 3 :

 3 Covariance and static error for INS sensors

	GPS	σ
	X	0.88 m
	Y	1.45 m
	Z	20.29 m
	V X	0.248 m/s
	V Y	0.146 m/s
	V Z	0.280 m/s

Table 4 :

 4 GPS variance

  • /s/ (Hz) = (0.0082 • 60) • / (hr) on the roll axis rotation rate indicates, that the angular error (incertainty) due to random walk is e.g. 0.063 deg after 1minute or 3.8 deg after 1hour. The first two column represents static bias and standard deviation (σ) and the last column Q represents spectral density of ω(t) coming from Allan variance analysis. In terms of visual comparison, only gyrometers histogram are shown in Fig.N, and match with a gaussian distribution. Note that it is the same shape for accelerometers histogram. Finally, GPS data variance are tabulated in Table4for speed and position.

Micro Electric Motors. Likewise a micro electric motor database was created from the commercially available off-the-shelf motors. Figure3shows the candidate electric motors. A first order electric model is used to characterize the behavior of the motors. Only the internal resistance (R m ), no-load current (I 0 ), and speed constant (K v ) are used for characterization. Table2presents the mass of each motor with their K v coefficient.4 Design ParametrizationOnce the required electronics was defined, the sizing of the aircraft has taken into consideration in a parametric way, which starts with the minimum possible dimension of the internal components and goes up to the user defined maximum limits. The optmization was done by using CDSGN program

Acknowledgements

The authors would like to thank to Matthieu Navarro, Michel Gorraz for their effort on developing the smallest Paparazzi Autopilot and Alexandre Bustico for helping on the sensor measurements.


Murat Bronz
email: murat.bronz@enac.fr

Gautier Hattenberger
email: gautier.hattenberger@enac.fr

Jean-Marc Moschetta
email: jean-marc.moschetta@isae.fr

Development of a Long Endurance Mini-UAV : ETERNITY

This study presents the effort given for the first prototype of a Long Endurance Mini UAV concept called Eternity. A multi-disciplinary conceptual aircraft design program called CDSGN is developed and used for the design of the Eternity. Unlike the traditional design methods that uses statistical data from the previous well-flown aircrafts, CDSGN analyses numerous aircraft candidates and simulates each candidate for the given mission definition and outputs the corresponding performance. The unique property of the presented design methodology comes from a computationally fast and physically accurate modelling of the aerodynamic characteristics of each candidate by using a modified version of a vortex lattice program called AVL from Mark Drela. Two types of configurations have been analysed for the Eternity design, conventional and flying-wing. A wide envelope of variable design parameters used for both configurations such as wing surface area, cruise speed, battery capacity, different airfoils, etc... Integration of solar cells, and the management of solar energy is also considered for every candidate. Only the wing span size has been limited to one meter. Additionally, the on-board avionics and payload weights and sizes are fixed for every candidate as they are independent of the design. Analyses by CDSGN concluded the dominance of the conventional configuration for the given long endurance mission performance both on solar and non-solar conditions. Optimum wing surface area and the on-board battery energy found interactively by a post-filtering program developed in-house. A custom airfoil family, transitioning along the span, have been designed specifically for the corresponding local Reynolds number for specific spanwise locations. A wind-tunnel campaign is performed with a full-scale model and first flight tests have been performed in order to show the feasibility of long endurance flights.

Introduction

Endurance performance enhancement of mini and micro UAVs has been a common interest. The endurance performance relies on many aspect of the design where it has to be concentrated in a more specific way than the traditional conceptual design methods such as Raymer [START_REF] Daniel | Aircraft Design : A Conceptual Approach[END_REF] and Roskam [START_REF] Roskam | Airplane Design[END_REF].

This study will present the design optimization of a one-meter UAV, called Eternity, which shows the feasibility of an electric powered Long-Endurance MAV concept. The particular interest is on the developed CDSGN conceptual aircraft design tool specialized on MAV scale. Each of the subsystems (e.g. aerodynamics,structure, propulsion,control,...) modelled seperately in order to see their final effect on the performance of the vehicle.

CDSGN : The conceptual design tool

CDSGN [START_REF] Bronz | A Contribution to the Design of Long-Endurance Mini Aerial Vehicles[END_REF] is a multi-disciplinary conceptual aircraft design program. The core aerodynamic calculations are done by using a modified version of AVL program [START_REF] Drela | Athena vortex lattice[END_REF], which uses vortex lattice method. Modifications include the airfoil viscous drag addition coming from an XFOIL [START_REF] Drela | An analysis and design system for low reynolds number airfoils[END_REF] database, fuselage skin friction drag with flat plate approximation [START_REF] Dr | Fluid-Dynamic Drag[END_REF], and stall information according to the maximum local lift coefficients of the particular airfoils at corresponding Reynolds number.

It uses a computationally fast, yet physically realistic method whose results have been verified with several wind-tunnel and flight tests [START_REF] Bronz | A Contribution to the Design of Long-Endurance Mini Aerial Vehicles[END_REF][START_REF] Bronz | Experimental comparison of 1-meter conventional, pusher and tractor flying wing aircraft configurations in s4 wind-tunnel[END_REF]. Figure 1 shows the brief flow chart of the program.

Design Envelope

In order to see the endurance and range performance limits of a one-meter aircraft, various configuration "cases" are investigated by CDSGN. Table 1 shows the main variables and their ranges used as an input on the program. Each of these variables is examined on different aircraft configurations with different airfoils, wing planforms, energy systems and stability characteristics. This facilitates the performance comparison of each case for the whole variable space. Comparisons will be explained in the related sections.

Configuration Selection

Although several aircraft configurations exist, in order to keep the design simple and parametric over the main variables shown in table 1, only conventional and flying wing configurations are considered. Each aircraft is generated automatically by CDSGN for each single main variable value. Figure 2 shows some of the automatically generated aircraft configurations as an example.

Figure 2: Examples of some automatic generated aircraft configurations by CDSGN.

The main objective of this design phase is to see the high-end performance capabilities of each configuration. Figure 3 shows the theoretical performance potentials of 1-meter conventional aircraft (on the left) within the initial specified design constants.The on-board energy, relatively represented by the circle radius, varies between 20 W h-210 W h corresponding to 0.1 kg -1.5 kg pure battery mass respectively. It can be seen that as the wing area increases, the saturation of maximum on-board energy reduces the maximum endurance and range performance. This was an expected behaviour; the performance of greater surface area configurations can easily be increased by increasing the maximum on-board energy limitation. However, the ease of operation (like hand launching) and certain UAV regulations favour keeping the total mass of the UAV under 2 kg [START_REF]Arreté du 21 mars 2007 relatif aux aéronefs non habités qui evoluent en vue directe de leurs opérateurs[END_REF].

On the right of Figure 3, the same type of performance plot for flying wing configuration is presented. In contrast to the conventional configuration, flying wing configuration tends to have better endurance and range performance for increased wing area. This is mainly due to having a lower maximum lift coefficient C Lmax in trimmed level flight conditions, in comparison to conventional configuration because of airfoil type and not having a separate horizontal tail.

International Micro Air Vehicle Conference and Flight Competition (IMAV2013)

17-20 September 2013, Toulouse, France

Effect of Solar Energy

In order to see the effect of using Solar Energy, the same type of analysis was done while taking into account the additional weight and limitations of solar cells, such as the maximum power point tracker board, connection cables, etc. Figure 4 shows the endurance and range performance of 1-meter conventional (on the left) and flying-wing (on the right) configuration with on-board battery and solar cells. For the conventional configuration, the increase of both endurance and range performance is clearly visible if it is compared with figure 3. Solar energy enhances the endurance up to 8 + hours, where it was around 6 hours before. The range increases 100 -150 km on average. Additionally, maximum endurance is evenly distributed over the various wing areas. The need for additional on-board energy is still visible for the bigger wing areas, but this time with the effect of solar cells covering the wing area, it is possible to use a smaller amount of on-board energy in order to be light. This ensures reduction of required flight power, which is maintained by the solar cells on the big wing surface. This is the reason why we see the configurations with smaller on-board energy having better endurance performance for the bigger wing surfaces, although the bigger wing creates more drag than the smaller wings at the same flight speed.

For the flying-wing configuration, the endurance performance seems to increase from 4.5 hours to 8 hours in certain cases. The cases with bigger wing areas, which carry small amounts of on-board energy, have slower cruise speeds and obtain more advantages from the solar energy. It should still be noted that these cases are too sensitive for flight conditions and environment such as a small amount of difference in the cruise flight speed or a decrease in solar irradiation, which end up shifting the performance values dramatically.

Conventional or Flying-Wing ?

Previous figures showed that both conventional and flying wing configurations are able to achieve 8 + hours of flight time with the help of solar energy. However, flying wing configuration is only able to reach that goal for a limited number of cases, while conventional configuration has well proven that the goal is reachable for a wide range of design variables. This is mainly due to the airfoil performance and the longitudinal trim issues that are limited with the flying wing configuration, especially at higher lift coefficients. This proves that the conventional configuration will be more robust for a real-

International Micro Air Vehicle Conference and Flight Competition (IMAV2013)

17-20 September 2013, Toulouse, France world application, where the conditions and environment change constantly. Additionally, non-solar performance of conventional configurations highly outperformed flying wing configurations on average. While making these conclusions, the whole envelope of design variables has taken into account for both configurations.Finally, the conventional configuration is selected as it outperformed the flying-wing configuration in these mission requirements.

Required Optimum Energy (On-Board)

The weight of the stored on-board energy is one of the key points in enhancing the platform performances, so it is important to select the optimum on-board energy that will give the best range and endurance performances for both solar and non-solar conditions. In order to see the effect of on-board energy on performance, the same analyses done previously are plotted in a different way. Figure 5 shows the endurance and range performance of different designs for non-solar configuration (on the left) and solar configuration (on the right) versus on-board energy. For non solar configuration, the effect of additional on-board energy is getting to a saturation around 150 W h , and above that value carrying more energy on board has no advantage for endurance performance. For the range performance, although there is a reduction on the increase, it keeps increasing with the additional energy. Caution should be given to the point that the given battery volume after a certain value will be impossible to carry on-board due to the size of the plane. On the right of figure 5, for solar configurations, it is assumed that 70% of the wing area is covered with solar cells. As seen, solar energy has a big influence on endurance performance, especially for the

International Micro Air Vehicle Conference and Flight Competition (IMAV2013)

17-20 September 2013, Toulouse, France configurations that have low on-board energy, thus lower weight. This result was already presented in the previous analysis, but the clear relation with the on-board energy has just been shown again. For solar configurations, an optimum on-board energy of 50 W h gives the best endurance performance. In these calculations, the specific energy of the battery was chosen as 190 W h/kg, which results in around 250 gr of battery weight for maximum endurance performance of a 1 m spanned aircraft in solar conditions. For the range performance, the behaviour is not different from the non-solar flight results. The more on-board energy that is carried on the aircraft, the better the range performance.

Design Philosophy

It was shown in the previous section that one configuration will not be capable of having the optimum performance in both endurance and range performance at solar and non-solar conditions. Designing the aircraft for only one particular condition will make it fail during almost all flight conditions different from the design point. The main focus of the design is to have a long-endurance aircraft, and the design philosophy is to obtain long-endurance performance in the operational cases as well, and not only on one particular condition.

Idea of Variable Configuration

It is possible to design different aircraft for different missions and conditions, however Eternity long endurance mini UAV should be as compact as possible. Therefore, concentrating on one wing planform which can be used with two different fuselages and various battery packages was found to be the most promising option. Finally, there will be different systems that will be appropriate for various missions and flight conditions. The battery pack is designed as a wing joiner simultaneously. The sizing of the joiner battery pack is also a key point, as it plays a big role on the wing planform and airfoil limitations. AMI cell1 (155mm × 60mm × 4mm), 4.5 Ah 3.7 V seemed to be the best fit for the wing joiner battery pack. Finally, a four cell battery pack is built in the joiner.

Design Details

With the battery cell type selected, the on-board energy was partially decided as a result. The wing joiner will consist of a four cell 4.5 Ah (66 W h) battery, which is suitable for long endurance mission in a sunny day. For overcast weather or a long range mission, it is preferable to use 2 packs (132 W h) of batteries. In order to fit an additional battery pack into the system, the fuselage has to be enlarged. This is not going to be beneficial for the long endurance version, as it will cause some additional drag. These points indicate that two separate fuselages are needed in order to optimize the system for both endurance and range missions using the same wing and tail parts. This decision also makes it possible to use different propulsion systems (propeller, motor, speed controller) optimized particularly for each mission, rather than compromising between them.

Figures 6 and7 shows the endurance performance for various wing areas and flight speeds for 66 W h and 132 W h on-board energy in both solar and non-solar conditions. The same performance evoluation is done for the range and taking into account the four conditions and additional operational requirements (such as take-off and landing), 0.14 m 2 of wing area performs the best.

Design of the Airfoil Family

In the conceptual design phase, suitable airfoils were selected as candidates by examining their maximum power coefficient (C 1.5 l /C d ) values. Then, the effect of different airfoils on endurance and range performance are investigated. Finally, by using the best existing airfoil as a starting geometry, a custom airfoil is designed by using XFOIL2 according to the main requirements of the design. The author would like to thank Philip Kolb for contributing to this very important phase with his knowledge and experience. PKMB500 airfoil is designed by Philip Kolb, which outperformed all previous airfoils existing in the selected database for our application.

The main reason why PKMB500 can give better results than its precedents relies on the fact that it is designed according to the Eternity's specifications. The wing area was designated 0.145 m 2 , and the battery weights were determined by the selected optimum capacities for solar and non-solar conditions. Since almost all of the main parameters of the design have been selected, it is easy to calculate the Re √ C L of the mean aerodynamic chord via equation 1 and later define all of the particular chord lengths proportionally.

Re C L = 1 µ 2ρW AR (1)
For the two configurations of Eternity (66 W h and 132 W h), calculated mean Re √ C L are 90k and 110k respectively. After obtaining the working regime, priority is given to the ((C 1.5 l /C d ) max ) value of the airfoil, as the endurance performance is the main objective for the design. SD7037 airfoil is selected as a reference because it is widely known and used in soaring competition gliders. By using the QDES routine of XFOIL, the laminar run on the top surface is extended by increasing the aftloading the airfoil. The transition ramp slope is slightly reduced and the pressure distribution at the bottom end is increased. This resulted an increase on the pitching moment increase on the final airfoil (∆C m = 0.025), which was acceptable as the conventional configuration was selected. The suction peak is smoothed, resulting in a better continuity of the flow and postponement of the transition. It is known that the bubble loss contributes significantly to the form drag on the low-Reynolds airfoils. One way of reducing the bubble size is to move it farther upstream by forcing the transition earlier. However, it will not be beneficial to completely prevent the occurrence of the bubble, because the laminar flow will be shortened dramatically. There is an optimum position for the bubble placement and size that corresponds to the minimum drag, and in low Reynolds airfoil design the control of the bubble becomes more critical compared to minimizing skin friction [START_REF] Drela | Low-reynolds number airfoil design for the mit daedalus prototype: A case study[END_REF].

In order to improve the airfoil performance along the span and decrease the viscous drag caused by inappropriate airfoil location, a series of airfoils are designed by referencing PKMB500.

PKMB803 is designed particularly for its own working regime, which is Re √ C L = 30000. Figure 8 shows the placement of the designed airfoils along the span, with a drag comparison of each with the root airfoil PKMB500 in their corresponding flow regime.

Having a smooth surface on the wing while transitioning from PKMB500 to PKMB803 on an elliptical planform created a dominant restriction on the thicknesses of PKMB601 and PKMB702. As the chord distribution has already been fixed, adding the root and the tip airfoil thickness values automatically defines the thickness distribution. With this thickness restriction, the only way found to slightly improve the performance of the middle airfoils was to compromise from the maximum lift coefficient and concentrate on the cruise coefficient regime corresponding to best endurance, which is around C L = 0.9 -1.0. Final dimensions of the two Eternity versions are shown in the figure 9, with their specifications in table 10.

Windtunnel Tests

In order to obtain the aerodynamic characteristics of the aircraft in a controlled environment, a wind-tunnel campaign was conducted at ISAE S4 wind-tunnel. An internal six component force and moment balance was used for the measurements.

Instead of manufacturing a separate wind-tunnel model, a fully equipped ready-to-fly version was manufactured and used for the tests. This was initiated in order to save some time, as one particular model takes one person about 3-4 weeks to build, including all the necessary equipment integration. Figure 11 shows the mounting of the aircraft in the wind-tunnel. The control surfaces and the motor of the aircraft were controlled by the on-board Paparazzi autopilot, which was connected to a computer via serial connection.

The experiments are conducted for different speeds in order to see the performance of the aircraft at each flight speed and evaluate the aerodynamic characteristics for a wider flight range. For each wind-tunnel speed the aircraft pitched from -6 to 12 degrees in order to obtain the performance 12 does not present the equilibrium condition performance. In order to have the equilibrium plot, elevator deflection has to be taken into account as well as at which it corresponds to zero pitching moment.

According to the previous analyses, the best endurance performance is achieved at C L = 0.9 -1.0, where the cruise speed corresponds to 12.5 -13.0 m/s, which is really close to stall speed of the big configuration that is 11.5 m/s. In real-life applications, for mini-micro UAV scales that are generally flying around 10 -25 m/s, sustaining airspeed accurately is really hard when the small momentum of the aircraft and the big ratio of the wind gusts to the flight speed are considered. Therefore, increasing the flight speed results in a reasonable safety margin. In the end of this decision, the cruise speed is selected as 14 m/s.

Flight Tests

Flight tests have been conducted by DelairTech3 UAV company under a partnership contract between ENAC, ISAE and TTT4 for the MIDLE project.

Close to 10 hours of flights have been recorded. The first prototype that has also been used in the wind-tunnel tests is given to the company for flights. Only one battery has been used, nevertheless the total weight of the aircraft (in flight ready state) was 930 g, just a bit heavier than the expected single battery configuration by 100 g. Additionally, the big fuselage is used as the thin fuselage was not build at that time. The longest flight that has been achieved was 2 hours and 7 minutes by using 54 W h of energy. The battery has given only 81% of its expected energy. The maximum glide ratio (L/D) and the power coefficient (C L 1.5 /C D ) is found to be 13.9 and 13.0 from the glide tests.

Considering all the differences of the test configuration compared to the design configurations, the flight test results act as a reference point rather than a direct comparison.

International Micro Air Vehicle Conference and Flight Competition (IMAV2013)

17-20 September 2013, Toulouse, France

Overall Design Conclusion

The main objective of this study was to show the feasibility of the Long Endurance Mini UAV concept.

Eternity shows a great perspective about what can be achieved with a one-meter aircraft. Finally, it can be concluded that with a one-meter aircraft that fits into an even smaller half meter carry-on luggage, 4 hours of flight are possible, and this can be enhanced up to 8 hours on a sunny day with the help of solar cells.

Figure 1 :

 1 Figure 1: Brief flowchart of CDSGN program.

Figure 3 :

 3 Figure 3: Performance plot of 1-meter conventional aircraft on the left and 1-meter flying wing aircraft on the right in various specifications with only battery. Circle radius represents the relative on-board energy being between 20 W h and 210 W h for min and max size.

Figure 4 :

 4 Figure 4: Performance plot of 1-meter conventional aircraft on the left and 1-meter flying wing aircraft on the right in various specifications with battery and solar cells. Circle radius represents the relative on-board energy being between 20 W h and 210 W h for min and max size.

Figure 5 :

 5 Figure 5: Effect of on-board energy on endurance and range performance, on the left with battery-only and on the right with solar energy addition.

Figure 6 :

 6 Figure 6: The effect of various flight speed and wing area on the endurance performance of non-solar and solar configurations with 66Wh of on-board energy.

Figure 7 :

 7 Figure 7: The effect of various flight speed and wing area on the endurance performance of non-solar and solar configurations with 132Wh of on-board energy.

Figure 8 :

 8 Figure 8: Placement of each airfoil is shown along the half span. The upper graphs shows the C l vs C d plots of each airfoils at their corresponding working regimes with the comparison to the root PKMB500 airfoil. The planform shown is given as demonstration and the real planform is shown later in figure 9.

Figure 10 :

 10 Figure 10: Geometrical specifications of the two Eternity configurations.

Figure 11 :

 11 Figure11: The mounting of Eternity on S4 Wind-tunnel, the payload compartment of the Eternity has been used for the mounting of the internal balance, and then the payload bay covered in order to protect the balance and to have the accurate fuselage shape.

Figure 12 :

 12 Figure 12: Effect of flight speed on the range and endurance performance characteristics.

Table 1 :

 1 The range of main variables used for the design envelope of Eternity.

	Main Variables	Range
	Aircraft Configuration Conventional or Flying-Wing
	Energy Source	Battery only or +Solar energy
	Wing Surface Area Flight Speed Battery Energy	0.05 -0.2 8 -20 20 -210	m 2 m/s W h

Table 2 :

 2 Table2shows the general specifications of off-the-shelf batteries. Lithium-Polymer batteries are used in all of the analyses, since Lithium-Sulfur batteries are almost impossible to obtain for our project because of their supplier policies(Sion Power[START_REF]Sion power, the rechargeable battery company[END_REF]). As a result of several discharge tests, the specific energy is taken as 190 W h/kg different from the general table. Typical battery specifications.

			Ni-Cd Ni-Mh Li-Po Li-S
	Specific Energy (W h/kg)	40	80	180	350
	Energy Density	(W h/l)	100	300	300	350
	Specific Power	(W/kg)	300	900	2800	600

Table 3 :

 3 Expected performance of the two Eternity configurations.

http://www.amicell.co.il/

xfoil.mit.edu

www.delair-tech.com

www.toulouse-tech-transfer.com