Diala Bitar
email: diala.bitar@femto-st.fr

Najib Kacem
email: najib.kacem@femto-st.fr

Noureddine Bouhaddi
email: noureddine.bouhaddi@femto-st.fr

MULTISTABILITY AND MODAL INTERACTIONS IN PERIODIC 2D COUPLED PENDULUMS ARRAY

The collective dynamics of an array of periodic two dimensional (2D) coupled pendulums under harmonic horizontal base excitation is investigated. The coupled differential equations governing the nonlinear vibrations of the considered system have been solved using an analytical-numerical solving procedure, based on the multiple scales method coupled with standing wave decomposition. It allows the identification of complex and wide variety of nonlinear phenomenon exhibited by the periodic nonlinear structure. The frequency responses for several coupled pendulums were calculated in order to analyze the stability, the modal interactions and the bifurcation topologies resulting from the collective dynamics of the coupled pendulums, while highlighting the large number of multimodal solutions for a small number of coupled pendulums. The complexity and the multivaludness of the responses were illustrated by a study of basins of attraction which display the large distribution of the multi-mode branches.

INTRODUCTION

Recent researches focused on the study of wave propagation in nonlinear periodic media seen the existence of complex wave phenomena such as multi-modal solutions, modal interactions, solitary waves and discrete breathers in addition to those exhibited by linear periodic structures. A large number of phys-ical systems in different scientific areas such as materials, optics, chemistry and acoustics can be modeled as one-dimensional (1D) periodic nonlinear chain. In addition, the localization phenomenon generated by the interaction between the nonlinearities and a discrete system has attracted a lot of attention in recent years, because of its important role in the qualification as well as quantification of systems operations. One of the most popular localization phenomena, that have attracted the interest of physicists, is the nonlinear energy localization, which are called intrinsic localized modes ILMs. For instance, Sivers et al. [START_REF] Sievers | Intrinsic localized modes in anharmonic crystals[END_REF] proved theoretically that ILMs could occur in pure anharmonic crystals. Recently, in the field of engineering, such as microelectromechanical systems (MEMS), ILMs were observed when micromechanical oscillators arrays were excited at high frequencies [START_REF] Sato | Study of intrinsic localized vibrational modes in micromechanical oscillator arrays[END_REF], and their characteristics were theoretically investigated [START_REF] Aj | Localization in mi-FIGURE 5[END_REF].

In mechanics, Several researches focused on ILMs in pendulums array, we can mention Ikeda et al. [START_REF] Ikeda | Intrinsic localized modes of principal parametric resonances in pendulum arrays subjected to vertical excitation[END_REF] that investigated the behavior of ILMs for an array of coupled pendulums subjected to horizontal or vertical excitations. In addition, modulation instability of the uniform mode in a pendulum array with harmonic coupling and horizontal, sinusoidal driving which at longer times makes possible the formation of stable ILMs has been observed [START_REF] Thakur | Driven intrinsic localized modes in a coupled pendulum array[END_REF]. Moreover, an experiment using a coupledpendulum chain has been reported to demonstrate various aspects of wave motion [START_REF] Liang | An undergraduate experiment of wave motion using a coupled-pendulum chain[END_REF]. Recently, The influence of adding external harmonic excitation on the intrinsic localized modes of coupled pendulums chains has been investigated [START_REF] Jallouli | Stabilization of solitons in coupled nonlinear pendulums with simultaneous external and parametric excitations[END_REF].

Taking into account the considerable work in the literature devoted to analyze the richness of the nonlinear dynamical effects found in an array of coupled pendulums [START_REF] Vyas | Dynamics of autoparametric vibration absorbers using multiple pendulums[END_REF], 2D systems are expected to present additional dynamics phenomena not realizable in the 1D case. Recently, the first analysis of a 2D array of coupled pendulums proposed a computational model to study the nonlinear dynamics of such system under parametric excitation [START_REF] Jallouli | Nonlinear dynamics of a 2d array of coupled pendulums under parametric excitation[END_REF]. Although the dynamics of coupled nonlinear pendulums was thoroughly investigated in the frequency and time-space domains, there is a real need to perform profound analysis of the collective dynamics of such systems in order to identify practical relations with the nonlinear energy localization phenomenon in terms of modal interactions and bifurcation topologies [START_REF] Bitar | Collective dynamics of periodic nonlinear oscillators under simultaneous parametric and external excitations[END_REF].

In this work, we investigate the collective dynamics of periodic 2D coupled pendulums array under external excitation. The coupled equations governing the nonlinear vibrations of the considered system have been solved using an analytical-numerical solving procedure, which allows a detailed study of complex model interactions of the periodic structure. The multiple scales method coupled with standing wave decomposition were used to transform the nonlinear system into a set of complex algebraic equations. The responses frequency for several coupled pendulums were plotted in the frequency domain in order to analyze the stability, the modal interactions and the bifurcation topologies resulting from the collective dynamics of the coupled pendulums, while highlighting the large number of multimodal solutions for a small number of coupled pendulums. The complexity and the multivaludness of the responses were illustrated by a detailed study of basins of attraction which display the large distribution of the multi-mode branches.

THEORETICAL ANALYSIS Problem Formulation

Figure 1 shows the discrete model of 2D coupled N × M -pendulums subjected to harmonic horizontal base excitation. The system is composed of N equidistant axles A n . Each axle has M equidistant pendulum "i, j" attached perpendicularly to the axle and designed as a uniform rod of length l and mass m fixed at the end. φ i, j is its rotational degree of freedom (dof) from its stable equilibrium position. Pendula "i, j"and "i + 1, j", "i, j" and "i, j + 1" are connected respectively by torsional and linear springs k x and k y . The periodic system is subjected to harmonic horizontal base excitation y = y 0 cos(Ωτ), with fixed boundary conditions φ i,M+1 = φ i,0 = φ N+1, j = φ 0, j = 0.

The kinetic energy T of the system can be expressed as Where (y i, j , z i, j ) represents the Cartesian coordinates of the pendulum {i, j}, given as

T = 1 2 N ∑ i=1 m( ẏ2 i, j + ż2 i, j ). (1)
y i, j = y + l sin φ i, j , z i, j = l cos φ i, j , (2)
and the dot denotes typically the derivative with respect to time. Substituting Eq. (2) into Eq. ( 1) the kinetic energy becomes

T = 1 2 N ∑ i=1 m l 2 φ 2 i, j + 2l cos(φ i, j ) ẏ φi, j + ẏ2 . (3)
The potential energy of such a problem can be written as:

V = N-1 ∑ i=1 M-1 ∑ j=1 1 2 k x (φ i+1, j -φ i, j ) 2 + (φ i-1, j -φ i, j ) 2 + 1 2 k y ∆ 2 j, j+1 + ∆ 2 j-1, j + mg l -l cos(φ i, j ) , (4)
where g represents the acceleration of gravity and ∆ j, j+1 the distance between pendulums {i, j} and {i, j + 1}, which can be written as

∆ j, j+1 = l(sin(φ i, j+1 ) -sin(φ i, j )).
(

) 5
The Lagrange equations of motion can be written in the following form

d dτ ( ∂ L ∂ φi, j ) - ∂ L ∂ φ i, j = Q i, j , (6)
with L the Lagrangian operator defined by L = T -V , Q i, j = αl φi, j represents the generalized forces applied to the pendulum {i, j} and α the viscous damping coefficient.

Substituting Equations ( 3) and ( 4) into the Lagrange Equation [START_REF] Liang | An undergraduate experiment of wave motion using a coupled-pendulum chain[END_REF], we obtain the equation of motion of each pendulum {i, j}

ml 2 φi, j + αl φi, j + k x (2φ i, j -φ i-1, j -φ i+1, j ) + mgl sin(φ i, j ) +k y l 2 cos(φ i, j )(2 sin(φ i, j ) -sin(φ i, j-1 ) -sin(φ i, j+1 )) = mly 0 Ω 2 cos(φ i, j ) cos(Ωτ) (7)
By expanding sin(φ i, j ) and cos(φ i, j ) into Taylor series up to the third and first order respectively, the equations ( 7) can be written as:

ml 2 φi, j + αl φi, j + k x (2φ i, j -φ i-1, j -φ i+1, j ) + mgl(φ i, j -1 6 φ 3 i, j ) +k y l 2 (2φ i, j -φ i, j-1 -φ i, j+1 ) -1 6 k y l 2 (2φ 3 i, j -φ 3 i, j-1 -φ 3 i, j+1 ) = mly 0 Ω 2 cos(Ωτ) (8)
In the following sections, and since we are interested in localized modes phenomenon, the third order expansion related to the linear coupling will be neglected.

Dimensionless equations

For convenience and equation simplicity, the following nondimensional variables are introduced:

θ i, j = φ i, j φ D , t = ω 0 τ, (9)
where φ D = my 0 Ω 2 αω 0 is the dynamic displacement of the associated linear system, Q = ω 0 ml α and ω 0 = g l represents the natural frequency of pendulum "i, j". After dropping Eqs. [START_REF] Jallouli | Nonlinear dynamics of a 2d array of coupled pendulums under parametric excitation[END_REF] into Eq. ( 8) and dividing by α m 2 l 2 y 0 ω 0 Ω 2 , we obtain the following nondimensional system of equations

θi, j + 1 Q θi, j + θ i, j + k x m g l (2θ i, j -θ i-1, j -θ i+1, j ) + k y l m g (2θ i, j -θ i, j-1 -θ i, j+1 ) -1 6 ( m y 0 Ω 2 α ω 0 ) 2 θ 3 i, j = α m l ω 0 cos( Ω ω 0 t) (10)
The aim of this paper is to present the influence of the modes localization on the collective dynamics of a 2D pendulums array. Therefore, a study on the natural frequency is required.

Linear study

Equation [START_REF] Bitar | Collective dynamics of periodic nonlinear oscillators under simultaneous parametric and external excitations[END_REF] describes a nonlinear Duffing system coupled by linear springs k x and k y in the x and y directions respectively and subjected to harmonic horizontal base excitation. To determine the natural frequencies and their eigenvectors Equation [START_REF] Bitar | Collective dynamics of periodic nonlinear oscillators under simultaneous parametric and external excitations[END_REF] can be written in matrix form as

M Θ + C Θ + K L Θ Linear part +F NL (θ ) Nonlinear part = F(t), (11)
with

Θ = [θ 1,1 , θ 1,2 , . . . , θ 1,M , . . . θ i,1 , θ i,2 , . . . , θ i,M , .
. . ] the N×M-dimensional variable vector, the excitation vector

F(t) = α m l ω 0 cos( Ω ω 0 t)[1, . . . , 1] T , M = I NM , C = 1 Q * diag(1, 1, . . . , 1
) , F NL (θ ) is the nonlinear stiffness vector and

K L =          K 1 -K 2 -K 2 K 1 -K 2 0 . . . . . . . . . 0 -K 2 K 1 -K 2 -K 2 K 1         
, where

K 1 =            1 + 2k x mgl + 2k y l mg - k y l mg - k y l mg 1 + 2k x mgl + 2k y l mg - k y l mg 0 . . . . . . . . . 0 - k y l mg 1 + 2k x mgl + 2k y l mg - k y l mg - k y l mg 1 + 2k x mgl + 2k y l mg            and K 2 = k x mgl I M .
It is convenient to express dimensionless eigenvalues ω I, j in the following form. For example,

N = 1, M = 1 : ω 1,1 = 1 + 2k x mgl + 2k y l mg N = 2, M = 1 : ω 1,1 = 1 + k x mgl + 2k y l mg ω 2,1 = 1 + 3k x mgl + 2k y l mg (12)
It is Convenient to express all natural frequencies relative to the same dimensional reference frequency which is 1, so that ω i, j = (1 + λ i ∆ x + λ j ∆ y ) (i = 1, . . . , N and j = 1, . . . , M), [START_REF] Lifshitz | Response of parametrically driven nonlinear coupled oscillators with application to micromechanical and nanomechanical resonator arrays[END_REF] where λ i and λ j are the coefficient depending on the mode i, j of vibration, ∆ x = k x mgl and ∆ y = k y l mg . We suppose that k x << mgl and k y << mg l . Consequently, ∆ x and ∆ y << 1 and:

ω i, j ≈ 1 + 1 2 λ i ∆ x + 1 2 λ j ∆ y (i = 1, . . . , N and j = 1, . . . , M). ( 14
)
This assumption leads to the creation of linear closed modes which permits to study the effect of the mode localization on the collective dynamics.

Semi-analytical approach

Several methods have been developed to solve similar nonlinear differential systems [START_REF] Bitar | Collective dynamics of periodic nonlinear oscillators under simultaneous parametric and external excitations[END_REF]. However, these methods are time consuming and not-suitable to solve a high number of nonlinear differential equations. In this paper, we use the multiple scales method [START_REF] Nayfeh | Introduction to Perturbation Techniques[END_REF] as an approximate analytical solution to solve nonlinear differential systems. The main advantages of the present approach are the capacity of handling weakly nonlinear systems, which permits to visualize all physical responses branches and their properties in terms of modal interactions and bifurcation topology transfer.

We introduce the parameters listed in Equation ( 15) and set the external frequency an amount εω 0 Ω D away from the resonant frequency, whereby they contribute to the equations of amplitude.

Q = εη, k x mgl = 1 2 εγ x , k y l mg = 1 2 εγ y , -1 6 ( m y 0 ω 2 α ω 0 ) 2 = ξ , α m l ω 0 = ε 3 2 f . (15)
We calculate the responses of the pendulums array to harmonic, horizontal excitation using secular perturbation theory, by expanding θ i, j (t) of a sum of standing-wave modes with slowly varying amplitude [START_REF] Lifshitz | Nonlinear Dynamics of Nanomechanical and Micromechanical Resonators[END_REF].

θ i, j (t) = ε 1/2 ∑ N r=1 ∑ N p=1 (A r,p (
T ) sin(nq r ) sin(md p )e iω r,p t + c.c.)

+ε 3/2 θ (1) i, j (t) + • • • , i = 1, ..., N and j = 1, ..., M (16)
where T = εt is a slow time variable, that authorizes the complex amplitude A r,p (T ) to vary slowly in time. Since we proposed fixed boundary conditions, the possible wave components q r , q p can be given as

q r = rπ N+1 , r = 1, ..., N d p = pπ M+1 , p = 1, ..., M (17)
After replacing the proposed solution ( 16) into the equation of motion, we can get at the order of ε 3 2 , N × M equations of the form:

θ (1) i, j + θ (1)
i, j = ∑ r ∑ p (r, p th secular term)e iω r,p t + other terms

We must eliminate the N × M secular terms so that θ i, j (1) remain finite which allow us to determine the equations for the slowly varying amplitudes A r,p (T ). To extract the equations for these amplitudes, we make use of the orthogonality of the modes, by multiplying the r, p th secular term by sin(nq r ) sin(md p ) and summing over n and m. We also express all normal frequencies relative to 1, so that:

ω r,p = 1 + εΩ r,p (18)
We find that the equation of the r, p th amplitude A r,p (T ) is given by: 2iω r,p dA r,p dT + iηω r,p A r,p + 2(γ x sin 2 ( q r 2 ) + γ y sin 2 (

d p 2 ))A r,p + 3 16 ξ ∑ i,k,l ∑ j,s,o A i, j A k,s A * l,o ∆ (1) ikl;r ∆ (1) jso;p e i(Ω i, j +Ω k,s -Ω l,o -Ω r,p )T = 2 (N+1)(M+1) f ∑ n ∑ m sin (nq r ) sin (md p )e i(Ω D -Ω r,p )T (19)
In solving Equations (19), we write the periodic steady state solution of the form A r,p = a r,p e i[Ω D -Ω r,p ]T (20)

Substituting Equation ( 20) into (19), we obtain the required equation for the complex amplitudes a r,p .

TABLE 1. DESIGN PARAMETERS FOR THE CORRESPONDING PERIODIC STRUCTURE DEPICTED IN FIGURE 1

Parameters

Design 1 Design 2

M (Kg) 25.10 -2 - g (m.s -2 ) 9.81 - l (m) 62.10 -3 - α (Kg.m.s -1 ) 10 -2 - y 0 (m.s -2 ) 5.10 -4 g - k x (N.m) 10 -5 5.10 -3 k y (N.m) k x /l 2 = 26.10 -4 10 -2
2iω r,p (Ω D -Ω r,p )a r,p + iηω r,p a r,p +2(γ x sin( q r 2 ) 2 + γ y sin(

d p 2 ) 2 )a r,p + 3 16 ξ ∑ i,k,l ∑ j,s,o a i, j a k,s a * l,o ∆ (1) ikl;r ∆ (1) jso;p = 2
(N+1)(M+1) f ∑ n ∑ m sin (nq r ) sin (md p ) (21)
This solving technique has been used by Lifshitz et al. [START_REF] Lifshitz | Response of parametrically driven nonlinear coupled oscillators with application to micromechanical and nanomechanical resonator arrays[END_REF] to calculate the response of nonlinear coupled microelectromechanical oscillators under parametric excitation and has been validated with respect to the experimental measurements of Buks and Roukes [START_REF] Buks | Electrically tunable collective response in a coupled micromechanical array[END_REF]. The complexity of the nonlinear algebraic System (21), limits the possibility to solve it analytically for a high number of coupled pendulums. We may use the appropriate time integration procedure in order to solve the differential System (19) numerically, which implements a version of the fourth-order Runge Kutta integration method.

RESULTS AND DISCUSSIONS

In order to study the collective dynamics of the 2D periodic structure, two different configurations were proposed in Table 1, which highlights the influence of linear coupling on the non-linear behavior in terms of frequency range and vibrational amplitude.

For the case of 2 × 2 coupled pendulums (N = 2, M = 2) and for the first design parameters listed in Table 1, the natural dimensionless Eigen frequencies of the associated linear system are very close to 1. In addition, as we considered k y = k x l 2 , the linear springs are expected to act the same way. Figures 3 show the response intensities as a function of the detuning parameter Ω D . The Single Mode (SM) solution branches corresponds to the Trivial solution of the System (19) with A 1,2 = A 2,1 = A 2,2 = 0. The multimodal solutions DM and QM are generated by the coupling between two and all modes respectively. These curves were ploted to underline: first the effect of modes localization on the high number of solution branches for a given frequency and second the identical contribution of the non-excited modes on the responses. Although, the advantages of modes localization, considering small coupling parameters does not affect the systems rotational degree of freedom (dof). Figure 4, displays the rotational dof of 2 × 2 coupled pendulums under harmonic horizontal base excitation as a function of frequency, which corresponds to a single forced Duffing oscillator. In addition, all multimodal solution branches are distributed between resonant and non-resonant branch, hence its complicated to identify their contribution.

We decide to increase linear coupling parameters in such a manner to preserve the localization phenomena. Therefore, we consider the second design parameters listed in Table 1, so that the natural dimensionless Eigen frequencies of the associated linear system are ω 1,1 = 1.0331, ω 1,2 = 1.0336, ω 2,1 = 1.0989 and ω 2,2 = 1.0994. Figure 5 shows the response intensities of 2 × 2 coupled pendulums. Unlike to the first case, linear springs k x and k y play different roles. In particular, the multimodal solution branches resulting from the interaction between A 1,1 and A 2,1 or A 2,2 present larger frequency range and their bifurcation point B 2 exceeds the one that corresponds to the non-resonant branch of the SM (B 1 ). In addition, in regards to the systems rotational dof, multimodal solution branches present additional features distinct from those presented by a single forced Duffing oscillator (larger bandwidth and higher vibrational amplitude).

Taking into account the high number of solution branches founded above, for a given frequency, we choose to explore the practical appearance of each attractor according to their modal contribution type. Figure 7, displays the distribution of the basins of attraction of the multimode solution branches in the Nyquist plane (Re(a 11 ), Im(a 11 )) of the case of 2 × 3 coupled pendulums under harmonic horizontal base excitation, for the second design parameters listed in Table 1. For the configuration ω = 11.63, and random initial conditions Re(a 12 ) = 0.047, Re(a 13 ) = 0.2447, Re(a 21 ) = 0.1868, Re(a 22 ) = -0.1315, Re(a 23 ) = 0.2802, Im(a 12 ) = -0.2037, Im(a 13 ) = -0.3110, Im(a 21 ) = -0.3165, Im(a 22 ) = 0.1256 and Im(a 23 ) = -0.4189. The solutions DM, T M, QM, 5 th M and 6 th M are related to the branches generated by the interaction between two, three, four, five or all modes respectively. This distribution highlights the robustness of the multimodal solutions and the high contribution of the modes into the responses.

CONCLUSION

The collective nonlinear dynamics of a 2D periodic array of coupled pendulums was modeled under harmonic horizontal base excitation. An analytico-numerical method suitable for nonlinear localized modes problems has been considered, based on a perturbation technique, combined with standing wave decomposition, transforming the normalized system into a set of coupled complex differential equations. The effects of modes localization on the collective nonlinear dynamics have been studied for the case of 2 × 2 coupled pendulums. The resulting complex system has been numerically solved for two different configurations, highlighting the impact of linear coupling parameters on the complexity and the multivaludness of the responses. The robustness of these additional branches was illustrated by a study of the basins of attraction for 2 × 3 coupled pendulums. In practice, this model can be adjusted to increase the number of multimodal solutions, their configuration and their vibrational amplitudes. A complete quantitative study is a veritable challenge that we will attempt in a future work.

FIGURE 1 .

 1 FIGURE 1. 2D ARRAY OF COUPLED PENDULUMS SUB-JECTED TO HORIZONTAL HARMONIC BASE EXCITATION.

FIGURE 2 .

 2 FIGURE 2. FRONT VIEW OF THE SYSTEM.

1

 1

FIGURE 3 .

 3 FIGURE 3. RESPONSE INTENSITY OF 2 × 2 COUPLED PENDULUMS AS A FUNCTION OF FREQUENCY, UNDER harmonic horizontal base excitation FOR THE FIRST SET OF DESIGN PARAMETERS LISTED IN TABLE 1.

FIGURE 6 .

 6 FIGURE 6. ROTATIONAL DOF AS A FUNCTION OF FRE-QUENCY OF 2×2 COUPLED PENDULUMS, UNDER harmonic horizontal base excitation FOR THE SECOND SET OF DESIGN PARAM-ETERS LISTED IN TABLE 1.

FIGURE 7 .

 7 FIGURE 7. DISTRIBUTION OF THE BASINS OF ATTRACTION FOR 2 × 3 COUPLED PENDULUMS FOR Ω D = -0.075 IN THE NYQUIST PLANE, FOR THE SECOND SET OF DESIGN PARAM-ETERS LISTED IN TABLE 1.

ACKNOWLEDGMENT

This project has been performed in cooperation with the Labex AC-TION program (contract ANR-11-LABX-01-01).


Stéphane Puechmorel

Daniel Delahaye

John Hansman

Jonathan Histon

AIR TRAFFIC COMPLEXITY MAP BASED ON NON LINEAR DYNAMICAL SYSTEMS

This paper presents a new air traffic complexity metric based on non-linear dynamical systems. Previous work has shown that the structure and organization of traffic are important factors in the perception of the complexity of an air traffic situation. The new metric captures these important aspects of complexity by identifying the organization of trajectories in a traffic pattern. This paper investigates only the features of this new metric without quantifying directly the connection between complexity and the metric. Authors of previous work in this area have proposed metrics that generally have not explicitly addressed the effects of organization in the traffic flow on complexity. In order to capture the effect of organization, the metric is based on a dynamical system which fits as closely as possible the observations given by the aircraft positions and speeds. Two approaches are presented. The first one is based on a linear dynamical system and produces an aggregate complexity metric. The second approach, uses a non-linear dynamical system model that fits the observations without error. This metric can be used to identify high (or low) complexity areas on a map, and, by capturing the organization properties of the traffic, captures some of the key factors involved in ATC complexity. A complexity map for the northwest area of France is shown as an example of the application of the model to real radar data. Such maps are an example of the usefulness of these methods for comparing the relative complexity of different regions of airspace.

INTRODUCTION

Currently, about 7500 flights occur everyday over France, a crossroad for European airspace. This traffic generates a large amount of controller workload. In order to safely handle this load, airspace is divided into sectors, managed by air traffic controllers. For several years, a constant increase of air traffic has induced more and more congestion in the control sectors. Two strategies can be used to reduce congestion. The first one consists in adapting the demand to the existing capacity (slot-route allocation, collaborative decision making etc...). The second one adapts the capacity to the demand (modification of the air network, new design of the sectorisation, new airports, new policies, etc.).

Currently, the capacity of a sector is measured as the number of aircraft flying across the sector during a given period of time. However, observations indicate that at times controllers will accept planes beyond the capacity threshold, while in other situations they refuse traffic before the maximum capacity is reached. This phenomenon indicates that the simple measure of traffic rate does not fully account for a controller's workload. In order to better quantify the congestion in a sector, this study synthesizes a traffic complexity metric that captures the effects of the traffic configuration and is thus thought to be more relevant than simply the number of aircraft in the sector.. Most complexity metrics that have been proposed to date have focused only on the speed vector distribution with the associated disorder metric capturing only some features of the traffic complexity. The objective of our work is to build a metric that captures the intrinsic organization of a set of trajectories in a 4D space (3D for the space and 1D for the time).

Such a metric has many applications in the air traffic management area. For instance, when a sector is designed [START_REF] Delahaye | Genetic algotithms for partitioning airspace[END_REF], a more sophisticated metric would help create airspace that achieves an appropriate distribution of controller workloads. Complexity metrics can also be applied to the traffic assignment [START_REF] Delahaye | Genetic algotithms for air traffic assignment[END_REF][START_REF] Delahaye | Aispace congestion smoothing by stochastic optimization[END_REF] problem for which an optimal time of departure and a route are searched for each aircraft in order to reduce the congestion in the airspace. A complexity metric may also be used to design new air networks, to support dynamic sectorization, and to define future ATM concepts (Free Flight). Complexity metrics may also be used to qualify and quantify the performance of Air Traffic service providers and enable a more objective consultation between airlines and providers.

The work presented in this paper is based on a dynamical systems modeling of an air traffic situation. A dynamical system describes the evolution of a given state vector. If such a vector is given by the position of aircraft  X =[ x,y,z ] T , a dynamical system associates a speed vector

 Ẋ = [ v x ,v y ,v z ]
T to each point in the airspace. The key idea is to find a dynamical system that models the observed aircraft trajectories. Based on this dynamical system modeling, a trajectory disorder metric can be computed.

The first part of the paper summarizes the previous related works. The second part presents a linear dynamical system modeling for which the complexity metric can be represented as a complex coordinate system. This system allows easy identification of any organization in the distribution of the speed vectors of aircraft in the air traffic situation. The third part introduces a non linear extension of the previous dynamical system modeling that can be applied to a set of radar observations to model any observed traffic. Such a non linear modeling can be used to produce maps of traffic complexity by identifying areas with high (or low) complexity. This non linear modeling can be extended with time and can be applied directly to trajectory segments instead of speed vectors. The fifth part presents results obtained from simulated traffic and gives an example of a traffic complexity map for the Northwest of France.

PREVIOUS RELATED WORK

Several reviews of complexity in ATC have been completed; these reviews and the literature in general, have recognized that airspace complexity is related to both the structure of the traffic and the geometry of the airspace [START_REF] Mogford | The Complexity Construct in Air Traffic Control: A Review and Synthesis of the Literature[END_REF][START_REF] Hilburn | Cognitive Complexity in Air Traffic Control: A Literature Review[END_REF]. However, the relationship between the organization of traffic generated by structure and complexity does not appear to have been captured in the many metrics that have been proposed to date.

Significant research interest in the concept of ATC complexity was generated by the "Free Flight" operational concept. Integral to Free Flight was the notion of dynamic density. Conceptually, dynamic density is a measure of ATC complexity that would be used to define situations that were so complex that centralized control was required [START_REF]Final report of rtca task force 3: Free flight implementation[END_REF]. [START_REF] Laudeman | Dynamic density : an air traffic management metric[END_REF] developed a metric called "Dynamic Density" that is based on the traffic characteristics of the airspace. The "Dynamic Density" is a weighted sum of the traffic density (number of aircraft), the number of heading changes (> 15 degrees), the number of speed changes (>0.02 Mach), the number of altitude changes (>750 ft), the number of aircraft with 3-D Euclidean distance between 0-25 nautical miles, the number of conflicts predicted in 25-40 nautical miles. These factors are summed together using weighting factors that were determined by showing different traffic scenarios to several controllers. B.Sridhar from NASA [START_REF] Sridhar | Airspace complexity and its application in air traffic management[END_REF], has developed a model to predict the evolution of such a metric in the near future. Efforts to define "Dynamic Density" have identified the importance of a wide range of potential complexity factors, including structural considerations.

Laudeman et al. from NASA

However, the instantaneous position and speeds of the traffic itself does not appear to be enough to describe the total complexity associated with an airspace. A few previous studies have attempted to include structural consideration in complexity metrics, but have done so only to a restricted degree. For example, the Wyndemere Corporation proposed a metric that included a term based on the relationship between aircraft headings and dominant geometric axis in a sector [START_REF] Inc | An evaluation of air traffic control complexity[END_REF]. The importance of including structural consideration has been explicitly identified in work at Eurocontrol. In a study to identify complexity factors using judgment analysis, "Airspace Design" was identified as the second most important factor behind traffic volume [START_REF] Kirwan | Investigating complexity factors in u.k. air traffic management[END_REF]. Histon,et,al. [9,[START_REF] Histon | Structural consideration and cognitive complexity in air traffic control[END_REF] investigated how this structure can be used to support structure-based abstractions that controllers appear to use to simplify traffic situations.

The previous models do not take into account the intrinsic traffic disorder thought to be a key driver of complexity. The first efforts addressing the effect of disorder can be found in [START_REF] Delahaye | Air traffic complexity : Towards intrinsic metrics[END_REF] which introduced two classes of metrics that measure the disorder of a traffic pattern. The first class is based on geometrical properties and proposes new metrics which are able to extract features on the traffic complexity such as proximity (measures the level of aggregation of aircraft in the airspace), convergence (the rate of change of the relative distance between aircraft) and sensitivity (this metric measures how sensitive the relative distance between aircraft is to a control maneuver). The second class is based on a dynamic system modeling of the air traffic and uses the topological entropy as a measure of disorder of the traffic pattern.

G.Aigoin has extended and refined the geometrical class by using a cluster based analysis [START_REF] Aigoin | Air traffic complexity modeling[END_REF]. Two aircraft are said to be in the same cluster if the product of their relative speed and their proximity (a function of the inverse of the relative distance) is above a threshold. For each cluster, a metric of relative dependence between aircraft is computed and the whole complexity of the cluster is then given by a weighted sum of the matrix norm. Those norms give an aggregated measure of the level of proximity of aircraft in clusters and the associated convergence. From the cluster matrix, it is also possible to compute the difficulty of a cluster (it measures how hard it is to solve a cluster). Multiple clusters can exist within a sector, and their interactions must also be taken into account. A measure of this interaction has been proposed by G.Aigoin [START_REF] Aigoin | Air traffic complexity modeling[END_REF]. This technique allows multiple metrics of complexity to be developed such as average complexity, maximum and minimum cluster complexities, and complexity speeds.

Another approach based on fractal dimension has been proposed by S.Mondoloni and D, Liang in [START_REF] Mondoloni | Airspace fractal dimension and applications[END_REF]. Fractal dimension is a metric comparing traffic configurations resulting from various operational concepts. It allows in particular to separate the complexity due to sectorization from the complexity due to traffic flow features. The dimension of geometrical figures is well-known: a line is of dimension 1, a rectangle of dimension 2, and so on. Fractal dimension is simply the extension of this concept to more complicated figures, whose dimension may not be an integer. The block count approach is a practical way of computing fractal dimensions: it consists in describing a given geometrical entity in a volume divided into blocks of linear dimension d and counting the number of blocks contained in the entity N. The fractal dimension D0 of the entity is thus :

log N log d ¿ Equation 1
The application of this concept to air route analysis consists in computing the fractal dimension of the geometrical figure composed of existing air routes. An analogy of air traffic with gas dynamics then shows a relation between fractal dimension and conflict rate (number of conflicts per hour for a given aircraft). Fractal dimension also provides information on the number of degrees of freedom used in the airspace: a higher fractal dimension indicates more degrees of freedom. This information is independent of sectorization and does not scale with traffic volume.

Fractal dimension is thought to be an aggregate metric for measuring the geometrical complexity of a traffic pattern.

Some new geometrical metrics have been developed in [START_REF] Delahaye | A new air traffic complexity metric based on dynamical system modelization[END_REF], which are able to capture the level of disorder or the level of organization for some traffic patterns. For instance, it is possible to create configurations of aircraft where the speed vectors are very different, even if there are no changes in the relatives distances between aircraft. An example would be if all aircraft were holding in a circle around a common point. To capture those features, the covariance and the Koenig metrics have been developed. The covariance metric identifies disorder or organization in movement associated with translation motion of aircraft in the situation. The Koenig metric identifies disorder or organization associated with rotational motion of the set of aircraft in the situation. Previous metrics of complexity do not appear to capture the impact on complexity of all the possible organizations of air traffic (high vs. low density, how vs. low convergence, translation organization, curl organization, etc.). In contrast, the topological entropy (Kolmogorov entropy) appears to be a metric that is able to capture most of these important drivers of complexity. The non-linear form is able to identify structure within the trajectories, for example flows of aircraft following the same path at the same speed. The next section will describe the linear and the nonlinear forms of this metric in detail.

LINEAR DYNAMICAL SYSTEM MODELING

Principle

This metric is based on the modeling of a set of trajectories by a linear dynamical system. This allows for the identification of different organizational structures of the aircraft speed vectors such as translation, rotation, divergence, convergence or a mix of them.

Site visits to air traffic facilities and a review of previously identified ATC complexity factors suggests that organization in the distribution of aircraft position and speeds can have an important effect on the perceived complexity of the traffic situation [START_REF] Histon | Introducing structural consideration into complexity metrics[END_REF]. In particular, situations where the relative distances between aircraft do not change over time are more predictable and would appear to be easier to control. We identify these situations as being composed of "fully organized traffic." In order to capture the complexity associated with a lack of such organization, the air traffic situation is modeled as a dynamic situation. The set of aircraft trajectories is modeled by a linear dynamical system, which is defined by the following equation :

 Ẋ =A⋅  X  B Equation 2
where r X is the state vector of the system :

 X = [ x y z ] Equation 3
This equation associates a vector speed to a position in the space coordinate r X . This represents a particular vector field. The static behavior of this vector field is given by the vector r B , and the linear mapping between the speed vector r & X and the position vector r X is given by the matrix A.

The eigenvalues of the matrix A describe and summarize the evolution of the system. The real parts of those eigenvalues are related with the convergence or the divergence property of the system in the direction of the eigenvector. When such eigenvalues have positive real parts, the system is in an expansion mode and when they are negative, the system is in a contraction mode. Furthermore, the imaginary part of the eigenvalues are related to the level of curl tendency of the system. Depending on those eigenvalues, a dynamical system can evolve in contraction, expansion, rotation or a combination of those three modes. Figure 1 shows four typical examples for which the A matrix and the associated eigenvalues are computed. The small squares are the initial positions of aircraft at a given time(this represents the observation given by a radar for instance with the associated speed vector). As it can be seen on the figures, the aircraft are initially located on a circle with diameter 2d. The positions of aircraft 2 and 3 are symmetric with respect to the x axis. The speed norms of aircraft are the same and have been fixed at v. In the four situations, the aircraft are evolving in a two dimensions space, so our state vector will be limited to those two dimensions:

 X = [ x y ] Equation 4
This restriction enhances the simplicity of the equations for those examples, but the method works the same way in the three dimensional case. The computation of the four matrices A and vectors r B is given by the following equations.

Situation 1 Situation 2 Situation 3 Situation 4

A= [ 0 0 0 0 ] A= v d [ -1 0 0 -1 ] A= v d [ 1 0 0 1 ] A= v d [ 0 1 -1 0 ]  B= [ v 0 ]  B= [ 0 0 ]  B= [ 0 0 ]  B= [ 0 0 ]  λ= v d [ 0 0 ]  λ= v d [ -1 -1 ]  λ= v d [ 1 1 ]  λ= v d [ +j -j ] Equation 5
Given the matrices in Equation 5, Equation 2 holds for all observations for all situations. For instance, the position r X of aircraft 1 in the fourth situation is given by:

 X = [ -d 0 ] Equation 6 So  X =A⋅  X  B= v d [ 0 1 -1 0 ] ⋅ [ -d 0 ]  [ 0 0 ] = [ 0 v ] =  V 1 Equation 7
which is the speed of the initial observation for this aircraft (1) in the situation 3.

For those four situations, the equation 2 fits exactly with the observations. This exact fitting is not always possible with a linear system but a LMS (Least Mean Square) procedure can be used to generate a best fit between the dynamical system and the observations. The eigenvalues associated to the matrices A are given by the vectors r λ and the corresponding loci are shown on the bottom part of figure 1. When the relative distances between aircraft remain unchanged with time (situation 1 and 4), it can be noticed that the real parts of the eigenvalues of matrix A are null. When the norms of the relative distances between aircraft diminish with time (situation 2) the real part of the eigenvalues are negatives; finally, when those relative distances increase with time (situation 3) the real part of the eigenvalues are positive. Thus, the absolute value of the real part identifies the speed at which the relative distances are expanding / contracting and can therefore be interpreted as an indication of the amount of organization present in the traffic situation. If the eigenvalues of matrix A of the associated dynamical system have a small real part (in absolute value), then the relative distances between aircraft are changing slowly with time and the situation will likely show a high degree of organization. This is illustrated in figure 2 where the small vertical strip illustrates the eigenvalues of situations that are examples of "organized" dynamical systems. In this coordinate system, the organized situations (in translation, curve, or both) are located around the imaginary axis (shadow area), which represents the eigenvalues with a small real part.Thus, the eigenvalues of the A matrix can be interpreted an indicator of the overall degree of organization of the linear dynamical system. The following section develops a method for computing these eigenvalues, or Kolmogorov entropy, for more complicated situations.

Kolmogorov Entropy Computation

A dynamical system can also be considered as a map T from the state space X to itself. Based on this mapping T, one can define the topological entropy (or Kolmogorov entropy) of the dynamical system which measures the level of mixing of X by T . This entropy is associated with the changes of the relative distances between points from X by T . The topological entropy of a dynamical system is then a disorder indicator of the distribution of the aircraft in the airspace. Unlike a statistical entropy that requires many samples to be reliable, the topological entroy is independent of the number of samples and thus can be computed for situations where there are few aircraft in a sector.

The Kolmogorov entropy is computed with the help of the eigenvalue decomposition of the matrix A :

A=P -1 ⋅D⋅P Equation 8
where D is the diagonal matrix of the eigenvalues. Based on the observations of the aircraft (positions and speed vectors), the dynamical system has to be adjusted to minimize errors. This fitting has been done with a Least Square Minimization (LMS) method [START_REF] Bjorck | Numerical methods for least squares problems[END_REF]. For each considered aircraft i, it is assumed that the position  X = [ x i ,y i ,z i ] and speed vector 

V = [ v x i ,v y i ,v z i ] T are
known. An error criterion between the dynamical system model and the observation is computed :

E= ∑ i=1 i=N ∥  V i - A⋅  X i   B  ∥ Equation 9
The minimization of this criterion with the LMS induces an optimum set of parameters:

A = P -1 ⋅ D ⋅ P and r B .
In order to avoid matrix ill-conditioning problems in the system identification, a singular value decomposition approach is used to perform the LMS regression of the A matrix and the r B vector. An immediate consequence is that when a singular value is not large enough, the associated dimension will not be taken into account in the computation of the Kolmogorov entropy of the system. This will ensure that the matrix A can always be determined for the dimension space where the traffic is evolving. The eigenvalues of the matrix A (equivalent to those of matrix D) are complex numbers. In order to produce a scalar metric, the Kolmogorov Entropy, the larger absolute value of the real part of those complex eigenvalues is computed and the largest value identified. The sign of this larger real part is related with the mode of the system (contraction or expansion). If those eigenvalues are drawn in the complex coordinates, they give information about the level of contraction/expansion of the system and also its curl tendency. Such a metric can identify any speed vector organization (translation, round about or both) because it is sensitive to relative distances between aircraft. The associated eigenvalue evolution is then computed and is represented in figures 4 and 5. The evolution of the real part of the larger eigenvalue for the situation is given in figure 4. The metric starts negative, indicating the situation is globally converging and after the "crossing" the metric becomes positive, indicating the situation is now globally diverging. If this eigenvalue is represented in the complex coordinate system (see figure 5), the locus of the conjugate eigenvalue begins from the left side (the start is symbolized by crosses) and moves toward the right side. The analysis of such loci may be used in order to compare different traffic situations.

Linear dynamical system modeling produces an aggregate metric for any traffic situation that is interpreted as an indication of the overall organization of the situation. In order to identify particular areas in the airspace with high (low) complexity, a non-linear dynamical system can be used.

NON LINEAR EXTENSION

Introduction

Evolution equations of the form :

 Ẋ t  =V   X  t ,ω  Equation 10
with V C 2 vector field depending on parameters ω , describe systems whose integral curves may fit exactly the observed trajectories (which is generally not possible with linear models). There are many classical ways of obtaining a class of parameterized vector fields which fulfill the fitting requirement. Vector splines control the smoothness of vector fields, and are consistent with the low acceleration guidance control laws used for civil aircraft maneuvers are based on low acceleration guidance laws. An approach similar to the linear case is used: first, the approximating model is computed; then Kolmogorov entropy is estimated using sampled integral curves of the underlying vector field.

Div-Curl Splines

Computing topological entropy for a given traffic situation requires interpolating a vector field given a limited set of samples (positions and speeds of aircraft at a given time). Vector spline interpolation seeks the minimum of a function of the form:

1 2 ∫ D ∥L  V   X  ∥ 2 d  X 1 2 ∑ i=1 i=N ∥  V   X i  - V i ∥ Equation 11
where r V is a vector field defined on a domain D ⊂ R n ; L is an elliptic differential operator controlling smoothness of the solution; and   X i ,  V i  i= 1 .. N are the interpolation data [8]. By introducing the adjoint operator L T , the optimal vector field can be shown to be a linear combination of the shifted version of the elementary solution kernel of the differential operator L T L. A special case is the so-called "div-curl" splines with the criterion :

∇ div {  V ∇ curl {  V ∫ R 2 α ∥¿  X ∥ 2 +β∥¿  X ∥ 2 d  X Equation 12
with α , β positive reals controlling the smoothness of the approximation by focusing on constant divergence or constant curl.

Dynamic Div-Curl Splines

Div-Curl splines work well for spatial vector field interpolation (wind or marine flows data for example). However, in our case, aircraft trajectories may intersect even if there is no conflict. In such situations, the classical div-curl interpolation cannot provide a solution since the crossing point will become a singular point in the approximating vector field. To deal with this situation, a family of spatio-temporal splines are introduced. The splines are obtained from the optimal solution of a dynamic interpolation problem. This problem is to find a time-dependent vector field t, {  X  V ¿  defined on [0, T ] × D and continuously differentiable up to order 2 in time that minimizes :

t, {  X t, {  X ∫ 0 T ∫ D ∥ ∂  V ¿ ∂t ∥ 2 +α∥L  V ¿ ∥ 2 d  X dt Equation 13
under the constraints:

 V  t i ,  X i  =  V i ,i= 1 .. N Equation 14
where α is a positive real. This value controls the relative importance of the vector field variation thru time over the discrepancy of r V as measured by the differential operator L. Taking α = 0 will yield to a constant vector field over time, while α → +∞ will focus on the differential part.

Solving the Associated Variational Problem

Let L be an elliptic differential operator of order p with constant coefficients. L T will denote the adjoint operator of L. We will assume in the following that the vector field r V has fixed value at 0 and T and that for all t ∈ [0, T ], r V (t, .) belongs to the Sobolev space H p . Furthermore, the mapping t → r V (t, .) must be at least two times continuously differentiable.

Let γ : [0, T ] → C p (R n , R n ) a variation of X (that is time-dependent H p vector field γ (0, .) = (T, .) = 0). Simple calculus shows that the gradient of :

t, {  X t, {  X ∫ 0 T ∫ D ∥ ∂  V ¿ ∂t ∥ 2 +α∥L  V ¿ ∥ 2 d  X dt Equation 15
with respect to γ , is the linear mapping:

γ  ∫ 0 T ∫ D 〈 γ, ∂ 2  V ∂t 2 -αL T L  V 〉 d  X dt Equation 16
The constraints on X require that:

γ  t i ,.  =0,i= 0 . . N Equation 17
One can find an elementary solution e of the partial differential operator:

∂ 2  V ∂ t 2 -αL T L  V Equation 18
and use its relocated representative e t i ,  X i to represent the evaluation operator at t i ,  X i  .

Combining this with the gradient previously obtained and using Lagrange multipliers (denoted λ i ,i=1 .. N ) shows that the optimal solution is of the form :

 V = ∑ i=1 N λ i ⋅e t i ,  X i   V 0 Equation 19
with r V 0 an element of the kernel of the operator:

∂ 2  V ∂ t 2 -αL T L  V Equation 20

Complexity Map Computation

The metric chosen for complexity computation relies on a measure of sensitivity to initial conditions of the underlying dynamical system called Lyapunov exponents. In order to determine the Lyapunov exponents, consider a point and trace its evolution as it is transported by the dynamical system. Let r X be fixed (initial point) and let φ be a point trajectory of the dynamical system given by:

φ t, {  X =  X ∫ 0 t  V⋅φ τ, {  X dτ ¿ ¿ Equation 21
Assume now that trajectory is disturbed by a small perturbation r ε , we have:

φ t, {  X ε =φ t, {  X ¿ +D  X φ t, {  X ¿ ⋅ ε +o∥ ε∥ ¿ Equation 22
where D  X φ t, {  X ¿ is the differential of the vector field at r X that satisfies:

D  X φ t, {  X  dt =D  X  V  φ t, {  X ¿ ⋅D  X φ t, {  X ¿ ¿ Equation 23
The Lyapunov exponents are the singular values of the matrix D  X φ t, {  X ¿ and can be thought of as local shear values for the dynamical system. When Lyapunov exponents are high, the trajectory of a point under the action of the dynamical system is very sensitive to initial conditions (or parameters on which the vector field may depend), and the future of the situation cannot be reliably predicted. On the other hand, small values of the Lyapunov exponents mean that the future is highly predictable, which, as described above, is thought to be an indicator of situations that controllers find less complex. Thus, the Lyapunov exponent map can be interpreted as a description of relative degree of organization of the underlying dynamical system. It identifies the places where the relative distances between aircraft do not change with time (low real value) and the ones where the relative distances change rapidly (high real value).

Let us now describe the practical procedure for computing complexity maps. First, the optimal dynamic div-curl approximation for the observed trajectories is computed, based on the defining equations. This step requires a solution to the linear system.

The second step computes the second derivatives matrix D  X φ t, {  X ¿ at each point of the grid  X for the dynamical system trajectory φ starting at r X . This is done by solving the differential equation ( 23) with a Runge-Kutta integrator. The complexity value at point r X is then obtained by averaging Lyapunov exponents over the time:

κ   X = 1 n ∑ i=1 i=n ∥D  X  V  φ t i ,  X ∥ 2
Equation 24

RESULTS

In order to validate the previous method, complexity maps have been computed for several scenarios. The complexity maps show the evolution of Lyapunov exponents in time and space (the level of predictability of the trajectories for a given 4D point (space and time)). Two situations are presented in this paper, but this method can be used for any traffic sample. The first experiment contained four artificial traffic situations with some common geometrical properties. The situations are presented in figure 5. The first one (5.1) consist in 8 aircraft flying parallel trajectories at the same speed. Those trajectories are well separated and the expected complexity must be near zero. The second situation (5.2) presents 8 aircraft converging at the same point with the same speed; the third one (5.3) shows 8 aircraft moving at the same speed , all with trajectories on a circle of constant radius. Finally, the fourth situation (5.4) is a random distribution of 8 trajectories with different speeds. The complexity in the airspace is computed in a 3D grid for which only the horizontal cuts are presented in the interface. This process enables the visualization of the complexity in any part of the airspace. The separation norm being different in the horizontal direction (ex 5 nmi) and the vertical direction (1000ft), the complexity has been computed with a non-isotrop distance, meaning that two aircraft separated by 5 nmi in the horizontal plane or 1000ft in the vertical direction have the same "distance" in the new non-isotropic coordinate system for which the unit are separation norm in any direction. As the method is based on deteriming a dynamic system that represents the traffic situation at any moment, a complexity map can be computed at any arbitrary moment in time. Figure 6 presents the situation after 5 minutes of flightbased on the trajectories shown in figure 5.

The black squares show the current positions of aircraft and moving directions are given by the associated little black line (these lines represent the current speed vectors of aircraft).

The results shown are consistent with the expectations from previous work. The first situation with 8 aircraft in parallel has a very small complexity value because small deviations will make little difference to the prediction of the future behavior of the aircraft situation (Kolmogorov entropy can also be interpreted as the level of predictability of a traffic situation). In the second situation, the metric shows a strong interaction of the trajectories at the central crossing point. This interaction was so strong that a different intensity scaling was used in order to reveal the peak at the centre. The third situation locates a potential interaction of the trajectories in the curl. Finally, in the fourth situation, the dark areas locate the positions where the aircraft are (or will be) in interaction: the darker the region the greater the interaction between aircraft in terms of rate of convergence.

The second experiment has been performed using radar traffic data with a large number of aircraft. The data was obtained from the "Coubron" radar site that covers the north west of France. The radar data used was collected over 3 hours on September 24, 2002. A software packaged was developed to produce a complexity map at one minute intervals. The maps can be combined in order to produce a movie.

Figure 8. Northwest France, September 24, 2002 Figure 7 presents one of those maps for a given time. The small squares are individual aircraft (black for the selected altitude level and shadow for the others) and the arrows give the distribution of the vector field for the selected altitude. The map shows the disorder in the aircraft trajectories located in a small layer around flight level 150 (from FL145 till FL165 (the size of this layer is a parameter of the software)). The traffic located in the dark areas is unpredictable as the relative distances between aircraft in those areas are chaotic with time. In contrast, the traffic in the white areas is very predictable (the relative distance between aircraft do not change a lot with time) which is associated with reduced complexity for the controller. For instance if a pair of aircraft are separated in a white area, that pair will stay separated in the near future.

A complexity map can be built for any airspace (e.g. a particular area, sector, TMA etc…). The complexity map is very useful for comparing different airspaces as it shows and quantifies regions of the highest complexity. For example, the complexity around Paris is clearly the highest on the map. For the same level of traffic (same number of aircraft), the complexity distribution of different areas can be compared, potentially supporting the quantification of relative differences between airspace. Furthermore, the complexity maps may also provide additional insights into how different examples of structure in the airspace affects the degree of organization of traffic and hence the impact on the perceived complexity of controlling that traffic..

CONCLUSION

We have presented in this paper a new air traffic complexity metric based on a non-linear vector field model of air traffic. Extending previous results on topological entropy, this method identifies organization in the traffic pattern, an important factor in assessing ATC complexity. This approach is based on trajectory segments instead of position and speed samples at a given time. The trajectory segments can take into account the past and the future (prediction) on the current aircraft position and thus is consistent with the controller's approach. Future work will expend the method to include intent information that is available to the controller throught the flight plan. The quasi-interpolation algorithm allows real-time processing on operational traffic even for large areas (Europe or US). Unlike linear models that produce mean complexity indicators, the nonlinear one may give local information, thus providing a way of displaying maps of complexity. In a future work, such a tool will be applied to a comparison of US and European airspace by producing complexity maps for a full day of traffic.

Figure 1 .

 1 Figure 1. Example of eigenvalue locations for several simple trajectory configurations.

Figure 2 .

 2 Figure 2. Impact of the eigenvalues on the dynamics of a system.

Figure 3 .

 3 Figure 3. Trajectories of 40 aircraft passing through a common area.

Figure 3 .Figure 4

 34 Figure 3. Evolution of the real part of the eigenvalue with the larger norm

Figure 5 .

 5 Figure 5. Four artifical traffic situations

Figure 6 .Figure [ 6 ]

 66 Figure 6. Complexity maps for the four artificial traffic situations. Intensity scaling for 6.2 has been adjusted in order to reveal peak in central region (6.1, 6.3, 6.4 have a common intensity scale)Figure[START_REF] Delahaye | A new air traffic complexity metric based on dynamical system modelization[END_REF] shows the complexity map for each situation when viewed as a horizontal cutting plane. The figure shows the associated complexity of the 4 situations in a grey scale. The results show that complexity is highest in the regions where aircraft trajectories interact with each other. The complexity in the airspace is computed in a 3D grid for which only the horizontal cuts are presented in the interface. This process enables the visualization of the complexity in any part of the airspace. The separation norm being different in the horizontal direction (ex 5 nmi) and the vertical direction (1000ft), the complexity has been computed with a non-isotrop distance, meaning that two aircraft separated by 5 nmi in the horizontal plane or 1000ft in the vertical direction have the same "distance" in the new non-isotropic coordinate system for which the unit are separation norm in any direction. As the method is based on deteriming a dynamic system that represents the traffic situation at any moment, a complexity map can be computed at any arbitrary moment in time. Figure6presents the situation after 5 minutes of flightbased on the trajectories shown in figure5.The black squares show the current positions of aircraft and moving directions are given by the associated little black line (these lines represent the current speed vectors of aircraft).