T Pottier
email: thomas.pottier@ensam.eu

H Louche
email: herve.louche@univ-montp2.fr

S Samper
email: serge.samper@univ-savoie.fr

H Favrelière

F Toussaint
email: franck.toussaint@univ-savoie.fr

P Vacher

A NEW FILTERING APPROACH DEDICATED TO HEAT SOURCES COMPUTATION FROM THERMAL FIELD MEASUREMENTS

This paper presents a new approach in heat source evaluation from thermal field measurements. A modal projection based on dynamics (Discrete Modal Decomposition) is used to estimate the spatial term of a heat diffusion problem. A numerical example is presented in order to validate this approach and compare it to a more classical spectral decompostion (based on strictly thermal considerations). Results show that the proposed projection basis not only provides closer assessement of the heat sources but is also more stable to measurement noise and side effects.

1.

INTRODUCTION

The study of calorific effects associated with material deformation has been widely investigated since the early works of Lord Kelvin. The thermoplasticty models now offer the ability to describe the heat generation during material loading from a phenomenological standpoint. Since the late 90s, experimenters have turned to the use of Infra-Red (IR) thermography in order to estimate, from temperature information, the heat sources related to plasticity, phase transformation and other dissipation related phenomenon [START_REF] Balandraud | Effet du couplage thermomécanique sur la propagation d'un front de changement de phase[END_REF][START_REF] Louche | Thermal observations associated with the Portevin-Le Chatelier effect in an Al-Mg alloy[END_REF][START_REF] Doudard | A probabilistic two-scale model for high cycle fatigue life predictions[END_REF].

Classically, the heat sources computation relies on performing both spatial and time derivation of the measured temperature fields. Under some classical assumptions [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localisation[END_REF], the heat diffusion equation is given by

ρC ∂θ ∂t -k ∂ 2 θ ∂x 2 + ∂ 2 θ ∂y 2 = w ch with θ(x, y,t) = T -T 0 , (1)
where ρ is the mass density, C the specific heat, k the conductivity coefficient, T (x, y,t) the current temperature field, T 0 (x, y) the initial temperature field and w ch the internal calorific heat sources. Such a definition allows the assessement of w ch from the time-derivative and the second space-derivative of θ(x, y,t) and the a priori knowledge of the material constants (ρ, C and k).

In practice, the time resolution of modern IR camera (up to 380 Hz) often allows the computation of a proper time derivation of θ. On the other hand, the space-derivative cannot be estimated as straight-forwardly. The measurement noise, intrinsic to IR-thermography, is amplified by space-derivation operators and massively impairs the calculation of the Laplace terms in Eq.(1).

Many numerical approaches have been developped in the last decade in order to estimate as accuratly as possible the spatial term of Eq.( 1). These regularisation methods can be sorted in 3 main categories:

• Filtering methods aimed at supressing the noise from the measurement using Gaussian [START_REF] Xiaoyan | Improving the peak wavelength detection accuracy of Sn-doped, H2-loaded FBG high temperature sensors by wavelet filter and Gaussian curve fitting[END_REF] or low pass filters based on linear, Fourier or wavelet definitions [START_REF] Yi | Source terms identification for the diffusion equation[END_REF][START_REF] Huon | Investigating the thermo-mechanical behaviour of cementitious materials using image processing techniques[END_REF].

• Inverse methods aimed at obtaining a quasi-solution through the use of one or several optimization algorithm [START_REF] Wong | A hybrid numerical/experimental technique for determining the heat dissipated during low cycle fatigue[END_REF][START_REF] Kaipio | Statistical inverse problems : Discretization, model reduction and inverse crimes[END_REF].

• Projection methods are among classical means in inverse heat transfer analysis. They consist of decomposing the measurement within a spectral basis built from eigen function of the Laplace operator [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localisation[END_REF][START_REF] Doudard | A probabilistic two-scale model for high cycle fatigue life predictions[END_REF]. These techniques have recently experienced various improvement through the use of wavelet decompositions [START_REF] Candau | Solving inverse heat diffusion problem with discrete wavelet transform[END_REF] or branch modes decomposition [START_REF] Neveu | Simulation de la conduction non linéaire en régime variable: décomposition sur les modes de branche[END_REF]. This latter approach allows the use of generalized boundary conditions which usually represent the main disadvantage of projection methods.

The present paper aims to propose the implementation of this latter approach and to compare the results obtained using a decomposition built from structural dynamics and another from thermodynamics considerations. For the sake of comparison, a numerical example is built and the heat sources are reconstructed using (i) the proposed modal basis and (ii) the classical spectral basis. Hence the influence the measurement noise is investigated.

SETUP OF PROJECTION BASES

Modal basis

This first projection basis is classically defined as the harmonic solution of a rectangular plate structural equation:

M • q + K • q = 0 with q = q(x, y) ∀k ∈ N, (2)
where M and K stand for the mass and the stiffness matrices respectively. Under such formalism, q(x, y) stands for the displacements which caracterize the modal shapes. Let's denote Q as the matrix whose columns are Q 0 , . . . , Q n-1 : the eigen vectors of Eq.( 2) and B D = (Q 0 , Q 1 , . . . , Q n-1 ) the associated modal basis. Thus, the set of modal coordinates α k resulting from the projection within a non-orthonormal basis is given as follows [START_REF] Goic | Multi-scale modal decomposition of primary form, waviness and roughness of surfaces[END_REF]

α k = α k,i ∈ R = Q T Q -1 Q T • q k (3)
The point of the present paper is to associate the displacement fields q k (x, y) to the measured temperatures θ k (x, y) through a projection operator as follows

θ k (x, y) = n-1 ∑ i=0 α k,i • Q i (x, y) + R n (x, y), ( 4
)
where n is the number of modes used for the projection (or truncation order). Even though B D is far from describing the actual underlaying physics, it may (in theory) lead to an exact solution if n tends toward infinity. In practice, the assumption of an infinite size is of course not relevant from a numerical standpoint, however the obtained basis exhibit several significant advantages. First of all, it naturally describes a space vector and thus allows a straight forward projection. In addition, the energy sorting of modes lead to an increasing complexity of shape (and frequency) when the truncation order n increase. Finally, it does not require the use of variable separation and therefore leads to a slow frequency increasing along i. In other words, the i th mode of B D has a smaller frequency than any other mode assessed assuming variable separation.

This latter point is of the utmost interest since the computation of the second order space-derivative (finite differences) is well known for significant divergence when high frequencies are involved.

Spectral Basis

Let us now consider the two dimensional Laplace equation in the rectangle

U = [0; a] × [0; b] × [t k ;t k+1 ],
with boundary conditions prescribed on the four edges of the rectangle. The Laplace problem is solved, using Dirichlet conditions over ∂U, as:

               ∆θ := ∂ 2 θ ∂x 2 + ∂ 2 θ ∂y 2 = ρC k ∂θ ∂t over U ∂θ ∂x (0; y,t) = ∂θ ∂x (a; y,t) = 0 over ∂U x ∂θ ∂y (x; 0,t) = ∂θ ∂y (x; b,t) = 0 over ∂U y , (5)
In order to sort the obtained basis in increasing order of frequencies, the modes are sorted by increasing products p × m. Double occurences are treated by placing the greater value of p ahead. One of the main consequence of the sorting operation is that it must be perfomed again if n (denoting the total number of modes or truncation order) varies. This disadvantage prevents from a single, once and for all calculation of B ∆ . Another disadvantage comes from the periodicity of the Q i vectors which requires high truncation order to discribe a non-periodic shape. However, this basis exhibits a major advantage: it is an hyperbolic solution of Eq.( 1) and is therefore close-to-physics since it is built from heat diffusion considerations.

NUMERICAL EXAMPLE

Description

The comparison of the two presented projection bases is performed by a numerical example. A square plate exhibiting two distinct heat sources is modeled in Abaqus. The finite element diffusion problem is solved numerically over a time period of t = 200 s. The resulting temperature fields are output to a matlab environment and interpolated over a 300×300 monotonic grid. A random white noise is added (standard deviation is denoted σ) and the projector Eq.( 4) is used to approximate the original thermal fields. Finally the space derivation of modes and a numerical time derivation of the obtained temperature fields (finite differences) allows the use of the heat diffusion equation Eq.( 1) and provides with an approximation of the imposed heat sources.

The chosen geometry is a 3D plane-parallel plate with dimensions (100; 100; 5). As shown in Fig. 1, two zones are chosen to be the active heat sources of the model. The first one is rectangular and host a pyramid-shaped distribution of heat. The second one is circular and exhibits a peak-shaped distribution of heat. The boundary conditions over the six sides of the model imposes a zero output heat flux.

Step #1

Diffusion Interpolation

Step #2

Step #3

Step #6

Step #5

Step #4

Comparison and results

The measurement noise generated by IR-camera capture may lead to significantly small signal-to-noise ratio if the recorded phenomenon involve small temperature rises. Indeed, the difficulties encountered in inverse heat transfer analysis mainly rely on the measurement noise which lead to a ill-posed problem which becomes very unstable to the double space derivation operator. Hence, the stability of the presented approach to noise is therefore highly relevant. Accordingly, a random white noise was artificially added to the numerical temperature fields with 30 increasing values of the σ. Fig. 2 shows the reconstructed heat sources profiles along the two lines that cross regions no.1 and no.2. Moreover, the global Root of mean square errors (RMS) are computed and plotted in Fig. 3. Two conclusions can be drawn from Fig. 2 and Fig. 3: firstly, as expected, the reconstruction error increases rapidly when the SNR decreases. Secondly, the modal projection presents an improvement from the spectral approach since the RMS error is smaller for every noise magnitude. Though, the gap between these two reconstructions becomes particularly visible for higly noised measurements (see Fig. 3), the relative improvement is actually about 10% for all the studied SNR.

This specific point constitues a major improvement of the proposed projection method since it allows a better heat source assessement when temperatures rises are small: e.g. small strain rates, low conductivity, early localization, vicinity of peaks... Another relevant point can be seen in Fig. 2-(c). In this case, the heat sources in region no.1 (pyramid-shaped) is actually very poorly estimated while the heat sources in region no.2 is far better approximated (especially with B D ). This emphasizes the shape influence on the reconstruction capability. It is well known that square-like shape leads to the use of high frequency modes and therefore a good approximation is disabled by a finite truncation order n. However, as proposed in [START_REF] Doudard | Identification of heat source fields from infrared thermography: Determination of 'self-heating' in a dual-phase steel by using a dog bone sample[END_REF], basis enrichement performed with tailor made sharp modes can be helpfull in such case.

CONCLUSIONS

The presented study shows that regardless of how far from physics the modal projection is, it allows significant improvements in heat source estimation. It appears to be less sensitive to noise and thus allows a more accurate caraterisation of calorific effects from temperature full-field measurements. It is also noticable that building the modes from structual dynamics allows to account various geometry such as disk, cylinder, hollowed plates, notched specimen, etc since a very simple FE calculation is required. Boundary conditions constitue a classical problem of thermal field projection methods. The presented approach allows to impose a wide range of them. Indeed the mechanical formulation of the problem leads to generate the modal basis from Finite Element Analysis and hence allows to impose either displacement, forces or velocity which all have thermal equivalences.

The spectral basis, is built from thermal considerations and thus is closer to physics. However, it appears to be less efficient in the presented numerical tests. The formulation of this basis tested in this paper is the classical one and not a enhanced one such as proposed by [START_REF] Doudard | A probabilistic two-scale model for high cycle fatigue life predictions[END_REF][START_REF] Doudard | Identification of heat source fields from infrared thermography: Determination of 'self-heating' in a dual-phase steel by using a dog bone sample[END_REF]. Indeed, it should be accounted that several flaws of such a decomposition have already be pointed out in many former studies and improvement procedure has been proposed and implemented [START_REF] Neveu | Simulation de la conduction non linéaire en régime variable: décomposition sur les modes de branche[END_REF].

Though, such a complex filtering approach is probably unnecessary when measurements are clear from noise and exhibit sufficient space resolution, it constitues a robust alternative in ill-posed inverse thermal problems. It must also be pointed that the presented numerical implementation leads to quick calculation since the modal basis is inverted once and for all. Processing a dozen images require less than a minute with an up-to-date computed.

5.

  Assuming a steady state case within [t k ;t k+1 ], it classically comes a harmonic solution related to the Laplace operator of Eq.(1) as θ k (x, y)

Figure 1 -

 1 Figure 1 -The six steps solving flowchart of the proposed validation example.

Figure 2 -

 2 Figure 2 -Heat sources (w ch ) reconstruction using modal and spectral projection for increasing noises (decreasing signal-to-noise ration : SNR). Graphs depicte the evolution of the reconstructed heat sources along two lines that cross the main heat sources (region no.1 and no.2 presented in Fig.1).

Figure 3 -

 3 Figure 3 -Comparison of the reconstruction error for increasing values of the signal-to-noise ratio (SNR).


Géraud Granger
email: granger@recherche.enac.fr-durand@cena.dgac.fr-alliot@dgac.fr

Nicolas Durand

Jean-Marc Alliot

Optimal Resolution of En Route Conflicts

Designing a Control Simulator involves many difficult problems: modeling conflict detection, trajectory uncertainties, solving conflict inside sectors, respecting military areas constraints, coordinating aircraft between sectors, etc. . . Moreover, the n-aircraft conflict resolution problem is highly combinatorial and cannot be optimally solved using classical mathematical optimization techniques. The set of admissible solutions is made of many unconnected subsets enclosing different local optima, but the subset enclosing the optimum cannot be found a priori.

In this paper, we present a conflict solver and its implementation in an Air Traffic simulator, with statistical results on real traffic in a French Sector. This solver, which takes into account real flight plans, solves every conflict inside a sector or over the French airspace on a loaded day.

Introduction

Whereas civil aviation authorities all around the world are engaging more and more controllers as the traffic increases, it seems that the available power of computers is not used and principles of control have not changed for the last 30 years: tools provided for Air Traffic Control (ATC) remain very basic, whereas aircraft are highly automated and optimized systems.

When building a new aircraft, engineers are now able to simulate the aircraft behavior before building a prototype. Many choices can thus be made before tying up expenses.

The same principle should be applied for designing future control systems or simply test minor changes on existing ones. The CATS/OPAS traffic simulator was designed to answer to this objective.

✁

The LOG is a common laboratory of the Centre d'Etudes de la Navigation Aérienne and the Ecole Nationale de l'Aviation Civile

The first versions of the simulator [START_REF] Durand | Optimal resolution of en route conflict[END_REF] were able to handle conflict detection and resolution in a Free-Route context: aircraft were supposed to fly direct routes from origin to destination, the airspace was not divided in sectors, conflicting aircraft were dynamically clustered before resolution. The simulator was able to test Free Flight hypotheses but not to simulate the controller's context and task.

A new modeling is introduced in this article to handle uncertainties on aircraft speeds, climbing and descending rates on flight plan routes.

The first part of the paper presents the state of the art for problem solvers and discusses the constraints hypothesis and goals chosen. Modeling is introduced in the second part. These two parts widely rehearse the article presented at the first USA/Europe ATM ✂ ✄ ☎ Seminar in Saclay [START_REF] Durand | Optimal resolution of en route conflict[END_REF], but were added for the sake of completeness and clarity. Part three details the maneuver modeling in the flight plan routes context and brings in the conflict solver. Some examples of resolution on real traffic and statistical results are given in part four.

Automatic conflict resolution 1.State of the art

Many people have studied the two aircraft problem (for example [START_REF] Prandini | Aircraft conflict prediction and resolution using brownian motion[END_REF]) which is certainly the most frequent but the easiest to solve. Projects, as for example AERA [START_REF] Niedringhaus | Automated En Route Air Traffic Control Algorithmic Specifications: Flight Plan Conflict Probe[END_REF][START_REF] Niedringhaus | A mathematical formulation for planning automated aircraft separation for AERA3[END_REF][START_REF] Niedringhaus | Automated planning function for AERA3: Manoeuver Option Manager[END_REF] or Free-R [START_REF] Hoffman Vu | Initial results of investigation into autonomous aircraft concept[END_REF] failed on the question: what happens if more than two, three, four . . . aircraft are involved in the same conflict or how can you guarantee this will never happen ?. Durand and Alliot showed [START_REF] Durand | Optimisation de Trajectoires pour la Résolution de Conflits en Route[END_REF] that in such cases, conflict resolution becomes a very complex and combinatorial mathematical problem. Feron [START_REF] Shewchun | Design and analysis of conflict resolution algorithms via positive semidefinite programming[END_REF][START_REF] Oh | Resolution of conflicts involving many aircraft via semidefinite programming[END_REF] gives a method to solve ✝ aircraft conflict but the modeling required to simplify the problem is completely unrealistic.

Other project failed on the modeling. ARC-2000 [K✆ 89, FMT93] optimized aircraft trajectories using ✞ dimensional cones and priority rules between air- craft. Optimum was not reached, and the system relied on the availability of FMS-4D for all aircraft, with no uncertainty on speeds. Zeghal [START_REF] Zeghal | Vers une thèorie de la coordination d'actions[END_REF], with reactive techniques for avoidance, gave a solution to the problem of automation which was robust to disturbance, but completely disregarded optimization. Furthermore, the modeling adopted implied a complete automation of both on board and ground systems and required speed regulation which could not be handled by human pilots and would probably be very difficult to apply to aircraft engines without damaging them.

A first approach to conflict resolution by stochastic optimization algorithms (genetic algorithms) was done by Alliot and Durand [START_REF] Alliot | Using genetic algorithms for solving ATC conflicts[END_REF]; more advanced results are presented in [START_REF] Daniel Delahaye | Genetic algorithms for partitioning airspace[END_REF][START_REF] Durand | Automatic aircraft conflict resolution using genetic algorithms[END_REF].

Specifications of the system

The main idea, guiding the design of the solver introduced in this paper, is to be as close as possible to the current ATC system: Constraints: the solver has to handle the following constraints: ✟ Conflict free trajectories must respect both air- craft and pilot performances. Considering the evolution of ATC toward automation [START_REF] Patrick Dujardin | Different paths to automation[END_REF], trajectories must remain simple for controllers to describe as well as for pilots to understand and follow. ✟ Trajectories must take into account uncertainties in aircraft speed due to winds, turbulence, unusual load, etc. Vertical speed uncertainties are particularly important. ✟ Flight plans must be respected as much as pos- sible. ✟ Maneuver orders must be given with an advance notice to the pilot. When a maneuver has begun, it must not be called into question.

Goals:

We want to achieve the following goals:

✟ find conflict free trajectories ✟ Simultaneously minimize different criteria :

1. the number of maneuver orders 2. the conflict resolution duration 3. the delay due to maneuvers ✟ compute these trajectories in real time.

General architecture of the system

We just sketch here the architecture of the simulator; each part will be detailed in the following sections. The system architecture is presented in figure 1 and 2. The system relies on three main processes P1, P2, and P3: ✟ P1 is the traffic simulator.

✟ P2 is in charge of conflict pair detection, clustering of pairs, and verification of new trajectories built by the solver.

✟ P3 is the problem solver.

P1 sends current aircraft positions and flight plans to process P2. Process P2 builds trajectories forecast for ✠ ✡ minutes, does conflict detection by pairs and transforms 1to-1 conflicts in n-aircraft conflict. Then, process P3 (the problem solver) solves in parallel each cluster, as aircraft in each cluster are independent from aircraft in the other clusters. The problem solver sends to P2 new orders and P2 builds new trajectories forecast based on these orders. Then P2 once again runs a conflict detection process to check that modified trajectories for aircraft do not interfere with aircraft in another cluster, or with new aircraft. If no interference is found, new flight orders are sent to P1. If there are interferences, interfering clusters are joined and the problem solver is used again on that (these) cluster(s). The process is iterated until no interference between clusters remains, or no new aircraft is concerned by modified trajectories. The new orders are sent back to the traffic simulator.

The above process is iterated and all trajectories are optimized each ☛ minutes. However, during the computation time, aircraft are flying and must know if they must change their route or not. ☛ should be large enough to compute a solution, send it to the pilot and let him time enough to begin the maneuver. Consequently, for each aircraft, at the beginning of the current optimization, trajectories are determined by the previous run of the problem solver and cannot be changed for the next ☛ minutes.

The Air Traffic simulator

One of the main goals of this project was to test the algorithms on real traffic. The Air Traffic Simulator takes as input flight plans given by companies and pilots: no pre-regulation is done neither on departure time nor on requested flight levels. Consequently, flight plans only have to be deposited ✠ ✡ minutes before take off.

The simulator uses a tabulated model for modeling aircraft performances: for a given aircraft type, it gives a vertical speed and a ground speed which depends on the aircraft attitude (whether it is climbing, leveled or descending). For example, a B747 leveled at FL-300 has a GS of ✞ ☞✌ kts. If it is climbing, its GS will be ✞ ✍✌ kts and its VS ✎✌ ✌✌ fts/mn. At FL-150, values would be respectively ✞ ✏✌, ✞ ✑✌ and ✎✍✌✌. Aircraft performances are in tabulated form describing ground speed, vertical speed, and fuel burn as a function of altitude, aircraft type and flight segment (cruise, climb or descent.) The main dataset for aircraft flight performance used is the base of aircraft data (BADA) performance summary tables derived from the total energy model of EUROCONTROL. 69 different aircraft types are described. Synonym aircraft are used to model the rest of the fleet. The Airbus A320 (EA32) is used as default aircraft.

All aircraft speeds are modified by a random value to take into account uncertainties on different factors (aircraft load, winds, etc. . . ) This value can be either computed once at aircraft activation and remains the same for all the flight, or can be modified anytime during the flight. The conflict detector and the conflict solver are impervious to the way this value is computed as long as it remains inside a given interval. Uncertainty modeling for conflict detection and resolution is discussed later in the article. Aircraft follow classical routes (from way-point to waypoint). The flight model is simple: an aircraft first climbs up to its RFL, then remains leveled till its top of descent, then descends to its destination.

t=0 t=1 t=2 t=3 t=5 t=6 t=4 ✒ ✒ ✒ ✒ ✒ ✒ ✒ ✒ ✒ ✒ ✒ ✒ ✒ ✒ ✒ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✔ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✕ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✖ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✗ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✘ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✙ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✚ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✛ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✜ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢ ✢
Aircraft fly with a time step that can be chosen at the start of the simulation. The time step is always chosen in order to guarantee that two aircraft face to face flying at 500 kts could not cross without being closer than one standard separation at at least one time step. For most of our simulation, we use a ✎✣ ✤ time step.

Conflict detection and clustering

Trajectory forecast and 1-to-1 conflict detection

As described above, the P2 process does trajectory prediction for ✠ ✡ minutes. This trajectory prediction is done again by a simulation on a slightly modified version of the Air Traffic simulator. But, as stated above, we assume that there is an error about the aircraft's future location because of ground speed prediction uncertainties1 . The uncertainties on climbing and descending rates are even more important. As the conflict free trajectory must be robust regarding these and many other uncertainties, an aircraft is represented by a point at the initial time. But the point becomes a line segment in the uncertainty direction (the speed direction here, see figure 3). The first point of the line "flies" at the maximum possible speed, and the last point at the minimum possible speed.

When changing direction on a beacon, the heading of the line segment's "fastest point" changes as described on figure 3.

To check the standard separation at time ✥, we compute the distance between the two line segments modeling the aircraft positions and compare it to the standard separation at each time step of the simulation.

In the vertical plane, we use a cylindrical modeling (figure 3). Each aircraft has a mean altitude, a maximal altitude and a minimal altitude. To check if two aircraft are in conflict, the minimal altitude of the higher aircraft is compared to the maximal altitude of the lower aircraft.

Let's take an example. A B747 is leaving its departing airport (altitude ✌) at ✥ ✦ ✌. ✎☞✧ kts, so max ground speed is ✎☞✧ ✩ ✎✪✑ ✦ ✑✏✧ kts. At ✏✫✌ fts, standard ground speed is ✎✍☞ kts, with a minimal ground speed of ✎✣✎ kts. So, the height of the segment grows much faster than the ✑✌★ factor for some aircraft. Duration ✠✡ can be changed, but must be at least equal to ✑ ✩☛. A good evaluation of ✠✡ is difficult. With a perfect trajectory prediction, the larger ✠✡, the better. However, this is not true as soon as uncertainties are included in the model. A large value of ✠✡ induces a large number of 1-to- 1 conflict, as sizes of segments modeling aircraft positions grow quickly with time. Therefore, the conflict solver can become saturated.

Clustering

After pair detection, P2 does a clustering which is a transitive closing on all pairs. Each equivalence class for the relation "is in conflict with", is a cluster.

For example, if aircraft ✬ ✭ ✮ are in conflict in the ✠✡ window, and if ✮ is also in conflict with ✯ in the same time window, then ✬ ✭ ✮ ✭ ✯ is the same cluster and will be solved globally by the conflict solver.

The conflict solver sends back to P2 maneuvers orders for solving conflicts. Then P2 computes new trajectories for all aircraft and checks if new interferences appear. For example, if the new trajectory given to aircraft ✮ to solve conflict with ✬ and ✯ interferes with cluster ☎ ✭ ✰ and with aircraft ✽ , then ✬ ✭ ✮ ✭ ✯ ✭ ☎ ✭ ✰ ✭ ✽ will be sent back to the problem solver as one conflict to solve. The process will always converge: in the worst case, P3 will have to solve a very large cluster including all aircraft present in the next ✠✡ minutes. However, this technique is usually efficient as a very large number of clusters can be solved very quickly in parallel.

t=0 t=1 t=2 t=3 t=5 t=6 t=4 B 0 1 B ✱ ✱ ✱ ✱ ✱ ✱ ✱ ✱ ✱ ✱ ✱ ✱ ✲ ✲ ✲ ✲ ✲ ✲ ✲ ✲ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✳ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✴ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✵ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✶ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✹ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✺ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼ ✼

The conflict solver

Maneuver modeling

In the horizontal plane, classical maneuvers given to aircraft are heading deviation. In the simulator, ✎✌, ✑✌ or ✏✌ degrees deviations will be allowed. The deviation starts on a virtual beacon created on the route (see figure ??). This beacon is defined by the position of the head of the segment at some time ✥✾. It ends on a second virtual beacon, posi- tion of the head of the segment at time ✥✿. An angle criteria is defined to find on which beacon the modified and initial routes should connect.

A maneuver will be determined by:

✟ ✥✾ which defines the first virtual beacon ✮✾. ✟ the deviation angle ❀. ✟ ✥✿ which defines the second virtual beacon ✮ ✿.

In the vertical plane, the aircraft trajectory is divided in ✞ periods (figure 5): ✟ Climbing period. In this period, aircraft can be lev-

❁ ❁ ❁ ❁ ❁ ❁ ❁ ❁ ❁ ❁ ❂ ❂ ❂ ❂ ❂ ❂ ❂ ❂ ❂ ❂ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❃ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❇ ❇ ❇ ❇ ❇ ❇ ❇ ❇ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❉ ❉ ❉ ❉ ❉ ❊ ❊ ❊ ❊ ❊ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ ❋ • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • Beacons Vertical plane

Horizontal plane

Vertical Manoeuvre Virtual Beacon Figure 6: Vertical maneuver during the climbing period.

✽ ❍ ✾ and starts again on a virtual beacon ✮ ✿ as stated on figure 6. ✽ ❍ ✾ and ✮ ✿ are defined by the position of the head of the uncertainty segment at time ✥✾ and ✥ ✿.

✟ Cruising period. When aircraft have reached their desired flight level, they may be moved to the nearest lower level to resolve a conflict. Aircraft starts descending when reaching a virtual beacon ✮ ✾ and starts climbing at ✮ ✿ (❀ ✦ ✌, ✮ ✾ and ✮ ✿ are defined by the position of the head of the uncertainty segment at time ✥✾ and ✥ ✿). An example of maneuver is represented on figure 7.

✟ End of Cruising period. When aircraft are about ✣✌ nautical miles from beginning their descent to destination, they may be moved to a lower level to resolve a conflict. Aircraft start descending on ✮ ✾ and are leveled at ✽ ❍ ✿ (❀ ✦ ✌) (see figure 8). ✮ ✾ and ✽ ❍ ✿ are defined by the position of the head of the uncertainty segment at time ✥✾ and ✥ ✿.

✟ Descending period. During this period no vertical maneuver is possible.

t=0 t=1 t=2 t=3 t=5 t=6 t=4 Vertical plane B 0 FL 0 B 1 ■ ■ ■ ■ ❏ ❏ ❏ ❏ ❑ ❑ ❑ ❑ ▲ ▲ ▲ ▲ ▼ ▼ ▼ ▼ ▼ ◆ ◆ ◆ ◆ ◆ ❖ ❖ ❖ ❖ ❖ ❖ ❖ ❖ ❖ ❖ ❖ ❖ ◗ ◗ ◗ ◗ ◗ ◗ ◗ ◗ ◗ ◗ ◗ ◗ | | | | | | | | | | | | ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ ❚ Beacons Horizontal plane
Vertical Manoeuvre Virtual Beacon ( ( ( ( ( ( ( ( ( ( ( ( ( ( ) No maneuver will be simultaneously done in the horizontal and vertical plane. This model has the great advantage of reducing the size of the problem. In order to solve conflict due to aircraft taking off or entering the airspace simultaneously at the same point, a variable of delay ✥❛ is introduced.

❯ ❯ ❯ ❯ ❯ ❯ ❱ ❱ ❱ ❱ ❱ ❱ ❲ ❲ ❲ ❲ ❲ ❲ ❲ ❲ ❲ ❲ ❲ ❲ ❲ ❲ ❲ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ( (
) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ❬ ❬ ❬ ❬ ❬ ❭ ❭ ❭ ❭ ❭ ( ( ( ( ( ( ( ( ( ( ) ) ) ) ) ) ) ) ) ) ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❴ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ Beacons Horizontal plane FL 1 Vertical Manoeuvre Virtual Beacon
For a conflict involving ✝ aircraft, the dimension of the search space is ✞ ✝. This will allow us to solve very difficult conflicts with many aircraft without investigating a large solution space.

Maneuver decision time

Because of uncertainties, a conflict that is detected early before it could occur may finally not happen. Consequently, deciding to move an aircraft in that case could sometimes be useless, and could even generate other conflicts that would not occur if no maneuver had been decided. This explains why controllers do not solve conflicts too early. With the turning point modeling, when there is no uncertainty, the earlier the maneuver is started, the lower the delay. However, if speed is not strictly maintained, the earlier the conflict is detected, the lower the probability it will actually happen. Thus, a compromise must be reached between the delay generated and the risk of conflict.

Choosing the model

Initially, aircraft are allowed to use their flight plan routes.

If we do not want to call into question previous maneuvers and be able to solve very large conflicts, we must try to start maneuvers as late as possible with respect to the aircraft constraints. This argument is enforced by the fact that we allow aircraft to have large uncertainties on their speeds2 .

For example, the first trajectory of figure 9, at ✥ ✦ ✌, cannot be modified before ✥ ✦ ☛. At the end of the first op- timization run, at ✥ ✦ ☛, the current position of the aircraft is updated. The maneuver that occurred between ✥ ✦ ☛ and ✥ ✦ ✑☛ is kept as a constraint for the second optimiza- tion run (on the example, no maneuver is decided). In the above example, we can see that the maneuver described on line ✑ (resulting from an optimization at ✥ ✦ ☛) is more pe- nalizing than the maneuver described on line ✏ (resulting from an optimization at ✥ ✦ ✑☛). This phenomenon occurs because of uncertainties. If uncertainties on speed are important, having a small ☛ will be very helpful to minimize the resolution costs in the real time situation. Pilots should only be given maneuver orders that will not be modified; if no conflict occurs, no order will be given.

The function to optimize

One of the principal algorithm design challenges is to define a suitable function to optimize. A multiple-criteria function is required that simultaneously attempts to:

✟ minimize the delay due to deviations imposed on air- craft. ✟ minimize the total number of resolution maneuvers required and the total number of aircraft that will be moved3 .

✟ minimize the maneuver duration so that aircraft are freed as soon as possible for maneuvers that may be necessary subsequently. ✟ enforce all separation constraints between aircraft.

Instead of considering a single scalar value that takes into account the different lengthenings of trajectories, the number of maneuvers and the conflicts between the aircraft, the contributions from each separate aircraft pair are maintained in a matrix F of size ✝ ✩ ✝ (where ✝ is the number of aircraft):

✟ If ❜ ❝ ❞ , ✽❡❢ ❣ measures the conflict between aircraft ❜ and ❞ in the optimization time window ✠✡. It is set to ✌ if no conflict occurs in this period and increases with the severity of the conflict. At each time step ✥, we compute ✯❤❢❡❢ ❣ as the difference (when positive) of the standard separation and the distance between the polygons ❜ and ❞ describing aircraft ❜ and ❞ position at time ✥. These values are added and give a measure of the conflict between ❜ and ❞ .

✽❡❢ ❣ ✦ ❤✐❤•❦ ❤❡• ♠ ♥ ❤♦✾ ♣ ✯❤❢❡❢ ❣ q
✟ If ❜ r ❞ , ✽❡❢ ❣ measures the efficiency of the resolution between aircraft ❜ and ❞ . It is set to ✌ if no conflict can happen between ❜ and ❞ after the optimization time window ✠✡. If a conflict may remain after this period, ✽❡❢ ❣ gives a bad mark to pairs of aircraft for which the difference of heading and speed are small (these conflicts are difficult to solve).

✟ ✽❡❢❡ (see equation 1,2,3) measures the takeoff (or en- tering) delay given to aircraft ❜ (✯ ❛ is a constant), the maneuver duration time ♣ ✥✿ s ✥✾ q and trajectory lengthening (✯ t is a constant depending on the ma- neuver angle ✤), and the number of maneuvers (✯ • is a constant multiplied by ✎ if a maneuver is supposed to become definitive and ✌ if not):

✽❡❢❡ ✦ ✯❛ ✥❛

(1)

✉ ✯t ♣ ✥✿ s ✥✾ q

(2)

✉ ✯• ✈ ♣ ✥✾ ✇ ✑ ☛ q ✄ ♣ ✥✿ r ✥✾ q① (3)
This matrix contains much more information than a scalar global value ✽ , and is useful in the optimization algorithm used. However, a global scalar value is required, and can be defined as follows:

② ♣ ❜✭❞ q ✭ ❜ ③ ✦ ❞ ✭ ✽❡❢ ❣ ③ ✦ ✌ ④ ✽ ✦ ✎ ✑ ✉ ⑤ ❡ ⑥ ♦ ❣ ✽❡❢ ❣ ⑦ ♣ ❜✭❞ q ✭ ❜ ③ ✦ ❞ ✭ ✽❡❢ ❣ ✦ ✌ ④ ✽ ✦ ✎ ✑ ✉ ✎ ✑ ✉ ⑤ ❡⑧ ❣ ✽❡❢ ❣

A global optimization problem

Use of local methods, such as gradient for example, is useless here, because these methods rely on the arbitrary choice of a starting point. Each connected component may contain one or several local optima, and we can easily understand that the choice of the starting point in one of these components cannot lead by a local method to an optimum in another component. We can thus expect only a local optimum.

Genetic Algorithms applied to conflict resolution

Genetic algorithms (GAs) are global stochastic optimization technics that mimic natural evolution. They were initially developed by John Holland [START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF] in the sixties. The subject of this paper is not GAs and the interested reader should read the appropriate literature on the subject [START_REF] Goldberg | Genetic Algorithms[END_REF]. The general principles are given on figure 10.

Genetic algorithms are a very powerful tool, because they do not require much information and are able to find many different optima that can be presented to a human operator. Moreover, we know much about the function to optimize and this information can be used to create adapted crossover [START_REF] Durand | Automatic aircraft conflict resolution using genetic algorithms[END_REF] and mutation operators, an other advantage of GAs over other optimization technics.

Genetic algorithms are very efficient for solving global combinatorial optimization problems but are not very efficient for solving local searches with a good precision. Consequently, in the last generation of the genetic algorithm, a local optimization method is used to improve the best solution of each chromosome class defined above: a simple hill-climbing algorithm is applied to the best chromosome at the end of the GA run.

Results

We present here examples of resolution that illustrate5 the performance of the algorithm. These examples were computed on a Pentium III ☞✏✏. In the following, the time window for prediction is fixed at ✎✑ minutes (✠ ✡ ✦ ✎✑ mn) and an optimization is computed every ✏ minutes (☛ ✦ ✏ mn).

No uncertainty is considered in the following examples. Implementation of the uncertainty modeling is still in progress for the moment. The final version of the paper will include uncertainty parameters.

Example of Two-Aircraft Conflict

In this first application (see figure 11), at 15:06:00 UT a conflict is detected between two aircraft numbered ✫✧✞✞ and ✫✫✞✍. Aircraft ✫ ✧✞✞ is descending from Flight Level ✏✣✌ to its destination and aircraft ✫✫✞✍ is flying at FL ✑✏✌. The conflict appears near the beacon. To solve the conflict, aircraft ✫✫✞✍ is turned left (30 degrees) on a point created ✎✞ nm before the beacon. A second turning point is created at a distance of 9 nm from the first one. After the second turning point aircraft is directed to the next beacon of its flight plan.

Complex conflict involving ⑩ aircraft

In this example, at 05:45:00 UT, ✏ aircraft are cruising at FL 370 and one is cruising at FL 350. Aircraft ✑✑☞✍ is first in conflict with aircraft ✎✌✣✧, then in conflict with aircraft ✍✑. It cannot be moved vertically because of aircraft ✑✏✍ cruising at FL 350. Three horizontal maneuvers (figure12) are given to aircraft ✑✑☞✍, ✍✑ and ✑✏✍.

A complete test

First simulations have been done with the new modeling without uncertainties. Simulations with uncertainties will be available within six month. A complete experiment done with unregulated flight plans of the 21th of May 1999 is described here. It involves ✧✣✞✌ aircraft over France. We only detect and solve conflicts above ✎✌✌✌✌ feet, as we are only interested in En Route conflicts. Aircraft entering Paris TMA control area are sequenced on the TMA entry points, but no control is done inside the TMA. Slots are given to aircraft in order to prevent simultaneous take off and landing.

When running this one day test with a very basic conflict detection algorithm (only actual conflicts are detected, with no uncertainty on speed) and with no conflict resolution, ✑✧✣✍ conflicts are detected with flight plan routes and ✎✍✌✏ conflicts are detected with direct routes. The mean flight duration is ✣ ✎ minutes with flight plan routes and ✞✫ minutes with direct routes.

When running the complete simulation with detection and resolution, fixing ☛ ✦ ✏ minutes and ✠ ✡ ✦ ✎✑ With direct routes, the number of maneuvers falls to ✎✫✫✧ and the number of maneuvered aircraft is ✎✏✧✏. The mean maneuver duration per maneuvered aircraft is ✎✍✤ which represents ✌ ✪✫✧★ of the flight duration. Maneuvered aircraft are delayed of ✎✑✤ on average. The global mean delay is ✑✤.

Further work:

Within the next six month, results with the uncertainty modeling will be available and added in the article.

The next step will include in the modeling the division of airspace into sectors. Clusters will have to be solved inside sectors and a coordination process should be defined.

The modeling should also take into account military zones.

Conclusion

The conflict solver introduced in this paper is a step toward simulation of en route control.

The goal of this work was to show that previous work done on automatic traffic resolution could help designing a controlled traffic simulator that could lead to the development of future ATC systems.

The simulator remains small in size (5000 lines of code in Caml language).

Figure 2 :

 2 Figure 1: General architecture

Figure 3 :

 3 Figure 3: Modeling of speed uncertainties.

Figure 4 :

 4 Figure 4: Horizontal maneuver modeling.

  Figure 5: Vertical maneuver modeling.

Figure 7 :

 7 Figure 7: Vertical maneuver during the cruising period.

Figure 8 :

 8 Figure 8: Vertical maneuver during the end of cruising period.

Figure 9 :

 9 Figure 9: The model and real time optimization.

Figure 10

 10 Figure 10: GA principle the value is smaller than ✿ ⑨ , maximizing the function min- imizes the remaining conflicts. When the value is larger than ✿ ⑨ , maximizing the function minimizes the possible re- maining conflicts after the optimization time window, the number of maneuvers, their duration, and the delays induced by maneuvers. When no conflict and no maneuver occurs, the function is equal to ✎.

  Its climb rate is ✎✍✌✌ fts/mn and its ground speed is ✎✧✣ kts. If we suppose that ground speed uncertainty is ✣★ and vertical speed uncertainty ✑✌★, maximal and minimal climb rate are ✎✍✌✌ ✩ ✎✪✑ ✦ ✑✎✫✌fts and ✎✍✌✌ ✩ ✌✪✍ ✦ ✎✞✞✌ fts/mn and ground speeds are respectively ✎✍✞ and ✎✫✫ kts. This means that ✎✣✤ later, the fastest and higher point has traveled ✌✪✧✫ Nm and ✣✞✌ fts while the slowest and lowest has only traveled ✌ ✪✫☞ Nm and ✏✫✌ fts. But this time, when

computing maximal and minimal speeds, the difference of altitude of both points must be taken into account. At ✣✞✌ fts, the tabulated model gives a standard ground speed of

Table 1 :

 1 minutes, the P2 process finds ✍✌✑✣ clusters of different size (see figure 1. There is no unsolved cluster and consequently no conflict remains. Cluster size Only ✎☞✧✌ aircraft are given ✑✣✣ ✎ maneuvers which represents ✎✪✏ maneuver per aircraft. The mean duration of a maneuver is ✑❶ ✝ ✎✞ ✤. Details on maneuvers are given in table2.The mean maneuver duration per maneuvered aircraft is ✞ ✣✤ which represents ✎✪✞ ✧★ of the flight duration. Maneu- vered aircraft are delayed of ✎☞✤ on average. The global mean delay is ✞✤.

	Cluster size Number
	2	6325
	3	1158
	4	331
	5	112
	6	52
	7	19
	8	15
	9	8
	10	1
	11	3
	12	0
	13	1

Table 2 :

 2 Maneuvers repartition.

Uncertainties on ground track will not be considered, as they do not increase with time and will be included in the standard separation

We do not plan to solve conflicts by speed modifications. Theoretical study shows that optimal En Route conflict resolution by speed modifications would require large anticipation time (anticipation time depends on different parameters such as angle of convergence, speed margins for each aircraft, standard separation etc; more details can be found in[START_REF] Durand | Optimisation de Trajectoires pour la Résolution de Conflits en Route[END_REF]). This is quite unrealistic due to aircraft speed uncertainties.

Thus, instead of sharing the global delay on all the aircraft, some aircraft will support a part of the delay and others will not.

The label gives the number of the aircraft, its heading, its flight level and its horizontal speed


Young-Hwan Kim
email: young-hwan@inria.fr

Roberto Konow

Diego Dujovne

Thierry Turletti

Walid Dabbous

Gonzalo Navarro

PcapWT: An Ecient Packet Extraction Tool for Large Volume Network Traces

Keywords: Network traces, Packet indexing, Packet extraction, Wavelet Tree, Multi-thread le I/O

Network packet tracing has been used for many dierent purposes during the last few decades, such as network software debugging, networking performance analysis, forensic investigation, and so on. Meanwhile, the size of packet traces becomes larger, as the speed of network rapidly increases. Thus, to handle huge amounts of traces, we need not only more hardware resources, but also ecient software tools. However, traditional tools are inecient at dealing with such big packet traces. In this paper, we propose pcapWT, an ecient packet extraction tool for large traces. PcapWT provides fast packet lookup by indexing an original trace using a Wavelet Tree structure. In addition, pcapWT supports multi-threading for avoiding synchronous I/O and blocking system calls used for le processing, and is particularly ecient on machines with SSD. PcapWT shows remarkable performance enhancements in comparison with traditional tools such as tcpdump and most recent tools such as pcapIndex in terms of index data size and packet extraction time. Our benchmark using large and complex traces shows that pcapWT reduces the index data size down below 1% of the volume of the original traces.

Moreover, packet extraction performance is 20% better than with pcapIndex. Furthermore, when a small amount of packets are retrieved, pcapWT is hundreds of times faster than tcpdump.

Introduction

The volume of network packet traces becomes larger and larger, and much more complex than before. The reason is that the speed and size of networks have increasingly expanded. Thus, much more H/W resources (e.g., massive storage space, stronger computing power, and so on) are necessary to deal with big traces including numerous packets and complex trac patterns. Furthermore, ecient tools are required to analyze the data eectively. However, traditional tools such as tcpdump [START_REF] Jacobson | Tcpdump[END_REF], tcpow [START_REF] Elson | Tcpow[END_REF], or tcpstat [START_REF] Herman | Tcpstat tool[END_REF] are inecient to handle very large packet traces.

Among current network trac analysis tasks such as protocol performance evaluation, network monitoring and forensic investigation, packet lookup is one of the most basic and important functions to investigate errors and to evaluate performance. In particular, packet lookup is a CPU-greedy task, especially when it deals with a huge packet trace including many complex trac patterns. However, most of the traditional tools use linear search algorithms, which are simple, but take a long processing time proportionally to the packet trace size. Basically, a packet trace is a very long array, so the time complexity of a linear search is comparatively higher than that of other algorithms, such as binary search, balanced tree (B-tree), and hashing. Moreover, most of available packet trace formats (e.g., pcap [START_REF] Jacobson | Tcpdump[END_REF]) do not include the number of matched packets and their locations. Consequently, the search delay can rapidly increase with the length of the trace and the complexity of the ltering query. In addition, traditional tools (e.g., tcpdump and wireshark [START_REF] Combs | [END_REF]) are inecient for iterative operations and possible future reuse, since they do not maintain historical data or reusable information.

In order to enhance packet trace analysis, a number of contributions have been published, such as fast packet ltering [START_REF] Deri | High-speed Dynamic Packet Filtering[END_REF], packet trace compression [START_REF] Holanda | Performance Analysis of A New Packet Trace Compressor Based on TCP Flow Clustering[END_REF][START_REF] Politopoulos | Evaluation of Compression of Remote Network Monitoring Data Streams[END_REF], and network statistics information extraction [START_REF] Estan | [END_REF]. Moreover, a few practical tools have also been recently proposed, such as PCAP-Index [START_REF]Pcap-index[END_REF], Compressed Pcap Packet Indexing Program (CPPIP) [START_REF]Cppip[END_REF], and pcapIndex [START_REF] Fusco | PcapIndex: An Index for Network Packet Traces with Legacy Compatibility[END_REF][START_REF]PcapIndex[END_REF]. Overall, they use extra data sets for fast search, which are extracted from the original trace.

However, these tools result in poor performance, as we show in Section 5.3. PCAP-Index running on top of a database requires too much space for handling the index data, and also its packet extraction procedure is consid-erably slow. CPPIP can save storage space by compressing original traces, but it can only support a few simple queries. PcapIndex is much faster than these tools, and is able to retrieve and extract matched packets from large size trace les (up to 4.2 GB). However, the index data size built by pcapIndex is abnormally increased when the number of packets is in the millions. Thus, there is a trade-o between the index data size and the packet extraction performance when processing a large packet trace.

In this paper, we propose pcapWT, a fast packet extraction tool for large network traces, and evaluate its performance in terms of index data size and packet extraction time. pcapWT adopts an advanced data structure named Wavelet Tree (WT) [START_REF] Grossi | High-order Entropy-compressed Text Indexes[END_REF], which enables a fast search and high compression ratio at the same time. Moreover, this tool supports multi-threading, which is able to enhance random le read and write performance over Solid State Drive (SSD). This paper is structured as follows. Section 2 presents the state of the art. Section 3 provides a background on Wavelet Tree, and Section 4 describes the design of pcapWT in detail. Section 5 evaluates the performance of pcapWT and compares its performance with other tools. Finally, Section 6 concludes the paper.

Related Work

Today most network trace analysis tools, such as tcpdump [START_REF] Jacobson | Tcpdump[END_REF] and wireshark [START_REF] Combs | [END_REF], run on a single thread, and their complexity increases linearly with the volume of original trace les. Possible solutions to improve their performance include using a higher processor clock speed, replacing the SSD by a Hard Disk Drive (HDD), or splitting the large packet trace into multiple pieces. For that reason, several tools have been recently published to enhance performance of packet extraction on large traces, such as CPPIP [START_REF]Cppip[END_REF], PCAP-Index [START_REF]Pcap-index[END_REF], and pcapIndex [START_REF] Fusco | PcapIndex: An Index for Network Packet Traces with Legacy Compatibility[END_REF].

CPPIP uses bgzip [START_REF] Li | The Sequence Alignment/Map Format and SAMtools[END_REF] (i.e., a block compression/decompression software) to reduce the volume of original packet trace les. This tool extracts an index data from the compressed trace le, and lters out matched packets directly from the compressed le. Although this tool is able to reduce the storage space for the original trace les, it can only support simple queries, such as packet number and received timestamp. In addition, CPPIP needs a signicant amount of space (about 7% of the original trace le) to store the index data.

PCAP-Index uses a database (SQLite Ver. 3) to build the index data and to perform packet lookup. Thus, this tool is more exible than other command-line based tools to express queries. However, as discussed in Section 5.3, the performance obtained is poor in terms of time needed to build index data, packet lookup time and index data size.

PcapIndex, having a similar name as PCAP-Index, is currently part of a commercial network monitoring solution [START_REF]PcapIndex[END_REF]. In order to reduce both index data size and packet lookup time, this tool adopts an encoding method using a bitmap compression technique based on a pre-dened codebook, named COMPressed Adaptive index [START_REF] Fusco | NET-FLi: On-the-y Compression, Archiving and Indexing of Streaming Network Trac[END_REF]. Thus, it obtains better performance than CPPIP and PCAP-Index, in terms of index data size and packet extraction time. In comparison with tcpdump, this tool can reduce packet extraction time up to 1,100 times. Moreover, its index data size only takes about 7 MB per GB of the original trace le. However, as discussed in Section 5.3, PCAP-Index is not as fast compared to what is mentioned in the paper [START_REF] Fusco | PcapIndex: An Index for Network Packet Traces with Legacy Compatibility[END_REF],

when it extracts a large amount of packets from a big trace le. In particular, the index data size rapidly increases when the volume of original trace le is greater than 4.2 GB.

Wavelet Tree

In this work, we present a compact index for pcap-traces that is built on top of compact data structures. To understand how these data structures are employed, we rst need to dene three basic operations: rank, select and access.

Given a sequence S containing n symbols from an alphabet of size σ (e.g., σ = 2 for binary sequences), rank b (S, i) counts the number of occurrences of symbol b until position i, select b (S, j) nds the position of the j-th occurrence of b in the sequence, and access(S, k) returns the symbol that is located at position k.

In the case of binary sequences, the queries described above can be solved in constant time O(1) by using a bit sequence representation presented by Munro [START_REF] Munro | Tables[END_REF]. Gonzales et al. [START_REF] González | Practical Implementation of Rank and Select Queries[END_REF] show that practical implementations (RG) are possible to achieve constant time for rank and access queries, and logarithmic time for select queries using 5 -37% extra space. Raman, Raman and Rao (RRR) [START_REF] Raman | Succinct Indexable Dictionaries with Applications to Encoding k-ary Trees, Prex Sums and Multisets[END_REF] presented another solution that is able to compress the bit sequence to H(S)+o(n), which corresponds to the empirical entropy plus some low order terms. This representation is able to compress the sequence up to 50% in practice while allowing to perform rank, select and access queries in constant time. Note however that performing these operations using RRR is signicantly slower than using the RG representation.

For general sequences (i.e., σ > 2), the Wavelet Tree [START_REF] Grossi | High-order Entropy-compressed Text Indexes[END_REF][START_REF] Navarro | Wavelet Trees for All[END_REF] data structure can be used. Its name is derived from an analogy with the wavelet transform for signals, which operates in a similar way. The wavelet tree is a perfect balanced binary tree, that stores a bit sequence in every node except the leaves. Every position of the bit sequences is marked with either a `0' or `1' depending on the value of the most signicant bit of the symbol at the position in the sequence; this can be seen as dividing the alphabet into halves. Starting from the root, the symbols that are marked with a `0' go to the left subtree, while the rest go to the right subtree. This segmentation of the values continues recursively at the subtrees with the next highest bit. The tree has σ leaves, and each level of the tree requires n bits, the height of the tree is log 2 σ, thus we need n⌈log 2 σ⌉ bits to represent the tree. Each sequence of bits at every level must be capable to answer access, rank and select queries, thus we can use the RG or RRR representation to handle the bit sequences. The wavelet tree is able to solve rank, select and access operations in O(log σ) time, that is to say, the execution time depends only on the size of the alphabet, not on the length of the sequence. Note that both the time to execute the operations and the size of the data structure heavily depends on the size of the alphabet σ. In the case that we use RRR to represent the bit sequences at each level, the wavelet tree uses nH(S) + o(n log σ) bits of space. If we want to compress the data structure even more, we can change the shape of the tree to the Human shape of the frequencies of symbols appearing in S, and maintaining the O(log σ) complexity for handling the operations. In practice, if we use the Human shaped wavelet tree using using RRR for the bitmaps, we can reduce the space of the data structure by 12% [START_REF] Claude | Practical Rank/Select Queries over Arbitrary Sequences[END_REF]. Figure 1 shows an example of a regular wavelet tree and a Human shaped wavelet tree built over a small sequence.

We will explain how to perform access operation by an example based on the regular wavelet tree from Figure 1. Let us assume that we want to perform access(S, 6), that is, return the symbol located at position 6. We start at position 6 at the root bit sequence (B root ) and ask if the corresponding bit is marked as `1' or `0'. If the bit is `0', we go to the left branch, if not, to the right. In this case B root [START_REF] Holanda | Performance Analysis of A New Packet Trace Compressor Based on TCP Flow Clustering[END_REF] = 1, so we go to the right branch of the tree. Now we have to map the original position (k = 6) to the corresponding position at the right branch. In other words, we want to know how many 1 ′ s were at the root bit sequence before position 6. We can easily do this in constant time by performing the operation supported by RG and RRR bit sequences: rank 1 (6, B root ), which returns the value 3. Using this value we can now go to the right branch of the root and execute the same procedure by setting k ′ = rank 1 (k, B root ) = 3. This process is done recursively until we reach a leaf. At the leaf we know the symbol at position 6 from the original sequence S. A very similar procedure is done to perform rank b (S, i) and select b (S, j). We refer the reader to the previous work from Navarro [START_REF] Navarro | Wavelet Trees for All[END_REF] and Gagie et al. [START_REF] Gagie | New Algorithms on Wavelet Trees and Applications to Information Retrieval[END_REF] for an extensive explanation of how these operations are performed and the virtues of wavelet trees with its wide range of possible applications. every item from the data set. The index data is created using part of packet headers or additional information possibly available.

The rst step is to extract index data from a packet trace le, such as packet oset, Ethernet source (Src) / destination (Dst) addresses, IPv4 Src/Dst addresses, Src/Dst port numbers, and Layer 4 (L4) protocol type. Here, the packet oset stands for the distance in bytes between two consecutive packets, which corresponds to the packet size in bytes. The reason why the oset is replaced by the packet size is that sequentially accumulated values consume a lot of memory. In addition, such large numbers are inecient for the WT compression, since they increase the required number of bits to represent them. As shown in Figure 2, all index data are stored into each array, and they include the same number of elements with the number of total packets in the source trace le.

However, as mentioned in Section 3, the size of the WT compressed data is highly aected by the number of bits per element and the number of elements. For instance, when an element set (e.g., packet sizes referring packet osets) of index data is S = {100, 200, 200, 300..., n}, the required number of bits per element is (⌈log 2 σ⌉), where σ is the maximum value among these input elements. Thus, for enhancing the compression performance, the elements have to be converted into sequential positive integer numbers (Σ = {1, 2, ..., σ}), and long arrays must be divided into multiple pieces (S = {S 1 , S 2 , ..., S N }). For example, as shown in Figure 3, IPv4 addresses are mapped into positive integer numbers in consecutive order. In addition, to increase eciency, we use a balanced tree (B-tree) to map the addresses with a positive integer. Consequently, a trace le accompanies a number of mapping tables, and the tables are provided as part of the index data.

In the last step for building index data, as shown in Figure 4, the WT builder creates one chunk of compressed index data every one million elements, and each chunk is stored in an independent le. This value has been empirically set to minimize the compressed index data size.

Meanwhile, the packet oset is not compressed by WT, because the large value of σ causes signicant degradation of the compression performance. As an alternative, we use an ecient array provided by the libcds package [22], which removes redundant bits.

Packet lookup and extraction

Figure 5 illustrates the packet lookup procedure. First, this procedure reconstructs the packet oset by accumulating the packet sizes. For example, if the packet sizes are P = {100, 200, 200, 300, ..., p n }, their osets are O = {0, 100, 300, 500, 800, ..., o n }. Then, it loads an index data related to the query. At this step, the index data is located in main memory, however it does not need to be decompressed to execute a lookup. The query syntax is as follows,

Q := Q type Q operand := Q 1 and/or Q 2
A basic query consists of a query type (Q type ) and a query operand (Q operand ).

The query type denes the kind of index data extracted from the traces, and pcapWT supports pkt_num, eth_src, eth_dst, ipv4_src, ipv4_dst, port_src, port_dst, and l4_proto. The query operand is the value that needs to be found on a data set of index data which is indicated by the query type. For example, the Q type is ipv4_src and the valid Q operand is 10.1.1.1. The query syntax also allows to perform intersection (and) or union (or) of dierent queries, thus we can construct more complex queries, such as Q 1 and Q 2 and ... or Q n . More precisely, the query operand is converted to a positive integer number, using the mapping table provided from the index data build process. As mentioned above, the reason is that all index data, except oset, is mapped to positive integers. If the mapping table does not have the same entry as the query operand, this indicates that there are no matched packets for the query in the trace.

Regarding basic functions, as described in Section 3, we can get the number of matched packets (m) containing the query operand by performing rank() of the WT. Next, we can gure out indexing numbers of matched packets by iterative executing of select(), for example Y = {2, 301, 303, ..., y m } as shown in Figure 6.

The matched packets are marked on the matching table, which is composed of a boolean array with the same number of entries as the number of packets included in the source pcap le. If there are multiple queries concatenated by and/or operators, those steps are repeated -from loading related index data-through updating the matching table. As shown in Figure 6, the matching table assigns `1' to its entries for the matched packets which are found by the rst query or queries located after the or operator. On the other hand, if the packets are found by queries located after the and operator, the table assigns `1' only to entries kept `1', otherwise `0'. Once the lookup procedure is completed, two lists are created where osets and packet sizes for the nal matched packets are marked `1' on the matching table.

The performance of packet extraction highly depends on the capability of storage devices. In particular, this function performs random read and sequential write from/to les. The random access performance is signicantly degraded in comparison with sequential access, in terms of read and write throughput. The main reason is related to the total idle time caused by multiplication of the average access delay (e.g., ranging from dozen to hundred milliseconds) and the number of random accesses. In addition, performance enhancement using cache for pre-reading consecutive data and a buer for full padded block writing are not allowed in random I/O operations, unlike the sequential ones. However, it is possible to reduce the access delay by replacing HDD with SSD. On the other hand, le I/O libraries are also deeply related to performance degradation. Especially, in the standard C/C++ library, le seeking functions such as lseek() and fseek(), take a lot of time at each call, and generate frequent idle operations.

Therefore, we need an ecient way to compensate the performance degradation. SSDs have a very short access delay (e.g., below 1 millisecond) compared to HDDs, thus we can reduce the idle time by distributing the I/O requests into multiple processes, without signicant delay overhead. Consequently, the proposed design uses multiple threads to handle multiple requests simultaneously, which are performed by OpenMP provided by gcc [23]. For more details, we provide numerical performance comparisons in Subsection 5.1.

As shown in Figure 7, a number of threads are used to balance the overall workload for the packet lookup and extraction tasks. For example, if the number of matched packets is m and the number of threads is w, each thread will retrieve m/w packets to an output le simultaneously. More precisely, each thread already knows the exact writing locations in the output le before the packet retrieval stage, by summing up the packet sizes that are assigned to each thread. Inside of each thread, the target packets are copied into user space memory from the source trace le, and then they are written into the output le. Then the oset position and the size of selected packet are written into the oset list and the size list, respectively.

Performance Evaluation

In this section, we evaluate the performance of pcapWT, in terms of index data building time and packet extraction time. PcapWT has been developed Figure 7: Retrieving the packets using multi-threaded le I/O. Multiple threads (m) evenly divide the amount of the nal matched packets (n), and each thread (T 1 , T 2 , ..., T m ) copies the assigned packets from the source to output pcap le, simultaneously. using gcc and libcds [22] on the Linux platform. In particular, pcapWT and libcds have been compiled with g++ (Ver. 4.4.7) with the highest optimization ag (-O3 ). In addition, we set the -fopenmp option of g++ for enabling OpenMP. The machine used for benchmarking consists of CentOS 6.4 (kernel Ver. 2.6.32 for 64 bits x86 system), Intel i7-2600k processor including 8 threads (4 cores), 4 GB of main memory, a 60 GB SSD (OCZ Agility 3, low-grade product), and a 500 GB HDD (Toshiba HDD2F22).

Performance benchmark of SSD and HDD

First we benchmark read and write operations on the storage devices that will be used in the experiments. As shown in Table 1, we measure the average read or write throughput with six types of tests: sequential read/write, 4KB random read/write, and 4KB random read/write using 64 threads.

As shown in Table 1, SSD and HDD exhibit similar capabilities for the sequential write. However SSD is 3.5 times faster than HDD for the sequential read. On the contrary, random read and write performance tested on both devices is degraded, compared to the sequential read/write operations. On the other hand, multi-threaded random write obtains almost the same performance as the sequential write on the SSD. Likewise, random read performance is considerably improved when using multi-threading. However, the performance drops below 1 MBps on the HDD because of the long seeking time of the HDD, which results in inecient multi-threading.

Basically, the packet extraction process is based on random reads and sequential writes from/to le. As we can observe in Table 1, multi-threaded read is 2.5 times faster than single threaded read on the SSD, even though the drive used in the benchmark is one of the most basic SSDs 1 . Thus, the random read capability substantially impacts performance of packet extraction on SSD. Nevertheless, other tools such as tcpdump and pcapIndex do not use multi-threading for random read.

Sample packet traces

PcapWT supports pcap [START_REF] Jacobson | Tcpdump[END_REF], one of the most popular packet trace formats.

Table 2 shows four sample pcap les used in the performance evaluation, such as simulated trac (ST) and real trac 1/2/3 (RT-1/2/3). ST was generated by the ns-3 network simulator [START_REF]Ns-3 Network Simulator[END_REF], and contains 15 million packets in a 4.2

In the same test made on one level higher product (OCZ-Agility 3, 240 GB) than the SSD used in this evaluation, the gap is greater than 10 times (Seq Read: 211.4 MBps, 4K Rand Read: 13.6 MBps, and 4K 64T Rand Read: 140.9 MBps). GB pcap le. RT-1 was generated by capturing packets in our local network, and also contains 15 million packets in a 3.6 GB pcap le. ST consists of 8 endpoints at the transport layer, whereas RT-1 has 2416 endpoints.

Meanwhile, those two samples are smaller than 4.2 GB because pcapIndex cannot support larger pcap les than the volume.

RT-2 and RT-3 are used to evaluate larger source trace les than 4.2 GB. Those two samples are generated by adding virtual packets to RT-1 in order to compare the results with RT-1 directly. Thus, RT-2 contains 115 million packets in a 13.4 GB pcap le, and the latter contains 200 million packets in a 21.6 GB pcap le.

As mentioned above, the packet extraction performance depends on the random read and write capabilities. Thus, we have manipulated those two samples in order to contain a larger amount of packets than ST and RT-1, instead of reducing their average packet sizes.

Other tools using packet indexing

We rst analyze the performances of other tools based on packet indexing, such as Compressed Pcap Packet Indexing Program (CPPIP) [START_REF]Cppip[END_REF], PCAP-Index [START_REF]Pcap-index[END_REF], and pcapIndex [START_REF]PcapIndex[END_REF].

These three tools use pre-built index data, but have dierent architectures. In particular, CPPIP adopts the bgzip data compression utility [START_REF] Li | The Sequence Alignment/Map Format and SAMtools[END_REF], PCAP-Index runs the SQLite Ver. 3 database, and pcapIndex uses the COM-PAX bitmap compression data structure [START_REF] Fusco | NET-FLi: On-the-y Compression, Archiving and Indexing of Streaming Network Trac[END_REF]. We compare their performance on the SSD described in Table 1, in terms of index data size and packet extraction time.

In order to build the index data from RT-1, pcapIndex and CPPIP need 27 MB (0.8% 2 ) and 289 MB (8.0%) of storage, respectively. However, PCAP-Index needs around 2.6 GB (72%), even though PCAP-index deals with a few simple elds (e.g., addresses, port numbers, protocol types, and so on). This is due to the database, which is not congured to use data compression. Moreover, PCAP-Index needs around 265 seconds to complete, whereas pcapIndex and CPPIP take only around 25 and 40 seconds, respectively. The reason is that PCAP-Index requires a lot of time to inject the index data into the database.

CPPIP and PcapIndex take around 20 and 30 seconds to extract the largest trac corresponding to 48% of the volume and 36% of the total number of the packets of RT-1, respectively. However, the available operations are not the same ones used as in the prior case, since this tool supports only simple ltering queries using a range of packet timestamps and striding (e.g., selecting every n-th packets). Moreover, PCAP-Index takes about 80 seconds (packet lookup time from the database about 30 seconds and le writing time about 50 seconds), which is signicantly slower compared to the other tools.

Even tcpdump takes less than 32 seconds in the same use case.

CPPIP aords the advantage of being able to extract packets directly from a compressed trace le. However, some eorts are needed to support more complex ltering queries and to reduce the index data size as much as that of pcapIndex. PCAP-Index is not satisfactory from many points of views, even though it uses a database. In contrast, pcapIndex is the only reliable one among those tools, since it provides not only comprehensive query operations, as well as remarkable performances in terms of index data size and packet extraction time. Therefore, we select pcapIndex to compare performance with pcapWT in the following.

Building index data

PcapWT aims to minimize the index data size, for example to compress it below 1% of the volume of original trace les. As shown in Table 3, pcapwt-index is successful for all sample les. More precisely, as mentioned in Section 3, WT supports two bitmap representation methods (i.e., RG and RRR) and its encoding unit size that can be set by the user. Figure 8 shows that RRR is slightly better at reducing a percentage of the volume of original trace le. index data size than RG. On the other hand, according to Figures 9 and10, RG is more ecient than RRR, in terms of index data building time and the packet lookup time. Especially, as shown in Figure 8, RG using a wide sampling provides mostly complementary performance between those features. Thus, we decided to use RG(50) as the default conguration.

In addition, for improving compression eciency and reducing search time at the same time, we use Human shaped WT which combines multi-layered wavelet trees and entropy compression of the bitmap sequences. Moreover, we minimize redundant bits (e.g., consecutive `0' bits) that can be generated in between multi-layered wavelet trees, using Human compression.

In contrast, pcapIndex obtains an abnormally large index data size and causes a segmentation fault when building the index data for trace les larger than 4.2 GB. In the case of RT-2, the index data size corresponds to 6.86% of the volume of the original trace le, and around 3.0 GB of main memory is required. Furthermore, this tool fails to build the index data for RT-3 because of lack of memory. Note that RT-2 and RT-3 have been generated by adding virtual packets to RT-1, in order to provide a direct comparison with the index data size and packet extraction performance of RT-1. More precisely, the virtual packets are equal except for the timestamp eld. Thus, the index data from those samples has to be minimized regardless of the number of packets. Nevertheless, the index data size produced by pcapIndex exponentially increases with the number of packets. Figure 10: A comparison of the packet lookup time for dierent sampling sizes. In case of ST, the number of the retrieval targets is 9.6 million packets (i.e., 63.8% in comparison to the total number of packets included in this sample.), and the volume of the target packets corresponds to 1.82 GB (i.e., 43.7% in comparison with the original volume). In case of RT-1/2/3, the amount of target packets are all same at 5.6 million packets and 0.97 GB, whereas the percentages of the number of the target packets are 37.5%, 3.8%, and 2.8%, and the percentages of the volumes of them correspond to 23.3%, 7.3%, and 4.5%, respectively. Figure 10, pcapWT takes less than 1 second to retrieve millions of packets.

The rest of time is consumed by le I/O, including le seeking, reading, and writing. Especially, as mentioned in Section 4.2, le seeking is one of the main reasons causing performance degradation of random le read and write operations. Thus, with pcapWT, we propose to apply multi-threaded le I/O to increase the throughput of random le read and write.

As we can observe in Figure 11, pcapIndex and pcapWT have almost the same packet extraction performance on both the SSD and the HDD when using a single thread. Meanwhile, when pcapWT extracts a large portion (88% of the total packets) from ST on the SSD, 32 threads can reduce the time down to 26%, compared to when using a single thread. However, multi-threading on HDD does not provide any gain. Even worse, the time is increased to 24% when using 32 threads. In this case, multi-threading generates frequent random accesses which increase delay and degrade performance on HDD.

Performance evaluation of packet extraction

To analyze the performance of packet retrieval, we use various queries, and extract a large amount of packets. As shown in Table 4, simple queries (from A through E) extract packets specied by a pair of the query type and its operand, and complex queries (from F through J) indicate combining simple queries by using the union operator. For example, query A extracts 4,030 packets from RT-1/2/3, corresponding to 0.03%, 0.003%, and 0.002% in comparison with the total numbers of packets, respectively. Likewise, query F (combined with A and B) extracts 525K (524,684) packets, and the total amount corresponds to 20.97%, 5.66%, and 3.49% compared to the volume of the original traces, respectively. In the same way, query I includes the former eight queries (from A through H), and contains over 12 million packets. On the other hand, query J searches packets belonging to queries B and E in common, and the extracted packets are same as the result of query B.

As shown in Figure 12, in experiments using RT-1, pcapIndex and pcapWT are faster than tcpdump in retrieving small amounts of packets, such as queries A, B, C, F, G, and J. On the contrary, the gap of packet extraction performance becomes smaller between tcpdump and those two tools using the packet indexing scheme. The reason is that the performance depends highly on the amount of packets extracted. Nevertheless, pcapWT (using 32 threads) is faster than pcapIndex in experiments performing complex queries, such as 21.8% at D, 20.3% at E, 20.3% at H, and 9.9% at I.

In experiments using larger samples (i.e., RT-2/3), we evaluate the performance of pcapWT and compare it with tcpdump, since pcapIndex cannot support such large traces. In these experiments, pcapWT uses a prompt mode which supports iterative operations once the index data is loaded. Thus, this mode can avoid unnecessary time loading the same index data repetitively.

As shown in Figure 13, tcpdump takes much longer (i.e., around 100 seconds in the case of RT-2 and 170 seconds in the case of RT-3 along the all queries) compared to the results of RT-1, even if it extracts the same amount of packets. In other words, the performance based on the sequential approach is signicantly impacted by the number of packets, regardless of the amount of packets retrieved. On the other hand, pcapWT using 32 threads is much faster than tcpdump, and the time increases proportionally to the amount of matched packets, unlike the sequential approach.

Conclusion

In network analysis, packet traces have always been important to analyze, since they record the complete information exchanged through networks. However, as network links, trac complexity, and applications throughput have signicantly increased, the analysis of huge amounts of those traces is becoming a challenging task, and traditional analysis tools are not ecient at dealing with such traces. Even though a few tools have been recently published for fast packet extraction, their performance is still not satisfactory in terms of index data size, and packet extraction delay.

In this paper, we propose a new tool designed to process eciently very large network traces, which allows not only enhancing packet extraction performance on large traces, but also reducing storage requirements. PcapWT uses a Wavelet Tree data structure and a multi-threading scheme for parallel le I/O. Our benchmarks including SSD show that pcapWT is about 10 to 20% faster than pcapIndex in the worst case scenario and about 200% in the

Figure 1 :

 1 Figure 1: On the left the regular wavelet tree, on the right the Human shaped wavelet tree, for the sequence S = 1511863875743288 , σ = 8 and n = 16.

Figure 2 :

 2 Figure 2: Extracting index data from a source packet trace le, such as packet size, Ethernet Src/Dst address, IPv4 Src/Dst address, Src/Dst port number, and Layer 4 protocol type.

Figure 3 :

 3 Figure 3: Mapping between the extracted index data and positive integers in consecutive order for reducing the index data size and packet lookup time.

Figure 4 :

 4 Figure 4: Building WT chunks from the mapped index data once every one million packets. This partitioning can reduce the size of index data.

Figure 5 :

 5 Figure 5: Searching matched packets from the index data.

Figure 6 :

 6 Figure 6: Listing up the nal matched packets marked by `1' at the matching point table.Then the oset position and the size of selected packet are written into the oset list and the size list, respectively.

Figure 8 :

 8 Figure 8: A comparison of the index data size by sampling methods (RG and RRR) and sizes (5, 20, 35, 50).

Figure 9 :

 9 Figure 9: A comparison of the building time of index data by sampling methods (RG and RRR) and sizes (5, 20, 35, 50).

Figure 11 :

 11 Figure 11: A performance comparison of packet extraction with various number of threads for parallel le I/O. In this test, pcapWT extracts two dierent size of tracs from ST: one small trac containing 1.8 million packets (12.2%) corresponding to 146 MB and one large trac containing 13.2 million packets (88.0%) corresponding to 3.0 GB. Note, the percentages indicate the ratio of the extracted packets compared to the total number of packets included in the sample.

Figure 12 :

 12 Figure 12: A performance comparison of packet extraction using RT-1.

Figure 13 :

 13 Figure 13: A performance comparison of packet extraction using RT-2/3.

Table 1 :

 1 Sequential (Seq), random (Rand), multi-threaded by 64 threads (64T) random I/O performance measurements of the SSD and the HDD used in the experimentation.

		SSD	HDD
		[MBps]	[MBps]
	Seq Read	198.8	56.2
	Seq Write	52.7	52.0
	4 KB Rand Read	17.7	0.4
	4 KB Rand Write	28.0	0.9
	4 KB 64T Rand Read	47.0	0.8
	4 KB 64T Rand Write	51.0	0.8
		[ms]	[ms]
	Read Access Time	0.18	17.39
	Write Access Time	0.31	4.27

Table 2 :

 2 Description of four sample pcap les: simulated trac (ST) and real trac 1/2/3/ (RT-1/2/3).

		ST	RT-1 RT-2 RT-3
	Number of packets [million]	15	15	115	200
	File size [GB]	4.2	3.6	13.4	21.6
	User data size [GB]	3.9	3.4	11.5	18.4
	Average packet size [Byte] 261.7 224.1 100.1	92.2
	Number of end-points	8 2416 2417 2417

Table 3 :

 3 A comparison of the index data size produced by pcapIndex and pcapWT with RG(50). The percentage indicates the ratio between the index data size and the volume of the original trace le.

		pcapIndex	pcapWT
		[%]	[%]
	ST	0.45	0.30
	RT-1	0.78	0.86
	RT-2	6.86	0.93
	RT-3	Error	0.96
	5.5. Performance evaluation of multi-threaded le I/O
	The packet lookup time required to nd matched items from index data
	takes a small portion of the overall packet extraction time. As shown in

Table 4 :

 4 Packet ltering queries and specications of RT-1/2/3. The percentage values indicate (1) the ratio of the number of matched packets and (2) the ratio of the size of matched packets compared to the packet number and the volume of the original trace le, respectively. Note that all numbers are rounded up. {G ∪ D} 7,317 48.8 4.88 3.66 2024 56.21 15.16 9.35 I = {H ∪ E} 12,225 81.5 8.15 6.11 3000 83.31 22.47 13.86 J = {B ∩ E}

	Query (IP Addr)	Matched packets [×10 3 ] [%] (1) (I) (II) (III)	Total packet size [%] (2) [MB] (I) (II) (III)
	A	4 0.03 0.00 0.00	0.3	0.01	0.00	0.00
	B	521	3.5 0.35 0.26	754 20.94	5.65	3.48
	C	1,161	7.7 0.77 0.58	298	8.28	2.23	1.38
	D	5,631 37.5 3.75 2.82	971 26.96	7.27	4.49
	E	5,428 36.2 3.62 2.71	1731 48.07 12.97	8.00
	F = {A ∪ B}	524	3.5 0.35 0.26	755 20.97	5.66	3.49
	G = {F ∪ C}	1,686 11.2 1.12 0.84	1053 29.24	7.89	4.87
	H = 521	3.5 0.35 0.26	754 20.94	5.65	3.48

The design of pcapWTIn this section, we describe how to build the index data, and explain the process for lookup querying and for extracting packets. The overall processes consist of six steps, as shown in Figures2 through 7.4.1. Building index dataIn this work, the generated index is a data set that improves the speed of packet lookup operations on network traces. The index data is used to quickly locate packets specied by user query without the need to search for

Acknowledgement

The authors would like to thank Dr. Renato Cerro for his great help to improve the quality of the paper. This study was supported by Anillo Project ACT-53, the Fondecyt project No. 11121475, CIRIC (INRIA-Chile) Project "Network Design", Project Semilla -UDP "ANDES" and "Anàlisis y diseño de algoritmos en redes de bajo consumo aplicado a condiciones extremas de los Andes": Programa de Cooperaciòn Cientìca Internacional CONICYT/MINCYT 2011.