Clement Peyronne
email: clement.peyronne@capgemini.com

Daniel Delahaye
email: delahaye@recherche.enac.fr

Marcel Mongeau
email: marcel.mongeau@math.univ-toulouse.fr

Laurent Lapasset
email: llapasse@capgemini.com

OPTIMIZING B-SPLINES USING GENETIC ALGORITHMS APPLIED TO AIR-TRAFFIC CONFLICT RESOLUTION

Keywords: Air-traffic conflict resolution, B-splines, Genetic algorithms

Conflict resolution has always been a sensitive matter in air-traffic management. Current European projects aim partial or total automation of air traffic control to deal with the constant growth of air traffic. Technological advances on flight management system allows us to consider an automatic conflict resolution using continuous trajectories. In this paper, we present a new methodology that, first, relies on B-splines to model trajectories, secondly models air-traffic conflict resolution as an optimization problem whose decision variables are the spline control points. Finally, we use genetic algorithms to tackle this optimization problem in order to generate optimal conflict-free situations.

INTRODUCTION

From the beginning, the most critical point of air traffic control was to ensure safety separation between airplanes. To achieve this goal a safety standard separation distance has been defined: 5 Nm (Nautical miles) horizontally and 1000 feet vertically. Air traffic controllers are responsible for ensuring the respect of these separation rules.

Context

Air traffic being constantly increasing, controllers in charge of an aviation sector must handle more and more flights. Up to now, decreasing the size of aviation sectors could compensate the growth of traffic but we are reaching the point where a decrease of the size of sectors is no longer efficient.

Nowadays, air traffic management has already used every available ressource to increase airspace capacity. However, from now to 2020, air traffic is expected to be multiplied by two or three. Consequently, air traffic management (ATM) will have to deal with this overload while ensuring at least equivalent standards of safety [START_REF] Alliot | L'ATM : 20 ans d'efforts et perspectives[END_REF]. The SESAR European project aims at finding solutions to this problem by providing a decision support to air traffic controllers in order to decrease their workload. Considering the technological advances on the aircraft flight management system (FMS), we will explore in this paper the possibility of a full automation generating continuous trajectories that new FMS can follow.

Previous Related Work

By the past, two kinds of method have become predominant in automatic air traffic conflict resolution for their good results : genetic algorithms (GA) [START_REF] Durand | Algorithmes génétiques et autres outils d'optimisation appliqués à la gestion du trafic aérien[END_REF] and navigation functions [START_REF] Zeghal | Vers une théorie de la coordination d'actions, application à la navigation aérienne[END_REF].

Genetic algorithms are part of what is known as natural computation or evolutionary methods. This optimization method is based on the evolution theory and uses concepts such as mutation, crossover and, of course, selection. Each possible solution of our optimization problem is encoded as a chromosome via a specific encoding (for example, number and duration of simple maneuvers). The algorithm creates randomly a first chromosome population. Each chromosome ability to solve the problem (fitness) is then evaluated. The best individuals (according to their fitness) are selected and mutation, crossover are applied to obtain a new population of chromosome (next generation). The user chooses the necessary generation number for the algorithm to converge (this choice is usually done empirically). The obtained trajectories' velocity stays within acceptable bounded range with respect to ATM operational constraints. However, these techniques have not been tested yet with curved trajectories.

Navigation functions use a different representation. Indeed, the airspace is considered as a potential field, and airplanes like particles navigating in it. Naturally, negative charges represent obstacles to the airplane (other aircrafts, congested areas) and positive charges, its destination. As a result, each airplane is attracted to its destination while being repulsed by obstacles which enables to the automatic generation of conflict-free trajectories. The major drawback of this method is that the obtained solution does not respect the ATM constraints such as limited speed, or trajectory smoothness. Evenmore, they can lead to major delays and overcosts as they tolerate any deviation from the business trajectories (BT : the straight line between the departure and arrival points).

In this paper, we present a new methodology that, first, relies on B-splines to model trajectories, secondly models air-traffic conflict resolution as an optimization problem whose decision variables are the spline control points. Finally, we use genetic algorithms to tackle this optimization problem in order to generate optimal conflict-free situations.

We present our trajectory model in Section 2. Sec-tion 3 presents the optimization method we chose to solve our optimization problem: genetic algorithms (GA). We detail in Section 4 the objective function (conflict detection) of our optimization problem. Numerical results are presented in Section 5. We draw conclusions in Section 6.

TRAJECTORY MODEL

B-spline was primarily introduced to find a curve interpolating a set of points of R 2 called control points. It was later extended to approximation, thereby avoiding the undesirable oscillation inherent to interpolation. In our study, we shall focus on this use of splines to approximate a set of control points. The control polygon, the piecewise linear curve linking the control points, completely defines the curve [START_REF] Rabut | Courbes, outils pour la CAO, approximation à une variable[END_REF].

Basically, B-splines are parameterized curves generalizing the Bezier curve concept. It is an efficient approximation tool which is constructed from polynomial pieces joined at certain parameter's values called knots, stored in a knot vector. In a very simplified way, if we consider a set of control points

(X i ,Y i ) = P i ∈ R 2 (i = 0, 1, ..., n)
, and a parameter u, we can define the B-spline as follows :

C(u) = (σ x (u), σ y (u)), u ∈ [a, b]
where σ x (u) and σ y (u) are the B-splines approximations of the couples (i, X i ) i=0:n and (i,

Y i ) i=0:n for u ∈ [a, b] .
Theoritically, the curve is calculated by multiplying the B-splines basis functions by the control points coordinates. For σ x :

σ x (u) = n+1 ∑ i=-1 X i B i (u) (1)
However, in practice, B i construction is made by recurrence. Consequently, we used a more numerically convenient way to calculate the B-spline values only using control points values. The B-spline and its derivative values at knot u i are defined by :

σ i = X i+1 +4X i +X i-1 6 σ ′ i = X i+1 -X i-1 2h σ ′′ i = X i+1 -2X i +X i-1 h 2 σ ′′′ i = X i+2 -3X i+1 +3X i -X i-1 h 3 (2)
Therefore, by using Taylor series, we obtain the Bspline value for all u ∈ [u i , i+1 ] :

σ x (u) = σ x i +(u-u i )(σ ′ x i +(u-u i )( σ ′′ x i 2 +(u-u i ) σ ′′′ x i 6
))

(3) We choose to rely on B-spline modeling trajectories because it is a very efficient tool for curve approximation in terms of both approximation quality and computational time. Moreover, B-splines feature interesting properties such as C 2 -continuity (crucial for modeling smooth aircraft trajectories), robustness and flexibility (if one control point is displaced, only a small part of the curve will be affected). The last good point in favor of B-splines is its compatibility with GA : the curve is completely determined by few control points only, which will be the optimization parameters of the GA. In our study, we use a fixed maximal number of control points between the departure and the arrival points.

OPTIMIZATION METHOD

When several aircraft are involved in a conflict, the conflict resolution problem has been shown to be NPhard [START_REF] Durand | Algorithmes génétiques et autres outils d'optimisation appliqués à la gestion du trafic aérien[END_REF]. Moreover, the optimization variables being the B-splines control-point location, we shall see that our objective function ( 4) is not differentiable with respect to these variables. Consequently, we must rely on black box (direct) optimization methods to address our problem. In this paper, we choose to use a stochastic global optimization method, genetic algorithms, to guide the control-point location. As we mentioned in the previous part, GA seem to be appropriate with B-splines. In this section, we explain briefly GA theory and then we shall detail how we adapt the different operators (selection, mutation and crossover) to our air traffic conflict resolution problem.

Basis of Genetic Algorithms

Evolutionary algorithms use techniques inspired by evolutionary biology such as inheritance, mutation, natural selection, and recombination (or crossover) to find approximate solutions to optimization problems [START_REF] Goldberg | Genetic Algorithms in Search Optimization and Machine Learning[END_REF], [START_REF] Koza | Genetic Programming[END_REF]). An individual, or solution to the problem to be solved, is represented by a list of parameters, called chromosome or genome. Initially several such individuals are randomly generated to form the first initial population (POP(k) with k = 0 on Figure 2). Then each individual is evaluated, and a value of fitness is returned by a fitness function (generally, the objective function of our optimization 2) that is different from the initial generation. This generational process is repeated until a termination condition has been reached. The next section presents the application of this algorithm to our specific problem.

Chromosome Encoding and Fitness Calculus

In this section we detail the link between our modeling method and genetic algorithms. Indeed, our optimization method is using B-splines to calculate the aircraft trajectories in order to determine each individuals fitness. To define this interaction, we first define an encoding.

A chromosome will represent the trajectories of N aircrafts using a matrix of N × N c where N c is the maximal number of control points we wish to consider. For example, the case where N = 20 airplanes and N c = 2 we encode the trajectories as follows : Let us now define precisely how we specify the control points P j i . We want a compromise between allowing the trajectory to deviate freely (in any directions) from the business trajectory (BT) and staying as close as possible from it. For that purpose, we define a fixed width band (D max ), depending on the BT length (D tot ), around the BT trajectory where the control points are allowed to lie.

A matter for the control point location is that, if there are all gathered in the same region, it can create oscillations on the resultant trajectory. Consequently, we distribute uniformly the control points along the BT.

To summarize, a control point can be defined by a single parameter representing a bandwidth percentage (± d D max % ).

Moreover, it is desirable to be able to have a different number of control points for a given chromosome. In order to satisfy this goal, we used a convention: percentage are generated within [-200%; 200%], if the percentage is between -100% and 100%, then we create a control point, if it belongs to [-200%; -100%[∪[100%; 200%], we do not. Consequently, our chromosome is a N × N c matrix of percentage. Several cases are presented in Figure 3 and4.

Crossover and Mutation

Let us now describe how we adapted the genetic operators we use in our GA.

Slicing crossover

The crossover operator aims at finding better solutions by mixing two good individuals of the previous generation. Therefore, we decide to keep trajectory consistencies by inverting only complete trajectories using a slicing crossover (we do not separate control points belonging to the same trajectory). Accordingly, the crossover consists in picking a plane number randomly and inverting the parent trajectories corresponding to the next aircrafts. An example is given in Figure 5 :

Mutation

The mutation operator aims at diversifying the genes in the population in order to explore as much as possible the problem space. Thus, mutation consists in choosing randomly one control point in the chromosome and to assign to it a new random number (using a uniform distribution) in the interval [-200%; 200%]. (see Figure 6). CONFLICT DETECTION METHOD

In order to evaluate each chromosome fitness, we decode it into N trajectory curves (one per aircraft) and evaluate two quantities. First, how many conflicts the situation engenders and secondly what is the total extra distance induced with respect to the BT. To calculate these quantities, we discretize the airspace into square cells of size half the standard separation (noted d). Our conflict detection is performed in two steps :

• First, for each airplane, we store the grid's cells through which the aircraft flies, the aircraft number (its label), the entry and exit times in and out each of the stored cell, • Then, we select each stored cell and we check whether any other airplane goes through any of the eight neighboring cells for other airplanes. If there are such airplanes, we check the time to see whether there is a conflict between these two airplanes. If so, we calculate the conflict duration.

Our conflict detection procedure send back the chromosome fitness to the GA. Here is the formula we use to calculate the fitness:

f (X) = -(CN + ( N ∑ i=1 NR i BT i -N)) ( 4
)
where CN is the number of conflict, NR i the length of the new route calculated by the algorithm for the aircraft i, while BT i the length of the business trajectory (straight line from the departure to the arrival point) and N the number of aircraft. High fitness corresponds to good individuals. Indeed, the lower are the number of conflict and the route lengthening, the better is the chromosome.

RESULTS AND PERSPECTIVES

In this section, we present several results we obtained using our methodology. First, we test our method on the roundabout test problem which is a very common test on automatic conflict-detection method. It consists in making a certain number of planes fly to the diametrically opposed point at a common speed (each point on the circle has an outgoing and an incoming trajectory).

We choose to fly N = 16 aircrafts for results' readability) equidistributed on a circle of 100Nm(= 185200m) radius To address this conflict resolution, we use the following parameters :

• The fitness' evolution with respect to generations is shown in Figure 8. The fitness is meant to increase to 0 as formula (4) shows. Moreover, when the fitness is in [-1, 0] the situation is conflict-free. Consequently, the algorithm stops only when the best individual fitness is in [-1, 0]. Although one can easily solve intuitively this academic problem due to its symmetry, our automatic implementation does not exploit any symmetry here. This result shows that our methodology is promising as the obtained conflict resolution is consistent with experts' experience.

We also test our methodology on a more realistic operational problem where N = 15 aircrafts are involved in a fuzzy convergence. In this situation, for each aircraft, we have its departure point, its heading and its speed. To calculate the arrival points, we make the 15 aircrafts fly for one hour at a constant speed, on a straight line following the initial direction (heading). We present the conflict-free situations our algorithm obtains is displayed in Figure 9. We can see on Figure 10 that on an operational application, our method converges far more rapidly in three generations (90 generations for the roundabout test problem). This results shows that our method can be very efficient from an operational point of view. Indeed, in operational problems, there are never more than four aircrafts involved in the same conflict. Consequently, our algorithm will be able to deal very efficiently with operational air traffic. For this situation, we used a different number of generations because of its extremely quick convergence (number of generation: 20). Moreover, we can see that few aircrafts are deviated, which is consistent with a controller behaviour, who would prefer to deviate largely one aircraft instead of deviating lightly several aircrafts to solve conflicts.

CONCLUSION

We have shown in this paper that the combination of B-splines and genetic algorithms can be a promising methodology for automatic conflict resolution in air traffic control. However, we have in mind several developments to improve our approach such as using sharing in GA (deals with equirepartition of the population on the different maximums) in our GA or implementing a self-adaptative GA (every parameter such as bandwidth, control points' number, etc will be considered within the chromosome encoding, as proper optimization variables). Furthermore, we plan to exploit our B-spline model of trajectory to address the conflict resolution problem with deterministic derivative-free optimization methods [START_REF] Conn | Introduction to derivative-free optimization[END_REF]. Indeed, despite the local aspect of these methods, they can also be adapted to global optimization.

Figure 1 :

 1 Figure 1: Principle of navigation functions between two aircrafts (left), between an aircraft and an obstacle (right) It has been demonstrated that navigation function methods ensure collision avoidance and connection between departure and destination.The major drawback of this method is that the obtained solution does not respect the ATM constraints such as limited speed, or trajectory smoothness. Evenmore, they can lead to major delays and overcosts as they tolerate any deviation from the business trajectories (BT : the straight line between the departure and arrival points).

Figure 2 :

 2 Figure 2: Genetic algorithm with tournament selection. The first step consists in the selection of the best individuals from population POP(k). Afterwards, recombination operators are applied in order to produce the POP(k + 1) population. problem). This initial population undergoes a selection process which identifies the most adapted individuals. The selection process which has been used in our experiments is a deterministic (λ, µ)-tournament selection. This selection begins by randomly selecting λ individuals from the current population POP(k) and keeps the µ bests (λ > µ). These two steps are repeated until a new intermediate population (POP i ) is completed. Following selection, one of the three following operators is applied : nothing, crossover, and mutation. The associated probability of application are respectively (1 -p c -p m )), p c and p m . Crossover results in two new child chromosomes, which are added to the next generation population. The chromosomes of the parents are mixed during crossover. These processes ultimately result in the next generation population of chromosomes (POP(k+1) on Figure2) that is different from the initial generation. This generational process is repeated until a termination condition has been reached. The next section presents the application of this algorithm to our specific problem.

  and arrival points are not encoded in the chromosome as they are not meant to move, they are not optimization variable.

Figure 3 :

 3 Figure 3: Chromosome encoding and corresponding Bspline with one percentage exceeding 100%

Figure 4 :

 4 Figure 4: Another chromosome encoding and its corresponding B-spline with both percentages within [-100%, -100%]

Figure 5 :

 5 Figure 5: Parent and children chromosomes

Figure 6 :

 6 Figure 6: Initial and mutated chromosome 4 OBJECTIVE FUNCTION : CONFLICT DETECTION METHOD

  Number of generation : 100 • Size of the population : 500 • Mutation probability : p m = 0.3 • Crossover probability : p c = 0.6 • Maximal number of control points : N c = 2 To compute the conflict resolution, we used a 2.53GHz Intel(R) Core(TM)2 Duo on a Windows Vista Operating System and we coded in Java. The obtained conflict-free situation is presented in meterscale (1Nm = 1852m) in Figure 7. This whole resolution computation required approximately 4 minutes.

Figure 7 :

 7 Figure 7: The roundabout configuration after resolution

Figure 8 :

 8 Figure 8: The fitness evolution with respect the generation

Figure 9 :

 9 Figure 9: The operational problem after resolution

Figure 10 :

 10 Figure 10: The fitness evolution with respect the generation

ACKNOWLEDGEMENTS

This work has been supported by French National Research Agency (ANR) through COSINUS program (project ID4CS n°ANR-09-COSI-005) and Capgemini through the first author PhD funding.


Marie-Laurence Mazure
email: mazure@imag.fr

Numerical Analysis Which spline spaces for design?

We recently determined the largest class of spaces of sufficient regularity which are suitable for design. We are now interested in connecting different such spaces to obtain the largest class of splines usable for design.

Introduction

To justify the framework presented in Section 3, we draw the reader's attention on two points. First, a "good" spline space is expected to permit the development of all classical design algorithms, among which any knot insertion algorithm. This is the reason why we naturally have to require the sectionspaces themselves to be "good for design". This explains why Section 2 is devoted to a brief presentation of the major results concerning Quasi Extended Chebyshev spaces, see [START_REF] Mazure | Which spaces for design[END_REF] and [START_REF] Mazure | Quasi-Extended Chebyshev spaces and weight functions[END_REF]. Second, the presence of connection matrices in our spline spaces is justified by the fact that requiring that left and right derivatives coincide at a given knot up to some order has no special meaning for parametric spline curves.

Quasi Extended Chebyshev spaces

Given a non-trivial real interval I and an integer n 1, let E ⊂ C n-1 (I) be an (n+1)-dimensional space. It is said to be of strict dimension (n+1) if it remains of dimension (n+1) by restriction to any non-trivial subinterval of I. We say that E is a Quasi Extended Chebsyhev space (for short, QEC-space) on I when any Hermite interpolation problem in (n + 1) data which is not a Taylor interpolation problem has a unique solution in E [START_REF] Mazure | Which spaces for design[END_REF]. By comparison, E is an (n + 1)-dimensional EC-space on I when any Hermite interpolation problem in (n + 1) data has a unique solution in E. This includes Taylor interpolation problems and therefore being an (n + 1)-dimensional EC-space on I requires E to be contained in C n (I). The interest of QEC-spaces lies in the following result: Theorem 2.1 Given n 2, let E ⊂ C n-1 (I) be of strict dimension (n + 1) and let it contain constants. The following two properties are then equivalent: (i) E is "good for design", in the sense that E possesses blossoms;

(ii) the space DE := {DF := F ′ | F ∈ E} is an n-dimensional QEC-space on I.

The reader interested in the exact definition of blossoms can refer to [START_REF] Mazure | Which spaces for design[END_REF]. Let us just recall that these powerful tools are defined in a geometrical way using intersection of osculating flats. Our terminology "good for design" is highly justified by the fact that, when (i) is satisfied, each function F ∈ E "blossoms" into a function f : I n → IR, called the blossom of F which satisfies the following properties: [n] standing for x repeated n times;

(B) 1 symmetry: f is symmetric on I n ; (B) 2 diagonal property: for all x ∈ I, f (x [n] ) = F (x), x
(B) 3 pseudoaffinity property: given any y 1 , . . . , y n-1 , a, b ∈ I, with a < b, there exists a strictly increasing function β(y 1 , . . . , y n-1 ; a, b; .) : I → IR (independent of F ) such that:

f (y 1 , . . . , y n-1 , x) = 1 -β(y 1 , . . . , y n-1 ; a, b; x) f (y 1 , . . . , y n-1 , a) + β(y 1 , . . . , y n-1 ; a, b; x)f (y 1 , . . . , y n-1 , b), x ∈ I. (1)
The latter three properties are crucial because they permit the development of all the classical geometric design algorithms in E ( e.g., de Casteljau algorithms), and they guarantee shape preserving proprties.

Example. Let us recall a simple procedure to build QEC-spaces on I [START_REF] Mazure | Which spaces for design[END_REF]. Take -any sequence of weight functions on I, i.e., any sequence (w 0 , . . . , w n-1 ) such that, for 0 i n -1, w i is positive and C n-1-i on I;

-any two-dimensional space U ⊂ C 0 (I) supposed to be a Chebyshev space on I (for short, C-space on I: any non-zero element in U vanishes at most once in I, not counting possible multiplicities) containing constants, that is, any space spanned by 1I, U where U is strictly monotone on I.

As is classical one can define differential operators L 0 , . . . , L n-1 on C n-1 (I) as follows:

L 0 F := F w 0 , L i F := 1 w i DL i-1 F, 1 i n -1. ( 2
)
Then the set of all functions F ∈ C n-1 (I) such that L n-1 F ∈ U is an (n + 1)-dimensional QEC-space on I, denoted QEC(w 0 , . . . , w n-1 ; U). For instance, given any positive numbers p, q > n -1, the linear space E p,q n spanned by the functions 1, x, . . . , x n-2 , (1x) p , x q is a QEC-space on [0, 1], see [START_REF] Costantini | Curve and surface construction using variable degree polynomial splines[END_REF], [START_REF] Costantini | On a class of weak Tchebycheff systems[END_REF], [START_REF] Kaklis | Convexity preserving polynomial splines of non-uniform degree[END_REF]. We recently established the following converse property, of which the proof strongly relies on Theorem 2.1 and on the pseudoaffinity property of blossoms.

Theorem 2.2 Let E be an (n + 1)-dimensional QEC-space on a closed bounded interval I. Then there exists infinitely many ways to write E as E = QEC(w 0 , . . . , w n-1 ; U).

The result

Given I = [a, b], a < b, the ingredients to build our spline space S are :

-a sequence of interior knots: a < t 1 < • • • < t q < b and an associated sequence of multiplicities m k , with 1 m k n for 1 k q;

-a sequence of section spaces E k , 0 k q: for each k, E k contains the constant function 1I k and DE k is an n-dimensional QEC-space on [t k , t k+1 ], where t 0 := a, t q+1 := b;

-a sequence of connection matrices M k , 1 k q: for each k, M k is a lower triangular square matrix of order (nm k ) with positive diagonal entries. The spline space S is then defined as the space of all continuous functions S : I → R such that (1) for k = 0, . . . , q, the restriction of S to [t k , t k+1 ] belongs to E k ;

(2) for k = 1, . . . , q, the following connection condition is fulfilled:

S ′ (t + k ), . . . , S (n-m k ) (t + k ) T = M k • S ′ (t - k ), . . . , S (n-m k ) (t - k ) T . ( 3
)
As a consequence of Theorems 2.1 and 2.2, for each k = 0, . . . , q one can choose a system (w k 1 , . . . , w k n-1 ) of weight functions on [t k , t k+1 ] and a two-dimensional C-space U k on [t k , t k+1 ] containing the constant function 1I k such that E k = QEC(1I k , w k 1 , . . . , w k n-1 ; U k ). For each such choice, when denoting by

L k 0 = Id, L k 1 , . . . , L k n-1 the differential operators on C n-1 ([t k , t k+1 ]) associated with (1I k , w k 1 , . . . , w k n-1
), the connection conditions (3) can be replaced by

L k 1 (t + k ), . . . , L k n-m k (t + k ) T = M k • L k-1 1 (t - k ), . . . , L k-1 n-m k (t - k ) T , 1 k q, (4)
where M k is in turn a lower triangular matrix of order (nm k ) with positive diagonal entries.

The result we announce is then the following:

Theorem 3.1
The following two statements are equivalent: (i) the spline space S is good for design, in the sense that it possesses blossoms;

(ii) for each k = 0, . . . , q, there exists a system (w k 1 , . . . , w k n-1 ) of weight functions on [t k , t k+1 ] and a two-dimensional C-space U k on [t k , t k+1 ] containing the constant function 1I k such that E k = QEC(1I k , w k 1 , . . . , w k n-1 ; U k ) and such that each matrix M k in (4) is the identity matrix of order (nm k ). Again, for the sake of simplicity we omit the precise definition of blossoms for splines, simply mentioning that, when (i) holds, each spline S ∈ S "blossoms" into a function s (the blossom of S), of which the natural domain of definition is a restricted set A n (K) of n-tuples said to be admissible w.r.t. the knotvector K := (ξ -n , . . . , ξ m+n+1 ) := (t 0

[n+1] , t 1 [m1] , . . . , t q [mq] , t q+1 [n+1] ), with m := q k=1 m k . To conclude the section let us stress that Theorem 3.1 presents a twofold interest:

• first, it provides us with a very simple recipe to build all good spline spaces;

• for given section spaces, it enables us to answer the following question: a sequence M k q k=1 of connection matrices being given, is the corresponding spline space suitable for design or not?

Sketch of the proof -Illustration

Theorem 3.1 extends to QEC-spaces a recent result obtained for EC-spaces [START_REF] Mazure | How to build all Chebyshevian spline spaces good for Geometric Design[END_REF]. Though guided by the ideas developed in [START_REF] Mazure | How to build all Chebyshevian spline spaces good for Geometric Design[END_REF], replacing EC-spaces by QEC-spaces is not a trivial adaptation of existing results, especially as we are working with different section spaces. Below we mention the two majors difficulties encountered in the present situation:

• When (i) holds, blossoms are obviously symmetric on A n (K) and they obviously satisfy the diagonal property. A major difficulty consists in proving that they also satisfy the crucial pseudoaffinity property (B) 3 , of course limited to A n (K). The tricky proof of (B) 3 in QEC-spaces [START_REF] Mazure | Which spaces for design[END_REF] involved difficult generalised convexity arguments. The spline framework will require that we go deeper into the arguments in question. We would like to stress that pseudoaffinity is THE property justifying that a spline space is considered good for design when it possesses blossoms. Indeed, it permits all the geometric design algorithms and leads to the important intermediate result stated below (see [START_REF] Mazure | Ready-to-blossom bases and the existence of geometrically continuous piecewise Chebyshevian Bsplines[END_REF] for the EC case):

Theorem 4.1 If (i) holds, then S possesses a quasi B-spline basis (see below) which is its optimal normalised totally positive basis. Conversely, if S and any spline space derived from S by knot insertion possess B-spline bases, then (i) holds.

Normalised totally positive bases are crucial for design [START_REF] Goodman | Total Positivity and the Shape of Curves, in Total Positivity and its Applications[END_REF], and optimality should simply be understood as "the best possible" [START_REF] Carnicer | Total positivity and optimal bases[END_REF]. For the links with blossoms, see [START_REF] Mazure | Blossoms and optimal bases[END_REF]. The presence of quasi B-spline basis is crucial too. Let us recall that it means a normalised sequence N ℓ ∈ S, -n ℓ m, ( i.e., m ℓ=-n N ℓ = 1I), each N ℓ being positive on the interior of its support [ξ ℓ , ξ ℓ+n+1 ], and satisfying some additional condition of zeros at the endpoints of its support. The term "quasi" simply refers to the fact that we are dealing with QEC-spaces, not with EC-spaces, and the count of zeros should take this into account.

• For simplicity we have assumed all interior mulitplicities to satisfy m k 1. Still, it is essential to understand what does m k = 0 would mean when dealing with QEC-space (see Example below).

Example. Let us give an elementary example. Given any real numbers p k , q k > n -1, 0 k q, here we consider a spline space S with kth section space E p k ,q k n up to the convenient affine change of variable.

1-To remain inside the framework described at the beginning of Section 3, it is sufficient to consider the case where all interior knots are simple. Given any positive numbers a k , 1 k q, the splines are then supposed to be C n-2 on I and to meet the additional requirement:

Theorem 3.1 guarantees that the spline space S is then good for design. In particular it therefore possesses a quasi B-spline basis which is its optimal normalised totally positive basis.

2-To illustrate multiplicities equal to 0, let us now assume that m k ∈ {0, 1} for all k = 1, . . . , q. Any interior knot t k of multiplicity m k = 0 is then allocated two additional positive numbers b k , c k . Along with (5) the splines in S are now supposed to satisfy: S (n-1) (t - k+1 ) = -b k S (n-1) (t + k-1 ) + c k S (n-1) (t - k ) for any k ∈ {1, . . . , q} such that m k = 0.

Again, this yields a spline space good for design.


Caroline Bernard-Michel

Sylvain Douté

Mathieu Fauvel

Laurent Gardes

Stéphane Girard

Support vectors machines regression for estimation of Mars surface physical properties

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Hyperspectral images are now a key tool for the analysis of remote planets [START_REF] Carlson | the NIMS Science, and Engineering Teams. Near infrared spectrometer experiment on Galileo[END_REF]. The very high resolution in the spectral domain allows a fine characterization of the physical properties of the scene. For instance, the OMEGA sensor acquires the spectral radiance coming from the planet in more than 200 contiguous spectral channels. A pixel of the image is represented by a spectrum/vector x ∈ R d , each component corresponds to a particular wavelength, d being the total number of wavelengths. Chemical composition, granularity, texture, and physical state are some of the parameters that characterize the morphology of spectra and thus the area of interest.

Deducing the physical parameters y from the observed spectra x is a central problem in classical physics, called an inverse problem. Since it generally cannot be solved analytically, optimization or statistical methods are necessary. Solving inverse problems requires an adequate understanding of the fundamental physics, i.e. a functional relation between x and y must be specified: x = g(y). Given g, different methods can be used to deduce the parameter from new observations. Current solutions to inverse problems can be divided into three main categories (for further details and comparisons, see [START_REF] Kimes | Inversion methods for physically-based models[END_REF]):

1. Optimization algorithms: These methods minimize an objective quality function that measures the fit between x and g(y). Inverse problems are often ill-posed, therefore estimations can be unstable and a regularization is needed. For instance, a prior distribution on model parameters can be used.

2. Look-up table (LUT) / k-nearest neighbors approach (k-NN): A large database (LUT) is generated by radiative transfer for many parameter values. Each observed spectrum is then compared with the LUT spectra in order to find the best match (the nearest neighbor), according to an objective function minimization. Parameters are then deduced from this best match.

Training approaches:

A functional relation y = f (x) between spectra and parameters is assumed, such as f -1 = g, and a LUT is used to estimate f . For instance, training approaches include neural network approaches. Inversion algorithms for hyperspectral images must be defined with the following constraints: a) Large datasets and various models require fast methodologies, b) High-dimensional data and "curse of dimensionality" with the associated sparsity issues require simple model with few parameters, c) Noisy spectra require robust algorithms. Optimization based approaches suffer from heavy computational cost and are thus not adapted to hyperspectral images. Nearest neighbors algorithms are faster, but the solution is unstable and noise-sensitive. Problems associated to training methods are the difficult interpretation of the estimated function f and the choice of the model. Moreover, neural networks are in general difficult to train for high dimensional data [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]Chapter 11].

The results presented in this paper follow previous research on hyperspectral image inversion [START_REF] Bernard-Michel | Inverting hyperspectral images with gaussian regularized sliced inverse regression[END_REF]. Support Vectors Machines (SVM) are proposed to estimate the functional f . This training approaches is robust with respect to the dimensionality and has good generalization abilities, see [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]Chapter 12] or [START_REF] Plaza | Recent advances in techniques for hyperspectral image processing[END_REF] for the classification of hyperspectral images. Using the kernel trick, the estimation of non-linear functions f is possible. Similarly to neural networks, the choice of the kernel function is not straightforward. In this article, several kernel functions were investigated, from linear to non-linear, leading to different levels of accuracy. For linear SVM, the connection with a dimension reduction technique, the Gaussian Regularized Sliced Inverse Regression (GRSIR) [START_REF] Bernard-Michel | Inverting hyperspectral images with gaussian regularized sliced inverse regression[END_REF], is discussed. SVM are compared with GRSIR and Partial Least Square regression both in terms of accuracy and computing time.

SVM and the different kernels are presented in Section 2. Experimental results are reported in Section 3 and discussed in Section 4.

Support Vectors Machines Regression

Regression Problem

SVM are a supervised methods for regression or estimation stemming from the machine learning theory. For inversion problems, the algorithm, which is called the ǫ-SVR, approximates the functional f : y = f (x) using solutions of the form

f (x) = n i=1 α i k(x, x i ) + b,
min 1 n n i=1 l f (x i ), y i + λ f 2 with l f (x), y = 0 if |f (x) -y| ≤ ǫ |f (x) -y| -ǫ otherwise.
The ǫ-SVR satisfies the sparsity constraint: Only some α i are non-null which corresponding sample x i are called "Support Vectors" (SVs). Some limitations come from the learning step: With a large training set, the processing time can be very long. Moreover, the problem is exacerbated when several optimizations for parameter selection are considered.

Kernel function

The choice of the kernel function is a crucial step with SVM. A kernel function is a similarity measure between two samples and corresponds to a dot product in some feature space. To be an acceptable kernel, the function should be positive semi-definite [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]Chapter 12]. In this work, several kernels were investigated:

• The linear kernel k(x, z) = x, z .

• The inhomogeneous polynomial kernel k(x, z) = x, z + q p , where q ≥ 0 induces a weight of each degree:

x, z + q p = p k=0 p k x, z p-k q k . • The Gaussian kernel k(x, z) = exp -γ x -z 2 . • The spectral kernel k(x, z) = exp -γα(x, z) 2 , α(x, z) = acos x,z
x z . The spectral kernel was first introduced in hyperspectral imagery for the classification purposes [START_REF] Fauvel | Evaluation of kernels for multiclass classification of hyperspectral remote sensing data[END_REF]. It is based on the angle between two spectra. It is a scale invariant kernel which is used in k-nn approaches or spectral unmixing [START_REF] Keshava | Distance metrics and band selection in hyperspectral processing with application to material identification and spectral librairies[END_REF].

Before running the algorithm, some hyperparameters need to be fitted: a) ǫ: Which controls the resolution of the estimation. Large values produce rough approximations while small values produce fine estimations. It can be set using some prior on the noise. b) λ: Which controls the smoothness of the solution. Large values imply nearly linear functions. c) Kernel parameters: γ for the Gaussian kernel, for instance. Cross-validation is used in this work to select the optimal parameters.

Experiments

Datasets

In this paper, real and simulated datasets are used. Real data have been collected during orbit 103 by the French imaging spectrometer OMEGA on board Mars Express Mission. A detailed analysis of this image by an expert led to a surface reflectance model. This model allows the generation of many synthetic spectra with the corresponding parameters: The proportions of CO 2 , water and dust; and the grain sizes of water and CO 2 . Centered multiGaussian noise has been added, its covariance matrix was determined experimentally from the real image. For the validation sake, separate training and testing datasets have been generated. The notations are the following: n (respectively n t ) is the number of spectra from the training data (respectively test data), x i ∈ R d , i ∈ 1, . . . , n denotes the spectra from the training data and y i (p) ∈ R, i ∈ 1, . . . , n is one of the 5 associated parameters (respectively x t j , y t j , j ∈ 1, . . . , n t ). In these experiments, n = 3584, n t = 3528 and the number of spectral bands is d = 184. Each parameter p = 1, . . . , 5 takes a finite number of values y(p) regularly distributed in a given interval of variation. The different realizations of the vector y i are generated by building all possible combinations of the individual component values. In the following, the index p is omitted in y and in its associated functional f :

f p (x) = y(p) is written f (x) = y.

Results

The quality of the estimation is assessed by computing the Normalized Root Mean Square Errors (NRMSE):

NRMSE = v u u t 1 nt n t X i=1 (ŷi -yi) 2 1 nt n t X i=1 (yi -y) 2 with y = 1 nt n t X i=1 yi (1)
where y i is the real value and ŷi the estimated one. It is close to zero when the predicted values are accurate and becomes larger when predictions are poor. SVM regression is compared with GRSIR and PLS. These methods are "training approaches", GRSIR and PLS both first reduce the dimension of the data using x and y, and then perform the estimation in the reduced space. PLS is a linear estimator while GRSIR is a piecewise linear one, for details see [START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF][START_REF] Bernard-Michel | A note on sliced inverse regression with regularizations[END_REF].

From Table 1, it is clear that the worst results 1 in terms of accuracy are obtained with linear algorithms: PLS and SVM with the linear kernel. Best accuracy are obtained with the SVM-Gaussian kernel, followed by SVM-inhomogeneous polynomial kernels with q ≥ 1. However, the last ones require more time to be trained. SVM with spectral or homogeneous polynomial kernel yield the lowest accuracy. Considering the training time, SVM lead to longer processing times, for every kernels. Inversions are performed on the real image using the functional f . Results are given in Fig. 1, for GRSIR and SVM-Gaussian kernel. The histogram of the 1 Results for GRSIR and PLS differ from [START_REF] Bernard-Michel | Inverting hyperspectral images with gaussian regularized sliced inverse regression[END_REF] because the selection of optimal parameters is done in a different way, cross-validation in this article, and the simulated datasets are different. estimated values is also reported and shows a shift between GRSIR and SVM estimates (astrophysical interpretations are given in [START_REF] Bernard-Michel | Retrieval of mars surface physical properties from omega hyperspectral images using regularized sliced inverse regression[END_REF]). However, since no ground data is available, the inversion accuracy of real image is difficult to appreciate.

Discussion

With SVM, the training quality can be assessed by looking at the Support Vectors (SVs). Selection of many SVs may indicate that the optimization problem is hard (but not necessarily badly solved). Unlike statistical methods such as GRSIR, interpretation of SVM's functional estimation is not possible, except for the linear kernel:

f (x) = n i α i x, x i +b = x, n i α i x i + b = x, w + b.
Surprisingly, the vector w is akin to the vector β estimated with GRSIR [START_REF] Bernard-Michel | Inverting hyperspectral images with gaussian regularized sliced inverse regression[END_REF], see Fig. 2.(c). In GRSIR, β is used to reduce the dimension and the estimation f ′ is done in the reduced space: f (x) = f ′ ( x, β ). It explains why GRSIR performs better than linear SVM: The estimation is done in the one-dimensional space spanned by β (respectively w) with a piecewise linear estimator (respectively linear estimator). The precise theoretical comparison is omitted here for the sake of conciseness. As a conclusion of this work, SVM with non-linear kernels produce accurate results for the inversion of hyperspectral images. Experiments on real data images confirm their robustness to the training set and to the high dimensionality of the data. However, the training time can be dramatically increased for some kernel, making the use of SVM in a practical situations sensitive. Current research is now oriented in the definition of kernel handling more efficiently the saturation in the physical model.

  where x i are samples from the training set, k is a kernel function and (α i ) i=1,...,n , b are the parameters of f which are found during the training process. The kernel k is used to produce non-linear functions. Given a training set, (x i , y i ) i=1,...,n ∈ R d × R, the training of an ǫ-SVR entails the following optimization problem:

Fig. 1 :

 1 Fig. 1: Proportion of dust estimated by (a) GRSIR, (b) SVM with a Gaussian kernel from the hyperspectral image observed from orbit 103 and (c) is the histogram of (a) and (b).

Fig. 2 .

 2 (a) shows the proportion of y corresponding to SVs (light gray) from the training set (dark gray). The fraction increases as y increases too. To understand why the problem becomes harder, the spectral variance of the simulated spectra as a function of the wavelength is reported in Fig.2.(b). It is computed as the trace of the covariance matrix of samples x for which y is constant. It can be seen that a high y value leads to: on the one hand a small variance of the simulated spectra, see Fig.2.(b), and, on the other, a high number of SVs, see Fig.2.(a). This phenomenon is explained by the properties of the physical model: For some values, the model saturates. Consequently when the saturation is reached, increasing the values of the parameters has little consequences on the simulated spectra (and the variance remains low) thus making the estimation more difficult, due to strong non-linearity between x and y. The main practical implication is in the generation of the training samples: Looking at the SVs after a first training provides information to possibly generate new samples for which the estimation is difficult. Finally, note that SVM and other training methods face difficulties in estimating f when saturation is reached.

Fig. 2 :

 2 Fig. 2: (a) Proportion of y (light gray) corresponding to SVs (Gaussian kernel) from the training set (dark gray), y was "Prop. of dust"; (b) Spectral variance of simulated x as function of y, y(p) was "Prop. of dust" ; (c) GRSIR axis β and SVM (linear kernel) normal vector w as function of the wavelength λ.

Table 1 :

 1 NRMSE and computing time for GRSIR, PLS and SVM with various kernels. "x-Pol" is p = x in the polynomial kernel. The power of the polynomial kernel was fixed to 9 for each parameter, after cross-validation. The bottom line of the table corresponds to the training time of parameter "Prop. of H 2 O" after the selection of optimal hyperparameters.

	Parameter	GRSIR	PLS	lin.	Gauss.	Spect.	SVM 0-Pol	0.5-Pol	1-Pol	2-Pol
	Prop. of H 2 O	0.28	0.32	0.31	0.14	0.25	0.24	0.17	0.14	0.13
	Prop. of CO 2	0.19	0.31	0.30	0.15	0.27	0.27	0.18	0.16	0.15
	Prop. of dust	0.11	0.22	0.22	0.09	0.19	0.19	0.11	0.10	0.10
	Grain size of H 2 O	0.34	0.39	0.39	0.15	0.34	0.33	0.23	0.19	0.18
	Grain size of CO 2	0.16	0.24	0.25	0.11	0.21	0.20	0.14	0.12	0.11
	CPU time (s)	0.16	0.66	3.57	10.30	5.89	5.98	10.20	60.30	478

* This work is supported by a contract with CNES through its Groupe Système Solaire Program and by INRIA and with the financial support of the "Agence Nationale de la Recherche" (French Research Agency) through its MDCO program ("Masse de Données et COnnaissances"). The Vahiné project was selected in 2007 under the reference ANR-07-MDCO-013.


S Khardi
email: salah.khardi@inrets.fr

F Nahayo

M Haddou

The Trust Region Sequential Quadratic Programming method applied to two-aircraft acoustic optimal control problem

Keywords: TRSQP algorithm, optimal control problem, aircraft noise levels, AMPL programming, KNITRO

The Trust Region Sequential Quadratic Programming method applied to two-aircraft acoustic optimal control problem

Salah Khardi, Fulgence Nahayo, Mounir Haddou

1. Introduction. Considering the current trend in the eld of air transport [START_REF] Peyrat-Armandy | Les avions de transport modernes et futurs[END_REF][START_REF] Martin | L'Analyse des Nuisances sonores autour des aéroports[END_REF][START_REF] Khardi | Reduction of commercial aircraft noise emission around airports. A new environmental challenge[END_REF], economic and environmental considerations related to the rising cost of oil and the need to preserve the environment impose more severe constraints on the next generation of aircraft [START_REF] Filippone | Comprehensive analysis of transport aircraft ight performance[END_REF]. In sight of reaching one of the 2020 Advisory Council for Aeronautics Reasearch in Europe 'ACARE' objectives [START_REF] Ventre | Les challenges environnementaux pour le transport aérien[END_REF], that consists of the reduction of the environmental pollution and noise impact, ACARE requires a 50% reduction of perceived noise for 2020. This goal represents a dicult scientic and engineering challenge as this requires aerodynamic models and mathematical optimization [START_REF] Zapolozhets | Predicted Flight Procedures for Minimum Noise Impact[END_REF][START_REF] Barriety | La voie de la science, Des outils pour optimiser la conception des avions en phase d'approche[END_REF]. Some work addressing this problem has been carried out. The majority of this work addresses the problem of minimization of aircraft noise around the airport by considering a single plane [START_REF] Abdallah | Minimisation des bruits des avions commerciaux sous contraintes physiques et aérodynamique[END_REF][START_REF] Abdallah | Optimization of operational aircraft parameters reducing noise emission[END_REF]. The other work concerns the stochastic conict detection for air trac management [START_REF] Blin | Stochastic conict detection for air trac management[END_REF], the dynamics of ight [START_REF] Boier | The Dynamics of Flight, The Equations. SUPAÉRO[END_REF] and the comprehensive analysis of transport aircraft ight performance [START_REF] Filippone | Comprehensive analysis of transport aircraft ight performance[END_REF].

The aim of this work is the development of a theoretical model of noise optimization while maintaining a reliable evolution of the ight procedures of two commercial aircraft on approach. These aircraft are supposed to land successively on one runway without conict [START_REF] Roux | Modèle de longueur de piste au décollage-atterrissage, Avions de transport civil[END_REF]. It is all about the evolution of ight dynamics and minimization of noise for two similar commercial aircraft landing taking into account the energy constraint. This model is a non-linear and non-convex optimal control. It is governed by a system of ordinary non-linear dierential equations [START_REF] Kokkinis | Classical and relaxed optimization methods for optimal control problems[END_REF]. The 3-D movement of the two planes is described by a system depending on ordinary non-linear dierential equations with mixed constraints. The function to be minimized is the integral describing the overall level of noise emitted by the two aircraft on approach and collected on the ground. We take into account constraints related to joint stability, performance and ight safety.

The problem of optimal control is described and solved by a Trust Region Sequential Quadratic Programming method 'TRSQP' [START_REF] Mathew | A feasible trust-region sequential quadratic programming algorithm[END_REF][START_REF] Bergounioux | Optimisation et Controle des systèmes linéaires, Cours et exercices corrigés[END_REF][START_REF] Helmut | SQP-methods for solving optimal control problems with control and state constraints: adjoint variables, sensitivity analysis and real-time control[END_REF][START_REF] Hait | Optimisation des Systèmes dynamiques hybrides, Application en génie industriel à l'élaboration des scenarios[END_REF][START_REF] Lemarchal | Numerical Optimization: Theorical and pratical Aspects[END_REF]. By using a merit function, a sequential quadratic programming method associated with global trust regions bypasses the non-convex problem. This method is established by following a tangent quadratic problem obtained from the optimality conditions of Karush-Khun-Tucker applied to the problem considering the objective function as the merit function.

The TRSQP methods are suggested as an option by a Nonlinear Interior point Trust Region Optimization solver 'KNITRO' [START_REF] Waltz | KNITRO: An integrated Package for nonlinear optimization[END_REF][START_REF] Plantenga | KNITRO user's Manual, Version 5.2[END_REF] under A Mathematical Programming Modeling Language 'AMPL' [START_REF] Gay | A modelling Language for Mathematical Programming[END_REF][START_REF] Laboratories | A modelling Language for Mathematical Programming[END_REF]. The global convergence properties are analyzed under dierent assumptions on the approximate Hessian. Additional assumptions on the feasibility perturbation technique are used to prove quadratic convergence to points satisfying second-order sucient conditions. Details of the two-aircraft ight dynamic, the noise levels, the constraints, the mathematical model of the two-aircraft acoustic optimal control problem and the trust region sequential quadratic programming method processing are presented in section 2, 3 and 4 while the numerical experiments are presented in the last section.

Mathematical Modelization.

The motion of each aircraft A i , i := 1, 2 is three dimensional analyzed with 3 frames: the landmark (O,

-→ X 1 , -→ Y 1 , -→ Z 1 ), the aircraft frame (G i , -→ X Gi , -→ Y Gi , -→ Z Gi
) and the aerodynamic one

(G i , -→ X ai , -→ Y ai , -→
Z ai ) where i := 1, 2 [START_REF] Blin | Stochastic conict detection for air trac management[END_REF]. The transition between these three frames is shown easily [START_REF] Boier | The Dynamics of Flight, The Equations. SUPAÉRO[END_REF]. In general, the equations of motion of each aircraft are summarized as:

∑ -→ F exti -dmi dt -→ V ai = mid - → V a i dt ∑ -→ M extG i = d dt [I Gi -→ Ω i ] d - → X o dt = d - → X 1 dt + -→ Ω 10 × -→ X (1)
The index i = 1, 2 reects the aircraft. In the system above, -→

F exti represents the external forces acting on the aircraft, m i (t) the mass of the aircraft, v i the airspeed of aircraft, -→ M extG i the moments of each aircraft, J(G i , A i ) the inertia matrix and Ω i the angular rotation of the aircraft. After transformations and simplications, the system takes the following explicit form:

                                                           V ai = 1 mi [-m i gsinγ ai -1 2 ρSV 2 ai C D + (cosα ai cosβ ai + sinβ ai + sinα ai cosβ ai )F xi -dmi dt u i -m i ∆A i u ] β ai = 1 miVa i [m i gcosγ ai sinµ ai + 1 2 ρSV 2 ai C yi + [-cosα ai sinβ ai + cosβ ai -sinα ai sinβ ai ]F yi -dmi dt v i -m i ∆A i v ] α ai = 1 miVa i cosβa i [m i gcosγ ai cosµ ai -1 2 ρSV 2 ai C Li + [-sinα ai + cosα ai ]F zi -dmi dt w i -m i ∆A i w ] ṗi = C AC-E 2 {r i q i (B -C) -Ep i q i + 1 2 ρSlV 2 ai C li + ∑ 2 j=1 F j [y b Mij cosβ mij sinα mij -z b Mij sinβ mij ]} + E AC-E 2 {p i q i (A -B) -Er i q i + 1 2 ρSlV 2 ai C ni + ∑ 2 j=1 F j [x b Mij sinβ mij -y b Mij cosβ mij cosα mij ]} qi = 1 B {-r i p i (A -C) -E(p 2 i -r 2 i ) + 1 2 ρSlV 2 ai C mi + ∑ 2 j=1 F j [z b Mij cosβ mij cosα mij -x b Mij cosβ mij sinα mij ]} ṙi = E AC-E 2 {r i q i (B -C) + Ep i q i + 1 2 ρSlV 2 ai C li + ∑ 2 j=1 F j [y b Mij cosβ mij sinα mij -z b Mij sinβ mij ]} + A AC-E 2 {p i q i (A -B) -Er i q i + 1 2 ρSlV 2 ai C ni + ∑ 2 j=1 F j [x b Mij sinβ mij -y b Mij cosβ mij cosα mij ]} ẊGi = V ai cosγ ai cosχ ai + u w ẎGi = V ai cosγ ai sinχ ai + v w ŻGi = -V ai sinγ ai + w w φi = p i + q i sinϕ i tanθ i + r i cosϕ i tanθ i θi = q i cosϕ i -r i sinϕ i ψi = sinφi cosθi q i + cosφi cosθi r i ṁi = -1 g dWi dt
(2) where j means the engine index, the expressions A = I xx , B = I yy , C = I zz , E = I xz are the inertia moments of the aircraft, ρ is the air density, S is the aircraft reference area, l is the aircraft reference length, g is the acceleration due to gravity,

C D = C D0 + kC 2 L is the drag coecient, C yi = C yβ β + C yp pl V + C yr rl V + C Y δ l δ l + C Y δn δ n is the lateral forces coecient, C Li = C Lα (α a -α a0 ) + C Lδm δ m + C LM M + C Lq q b a l V is the lift coecient, C li = C lβ β + C lp pl V + C lr rl V + C lδ l δ l + C lδn δ n is the rolling moment coecient, C mi = C m0 + C mα (α -α 0 ) + C mδm δ m is the pitching moment coecient, C ni = C nβ β +C np pl V +C nr rl V +C nδ l δ l +C nδ n δ n is the yawing moment coecient, (x b M ij , x b M ij , x b M ij )
is the position of the engine in the body frame,

F = (F xi , F yi , F zi ) is the propulsive force, V ai = (u i , v i , w i ) is the aerodynamic speed, (∆A i u , ∆A i v , ∆A i w )
is the complementary acceleration, (u w , v w , w w ) is the wind velocity, β mij is the yaw setting of the engine and α mij is the pitch setting of the engine. The expressions α ai (t), β ai (t), θ i (t), ψ i (t), ϕ i (t), V ai (t), X Gi (t), Y Gi (t), Z Gi (t), p i (t), q i (t), r i (t), m i (t) are respectively the attack angle, the aerodynamic sideslip angle, the inclination angle, the cup, the roll angle, the airspeed, the position vectors, the roll velocity of the aircraft relative to the earth, the pitch velocity of the aircraft relative to the earth, the yaw velocity of the aircraft relative to the earth and the aircraft mass.

Transforming the system (2) in state function, one has:

dy i (t) dt = f i (y i (t), u i (t)), i = 1, 2 (3)
where the state vector is:

y i (t) : [t 0 , t f ] -→ R 13 y i (t) = (α ai (t), β ai (t), θ ai (t), ψ ai (t), ϕ i (t), V ai (t), X Gi (t), Y Gi (t), Z Gi (t), p i (t), q i (t), r i (t), m i (t)) (4)
The control vector is

u i (t) : [t 0 , t f ] -→ R 4 t -→ u i (t) = (δ li (t), δ mi(t) , δ ni (t), δ xi (t)) (5)
where the expressions δ li (t), δ mi(t) , δ ni (t), δ xi (t) are respectively the roll control, the pitch control, the yaw control and the thrust one. The dynamics relationship can be written as:

ẏi (t) = f i (y i (t), u i (t), t), ∀t ∈ [0, T ], y i (0) = y i0 (6)
The angles γ ai (t), χ ai (t), µ ai (t) corresponding respectively to the aerodynamic climb angle (air-path inclination angle), the aerodynamic azimuth (air-path track angle) and the air-path bank angle (aerodynamic bank angle) are not taken as state in this model.

To simplify the model, the atmosphere standards conditions are considered. The engine angles, the complementary acceleration and the aerodynamic sideslip angle are negligible because the wind is constant and there is no engine failure. With some complex mathematical transformations, the dynamic system (2) becomes:

                                       V ai = 1 mi [-m i gsinγ ai -1 2 ρSV 2 ai C D + (cosα ai + sinα ai )F xi -dmi dt u i ] α ai = 1 miVa i cosβa i [m i gcosγ ai cosµ ai -1 2 ρSV 2 ai C Li + [cosα ai -sinα ai ]F zi -dmi dt w i ] ṗi = C AC-E 2 {r i q i (B -C) -Ep i q i + 1 2 ρSlV 2 ai C li } + E AC-E 2 {p i q i (A -B) -Er i q i + 1 2 ρSlV 2 ai C ni ṙi = E AC-E 2 {r i q i (B -C) + Ep i q i + 1 2 ρSlV 2 ai C li + A AC-E 2 {p i q i (A -B) -Er i q i + 1 2 ρSlV 2 ai C ni } ẊGi = V ai cosγ ai cosχ ai ẎGi = V ai cosγ ai sinχ ai ŻGi = -V ai sinγ ai φi = p i + q i sinϕ i tanθ i + r i cosϕ i tanθ i θi = q i cosϕ i -r i sinϕ i ψi = sinφi cosθi q i + cosφi cosθi r i ṁi = -1 g dWi dt (7)
By the combination of this system with the aircraft control, one has the two-aircraft dynamic ight model as shown in [START_REF] Zapolozhets | Predicted Flight Procedures for Minimum Noise Impact[END_REF]. The state vector is

y i (t) : [t 0 , t f ] -→ R 12 y i (t) = (α ai (t), θ ai (t), ψ ai (t), ϕ ai (t), V ai (t), X Gi (t), Y Gi (t), Z Gi (t), p i (t), q i (t), r i (t), m i (t)) (8)
This will be added to the cost function and constraint function for the aircraft optimal control problem as shown in the following paragraphs.

The ob jective function model. Many noise indices (LAmax, L eq , EP N L, ...) exist and are presented in the open literature [START_REF] Abdallah | Minimisation des bruits des avions commerciaux sous contraintes physiques et aérodynamique[END_REF][START_REF] Martin | Noise monitoring in the vicnity of airports[END_REF]. In this paper, the Sound Exposure Level 'SEL' is considered. It provides a comprehensive way to describe noise events in modeling and analyzing noise environment impacts [START_REF] Harris | How do we Describe Aircraft Noise[END_REF]:

SEL = 10log [ 1 t o ∫ t ′ 10 0.1L A1,dt (t) dt ] (9)
where t o is the time reference taken equal to 1 s and t ′ the noise event interval. [t 10 , t 1f ] and [t 20 , t 2f ] are the respective approach intervals for the rst and the second aircraft, the objective function is calculated as:

SEL 1 = 10log [ 1 to ∫ t20 t10 10 0.1L A1,dt (t) dt ] , t ∈ [t 10 , t 20 ] SEL 12 = SEL 11 ⊕ SEL 21 = 10 log[ 1 to ∫ t 1f t20 10 0.1L A1,dt (t) dt + 1 to ∫ t 1f t20 10 0.1L A2,dt (t) dt], t ∈ [t 20 , t 1f ] SEL 2 = 10 log [ 1 to ∫ t 2f t 1f 10 0.1L A2,dt (t) dt ] , t ∈ [t 1f , t 2f ] SEL G = (t20-t10)SEL1⊕(t 1f -t20)SEL12⊕(t 2f -t 1f )SEL2 t 2f -t10 = 10 log{ 1 t 2f -t10 [(t 20 -t 10 ) ∫ t20 t10 10 0.1LA1(t) dt + (t 1f -t 20 ) ∫ t 1f t20 10 0.1LA1(t) dt + (t 1f -t 20 ) ∫ t 1f t20 10 0.1LA2(t) dt +(t 2f -t 1f ) ∫ t 2f t 1f 10 0.1LA2(t) dt, ]}, t ∈ [t 10 , t 2f ] (10
) where SEL G is the cumulated two-aircraft noise and the operator ⊕ means the convoluted addition of noises. Expressions L A1 (t), L A2 (t) are equivalent and reect the aircraft jet noise given by the formula [START_REF] Abdallah | Minimisation des bruits des avions commerciaux sous contraintes physiques et aérodynamique[END_REF][START_REF] Zola | An improved prediction method for noise generated by conventional prole coaxial jets[END_REF]:

L A1 (t) = 141 + 10 log ( ρ 1 ρ ) w + 10 log ( V e c ) 7.5 + 10 log s 1 + 3 log ( 2s 1 πd 2 1 + 0.5 ) + 5 log τ 1 τ 2 + 10 log      ( 1 - v 2 v 1 ) me + 1.2 ( 1 + s 2 v 2 2 s 1 v 2 1 ) 4 ( 1 + s 2 s 1 ) 3      -20 log R + ∆V + 10 log [ ( ρ ρ ISA ) 2 ( c c ISA ) 4 ]
where v 1 is the jet speed at the entrance of the nozzle, v 2 the jet speed at the nozzle exit, τ 1 the inlet temperature of the nozzle, τ 2 the temperature at the nozzle exit, ρ the density of air, ρ 1 the atmospheric density at the entrance of the nozzle, ρ ISA the atmospheric density at ground, s 1 the entrance area of the nozzle hydraulic engine, s 2 the emitting surface of the nozzle hydraulic engine, d 1 the inlet diameter of the nozzle hydraulic engine,

V e = v 1 [1 -(V /v 1 ) cos(α p )] 2/3
the eective speed (α p is the angle between the axis of the motor and the axis of the aircraft), R the source observer distance, w the exponent variable dened by: w = 3(V e /c) 3.5 0.6 + (V e /c) 3.5 -1, c the sound velocity (m/s), m the exhibiting variable depending on the type of aircraft: me = 1.1

√ s 2 s 1 ; s 2 s 1 < 29.7, me = 6.0; s 2 s 1 ≥ 29.7, the term ∆V = -15log(C D (M c , θ)) -10log(1 -M cosθ), means the Doppler convection when C D (M c , θ) = [(1 + M c cosθ) 2 + 0.04M 2 c ], M the aircraft Mac Number, M c the convection Mac Number: M c = 0.62(v 1 -V cos(α p )) c, θ is the Beam angle.
Formula above leads to the objective function

J G12 (y(t), u(t), i = 1, 2) = ∫ t ′ g(y(t), u(t), t, i = 1, 2)dt.
Constraints.

The two-aircraft optimal control problem which minimizes noise and provides optimal trajectory must be done in realistic ight domains. Operational procedures are performed while respecting dynamic limits related to the safety of ight and operational modes of the aircraft. This concerns aircraft airspeed, altitude, attack angle, aerodynamic sideslip angle, inclination angle, roll angle, throttle control position, roll control position, pitch control position, yaw control position, energy constraint, vertical, longitudinal and lateral aircrafts separation, roll velocity of aircraft relative to the earth, pitch velocity of aircraft relative to the earth, yaw velocity of aircraft relative to the earth and the aircraft mass.

1. The vertical separation given by Z G12 = Z G2 -Z G1 where Z G1 , Z G2 are respectively the altitude of the rst and the second aircraft and Z G12 the altitude separation.

2. The horizontal separation X G12 = X G1 -X G2 [START_REF] Dgac | Mémento à l'usage des utilisateurs des procédures d'approche et de départ aux instruments[END_REF][START_REF] Sors | Séparation et contrôle aérien[END_REF][START_REF] Dgac | Méthodes et minimums de séparations des aéronefs aux procédures[END_REF] where X G1 , X G2 are horizontal positions of the rst and the second aircraft and their separation distance.

3. The aircraft speed V ai must be bounded as follows 1.3V s ≤ V ai ≤ V f where V s is the stall speed,V f is the maximum speed and V o the minimum speed of the aircraft A i [START_REF] Dominique | Cisaillement de vent ou Windshear[END_REF][START_REF] Roux | Pour une approche analytique de la dynamique du[END_REF], the roll velocity of the aircraft relative to the earth p i ∈ [p 0 , p f ], the pitch velocity of the aircraft relative to the earth q i ∈ [q 0 , q f ] and the yaw velocity of the aircraft relative to the earth r i ∈ [r 0 , r f ] .

4. On the approach, the ICAO standards and aircraft manufacturers require ight angle evolution as follows: attack angle

α ai ∈ [α o , α f ], the inclination angle θ i ∈ [θ 0 , θ f
] and the roll angle

ϕ i ∈ [ϕ o , ϕ f ].
5. The aircraft control δ(t) = (δ li (t), δ mi(t) , δ ni (t), δ xi (t)) keeps still between the position δ l0 and δ lf for the roll control, δ m0 and δ mf for the pitch control, δ n0 and δ nf for the yaw control and δ xo and δ xf for the thrust.

The mass m i of the aircraft

A i is variable: m o < m i < m f , i = 1, 2.
This constraint results in energy consumption of the aircraft [START_REF] Boier | Dynamique de vol de l'avion. SupAéro, Départements des Aéronefs[END_REF][START_REF] Ifrance | Fiches techniques[END_REF].

On the whole, the constraints come together under the relationship:

C(y i (t), u i (t), t) ≤ 0 (11)
where

C(t) : R 12 × R 4 × R -→ R 17 , (y i (t), u i (t), t) -→ C(y i (t), u i (t), t), with C i (t) = (α i (t), θ i (t), ψ i (t), ϕ i (t), V ai (t), X Gi (t), Y Gi (t), Z Gi (t), p(t), q(t), r(t), δ li (t), δ mi (t), δ ni (t), δ xi (t), m i (t), t).
The following values reect the digital applications considered for the two-aircraft [START_REF] Abdallah | Minimisation des bruits des avions commerciaux sous contraintes physiques et aérodynamique[END_REF][START_REF] Boier | The Dynamics of Flight, The Equations. SUPAÉRO[END_REF][START_REF] Roux | Pour une approche analytique de la dynamique du[END_REF][START_REF] Boier | Dynamique de vol de l'avion. SupAéro, Départements des Aéronefs[END_REF].

V a1 = V a2 = 180 m/s V a1 = V a2 = 73.45 m/s
The A1 Aircraft altitude

Z 1 = 35 × 10 2 m Z 1 = 0 m
The A2 Aircraft altitude

Z 2 = 41 × 10 2 m Z 2 = 0 m
The aircraft roll control

δ l1 = δ l2 = 0.0174 δ l1 = δ l2 = -0.0174
The pitch control

δ m1 = δ m2 = 0.087 δ m1 = δ m2 = 0
The yaw control

δ n1 = δ n2 = 0.314 δ n1 = δ n2 = -0.035
The thrust control

δ x1 = δ x2 = 0.6 δ x1 = δ x2 = 0.2
The attack angle

α a1 = α a2 = 12 • α a1 = α a2 = 2 •
The inclination angle

θ a1 = θ a2 = 7 • θ a1 = θ a2 = -7 •
The air-path inclination angle

γ a1 = γ a2 = 0 • γ a1 = γ a2 = -5 •
The aerodynamic bank angle

µ a1 = µ a2 = 3 • µ a1 = µ a2 = -2 •
The air-path azimuth angle

χ a1 = χ a2 = 5 • χ a1 = χ a2 = -5 •
The roll angle

ϕ a1 = ϕ a2 = 1 • ϕ a1 = ϕ a2 = -1 •
The cup

ψ a1 = ψ a2 = 3 • ψ a1 = ψ a2 = -3 •
The limits of time

t 1f = 600 s,t 2f = 645 s t 10 = 0 s, t 20 = 45 s
The mass of the A1 Aircraft The Aircraft vertical separation

Z 12 = 2 × 10 3 f t ≃ 6 × 10 2 m
The Aircraft longitudinal separation

X G12 = 5 N M ≃ 9 × 10 3 m
The Aircraft roll velocity relative to the earth

p 1 = p 2 = 1 • s -1 p 1 = p 2 = -1 • s -1
The Aircraft pitch velocity relative to the earth q 1 = q 2 = 3.6 • s -1

q 1 = q 2 = 3 • s -1
The Aircraft yaw velocity relative to the earth

r 1 = r 2 = 12 • s -1 r 1 = r 2 = -12 • s -1
The two-aircraft acoustic optimal control problem.

The combination of the aircraft dynamic equation ( 3) and ( 7), the aircraft objective function from equations [START_REF] Blin | Stochastic conict detection for air trac management[END_REF] and the the aircraft constraints function [START_REF] Boier | The Dynamics of Flight, The Equations. SUPAÉRO[END_REF], the two-aircraft acoustic optimal control problem is given as follows:

           min (y,u)∈Y×U J G12 (y(t), u(t)) = ∫ t 1f t10 g 1 (y 1 (t), u 1 (t), t)dt+ ∫ t 1f t20 g 12 (y 1 (t), u 1 (t), y 2 (t), u 2 (t), t)dt + ∫ t 2f t20 g 2 (y 2 (t), u 2 (t), t)dt + ϕ(y(t f )) ẏ(t) = f (u(t), y(t)), u(t) = (u 1 (t), u 2 (t)) y(t) = (y 1 (t), y 2 (t)), ∀t ∈ [t 10 , t 2f ], t 10 = 0, y(0) = y 0 , u(0) = u 0 C(y(t), u(t)) ≤ 0 (12)
where g 12 shows the aircraft coupling noise function and J G12 is the SEL of the two A300-aircraft.

3. The numerical processing. The problem as dened in the relation ( 12) is an optimal control problem with instantaneous constraints. We aim to solve this problem with the Trust Region Sequential Quadratic Programming method. Applying SQP methods [START_REF] Wright | Numerical Optimization[END_REF] , we write the system (12) as:

       min J G12 (x), x = (y(t), u(t)) ẏ = f (x) n j (x) = 0, j ∈ Ξ n j (x) ≥ 0, j ∈ Γ (13)
where the expressions Ξ and Γ are the sets of equality and inequality indices. The function J G12 (x), f (x), n(x) must be twice continuously dierentiable. The Lagrangian of the system ( 13) is dened by the function L(x, λ) = J GP 12 (x) + λ T [b( ẏ, x) + n(x)] where the vector λ is the Lagrange multiplier and b( ẏ, x) = ẏ -f (x) = 0. Considering the feasible points of (12), one transforms the system (13) into a quadratic problem. A SQP method solves a succession of quadratic problems. The mathematical formulation of sub-problems obtained at the k-th step ∆x k is the following:

       min ∆x k [K G12 (x k )] = ∇ T J G12 (x k )∆x k + 1 2 (∆x k ) T H k ∆x k ∇ T b( ẏk , x k )∆x k + b( ẏk , x k ) = 0 ∇ T n Ξ (x k )∆x k + n Ξ (x k ) = 0 ∇ T n Γ (x k )∆x k + n Γ (x k ) ≥ 0 (14)
The vector ∆x k is a primal-dual descent direction, H k = ∇ 2 L(x k , λ k ) is the Hessian matrix of Lagrangian L from system ( 13) and K G12 (x k ) the quadratic model. The estimation of gradients is, in principle, calculated by nite dierences or the calculation of the adjoint systems for problems with many parameters and nally by the sensitivity analysis. This last technique is very eective in the case of a large number of variables with few parameters [START_REF] Helmut | SQP-methods for solving optimal control problems with control and state constraints: adjoint variables, sensitivity analysis and real-time control[END_REF][START_REF] Hait | Optimisation des Systèmes dynamiques hybrides, Application en génie industriel à l'élaboration des scenarios[END_REF]. The SQP method is a qualied local method. Its convergence is quadratic if the rst iterate is close to a solution ỹ satisfying the sucient optimality conditions [START_REF] Walter | SNOPT: An SQP Algorithme for large-scale constrained optimization[END_REF][START_REF] Gilbert | SQPpro: A solver of nonlinear optimization problems, using an SQP approach[END_REF][START_REF] Bricker | SQP: Sequential quadratic programming[END_REF]. This algorithm above must be transformed because the two-Aircraft problem is non-convex. For improving the robustness and global convergence behavior of this SQP algorithm, it must be added with the trust radius of this form:

||D∆x k || p ≤ ∆, p ∈ [1, ∞] (15)
where D is uniformly bounded. The relations ( 14) and ( 15) form a quadratic program when p = ∞. So, the trustregion constraint is restated as -∆e ≤ Dx ≤ ∆e, e = (1, 1, 1, ..., 1) T . If p = 2, one has the quadratic constraint

∆x T k D T D∆x k ≤ ∆ 2 .
In the following, we develop the convergence theory for any choice of p just to show the equivalence between the ||.|| p and ||.|| 2 . By the combination of some relation of ( 13) and the relation ( 14), all the components of the step are controlled by the trust region. The two-aircraft problem takes the following form

           min ∆x k [K G12 (x k )] = ∇ T J G12 (x k )∆x k + 1 2 (∆x k ) T H k ∆x k ∇ T b( ẏk , x k )∆x k + b( ẏk , x k ) = 0 ∇ T n Ξ (x k )∆x k + n Ξ (x k ) = 0 ∇ T n Γ (x k )∆x k + n Γ (x k ) ≥ 0 ||D∆x k || p ≤ ∆, p ∈ [1, ∞] (16)
In some situations, all of the components of the step are not controlled by the trust region because of some hypotheses on D. There is an other alternative which allows the practical SQP methods by using the merit function or the penalty function to measure the worth of each point x.

Several approaches like Byrd-Omojokun and Vardi approaches exist to solve the system (13) [START_REF] Xavier-Jonsson | Méthodes des points intérieurs et de régions de conence en optimisation non-linéaire et application à la conception des verres ophtamiques progressifs[END_REF]. It can also be solved with the KNITRO, the SNOPT and other methods [START_REF] Ouriemchi | Résolution de problèmes non-linéaires par les méthods des points intérieurs[END_REF]. In the latter case, we have an ordinary dierential system of non-linear and non-convex equations. The uniqueness of the solution of the quadratic sub-problem is not guaranteed. It therefore combines the algorithm with a merit function for judging the quality of the displacement. The merit function can therefore oer a way to measure all progress of iterations to the optimum while weighing the importance of constraints on the objective function. It is chosen in l 2 norm particularly the increased Lagrangian L I because of its smooth character. So, in the equation above, one replaces L by L I . Thus, this transforms the SQP algorithm in sequential quadratic programming with trust region globalization 'TRSQP'. Its principle is that each new iteration must decrease the merit function of the problem for an eligible trust radius. Otherwise, we reduce the trust radius ∆x K for computing the new displacement. A descent direction is acceptable if its reduction is emotionally positive. The advantages of the method are that the merit function will circumvent the non-convexity of the problem. This approach shows that only one point is sucient to start the whole iterative process [START_REF] Gilbert | Eléments d'optimisation diérentiable -Théorie et Algorithmes[END_REF][START_REF] Nahayo | Les méthodes numériques appliquées en optimisation non-linéaire et en commande optimale[END_REF][START_REF] Khardi | Mathematical Model for Advanced CDA and Takeo Procedures Minimizing Aircraft Environmental Impact[END_REF].

Meanwhile, we use an algorithm called feasibility perturbed SQP in which all iterates x k are feasible and the merit function is the cost function. Let us consider the perturbation ∆x k of the step ∆x k such that 1. The relation

x + ∆x k ∈ F ( 17
)
where F is the set of feasible points for (12),

The asymptotic exactness relation

||∆x -∆x k || 2 ≤ ϕ(||∆x k || 2 )||∆x k || 2 (18)
is satised where ϕ : R + -→ R + with ϕ(0) = 0.

These two conditions are used to prove the convergence of the algorithm and the eectiveness of this method. The advantages gained by maintaining feasible iterates for this method are:

• An inequality constraint n j is active at point x = (y * , u * ) if n j (x) = 0. Γ(x) = Γ * is the set of indices j corresponding to active constraints in x,

Γ + * = {j ∈ Γ * |(λ * Γ ) j > 0} Γ 0 * = {j ∈ Γ * |(λ * Γ ) j = 0} (19)
where the constraints of index Γ + * are highly active and those of Γ 0 * weakly active. • An element x ∈ Γ * veries the condition of qualifying for the constraints n if the gradients of active constraint ∇n Ξ (x), ∇n Γ (x) are linearly independent. This means that the Jacobian matrix of active constraints in x is full.

• An element x ∈ Γ * satises the qualication condition of Mangasarian-Fromowitz for constraints n in x if there exists a direction d such that

∇n Ξ (x) T d = 0∇n j (x) T d < 0∀j ∈ Γ(x) (20)
where the gradients {∇n(x)} are linearly independent.

The Karush-Kuhn-Tucker optimality conditions are obtained by considering that J, n functions of C 1 class and

x a solution of the problem ( 12) which satises a constraints qualication condition. So,there exists λ * such that:

∇ y L(x, λ * ) = 0, n Ξ (x) = 0, n Γ (x) ≤ 0, λ * Γ ≥ 0, λ * Γ n Γ (x) = 0 (21)
These equations are called the conditions of Karush-Kuhn-Tucker(KKT). The rst equation reects the optimality, the second and third the feasibility conditions. The others reect the additional conditions and Lagrange multipliers corresponding to inactive constraints nj(x) are zero. The couple (x, λ * ) such that the KKT conditions are satised is called primal-dual solution of [START_REF] Waltz | KNITRO: An integrated Package for nonlinear optimization[END_REF]. So, x is called a stationary point.

For the necessary optimality conditions of second order [START_REF] Bergounioux | Optimisation et Controle des systèmes linéaires, Cours et exercices corrigés[END_REF], taking x a local solution of [START_REF] Waltz | KNITRO: An integrated Package for nonlinear optimization[END_REF] and satisfying a qualication condition, then there exist multipliers (λ * ) such that the KKT conditions are veried . So we have

∇ 2 xx L(x, λ * )d.d > 0∀h ∈ C * where C * is a critical cone dened by C * = {h ∈ Y × U : ∇n j (x).h = 0 ∀j ∈ Ξ ∪ Γ + * , ∇n j (x).h ≤ 0∀j ∈ Γ 0 * }. The elements of C * are called critical directions.
For the sucient optimality conditions of second order [START_REF] Bergounioux | Optimisation et Controle des systèmes linéaires, Cours et exercices corrigés[END_REF], suppose that there exists (λ * ) which satisfy the KKT conditions and such that ∇ 2 xx L(x, λ * )d.d > 0∀h ∈ C * \{0}. So x is a local minimum of [START_REF] Roux | Modèle de longueur de piste au décollage-atterrissage, Avions de transport civil[END_REF]. [START_REF] Filippone | Comprehensive analysis of transport aircraft ight performance[END_REF] The TRSQP algorithm and convergence analysis. Assume that for a given SQP step ∆x k and its perturbation ∆x k , the ratio to predict decrease is

r k = J G12 (x k ) -J G12 (x k + ∆x k ) -K G12 ( ∆x k ) (22)
The two-aircraft acoustic optimal control TRSQP algorithm is written as:

1. Let x 0 a given starting point, ∆ ≥ 1 the trust region upper bound, ∆ 0 ∈ (0, ∆) an initial radius, ϵ ∈ [ϵ 0 , ϵ f ) and p ∈ [1, ∞]
2. Calculate ∆x k by solving the system

           min ∆x k [K G12 (x k )] = ∇ T J G12 (x k )∆x k + 1 2 (∆x k ) T H k ∆x k ∇ T b( ẏk , x k )∆x k + b( ẏk , x k ) = 0 ∇ T n Ξ (x k )∆x k + n Ξ (x k ) = 0 ∇ T n Γ (x k )∆x k + n Γ (x k ) ≥ 0 ||D∆x k || p ≤ ∆, p ∈ [1, ∞]
Seek also ∆x k by using the system

x + ∆x k ∈ F ||∆x -∆x k || 2 ≤ ϕ(||∆x k || 2 )||∆x k || 2
3. If no such for the perturbed counterpart ∆x k is found, the following aectations are considered.

∆x k+1 ← ( 1 2 )||D k ∆x k || p x k+1 ← x k ; D k+1 ← D k ; 4. Otherwise, calculate r k = JG12(x k )-JG12(x k + ∆x k ) -KG12( ∆x k ) ; if r k ≤ ϵ f , ∆ k+1 ← ( 1 2 )||D k ∆x k || p ; else if r k > a 0 × ϵ 0 and ||D k ∆x k || p = ∆ k ∆ k+1 ← min(2∆ k , ∆); else ∆ k+1 ← ∆ k ; 5. If r k > ϵ x k+1 ← x k + ∆x k ; Choose the new matrix D k+1 ; else x k+1 ← x k ; D k+1 ← D k ;
6. end.

At each major iteration a positive denite quasi-Newton approximation of the Hessian of the Lagrangian function, H, is calculated using the BFGS method, where λ i , i = 1, ..., m, is an estimate of the Lagrange multipliers.

H k+1 = H k + q k q T k q T k s k - H T k s T k s k H k s T k H k s k
where

s k = x k+1 -x k , q k = (∇J G12 (x k+1 + ∑ n j=1 λ j .∇n(x k+1 ) + b(x k+1 )) -(∇J G12 (x k + ∑ n j=1 λ j .∇n(x k ) + b(x k ))
A positive denite Hessian is maintained providing q T k s k is positive at each update and that H is initialized with a positive denite matrix. This algorithm is implemented by AMPL language programming and the KNITRO solver [START_REF] Waltz | KNITRO: An integrated Package for nonlinear optimization[END_REF][START_REF] Plantenga | KNITRO user's Manual, Version 5.2[END_REF].

Analysis of the algorithm and its convergence.

Let us dene the set F 0 as follows:

F 0 = {x|∇ T b( ẏ, x)∆x + b( ẏ, x) = 0, ∇ T n Ξ (x)∆x + n Ξ (x) = 0, ∇ T n Γ (x)∆x + n Γ (x) ≥ 0, J G12 (x) ≤ J G12 (x 0 )} ∈ F The trust-region bound ||D∆x k || p ≤ ∆, p ∈ [1, ∞]
species the following assumption.

1. There exists a constant β such that for all points x ∈ F 0 and all matrix D used in the algorithm, we have for any ∆x satisfying the following equations

∇ T b( ẏ, x)∆x + b( ẏ, x) = 0, ∇ T n Ξ (x)∆x + n Ξ (x) = 0, ∇ T n Γ (x)∆x + n Γ (x) ≥ 0 that β -1 ||∆x|| 2 ≤ ||D∆x|| p ≤ β||∆x|| 2 (23)
2. The level set F 0 is bounded and the functions J G12 , b, η are twice continuously dierentiable in an open neighborhood M(F 0 ) of this set.

Under certain assumptions as shown in [START_REF] Mathew | A feasible trust-region sequential quadratic programming algorithm[END_REF], this algorithm is well dened.

In this paragraph, one wants to prove that the algorithm has a convergence to stationary point of [START_REF] Kokkinis | Classical and relaxed optimization methods for optimal control problems[END_REF]. If we consider that all assumptions hold for each feasible point x for [START_REF] Roux | Modèle de longueur de piste au décollage-atterrissage, Avions de transport civil[END_REF], the Mangasarian-Fromowitz are satised for constraints. After all, the KKT optimality conditions are specied and that shows that there is at least a local convergence. With other added conditions as shown in [START_REF] Mathew | A feasible trust-region sequential quadratic programming algorithm[END_REF], the global convergence is held.

5. Numerical Results. For this result, one considers the 20 km radius for the rst aircraft around the airport when its maximum altitude is 1 km and 1.5 km for the second on 29 km radius. Figures 1, 2, 3 and4 show noise levels around the airport when the optimization is applied or not. This explains the importance of optimizing the acoustics in aircraft approach and the gains brought by this model when compared to what is done daily. As specied by this gure, the noise level increases and is maximum when the observation point lies below the aircraft. After all, noise level decreases gradually as the aircraft departs increasingly from the observation point. The noise levels increase comparing with the result of the rst gure because aircraft altitudes decrease. As specied by this gure, the noise level increases and is maximum when the observation point lies below the aircraft. After all, noise level decreases gradually as the aircraft departs increasingly from the observation point.The noise levels increase comparing with the result of the second gure because aircraft altitudes decrease. As specied by this gure, the noise level increases and is maximum when the observation point lies below the aircraft. After all, noise level decreases gradually as the aircraft departs increasingly from the observation point.

The general remark on gure 1, 2, 3 and 4 is that all noise levels evolution keeps the same behavior. By comparison, this result is close to standard values of jet noise on approach as shown by Harvey [START_REF] Hubbard | Noise sources and[END_REF]. Figure 5 shows the evolution of the trajectories which reect a continuous descent. The aircrafts' landing procedures are suciently separated. It is obvious that each aircraft follows its optimal trajectory when considering the separation distance. The same constraints on speed are considered, allowing a subsequent landing on the same track. Thus, as recommended by ICAO, the security conditions are met and ight procedures are good as shown by the following results. The maximum altitude considered are 3500 m and 4100 m for the rst and second aircraft. The duration approach considered is 600 s for the rst aircraft and 645 s for the second. Our suggestion is an optimal trajectory for the two-aircraft , not necessarily optimal for each of them , but resulting from the two trajectories combination. 6 shows the principle angles evolution as recommended by ICAO during aircraft landing procedures. As specied by this result, the aircraft roll angles oscillate around zero, the ight-path angles are negative and keep the recommended position for aircraft landing procedures. This is the same for the attack angles. Figure 7 shows the aircraft speed evolution during aircraft landing procedures. For the rst, the aircraft speed decreases and keeps a constant position, increases and keeps again a constant position till the end of the aircraft landing. This evolution remains the same for the speed of the second aircraft. The maximum value is 180 m/s when the minimum is around 124 m/s. Figure 8 shows that the nesse evolution is bang-bang for both aircraft. That means that the control lies constant on one bound of the variation interval. The throttle position keeps constant during all the landing procedures. This is the same for the two-aircraft roll velocity relative to the earth p 1 , p 2 , the two-aircraft pitch velocity relative to the earth q 1 , q 2 and the two-aircraft yaw velocity relative to the earth r 1 , r 2 , the aircraft mass, its roll control, pitch and yaw control and all variation reect the limiters conditions as shown in table 1.

This solution is achievable with KNITRO through the optimality conditions: nal feasibility error (abs/rel) = 3.29e-15 / 8.51e-18, nal optimality error (abs/rel) = 1.04e-13 / 1.04e-13, total program time (secs) = 17738,36914 (14893,13 CPU time), time spent in evaluations (secs) = 9815,84082). 6. Conclusion. We developed a numerical solution of an optimal control problem in the case of two-aircraft on approach. Theoretical considerations and practices of the feasible TRSQP algorithm are used for the establishment of a non-linear program, implementing the considered problem. The algorithm minimizes a sequence of merit function using a sub-problem of the quadratic problem at each step for all active constraints to generate a search trust direction for all primal and dual variables. An optimal solution to the discretized problem is found through a local convergence. The results show a reduction of noise at reception points during the approach of the two-aircraft. The obtained trajectories exhibit optimal characteristics and are acoustically eective. Some added conditions are necessary to prove the global convergence of the considered algorithm. Further researches are needed to complete the problem processing. We suggest one optimal trajectory for all the landing aircraft, although it may not necessarily be the optimal one for each aircraft.

m 10 ≃ 1 . 1 ×

 1011 10 5 kg, m 1f ≃ 1.09055 × 10 5 kg, The mass of the A2 Aircraft m 20 ≃ 1.10071 × 10 5 kg m 2f ≃ 1.09126 × 10 5 kg The A300 inertia moments [31] A = 5.555 × 10 6 kg m 2 B = 9.72 × 10 6 kg m 2 C = 14.51 × 10 6 kg m 2 E = -3.3 × 10 4 kg m 2

  Hypothesis:The observation points are taken on the ground under the path of the aircraft and are independent of each other.The observation positions are: (-20000 m, -20000 m, 0 m), (-19800 m, -19800 m, 0 m), ..., (0 m, 0 m, 0 m), for a space step of 200 m for x and y. The touching point on the ground is (0 m, 0 m, 0 m) while the temporal separation of aircraft varies from 45 s to 90 s . At each point of view, it is a vector of N noise levels as shown in the discretization. It is very important to consider the maximum value among the N values, which value corresponded to the shortest distance between the noise source and the observation point.

20 Figure 1 :

 201 Figure 1: Aircraft noise at the indicated reception point

Figure 2 :Figure 2

 22 Figure 2: Aircraft noise at the indicated reception point

Figure 3 :

 3 Figure 3: Aircraft noise at the indicated reception point

Figure 4 :Figure 4

 44 Figure 4: Aircraft noise at the indicated reception point

Figure 5 :

 5 Figure 5: Aircraft optimal ight paths

Figure 6 :

 6 Figure 6: Flight-path angles of the aircraft

Figure

  Figure6shows the principle angles evolution as recommended by ICAO during aircraft landing procedures. As specied by this result, the aircraft roll angles oscillate around zero, the ight-path angles are negative and keep the recommended position for aircraft landing procedures. This is the same for the attack angles.

Figure 7 :

 7 Figure 7: Aircraft speeds

Figure 8 :

 8 Figure 8: Aircraft nesse

Table of

 of

		limit digital values for the two-aircraft in approach phase
	Constraint denomination	maximum value	minimum value
	The Aircraft speed

• The trust region restriction[START_REF] Bergounioux | Optimisation et Controle des systèmes linéaires, Cours et exercices corrigés[END_REF] is added to the SQP problem[START_REF] Mathew | A feasible trust-region sequential quadratic programming algorithm[END_REF] without concern that it will yield an infeasible subproblem.• The objective function J G12 is itself used as a merit function in deciding whether to take a step.• If the algorithm is terminated early, we will be able to use the latest iterate x k as a feasible suboptimal point, which in many applications is far preferable to an infeasible suboptimum.Here are some considerations that are needed for the KKT optimality conditions.

Acknowledgments. This work is supported by the Agence Universitaire de la Francophonie-Région Afrique Centrale. Thanks to The International Joint Research and International Publication Project (Universitas Indonesia -Indonesia and INRETS -France), Memorandum of Understanding between Universitas Indonesia and INRETS, the French Embassy in Indonesia, and DIKTI (Indonesia) for their supports.