F Nahayo

S Khardi
email: salah.khardi@ifsttar.fr

M Haddou

Optimal control of two-commercial aircraft dynamic system during approach. The noise levels minimization

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

This paper aims to reduce noise levels of two-aircraft landing simultaneously on approach. Constraints related to stability, performance and flight safety are taken into account. The problem of optimal control is described and solved by a Sequential Quadratic Programming numerical method 'SQP' when globalized by the trust region method. By using a merit function, a sequential quadratic programming method associated with global trust regions bypasses the non-convex problem. This method used a nonlinear interior point trust region optimization solver under AMPL. Among several possible solutions, it is shown that there is an optimal trajectory leading to a reduction of noise levels on approach.

Introduction

The aim of this work is the development of a theoretical model of noise optimization while maintaining a reliable evolution of the flight procedures of two commercial aircraft on approach. These aircraft are supposed to land successively on one runway without conflict [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 flight 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 differential equations [START_REF] Chryssoverghi | Classical and relaxed optimization methods for optimal control problems International Mathematical Forum[END_REF]. The 3-D movement of the two planes is described by a system depending on ordinary non-linear differential 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 flight safety.

The problem of optimal control is described and solved by a Trust Region Sequential Quadratic Programming method 'TRSQP' [START_REF] Avdoshin | Optimisation des Systèmes dynamiques hybrides[END_REF][START_REF] Bergounioux | Optimisation et Controle des systèmes linéaires[END_REF][START_REF] Bonnans | Numerical Optimization: Theorical and pratical Aspects[END_REF][START_REF] Christof | SQP-methods for solving optimal control problems with control and state constraints: adjoint variables, sensitivity analysis and real-time control[END_REF][START_REF] Mathew | A feasible trust-region sequential quadratic programming algorithm[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] Byrd | KNITRO: An integrated Package for nonlinear optimization[END_REF] under A Mathematical Programming Modeling Language 'AMPL' [START_REF] Fourer | A modelling Language for Mathematical Programming[END_REF][START_REF] Laboratories | AMPL, A modelling Language for Mathematical Programming[END_REF]. The global convergence properties are analyzed under different assumptions on the approximate Hessian. Additional assumptions on the feasibility perturbation technique are used to prove quadratic convergence to points satisfying second-order sufficient conditions.

Details of the two-aircraft flight 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 conflict detection for air traffic management[END_REF]. The transition between these three frames is shown easily [START_REF] Boiffier | The Dynamics of Flight, The Equations, SUPA ÉRO[END_REF]. In general, the equations of motion of each aircraft are summarized as:

-→

F ext i -dm i dt -→ V a i = m i d -→ V a i dt -→ M ext G i = d dt [I G i -→ Ω i ] d -→ X o dt = d -→ X 1 dt + -→ Ω 10 × -→ X (1)
The index i = 1, 2 reflects the aircraft. In the system above, -→ F ext i represents the external forces acting on the aircraft, m i (t) the mass of the aircraft, v i the airspeed of aircraft,

-→ M ext G i the moments of each aircraft, J(G i , A i ) the inertia matrix, d -→ X o dt
is the derivative with respect to time of the vector

X in vehicle-carried normal Earth frame R O , d -→ X 1
dt is the derivative with respect to time of the vector X in frame R 1 , Ω i the angular rotation of the aircraft and Ω 10 is the angular velocity of the frame R 1 relative to the frame R O . After transformations and simplifications, the system (1) becomes:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ V a i = 1 m i [-m i gsinγ a i -1 2 ρSV 2 a i C D +(cosα a i cosβ a i + sinβ a i + sinα a i cosβ a i )F x i -dm i dt u i -m i ΔA i u ] β a i = 1 m i Va i [m i gcosγ a i sinμ a i + 1 2 ρSV 2 a i C y i +[-cosα a i sinβ a i + cosβ a i -sinα a i sinβ a i ]F y i -dm i dt v i -m i ΔA i v ] αa i = 1 m i Va i cosβa i [m i gcosγ a i cosμ a i -1 2 ρSV 2 a i C L i + [-sinα a i + cosα a i ]F z i -dm i 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 a i C l i + 2 j=1 F j [y b M i j cosβ m i j sinα m i j -z b M i j sinβ m i j ]} + E AC-E 2 {p i q i (A -B) -Er i q i + 1 2 ρSlV 2 a i C ni + 2 j=1 F j [x b M i j sinβ m i j -y b M i j cosβ m i j cosα m i j ]} qi = 1 B {-r i p i (A -C) -E(p 2 i -r 2 i ) + 1 2 ρSlV 2 a i C mi + 2 j=1 F j [z b M i j cosβ m i j cosα m i j -x b M i j cosβ m i j sinα m i j ]} ṙi = E AC-E 2 {r i q i (B -C) + Ep i q i + 1 2 ρSlV 2 a i C l i + 2 j=1 F j [y b M i j cosβ m i j sinα m i j -z b M i j sinβ m i j ]} + A AC-E 2 {p i q i (A -B) -Er i q i + 1 2 ρSlV 2 a i C ni + 2 j=1 F j [x b M i j sinβ m i j -y b M i j cosβ m i j cosα m i j ]} ẊG i = V a i cosγ a i cosχ a i + u w ẎG i = V a i cosγ a i sinχ a i + v w ŻG i = -V a i sinγ a i + 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 = -2.01 × 10 -5 (Φ-μ-K ηc ) √ Θ √ 5ηn(1+η tf λ) G+0.2M 2 η d η tf λ-(1-λ)M 1 2 ρSV 2 ai C D (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 coefficient, C yi = C yβ β + C yp pl V + C yr rl V + C Y δ l δ l +C Y δn δ n is the lateral forces coefficient, C Li = C Lα (α a -α a0 )+C Lδm δ m + C LM M + C Lq q b a l V is the lift coefficient, C li = C lβ β + C lp pl V + C lr rl V + C lδ l δ l + C lδn δ n is the rolling moment coefficient, C mi = C m0 + C mα (α -α 0 ) + C mδm δ m is the pitching moment coefficient, C ni = C nβ β +C np pl V +C nr rl V +C nδ l δ l +C nδn δ n is the yawing moment coefficient, (x b Mij , x b Mij , x b Mij ) 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 mass change is reflected in the aircraft fuel consumption as described by E. Torenbeek [START_REF] Roux | Pour une approche analytique de la dynamique du[END_REF] where the specific consumption is

C SR i = 2.01 × 10 -5 (Φ -μ -K i ηc ) √ Θ 5η n (1 + η tf i λ) G i + 0.2M 2 i η d i η tf i λ -(1 -λ)M i
with the generator function G:

G i = (Φ -K i ηc )(1 - 1.01 η ν-1 ν i (K i +μ i )(1- K i Φηcη t ) ) K i = μ i ( ν-1 ν c -1) μ i = 1 + ν-1 2 M 2 i
The Nomenclature of engine performance variables are given by G the gas generator power function, G0 the gas generator power function (static, sea level), K i the temperature function of compression process, M i the flight Mach number, T4 the turbine Entry total Temperature, T0 the ambient temperature at sea level, T the flight temperature, while the nomenclature of engines yields is η c = 0.85 the isentropic compressor efficiency,

η d i = 1 -1.3( 0.05 Re 1 5
) 2 ( 0.5 M i ) 2 L D , the isentropic fan intake duct efficiency, L the duct length, D the inlet diameter, Re the reynolds number at the entrance of the nozzle, η f i = 0.86 -

3.13 × 10 -2 M i the isentropic fan efficiency, η i = 1+η d i γ-1 2 M 2 i 1+ γ-1 2 M 2 i
the gas Generator intake stagnation pressure ratio, η n = 0.97 the isentropic efficiency of expansion process in nozzle, η t = 0.88 the isentropic turbine efficiency η tf i = η t η f i , c the overall pressure ratio (compressor), ν the ratio of specific heats ν = 1.4, λ the bypass ratio, μ i the ratio of stagnation to static temperature of ambient air, Φ the nondimensional turbine entry temperature Φ = T 4 T and Θ the relative ambient temperature Θ = T T 0 . The expressions

α ai (t), β ai (t), θ i (t), ψ i (t), φ i (t), V a i (t), X G i (t), Y G i (t), Z G i (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 a i (t), X G i (t), Y G i (t), Z G i (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) = (δ l i (t), δ m i (t) , δ n i (t), δ x i (t)) (5)
where the expressions δ l i (t), δ m i (t) , δ n i (t), δ x i (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 γ a i (t), χ a i (t), μ a i (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 a i = 1 m i [-m i gsinγ a i -1 2 ρSV 2 a i C D + (cosα a i + sinα a i )F x i -ṁi u i ] αa i = 1 m i Va i cosβa i [m i gcosγ a i cosμ a i -1 2 ρSV 2 a i C L i +[cosα a i -sinα a i ]F z i -ṁi w i ] ṗi = C AC-E 2 {r i q i (B -C) -Ep i q i + 1 2 ρSlV 2 a i C l i } + E AC-E 2 {p i q i (A -B) -Er i q i + 1 2 ρSlV 2 a i C ni qi = 1 B {-r i p i (A -C) -E(p 2 i -r 2 i ) + 1 2 ρS i lV 2 a i C mi }, ṙi = E AC-E 2 {r i q i (B -C) + Ep i q i + 1 2 ρSlV 2 a i C l i + A AC-E 2 {p i q i (A -B) -Er i q i + 1 2 ρSlV 2 a i C ni } ẊG i = V a i cosγ a i cosχ a i ẎG i = V a i cosγ a i sinχ a i ŻG i = -V a i sinγ a i φ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 = -2.01 × 10 -5 (Φ-μ-K ηc ) √ Θ √ 5ηn(1+η tf λ G+0.2M 2 η d η tf λ-(1-λ)M 1 2 ρSV 2 ai C D (7)
By the combination of this system with the aircraft control, one has the twoaircraft dynamic flight model as shown in [START_REF] Blin | Stochastic conflict detection for air traffic management[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 a i (t), X G i (t), Y G i (t), Z G i (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 objective function model. In this paper, the considered cost function is the Sound Exposure Level 'SEL' [START_REF] Abdallah | Minimisation des bruits des avions commerciaux sous contraintes physiques et aérodynamique[END_REF][START_REF] Harris | How do we Describe Aircraft Noise?[END_REF][START_REF] Martin | Noise monitoring in the vicnity of airports,DSNA-DTI[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 first and the second aircraft, the objective function is calculated as:

SEL 1 = 10log 1 to t 20 t 10 10 0.1L A1,dt (t) dt , t ∈ [t 10 , t 20 ] SEL 12 = SEL 11 ⊕ SEL 21 = 10 log[ 1 to t 1f t 20 10 0.1L A1,dt (t) dt + 1 to t 1f t 20 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 = (t 20 -t 10 )SEL 1 ⊕(t 1f -t 20) SEL 12 ⊕(t 2f -t 1f )SEL2 t 2f -t 10 = 10 log{ 1 t 2f -t 10 [(t 20 -t 10 ) t 20 t 10 10 0.1L A1 (t) dt +(t 1f -t 20 ) t 1f t 20 10 0.1L A1 (t) dt + (t 1f -t 20 ) t 1f t 20 10 0.1L A2 (t) dt +(t 2f -t 1f ) t 2f t 1f 10 0.1L A2 (t) dt, ]}, t ∈ [t 10 , t 2f ] ( 10
)
where SEL G is the cumulated two-aircraft noise and the operator ⊕ means the acoustic adding. Expressions L A1 (t), L A2 (t) are equivalent and reflect 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] Groesbeck | An improved prediction method for noise generated by conventional profile 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 effective 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 defined by: w = 3(V e /c) 3.5 0.6 + (V e /c)

ΔV = -15log(C D (M c , θ)) -10log(1 -Mcosθ), 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 considered constraints concern aircraft flight speeds and altitudes, flight angles and control positions, energy constraint, aircraft separation, flight velocities of aircraft relative to the earth and the aircraft mass:

1. The vertical separation given by

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

The horizontal separation

X G 12 = X G 1 -X G 2
3. The aircraft speed V a i must be bounded as follows 1.3V s ≤ V a i ≤ V if where V s is the stall speed,V if is the maximum speed and V io = 1.3V s 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 i0 , p if ], the pitch velocity of the aircraft relative to the earth q i ∈ [q i0 , q if ] and the yaw velocity of the aircraft relative to the earth

r i ∈ [r i0 , r if ] .
4. On the approach, the ICAO standards and aircraft manufacturers require flight angle evolution as follows: attack angle

α a i ∈ [α io , α if ], the inclination angle θ i ∈ [θ i0 , θ if
] and the roll angle

φ i ∈ [φ io , φ if ]. 5. The aircraft control δ(t) = (δ l i (t), δ m i (t) , δ n i (t), δ x i (t)
) keeps still between the position δ li0 and δ lif for the roll control, δ mi0 and δ mif for the pitch control, δ ni0 and δ nif for the yaw control and δ xi0 and δ xif for the thrust.

The mass m i of the aircraft

A i is variable: m i0 < m i < m if , i = 1, 2.
This constraint results in energy consumption of the aircraft [START_REF] Boiffier | 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:

k 1i (y i (t), u i (t)) ≤ 0 k 2i (y i (t), u i (t)) ≥ 0 ( 11
)
where

k(t) :R 12 × R 4 -→ R 16 , (y i (t), u i (t)) -→ k i (y i (t), u i (t)) k 1i (y i (t), u i (t)) = (α i (t) -α if , θ i (t) -θ if , ψ i (t) -ψ if , φ i (t) -φ if , V a i (t) -V aif , X G i (t) -X Gif , Y G i (t) -Y Gif , Z G i (t) -Z Gif , p i (t) -p if , q i (t) -q if , r i (t) -r if , δ l i (t) -δ l i f , δ m i (t) -δ mif , δ n i (t) -δ nif , δ x i (t) -δ xif , m i (t) -m if ) k 2i (y i (t), u i (t)) = (α i (t) -α i0 , θ i (t) -θ i0 , ψ i (t) -ψ i0 , φ i (t) -φ i0 , V a i (t) -V ai0 , X G i (t) -X Gi0 , Y G i (t) -Y Gi0 , Z G i (t) -Z Gi0 , p i (t) -p i0 , q i (t) -q i0 , r i (t) -r i0 , δ l i (t) -δ l i 0 , δ m i (t) -δ mi0 , δ n i (t) -δ ni0 , δ x i (t) -δ xi0 , m i (t) -m i0 ).
The following values reflect the digital applications considered for the twoaircraft [START_REF] Abdallah | Minimisation des bruits des avions commerciaux sous contraintes physiques et aérodynamique[END_REF][START_REF] Boiffier | The Dynamics of Flight, The Equations, SUPA ÉRO[END_REF][START_REF] Boiffier | Dynamique de vol de l'avion, SupAéro, Départements des Aéronefs[END_REF][START_REF] Roux | Pour une approche analytique de la dynamique du[END_REF].

V a1f = V a2f = 200 m/s V a10 = V a20 = 73.45 m/s The A1 Aircraft altitude Z G1f = 35 × 10 2 m Z G10 = 0 m The A2 Aircraft altitude Z G2f = 41 × 10 2 m Z G20 = 0 m The aircraft roll control δ l1f = δ l2f = 0.0174 δ l10 = δ l20 = -0.0174 The pitch control δ m1f = δ m2f = 0.087 δ m10 = δ m20 = 0 The yaw control δ n1f = δ n2f = 0.314 δ n10 = δ n20 = -0.035 The thrust control δ x1f = δ x2f = 0.6 δ x10 = δ x20 = 0.2 The attack angle α a1f = α a2f = 12 • α a10 = α a20 = 2 • The inclination angle θ a1f = θ a2f = 7 • θ a10 = θ a20 = -7 • The air-path inclination angle γ a1f = γ a2f = 0 • γ a10 = γ a20 = -5 • The aerodynamic bank angle μ a1f = μ a2f = 3 • μ a10 = μ a20 = -2 • The air-path azimuth angle χ a1f = χ a2f = 5 • χ a10 = χ a20 = -5 • The roll angle φ a1f = φ a2f = 1 • φ a10 = φ a20 = -1 • The cup ψ a1f = ψ a2f = 3 • ψ a10 = ψ a20 = -3
p 1f = p 2f = 1 • s -1 p 10 = p 20 = -1 • s -1
The Aircraft pitch velocity relative to the earth

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

r 1f = r 2f = 12 • s -1 r 10 = r 20 = -12 • s -1
Table 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] Bonnans | Numerical Optimization: Theorical and pratical Aspects[END_REF] and the the aircraft flight constraints [START_REF] Christof | SQP-methods for solving optimal control problems with control and state constraints: adjoint variables, sensitivity analysis and real-time control[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 k 1i (y i (t), u i (t)) ≤ 0 k 2i (y i (t), u i (t)) ≥ 0 ( 12
) where g 12 shows the aircraft coupling noise function and J G12 is the SEL of the two A300-aircraft.

The numerical processing

The problem as defined 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] Bonnans | Numerical Optimization: Theorical and pratical Aspects[END_REF][START_REF] Wright | Numerical Optimization, Second edition[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 differentiable. The Lagrangian of the system (13) is defined 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 finite differences or the calculation of the adjoint systems for problems with many parameters and finally by the sensitivity analysis. This last technique is very effective in the case of a large number of variables with few parameters [START_REF] Avdoshin | Optimisation des Systèmes dynamiques hybrides[END_REF][START_REF] Christof | SQP-methods for solving optimal control problems with control and state constraints: adjoint variables, sensitivity analysis and real-time control[END_REF]. The SQP method is a qualified local method. Its convergence is quadratic if the first iterate is close to a solution ỹ satisfying the sufficient optimality conditions [START_REF] Bricker | SQP: Sequential quadratic programming[END_REF][START_REF] Gilbert | SQPpro: A solver of nonlinear optimization problems, using an SQP approach[END_REF][START_REF] Gill | SNOPT: An SQP Algorithme for large-scale constrained optimizatio[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 trust-region 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 confience 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. Théorie et algorithmes[END_REF]. In the latter case, we have an ordinary differential 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 offer 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 sufficient to start the whole iterative process [START_REF] Gilbert | Eléments d'optimisation différentiable -Théorie et Algorithmes[END_REF][START_REF] Khardi | Mathematical Model for Advanced CDA and Takeoff Procedures Minimizing Aircraft Environmental Impact[END_REF][START_REF] Nahayo | Les méthodes numériques appliquées en optimisation non-linéaire et en commande optimale[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), 2. The asymptotic exactness relation

||Δx -Δx k || 2 ≤ φ(||Δx k || 2 )||Δx k || 2 (18)
is satisfied where φ : R + -→ R + with φ(0) = 0.

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

• The trust region restriction ( 15) is added to the SQP problem ( 14) 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.

• 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 ∈ Γ * verifies 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 ∈ Γ * satisfies the qualification 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 ( 12) which satisfies a constraints qualification 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 first equation reflects the optimality, the second and third the feasibility conditions. The others reflect the additional conditions and Lagrange multipliers corresponding to inactive constraints nj(x) are zero. The couple (x, λ * ) such that the KKT conditions are satisfied is called primal-dual solution of [START_REF] Gilbert | Eléments d'optimisation différentiable -Théorie et Algorithmes[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[END_REF], taking x a local solution of ( 19) and satisfying a qualification condition, then there exist multipliers (λ * ) such that the KKT conditions are verified . So we have

∇ 2 xx L(x, λ * )d.d > 0∀h ∈ C * where C * is a critical cone defined 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 sufficient optimality conditions of second order [START_REF] Bergounioux | Optimisation et Controle des systèmes linéaires[END_REF], suppose that there exists (λ * ) which satisfy the KKT conditions and such that [START_REF] Chryssoverghi | Classical and relaxed optimization methods for optimal control problems International Mathematical Forum[END_REF].

∇ 2 xx L(x, λ * )d.d > 0∀h ∈ C * \{0}. So x is a local minimum of

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 affectations 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 = J G12 (x k )-J G12 (x k + Δx k ) -K G12 ( Δ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 definite 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 definite Hessian is maintained providing q T k s k is positive at each update and that H is initialized with a positive definite matrix. This algorithm is implemented by AMPL language programming and the KNITRO solver [START_REF] Byrd | KNITRO: An integrated Package for nonlinear optimization[END_REF].

Analysis of the algorithm and its convergence. Let us define 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, ∞]
specifies 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 differentiable 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 defined.

In this paragraph, one wants to prove that the algorithm has a convergence to stationary point of [START_REF] Dgac | Mémento à l'usage des utilisateurs des procédures d'approche et de départ aux instruments[END_REF]. If we consider that all assumptions hold for each feasible point x for [START_REF] Chryssoverghi | Classical and relaxed optimization methods for optimal control problems International Mathematical Forum[END_REF], the Mangasarian-Fromowitz are satisfied for constraints. After all, the KKT optimality conditions are specified 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.

Numerical Results

The observation points are taken on the ground under the flight path and are independent of each other. Numerical processing is implemented by AMPL and KNITRO solver. KNITRO output optimality conditions for the obtained solution is achieved as follow: Multistart stopping, found local optimal solution.

MULTISTART: Best locally optimal point is returned. EXIT: Locally optimal solution found. for ONL 10 . In this figure, the legend ONL means optimal noise level. As specified, noise level increases (till 550 sec) and is maximum when the observation point lies below the aircraft. Noise levels decrease gradually as the aircraft moves away from the observation point. This is confirmed by Khardi analysis [START_REF] Khardi | Reduction of commercial aircraft noise emission around airports. A new environmental challenge[END_REF]. By comparison, this result is also close to standard values of jet noise on approach as shown by Harvey [START_REF] Abdallah | Optimization of operational aircraft parameters reducing noise emission[END_REF][START_REF] Hubbard | Aeroacoustics of flight vehicles, Theory and Practices[END_REF]. To conclude, numerical calculations carried out in this paper are efficient and fitted with experimental and theoretical researches related to acoustical developments. Figure 2 shows the trajectories which reflect a path in one level flight followed by a continuous descent till the aircraft touch point. The aircrafts' landing procedures are sufficiently separated. It is obvious that each aircraft follows its optimal trajectory when considering the separation distance. Constraints on speeds described in the previous table are considered, allowing a subsequent landing on the same track. Thus, as recommended by ICAO, the security conditions are met and flight procedures are good as shown by the presented results. The maximum altitudes considered are 3500 m and 4100 m for the first and the second aircraft. The duration approach is 600 s for the first aircraft and 690 s for the second. This figure shows that after some time, we have obtained the same optimal trajectory for the two-aircraft even the procedures are different. This shows the aircraft trajectory resulting from the two trajectories combination. This figure also shows aircraft speed evolution during landing. For the first, the aircraft speed decreases from 200 m/s to 140 m/s and keeps a constant position till the end of the aircraft landing. This evolution remains the same for the speed of the second aircraft. Processing calculation provided that the aircraft throttle position is kept constant (0.6) during the landing procedures. The two-aircraft roll velocity p 1 , p 2 , pitch velocity q 1 , q 2 and yaw velocity r 1 , r 2 , both related to earth frame, are obtained and they have a constant behavior.

Conclusion

We have developed a mathematical model in the case of two approaching aircraft landing in succession on the same track. An algorithm for solving the optimal control model has been developed. Theoretical considerations and practices of the feasible TRSQP algorithm are used by KNITRO 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 effective. Some added conditions are necessary to prove the global convergence of the considered algorithm. We found that the aircraft optimal trajectories coincide for a large portion of the flight as soon as the continuous descent is initiated. Further researches are needed to complete the problem processing.

ACKNOWLEDGEMENTS. This work is supported by the Agence Universitaire de la Francophonie-Région Afrique Centrale. This is also supported by the International Join Research and International Publication Project (Universitas Indonesia -Indonesia and IFSTTAR -France), Memorandum of Understanding between Universitas Indonesia and IFSTTAR, the French Embassy in Indonesia, and DIKTI (Indonesia).

  Final Statistics: Final objective value = 4.91134561926630e+01 Final feasibility error (abs / rel) = 5.79e-09 / 1.96e-10 Final optimality error (abs / rel) = 2.58e-07 / 2.58e-07 Number of iterations = 72 Number of CG iterations = 3 Number of function evaluations = 85 Number of gradient evaluations = 73 Number of Hessian evaluations = 72 Total program time (secs) = 258.19754 ( 258.104 CPU time) Time spent in evaluations (secs) = 240.00644 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 touch point on the ground is (0 m, 0 m, 0 m) while the temporal separation of aircraft is 90 s. At each point, it is a vector of N noise levels as shown in the discretization process. It is very important to consider the maximum value among the N values, which value corresponds to the shortest distance between the noise source and the observation point.

Figure 1 :

 1 Figure 1: Aircraft noise at the indicated reception point

Figure 2 :

 2 Figure 2: Aircraft optimal flight paths and speeds

Figure 3 :

 3 Figure 3: Flight-path angles of the aircraft

Figure 4 :

 4 Figure 4: Aircraft finesse

  The mass of the A1 Aircraft m 10 1.1 × 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 [8] 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 The Aircraft vertical separation Z 12 = 2 × 10 3 ft 6 × 10 2 m The Aircraft longitudinal separation X G 12 = 5 NM 9 × 10 3 m The Aircraft roll velocity relative to the earth

•

The limits of time t 1f = 600 s,t 2f = 645 s t 10 = 0 s, t 20 = 45 s


S Khardi
email: salah.khardi@ifsttar.fr

Aircraft Flight Path Optimization The Hamilton-Jacobi-Bellman considerations

Keywords: Hamilton-Jacobi-Bellman equation, Dynamic programming, Aircraft, Flight path optimization, environment

Flight path optimization is designed for minimizing aircraft noise, fuel consumption and air pollution around airports. This paper gives theoretical considerations and algorithms solving the Hamilton-Jacobi-Bellman equation (HJB) of aircraft trajectory optimization. Comparisons with direct and indirect methods are carried out. The OCP problem is transformed into new equalities-constrained as a viscosity problem. This constitutes an original dynamical system extension where subsystems are linked to the original dynamics via algebraic coupling equations. A feedback control method using dynamic programming has been developed. Comparisons show its fast computing times. It provides the best optimized ight paths which could be more suitable for CDA approach applicability. A two-segment approach is provided by HJB method which also favors fuel consumption saving. This improved CDA approach could benet both airlines and communities. Because of the processing speed and eciency of the HJB method, it can be better interfaced with the in-ight management system respecting airspace system regulation constraints.

Introduction

Due to the increase of air trac, populations living near airports and the environment are impacted by commercial aircraft. This is considered to be one of the most environmental concerns aecting people and the physical environment [START_REF]Dose-Eect European Commission. Position paper on dose response relationships between transportation noise and annoyance[END_REF][START_REF] Lambert | Perception and attitudes to transportation noise in France: A national survey[END_REF]. Technology development, airspace management, operational improvement and system eciency should be considered as an environmental innovation. There is no justication that air transport will not continue to progress without improving its environmental impacts [START_REF] Kurniawan | Comparison of methodologies estimating emissions of aircraft pollutants. Environmental impact assessment around airports[END_REF]. This is because all types of procedures are not optimized but rather generic in nature, that new ight path development, associated to new aircraft design and engines, is a solution which should contribute to a decrease in aircraft annoyances. This development cannot be carried out without improvement of the scientic knowledge in this eld, in particular the contribution of modeling. The latter consists in developing ecient data processing tools allowing in-ight diagnosis and control in real-time taking into account the FMS (ight management system) and the AMS (airspace management system) updates. In this paper, we have suggested a dynamic optimization method solving a model governed by an ODE system [START_REF] Bellman | Dynamic Programming[END_REF][START_REF] Pontryagin | Mathematical Theory of Optimal Processes[END_REF][START_REF] Bryson | Applied Optimal Control[END_REF]. The cost function of this model describes aircraft noise and fuel consumption [START_REF] Berman | Some analytical consideration in jet noise prediction. mechanics of Sound Generation in Flows[END_REF][START_REF] Stone | An improved prediction method for noise generated by conventional prol coaxial jets[END_REF][START_REF] Tester | Engine noise source breakdown. Theory, simulation and results[END_REF]10,[START_REF] Moris | Acoustic analogy and alternative theories for jet noise prediction[END_REF]. The ODE depends on the ight dynamics of the aircraft and considers ight safety and stability requirements (constraints and extreme conditions). Numerical methods solving control problems fall in many categories [START_REF] Osher | High-order Essentially Nonoscillatory Schemes for Hamilton Jacobi Equations[END_REF][START_REF] Krener | Optimal model matching controllers for linear and nonlinear systems[END_REF][START_REF] Wright | Interior-point methods for optimal control of discrete-time systems[END_REF][START_REF] Wright | Primal-dual interior-point methods[END_REF][START_REF] Bonnans | Optimisation Numérique. Aspects théoriques. Mathématiques et Applications[END_REF][START_REF] Bonnans | Large-scale Direct Optimal Control Applied to a Re-Entry Problem[END_REF][START_REF] Cormen | Introduction to Algorithms[END_REF][START_REF] Giegerich | A Discipline of Dynamic Programming over Sequence Data[END_REF][START_REF] Weiser | Interior-point methods in function space[END_REF][START_REF] Berend | An Interior-Point Approach to trajectory Optimization[END_REF] among which it is necessary to choose, improve or develop a new method. In this context, this paper gives theoretical considerations and algorithms solving Hamilton-Jacobi-Bellman equation "HJB" [START_REF] Barnett | Introduction to Mathematical Control Theory[END_REF][START_REF] Allaire | Numerical Analysis and Optimization[END_REF] because the problem is transformed into new equalities-constrained as a viscosity problem. This is an original system extension where subsystems are linked to the dynamics via algebraic coupling equations [START_REF] Bryson | Dynamic Optimization[END_REF][START_REF] Graichen | Transformation of optimal control problems with a state constraint avoiding interior boundary conditions[END_REF]. Among the existing methods, solving the HJB problem, a feedback control method using dynamic programming has been developed. The latter is a method used for solving complex problems by breaking them down into simpler subproblems. To solve the given problem, it solves subproblems, then combines solutions of them to reach a global solution. Because subproblems are generally the same, it seeks to solve each subproblem only once reducing the number of the total computations, in particular when the subproblem number is exponentially large. Comparisons have been performed between HJB, direct and indirect methods [START_REF] Stryk | Direct and indirect methods for trajectory optimization[END_REF][START_REF] Stryk | Numerical solution of optimal control problems by direct collocation[END_REF][START_REF] Chaplais | Inversion in indirect optimal control of multivariable systems[END_REF][START_REF] Graichen | Incorporating a class of constraints into the dynamics of optimal control problems[END_REF] stressing the computing times with the aim of nding the best aircraft approach also favoring fuel consumption saving. Indeed, the direct approach has been used reducing the OCP to a nite-dimensional nonlinear program which is solved by a standard nonlinear programming solver. Algorithms are adapted and modied versus constraints, limits of ight dynamic parameters, and location points on the ground. For the indirect approach, optimality conditions given by Pontryagin's principle, have been discretized. An AMPL model (A Modeling Language for Mathematical Programming) [START_REF] Ampl | A Modeling Language for Mathematical Programming[END_REF] combined with NLP solver [START_REF] Gill | A large-scale smooth optimization problems having linear or nonlinear objectives and constraints[END_REF][START_REF] Waltz | KINITRO User's Manual Version 4[END_REF][START_REF] Byrd | KNITRO An Integrated Package for Nonlinear Optimization[END_REF] has been performed for processing. In-depth details have been described in previous papers [START_REF] Abdallah | Optimization of Operational aircraft parameters reducing noise emission[END_REF][START_REF] Khardi | Optimal approach minimizing aircraft noise and fuel consumption[END_REF][START_REF] Houacine | A Gauss Pseudospectral method for Less Noise and Fuel Consumption of Aircraft Operations[END_REF][START_REF] Abdallah | Optimization of operational aircraft parameters reducing noise emission[END_REF][START_REF] Khardi | The Trust Region Esquential Quadratic Programming method applied to two-aircraft acoustic optimal control problem[END_REF]. Technically, we analyze the processing speed and algorithm eciency and their ability to be interfaced with the in-ight management system respecting airspace system regulation constraints. This has to be tted with the necessity to compensate both the growth in air trac and the encroachment of airport-neighboring communities. This paper presents in the rst two sections an introduction and the optimal control problem of aircraft trajectory minimization (ight dynamics, constraints, and aircraft noise model). For comparison, the third section gives applied methods (indirect, direct, Hamilton-Jacobi-Bellman). The last two sections show numerical results followed by the conclusion.

Optimal Control Problem

We present in this section a summary of the optimal control problem that will be solved and methods compared [START_REF] Abdallah | Optimization of Operational aircraft parameters reducing noise emission[END_REF][START_REF] Khardi | Optimal approach minimizing aircraft noise and fuel consumption[END_REF][START_REF] Abdallah | Optimization of operational aircraft parameters reducing noise emission[END_REF][START_REF] Khardi | The Trust Region Esquential Quadratic Programming method applied to two-aircraft acoustic optimal control problem[END_REF]. The system of dierential equations commonly employed in aircraft trajectory analysis is the following six-dimension system derived at the center of mass of the aircraft [START_REF] Abdallah | Optimization of Operational aircraft parameters reducing noise emission[END_REF][START_REF] Khardi | Optimal approach minimizing aircraft noise and fuel consumption[END_REF][START_REF] Abdallah | Optimization of operational aircraft parameters reducing noise emission[END_REF][START_REF] Mccormick | Aerodynamics, Aeronautics and Flight Mechanics[END_REF][START_REF] Etkin | Dynamics of Flight[END_REF][START_REF] Boier | The Dynamics of Flight[END_REF][START_REF] Abdallah | Minimisation des bruits des avions commerciaux sous contraintes physiques et aérodynamiques[END_REF]:

(ED)                                              V = g ( T cos α -D mg -sin γ ) γ = 1 mV ((T sin α + L) cos µ -mg cos γ) χ = (T sin α + L) sin µ mV cos γ ẋ = V cos γ cos χ ẏ = V cos γ sin χ ḣ = V sin γ
where V, γ, χ, α and µ are respectively the speed, the angle of descent, the yaw angle, the angle of attack and the roll angle. (x, y, h) is the position of the aircraft. The variables T, D, L, m and g are respectively the engine thrust, the drag force, the lift force, the aircraft mass and the aircraft weight acceleration given in previous paper [START_REF] Abdallah | Optimization of Operational aircraft parameters reducing noise emission[END_REF][START_REF] Khardi | Optimal approach minimizing aircraft noise and fuel consumption[END_REF][START_REF] Khardi | The Trust Region Esquential Quadratic Programming method applied to two-aircraft acoustic optimal control problem[END_REF][START_REF] Mccormick | Aerodynamics, Aeronautics and Flight Mechanics[END_REF][START_REF] Etkin | Dynamics of Flight[END_REF][START_REF] Boier | The Dynamics of Flight[END_REF][START_REF] Abdallah | Minimisation des bruits des avions commerciaux sous contraintes physiques et aérodynamiques[END_REF][START_REF] Mattingly | Elements of Gas Turbine Propulsion[END_REF]. ED can be written in the following matrix form:

ż(t) = f (z(t), u(t))
where

z : [t 0 , t f ] -→ IR 6 t -→ z(t) = [V (t), γ(t), χ(t), x(t), y(t), h(t)] are the state variables u : [t 0 , t f ] -→ IR 3 t -→ u(t) = [α(t), δ x (t), µ(t)]
are the control variables and t 0 and t f are the initial and nal times.

Along the trajectory, we have some safety requirements and comfort constraints. For that, we have to respect parameter limits related to the safety of ight and the operational modes of the aircraft.

Aircraft modeling continues to meet the increased demands associated with aviation and airport expansion. Aircraft noise footprints are commonly used for forecasting the impact of new developments, quantifying the noise trends around airports and evaluating new tools. Thus, aircraft models have become more sophisticated and their validation complex. A number of them are entirely based on empirical data [START_REF] Ollerhead | The CAA Aircraft Noise Countour Model[END_REF]. Because of this complexity [START_REF] Thodes | Aircraft Noise Model validation[END_REF], such models are not characterized by a given analytical form describing noise at reception points on the ground. This paper uses the basic principles of aircraft noise modeling. The cost function may be chosen as any of the usual aircraft noise indices, which describes the eective noise level of the aircraft noise event [START_REF] Pietrzko | Prediction of A-Weighted aircrfat noise based on measured directivity patterns[END_REF][START_REF] Zaporozhets | Predicted Flight Procedures for Minimum Noise Impact[END_REF][START_REF] Zaporozhets | Aircraft Noise Modelling for Environmental Assessment Around Airports[END_REF][START_REF] Clemente | Sensitivity of the FAA Integrated Noise Model to input parameters[END_REF]. This study is limited to minimize the index L eq,∆T using a semi-empirical model of jet noise [START_REF] Berman | Some analytical consideration in jet noise prediction. mechanics of Sound Generation in Flows[END_REF][START_REF] Stone | An improved prediction method for noise generated by conventional prol coaxial jets[END_REF]10,[START_REF] Moris | Acoustic analogy and alternative theories for jet noise prediction[END_REF][START_REF] Norum | Simulated High-Speed Flight Eects on Supersonic Jet Noise[END_REF][START_REF] Tam | Computational Aeroacoustics. Issues and methods[END_REF][START_REF] Tam | Noise from High Speed Jets[END_REF][START_REF] Zaporozhets | Optimisation of aircraft ight trajectories as a basis for noise abatement procedures around the airports[END_REF][START_REF] Filippone | Comprehensive analysis of transport aircrfat ight performance[END_REF]. The cost function is expressed in the form:

J : C 1 ([t 0 , t f ], IR 6 ) × C 1 ([t 0 , t f ], IR 4 ) -→ IR J(X(t), U (t)) = ∫ t f t 0 (ℓ(X(t), U (t)) + Φ(X(t), t f -t 0 ))dt
where J is the criterion which optimizes noise levels and fuel consumption.

The cost function can be written in the following integral function form:

J : C 1 ([t 0 , t f ], IR 6 ) × C 1 ([t 0 , t f ], IR 3 ) -→ IR J(z(t), u(t)) = ∫ t f t 0 (ℓ(z(t), u(t)) + ϕ(z(t), (t f -t 0 )))dt
where J is the criterion to be optimized. Finding an optimal trajectory can be stated as an optimal control problem as follows (t 0 = 0):

(OCP )                        min J(z, u) = ∫ t f 0 (ℓ(z(t), u(t)) + ϕ(z(t), t f ))dt ż(t) = f (z(t), u(t)), ∀t ∈ [0, t f ] z I 1 (0) = c 1 , z I 2 (t f ) = c 2 a ≤ C(z(t), u(t)) ≤ b
where J : IR n+m → IR, f : IR n+m → IR n and C : IR n+m → IR q correspond respectively to the cost function, the dynamic of the problem possessing a unique state trajectory and the constraints. The second equation giving the trajectory is a nonlinear system with states in IR n . At t f , terminal conditions are imposed. t f can be xed or unspecied. In this paper terminal conditions, describing the boundary conditions, are specied and their values are given. t f is also xed. In the next section, we present dierent methods solving the presented (OCP ) (optimal control problem) problem. There is no practical theoretical limitations to using those methods that cannot be guaranteed to provide a global solution. We assume that the OCP has an optimal trajectory solution with the optimal cost. Subsequently, to reduce notations of l and ϕ, together, are called (replaced by) l.

Applied resolution methods

In this paper, we have applied three dierent approaches solving the OCP problems: direct, indirect, and dynamic approaches [START_REF] Osher | High-order Essentially Nonoscillatory Schemes for Hamilton Jacobi Equations[END_REF][START_REF] Krener | Optimal model matching controllers for linear and nonlinear systems[END_REF][START_REF] Cormen | Introduction to Algorithms[END_REF][START_REF] Dreyfus | The art and theory of dynamic programming[END_REF][START_REF] Iofe | Theory of Extremal Problems[END_REF][START_REF] Kamien | Dynamic Optimization[END_REF][START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF]. The direct method discretized the OCP for obtaining a nite-dimensional parameter optimization problem and solving the resulting nonlinear programming problem [START_REF] Hargraves | Direct trajectory optimization using nonlinear programming and collocation[END_REF][START_REF] Seywald | Trajectory optimization based on dierential inclusion[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Graichen | Handling constraintes in optimal control with saturation functions and system extension[END_REF]. It is well appropriated because of the domain of convergence and the ecient handling of constraints and the dened limits. It is opposed to the indirect approach based on Pontryagin's principle [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF][START_REF] Kirk | Optimal Control Theory, An Introduction[END_REF][START_REF] Gramkelidze | Principles of optimal control theory[END_REF][START_REF] Ojasiewicz | Necessary conditions for a nonlinear control system[END_REF][START_REF] Geering | Optimal Control with Engineering Applications[END_REF][START_REF] Ross | A Primer on Pontryagin's Principle in Optimal Control[END_REF] [START_REF] Dreyfus | The art and theory of dynamic programming[END_REF][START_REF] Iofe | Theory of Extremal Problems[END_REF][START_REF] Kamien | Dynamic Optimization[END_REF][START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF] than can be used to nd the optimal state, costate and the control variables which is focusing on the optimal function value [START_REF] Seywald | Method for automatic costate calculation[END_REF]. Dynamic programming is a method for solving complex problems by breaking them down into simpler subproblems. To solve a given problem, it solves dierent parts of the problem called subproblems, then combines solutions of the subproblems to reach a global solution. Because of subproblems are generally the same, it seeks to solve each subproblem only once reducing the number of the total computations, in particular when subproblem number is exponentially large.

The rst-order partial dierential equation is derived using the Hamilton-Jacobi-Bellman equation which uses the principle of Optimality of bellman.

The optimal value of the control vectors depending on the the date, the state and the parameters of the control problems is obtained. This way is called the feedback or the closed-loop form of the control. The Open-loop, where the form of the optimal control vector ca be obtained by the necessary conditions of optimal control theory can be used. The optimal value of the control vector can be obtained as a function of the independent variable time, the dierent used parameters and initial/nal conditions of the cost function and the state vector. The optimal solution can be given for any period and for any possible state.

Indirect method

We set H : IR n × IR m × IR n × IR q × IR q -→ IR the hamiltonian function of the problem (OCP ):

H(z, u, p, λ, µ) = ℓ(z, u) + ϕ(z, t f ) + p t f (z, u) + λ t (C(z, u) -a) + µ t (b -C(z, u))
where λ, µ are the multiplicators associated to the constraints and p is the costate vector.

We describe now the optimality conditions (OC) for the (OCP ) problem:

(OC)            ż(t) = f (z(t), u(t)) ṗ(t) = -H z (z(t), u(t), p(t), λ(t), µ(t)) u(t) = Argmin w H(z(t), w, p(t), λ(t), µ(t)) 0 = λ.(C(z(t), u(t) -a), λ ≥ 0 0 = µ.(b -C(z(t), u(t)), µ ≤ 0
In this paper, we have used the interior point method [START_REF] Wright | Interior-point methods for optimal control of discrete-time systems[END_REF][START_REF] Weiser | Interior-point methods in function space[END_REF][START_REF] Berend | An Interior-Point Approach to trajectory Optimization[END_REF] discretizing the optimality conditions of the system. The method which solved the (OC) problem is described below. We explain the transformation of the (OC) problem into a sequence of problems. We also show that the solution of the optimality conditions is a solution of the (OCP ) problem: discretization used an Euler scheme and the resolution the Newton method [START_REF] Bally | The law of the Euler scheme for stochastic dierential equations, II. Convergence rate of the density[END_REF].

By perturbing the last two equations (the complementary conditions) by a positive parameter ε we obtain the following system:

(OC ε )            ż(t) = f (z(t), u(t)) ṗ(t) = -H z (z(t), u(t), p(t), λ(t), µ(t)) u(t) = Argmin w H(z(t), w, p(t), λ(t), µ(t)) 1ε = λ.(C(z(t), u(t)) -a), λ ≥ 0 -1ε = µ.(b -C(z(t), u(t))), µ ≤ 0
The previous system can be interpreted as the optimality conditions for the following problem:

(P ε )        min ∫ t f 0 (ℓ ε (z(t), u(t)) + ϕ ε (z(t), t f ))dt ż(t) = f (z(t), u(t)), t ∈ [0, t f ]
where ℓ ε is the barrier logarithmic of (P ε ), dened by:

ℓ ε (z, u) + ϕ ε (z, t f ) = ℓ(z, u) + ϕ(z, t f ) -ε ∑ i [log(C i (z, u) -a i ) + log(b i -C i (z, u))] -εD(z)
To solve (OC), we have to solve a sequence of problems (OC ε ) by tending ε to zero. When ε decreases to 0, the solution of optimal conditions (OC ε ) is a solution of (OC). To compute the solution of the continuous optimal conditions, we rst discretized them. We obtained a set of non-linear equations, which has to be solved for the discretized control, state and costate vectors using a Newton method [START_REF] Berend | An Interior-Point Approach to trajectory Optimization[END_REF][START_REF] Dennis | Least change secant updates for quasi Newton methods[END_REF][START_REF] Hildebrand | Introduction to numerical analysis[END_REF]. For the discretization, we have chosen an Euler schema [START_REF] Hargraves | Direct trajectory optimization using nonlinear programming and collocation[END_REF][START_REF] Bally | The law of the Euler scheme for stochastic dierential equations, II. Convergence rate of the density[END_REF][START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic dierential equations[END_REF][START_REF] Cattiaux | Calcul stochastique et opérateurs dégénérés du second ordre, II. Problème de Dirichlet[END_REF] providing for (OC ε ) the following system:

           z k+1 = z k + hf (u k , z k ), k = 0, . . . , N -1 p k+1 = p k -hH z (z k , u k , p k , λ k , µ k ), k = 0, . . . , N -1 0 = H u (u k , z k , p k , λ k , µ k ), k = 0, . . . , N 1ε = λ k .(C(z k , u k ) -a), λ k ≥ 0, k = 0, . . . , N -1ε = µ k .(b -C(z k , u k )), µ k ≤ 0, k = 0, . . . , N
We have obtained a set of equations to be solved under the boundary constraints corresponding to the multiplicators:

(N ε )    F ε (X) = 0 λ k ≥ 0 µ k ≤ 0
where F ε is the set of optimal conditions, and X = (z k , u k , p k , λ k , µ k ) the variable vector. OCP is then successfully solved for decreasing ε with a nongrowth of the cost.

Direct method

To solve (OC) problems, many methods exist in the open literature [START_REF] Mastov | On a new principle of supersposition for optimization problems[END_REF][START_REF] Teo | A unied computational approach to optimal control problems[END_REF][START_REF] Zhou | User's Guide for FFSQP Version3.7. A Fortran code for solving constrained optimization problems, generating iterates satisfying all inequality and linear constraints[END_REF][START_REF] Huang | Solving HJBE by a modied method of characteristics[END_REF]. In this section, we have used a direct optimal control technique. We discretize the control and the state for reducing the dimension of the optimal control problem. Then, we solve the resulting nonlinear programming problem using a standard NLP solver. The paragraph below gives discretization steps which used Euler scheme where the continuous set of the obtained equations is replaced by a discretized control problem which is solved thereafter.

To solve (OC), we have used in this section a direct optimal control technique.

We discretize the control and the state for reducing the dimension of the optimal control problem. Then, we solve the resulting nonlinear programming problem using a standard NLP solver. We use an equidistant discretization of the time interval as:

t k = t 0 + kh, k = 0, ..., N and h = t f -t 0 N
Then we consider that u(.) is parameterized as a piecewise constant function:

u(t) := u k for t ∈ [t k-1 , t k [
and we use an Euler scheme to discretize the dynamic:

z k+1 = z k + hf (z k , u k ), k = 0, . . . , N -1
The new cost function can be written as:

N ∑ k=0 (ℓ(z k , u k ) + ϕ(z k ))
The continuous set of equations is replaced by the following discretized control problem:

(N LP )              min (z k ,u k ) N ∑ k=0 (ℓ(z k , u k ) + ϕ(z k )) z k+1 = z k + hf (z k , u k ), k = 0, . . . , N -1 z 0 I 1 = c 1 , z N I 2 = c 2 a ≤ C(z k , u k ) ≤ b, k = 0, . . . , N

Hamilton-Jacobi-Bellman method

The main idea behind this section is how to reduce an innite-period optimization problem to a two-period or some-period optimization problem. Diculties appear when the optimization problem is continuous. Two approaches to dynamic optimization: the Pontryagin approach which is Hamiltonian and the Bellman approach [START_REF] Dreyfus | The art and theory of dynamic programming[END_REF][START_REF] Iofe | Theory of Extremal Problems[END_REF][START_REF] Kamien | Dynamic Optimization[END_REF][START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF]. Dynamic programming method solves a complex problem by dividing it into simpler subproblems solving dierent parts of the problem reducing the processing steps in particular for large dimensions. The global solution is reached by combining solutions of the subproblems. When used, the method is faster compared to other existing methods. Top-down or bottom-up dynamic programming exist. The rst record the stages of calculation which are used thereafter, and the second reformulate the problem by a recursive series which calculations are where the processing are made easy. On the one hand, a candidate trajectory called the candidate solution is an element of a set of possible solutions for the given problem. Usually, a candidate solution could not be the best solution of the considered problem. It is the solution satisfying constraints. On the other hand the solution must belong the feasible region (solution space). The existence of the the trajectory solution is conrmed by Weierstrauss theorem which states that in a compact set the maximum and minimum values are reached for continuous or semi-continuous functions. Necessary conditions for optimality have been conrmed. Indeed, optima of the suggested inequality-constrained problem is instead found by the Lagrange multiplier method. A system of inequalities or the Karush-Kuhn-Tucker conditions, calculating the optimum, are met.

Sucient conditions for optimality are met. First derivative tests identify the optima without dierences between the minimum or the maximum. We assume that the cost function is twice dierentiable. By controlling the second derivative or the Hessian matrix in the transformed unconstrained problem or the matrix of second derivatives of the cost function and the constraints / the bordered Hessian, we can easily distinguish maxima and minima from stationary points. It should be remembered that methods commonly used to evaluate Hessians (or approximate Hessians) are Newton's method or Sequential quadratic programming. The second is particulary used for small-medium scale constrained problems. Some versions exist and can handle large-dimensional problems. A bordered Hessian is performed for the second derivative test in ight path constrained optimization problem. This could be in relationship with the suggested further research described in the last section making possible the generalization method. If we consider m constraints, 0 of the above matrix is a m*m block of zeroes, m border rows at the top and m border columns at the left; positive denite and negative denite can not apply here since a bordered Hessian can not be denite, we have z'Hz=0 if the vector z has non-zero as its rst element followed by zeroes. The second derivative test consists of restrictions of the determinants of a set of n-m sub-matrices of the bordered Hessian. A constraint reduction of the problem by one with n-m free variables.

If we consider:

f : R n -→ R m B(x 1 , x 2 , . . . , x n ) = (f 1 , f 2 , f 3 , . . . , f n )
the array of the second partial derivatives is not a nxn matrix but a tensor of order 3 or a multi-dimensional array mxnxn which can be reduced to a usual Hessian matrix with m=1. A Riemannian manifold and its Levi-Civita connection could be used for considerations using a Hessian tensor.

Dynamic methods approximate continuous systems to discrete systems leading to recurrence relations making easier the processing. One approach that should be considered as alternative transforms the OCP system to a non-linear partial dierential equations, often called the Hamilton-Jacobi-Bellman equation [START_REF] Bonnans | Using logarithmic penalties in the shooting algorithm for optimal control problems, Optimal Control Applications and Methods[END_REF][START_REF] Graichen | Feedforward control design for nite-time transition problems of nonlinear systems with input and output constraints[END_REF][START_REF] Graichen | Feedforward control design for nite-time transition problems of nonlinear systems with input and output constraints[END_REF].

The following state equation can be written as:

ż(t) = f (z(t), u(t)), ∀t ∈ [0, t f ]
It minimizes:

J(z, u) = ∫ t f 0 (ℓ(z(t), u(t)) + ϕ(z(t), t f ))dt
J(z,u) can be transformed as:

J = κ(z(t f ), t f ) + ∫ t f 0 ℓ(z(τ ), u(τ ))dτ
where B and κ are given functions, and t f is xed. When we use imbedding principle, we can rst give the following weak modication of the suggested problem at the moment t (less than t f ) at any acceptable state z(t) to make more large:

J(z(t), u(τ ) t≤τ ≤t f ) = κ(z(t f ), t f ) + ∫ t f t ℓ(z(τ ), u(τ ))dτ
Optimization process depends on numerical values of the state at t moments and the optimal control history in the considered time intervals. For a performed acceptable state for all t ≤ t f , we can assess controls minimizing the cost function J. This minimum can be written as:

J * (z(t), t) = min [u(τ ),t≤τ ≤t f ] {κ(z(t f ), t f ) + ∫ t f t ℓ(z(τ ), u(τ ))dτ }
Dividing the main interval in innitesimal intervals, we can write:

J * (z(t), t) = min [u(τ ),t≤τ ≤t f ] {κ(z(t f ), t f ) + ∫ t+∆t t ℓ(z(τ ), u(τ ))dτ + ∫ t f t+∆t ℓ(z(τ ), u(τ ))dτ }
Application of the optimality principle gives:

J * (z(t), t) = min [u(τ ),t≤τ ≤t+∆t] {J * (z(t+∆t), t+∆t), t)+ ∫ t+∆t t ℓ(z(τ ), u(τ ))dτ }
We expressed the cost function versus the minimum cost function for the interval t + ∆t ≤ τ ≤ t f where initial state is z(t + ∆t). Good mathematical conditions are lled, the second partial derivatives exist and are limited. The Taylor series development of J * gives:

J * (z(t), t) = min [u(τ ),t≤τ ≤t+∆t] { ∫ t+∆t t ℓ(z(τ ), u(τ ))dτ + J * (z(t), t), t) + [ ∂J * ∂t (z(t), t), t)]∆t + [ ∂J * ∂t (z(t), t), t)] T [z(t + ∆t) -z(t)]

+ higher orders}

For innitesimal ∆t, as small as we can use, depending on the computer facilities during the processing step, we can write:

J * (z(t), t) = min u(t) {l(z(t), u(t)∆t) + J * (z(t), t) + J * t (z(t), t)∆t + J * T z (z(t), t)[f (z(t), u(t), t)]∆t + o(∆t)}
where o(∆t) ∼ = o(∆t) p higher orders of (∆t) 2 following from integral approximations and the stop order of terms of Taylor series (dichotomy of Taylor series). Simplication f the last equation and tending ∆t -→ 0, we obtain:

J * t (z(t), t) + min u(t) l(z(t), u(t), t) + J * T z (z(t), t)[f (z(t)
, u(t), t)] = 0 t = t f provides the limit of this partial dierential equation and then:

J * (z(t f ), t f ) = κ(z(t f ), t f )
Writing the Hamiltonian as: We nally obtain the following Hamilton-Jacobi equation which can be called Hamilton-Jacobi-Bellman equation because it is based on the continuous time recurrence equation of Bellman [START_REF] Wright | Interior methods for constrained optimization[END_REF][START_REF] Pesch | A practical guide to the solution of real-life optimal control problems[END_REF][START_REF] Betts | Practical Methods for Optimal Control Using Nonlinear Programming[END_REF][START_REF] Bonnans | Using logarithmic penalties in the shooting algorithm for optimal control problems, Optimal Control Applications and Methods[END_REF][START_REF] Graichen | Feedforward control design for nite-time transition problems of nonlinear systems with input and output constraints[END_REF][START_REF] Graichen | Feedforward control design for nite-time transition problems of nonlinear systems with input and output constraints[END_REF]:

J * t (z(t), t) + H(z(t), u * (z(t), J * z , t), J * z , t) = 0
The necessary optimality condition has then been lled because the cost function J * t (z(t), t) satises the Hamilton-Jacobi-Bellman equation. We notice for the boundary conditions that some variables can be kept free. The Hamilton-Jacobi-Bellman equation will always function and can be solved. We can generally conclude that if the function satises the Hamilton-Jacobi-Bellman equation then it is the minimum cost function. This is proven and conrmed [START_REF] Athans | Optimal Control, An introduction to the theory and its applications[END_REF][START_REF] Kirk | Optimal Control Theory. An introduction[END_REF]. Solutions are obtained using appropriate solvers. Discrete approximations of the given continuous OCP problem are established and solutions obtained using recurrence relations. Exact solutions of the discrete approximation of the the Hamilton-Jacobi-Bellman equation are obtained in the state-time space regions. The latter cannot be known since the beginning of the recurrence processing solving the Hamilton-Jacobi-Bellman equation. The question is to know if we obtain an exact solution of the discretized problem or an approximate solution of the exact optimization equation. The major assumption are: the state and the control variables are constrained, the nal time t f is xed and z(t f ) is free.

In practice, it is important to solve the equation numerically, if an analytical solution is not possible. The equation can be solved explicitly. In general, it is dicult to calculate the solution. Numerical method based on viscosity solutions to the Hamilton-Jacobi-Bellman equation is performed in this section [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF]. First order of the Hamilton-Jacobi-Bellman equation is perturbed by an added diusion term; a singular perturbation parameter is used. Time and variables are discretized.

There are many methods existing in the open literature solving the Hamilton-Jacobi-Bellman equation which can be considered as ecient. We have choose a feedback control method for the computation which solve the previous Hamilton-Jacobi-Bellman equation using the dynamic programming. The Hamilton-Jacobi-Bellman equation can be expressed as:

- ∂L ∂t + sup u∈7 [-∇ z L.f (z, u, t) -l(z(t), u(t), t)] = 0
We consider a nite region Ξ for z included R n+m , 7 the control set, and the initial condition v(t f , z(t f )) = κ(z(z(t f ), t f ) or the value function or the boundary condition, and -∇ z L the gradient of L related to z. L and u are unknown.

To assess gradients or approximate gradients (or even subgradients), the nite dierence methods can be used: Quasi-Newton methods, Conjugate gradient methods, Interior point methods (a large class of methods for constrained optimization. Some interior-point methods use only (sub)gradient information, and others of which require the evaluation of Hessians), Gradient descent, Bundle method of descent, Ellipsoid method, Reduced gradient method (Frank-Wolfe), ... Finding a feedback control for the OCP is equal to solve the initial value problem. Many works are performed showing the existence and the unicity of the solution of this problem well know as a viscosity solution [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF][START_REF] Douglas | Numerical methods for convection-dominated diusion problems based on combining the method of characteristics with nite element or nite dierence procedures[END_REF][START_REF] Zhou | Verication theorems within the framework of viscosity solutions[END_REF].

The extreme value theorem of Weierstrass or Weierstrauss theorem states that a continuous real-valued function on a compact set attains its maximum and minimum value. More generally, a lower semi-continuous function on a compact set attains its minimum; an upper semi-continuous function on a compact set attains its maximum.

Numerical approximations are need for initial value problem for solving the original OCP, which are dicult to nd in the open literature in particular ecient algorithms. If we disturb the latter equation by a diusion term ε∇ 2 L where ε << 1 a singular perturbating parameter, we can write the viscosity approximation of previous equation whose principal mathematical properties are shown by Zhou [START_REF] Zhou | Verication theorems within the framework of viscosity solutions[END_REF]:

ε∇ 2 L -∂L ∂t + sup u∈7 [-∇ z L.f(z, u, t) -l(z(t), u(t), t)] = 0
This initial value problem can be solved explicitly by a time stepping schemes.

t f is xed and no information is available on L. After domain extension of Ξ which is necessary for nding the exact solution, we can write:

ε∇ 2 L -∂L ∂t -∇ z L.f(z, u * , t) -l(z(t), u * (t), t) = 0 with u * = arg sup u∈7 [-∇ z L.f(z, u, t) -l(z, u, t)]
t f could not be exactly reach and the extension domain is a transitional region where the solution satises the articial boundary conditions of the solutions of the Hamilton-Jacobi-Bellman equation. Finally, we obtain the following equations:

l(z(t), u * (t), t) = ε∇ 2 L -{ ∂L ∂t + ∇ z L.f(z, u * , t)} u * = arg sup u∈7 [-∇ z L.f (z, u, t) -l(z, u, t)]
To solve, by approximation, the previous equations (having the form of the convection-diusion equation), we initially simplify the writing form of those equations:

l(z(t), u * (z, t n ), t n ) = ε∇ 2 L -{ L(x, t n-1 ) -L(x, t n ) ∆t u * (z, t n ) = arg sup u∈7 [-∇ z L.f (z, u(z, t n ), t n ) -l(z, u(z, t n ), t n )]
the approximation of the solution at t n = 1 -n∆t (∆t > 0) with x = x + f (t, z, u(z, t n ))∆t, and L t + f.∇L is an operator acting terms of the given equations which can express the dierentiation in the characteristic direction ζ = ζ(x) so that we can write:

l(z(t), u * (t), t) = ε∇ 2 L -{ ∂L ∂t + ∇ z L.f (z, u * , t)} l(z(t), u * (t), t) = ε∇ 2 L -[1 + |f (z, u, t)| 2 ] 1 2

∂L ∂ζ(z)

with:

∂ ∂ζ(x) ∼ = [ 1 [1 + |f (z, u, t)| 2 ] 0.5 ][ ∂ ∂t + f (z, u * , t).∇]
The problem is the assessment of L(x, t n-1 ). This achieved by an extrapolation of u(z, t n ) during processing. This discretization in time associated with the new form of those equations is called the continuous in space MMOC (the modied methods of characteristics) procedure of the original coupled and nonlinear system of equations. The idea behind this method [START_REF] Cheng | The Modied Method of Characteristics for Compressible Flow in Porous Media / Numerical Treatment of Multiphase Flows in Porous Media[END_REF][START_REF] Polyanin | Handbook of First Order Partial Dierential Equations[END_REF][START_REF] Scott | The Method of Characteristics with applications to Conservation Laws[END_REF] is to use a small extension of the domain. This method of characteristics usually solves partial dierential equations. It is generally performed for rst-order equations. The method can also be applied for any hyperbolic partial dierential equation. Its power allows, inter alia, a reduction a partial dierential equation to a family of ordinary dierential equations where the solution can be obtained when initial data are given.

MMOC development by Huang et al. and by Cheng and Wang, provided the approximation of the characteristic derivative as:

[[1 + |f (z, u, t)| 2 ] 0.5 ] ∂L ∂ζ(x) ∼ = [[1 + |f (z, u, t)| 2 ] 0.5 ] L(x, t n-1 ) -L(x, t n ) [|x -x| 2 + (∆t) 2 ] 0.5
and

[[1 + |f (z, u, t)| 2 ] 0.5 ] L(x, t n-1 ) -L(x, t n ) [|x -x| 2 + (∆t) 2 ] 0.5 = L(x, t n-1 ) -L(x, t n ) ∆t
They suggested the use of the operator splitting technique providing calculation of variables of the rst equation from the second one approximating ∇L and ∇ 2 L. For any (i, j, n):

x i,j = (-1 -a + iκ, -1 -a + jκ) For a function σ, σ n i,j = σ(z i,j , t n ) a. κ is given b. t n = 1 -n∆t c. ∆t = 1 N d. zi,j = z + f n i,j ∆t e. wn-1 i,j = w(z i,j )
f. M n i,j and O n i,j are approximated values of the solution (L,u)

The denite discretization of the last Hamilton-Jacobi-Bellman equation given by the previous nonlinear system is:

l n i,j = εδ 2 M n i,j - wn-1 i,j -w n i,j ∆t z n i,j = arg sup y n ∈7 [-δw n i,j .f n i,j -l n i,j ]
Terms l n i,j and M n-1 i,j

made the system strongly coupled. It is dissociated by replacing f n i,j by f n-1 i,j and l n i,j by l n-1 i,j . Now, we can write the dissociated terms, during numerical processing, using the following algorithm steps [START_REF] Osher | High-order Essentially Nonoscillatory Schemes for Hamilton Jacobi Equations[END_REF][START_REF] Krener | Optimal model matching controllers for linear and nonlinear systems[END_REF][START_REF] Cormen | Introduction to Algorithms[END_REF][START_REF] Giegerich | A Discipline of Dynamic Programming over Sequence Data[END_REF]

: 1. Give M 0 i,j = κ i,j . 2. Calculate O 0 i,j = arg sup y 0 ∈7 [-δM 0 i,j .f 0 i,j -l 0 i,j )].
3. n ∈ [1, N ] for a given total processing number N, calculate M n i,j from:

l n-1 i,j = εδ 2 M n i,j - Mn-1 i,j -M n i,j ∆t 4. Calculate O n i,j = arg sup y n ∈7 [-δM n i,j .f n i,j -l n i,j )].
The processing stages are carried out led with an high calculation rate. The main remarks are given below.

* Discretization scheme given previously is coherent, consistent and stable.

It is a copy of Euler's method. * The matrix of the system is symmetric and positive with 1 ∆t ≫ ε.

* High dimensions could be easily performed with those latter properties.

* The OCP dimension is then reduced to a nite-dimensional nonlinear program.

* The large nonlinear program is solved by a standard a robust NLP solver according to the discretized variables.

* To solve the obtained NLP problem, we developed an AMPL model.

* The viscosity coecient ε = 10 -12 is used with this value in processing steps. * M = 81, a and b can be choose in a random way. They are dimensionindependent. Domain extension must gradually increase by steps of 10 -3 . * ∆ t=2.5 10 -3 sec is the used time step. We have also performed calculation with a time step of ∆t = 0.5 sec because radar data are updated every 0.5 sec. Our algorithms could be then interfaced with in ight and on the ground radars controlling aircraft ight paths. * Asymptotic convergence has been reached because decoupling variables as described.

To conclude, the starting OCP problem is described by a Hamilton-Jacobi-Bellman equation which is transformed to being a convection-diusion equation. The modied method of characteristics approximating the solution is used. It is then solved the problem in time and a nite-dierence in the state space with a high accuracy.

In addition, Lie-Ovsyannikov innitesimal approach applied for reduction of the corresponding Bellman equation has been described rst by Garaev [START_REF] Garaev | On one corollary from Noether theorem for two-dimensional variation problems of Mayer type[END_REF][START_REF] Bardi | A remark on the Bellman principle of optimality[END_REF] and the Noether theory of invariant variation problems. It could be suggested for the problem of optimum control [START_REF] Djordje | Noether theorem for optimum control systems[END_REF][START_REF] Ovsyannikov | Group Analysis of Dierential Equations[END_REF][START_REF] Olver | Applications of Lie groups to dierential equations[END_REF][START_REF] Garaev | Theory of invariant variation problems in the problem of optimization of dynamic problems with control[END_REF]. In the open literature, it is not usually considered. This approach is new in the area of optimal control problem related to aircraft annoyances minimization. The Bellman equation corollary can be obtained in the form of a linear partial dierential equation. The use of the equation simplies the construction of synthesizing controls.

Numerical results

We consider an aircraft landing by xing initial and nal ight conditions:

180 m/s = V min ≤ V ≤ V max = f ree 0.2 = δ x min ≤ δ x ≤ δ xmax = f ree -10 • = γ min ≤ γ ≤ γ max = f ree 0 • = α min ≤ α ≤ α max = +20 • -5 • = χ min ≤ χ ≤ χ max = +5 • -5 • = µ min ≤ µ ≤ µ max = +5 • -5 • = ϕ min ≤ ϕ ≤ ϕ max = +5 • -60 km = x min ≤ x ≤ x max = 0 km -10 km = y min ≤ y ≤ y max = +10 km 3500 m = h min ≤ h ≤ h max = f ree t min = 0 ≤ t ≤ t max = +10 min
V min represents the aircraft velocity (stall velocity). As shown, some of these parameters are kept free. Once the processing steps and calculation eciency are conrmed, their limit values are found and given. These inequalities are represented by:

a ≤ C(z(t), u(t)) ≤ b
where a and b are constant vectors. The used data in this optimization model are from an Airbus A300. The three-dimensional analysis is useful in enhancing the reliability of the optimization model applied in automatic detection of aircraft noise and in the aircraft noise features. We consider R the distance aircraft-observer:

R = (x -x obs ) 2 + (y -y obs ) 2 + h 2
where (x obs , y obs , 0) is the coordinates of the observer on the ground. OCP is discretized along its state z = (V, γ, χ, x, y, h) and control u = (α, δ x , µ) variables.

P1 to P12 are the considered observers on the ground for which noise levels have to be calculated: Location points P2 and P3 are symmetrical compared to P1 and are regarded as side points. P5-P6, P8-P9 and P11-P12 are considered as side points compared respectively to P4, P7 and P10. P1, P4, P7 and P10 are under the ight path on the ground. We minimize noise levels in the cost function previously described. The problem to solve is written as follows :

(OCP ) 3            min ϑ ϑ >= J obs j ż(t) = f (z(t), u(t)) z I 1 (0) = c 1 , z I 2 (t f ) = c 2 a ≤ C(z(t), u(t)) ≤ b (1)
where J obs j are noise levels corresponding to j xed observers. For several observers, the method is applied; SNOPT found a solution with a very high accuracy.

The discretization parameter is N = 100 points because the solution stability.

Results will conrm this state. To solve the N ε and NLP problems, we have used the AMPL model [START_REF] Ampl | A Modeling Language for Mathematical Programming[END_REF] and a robust solver SNOPT [START_REF] Khardi | Optimal approach minimizing aircraft noise and fuel consumption[END_REF][START_REF] Houacine | A Gauss Pseudospectral method for Less Noise and Fuel Consumption of Aircraft Operations[END_REF][START_REF] Gill | A large-scale smooth optimization problems having linear or nonlinear objectives and constraints[END_REF][START_REF] Waltz | KINITRO User's Manual Version 4[END_REF][START_REF] Byrd | KNITRO An Integrated Package for Nonlinear Optimization[END_REF].

They have been chosen after numerous comparisons among other standard solvers available on the NEOS optimization platform. We have used the callby-need mechanism which memorized automatically the result of the cost function in order to speed up call-by-name evaluation ().

The (OCP ) is transformed with the direct method into a NLP problem. The algorithm is adapted, it changes rules, and initializes points. The objective function has been minimized using 897 variables, 503 constraints, 500 nonlinear equalities, and 3 inequalities. Concerning the indirect method, a sequence of N ε problems (tending ε to zero). We initialize the problem N ε by centering the state and the control. Then, we initialize the Lagrange multiplicators as follows:

λ = ε(C(z, u) -a) -1 , µ = ε(b -C(z, u)) -1
For the implementation of the penalty parameter ε and computation, we used the following strategy [START_REF] Dussault | On the superlinear convergence order of the logarithmic barrier algorithm[END_REF]: For each iteration of the interior point method, the algorithm found an optimal solution. First, it should be remembered that direct and indirect methods (DIM) provided the same optimal trajectory and the same throttle setting δ x . The solution trajectory (optimal trajectory) and the control δ x are shown in gure 1. The optimization processing found a constant throttle setting δ x which corresponds to a stabilization ight or a constant ight level for the three applied methods (DIM and JHB). δ x is bang-bang between its bounds, in particular for DIM methods where its increase is made in only once to 0.6 lasting slightly more with JHB that DIM. The altitude H, having a predominant role in the noise level behavior called the cost function, decreased with three soft slopes for DIM and two for JHB. These provided three and two constant ight segments in favor of the JHB method. One trajectory stage is observed for JBH before alignment on the runway with a slope of 3 degrees accompanied by a reduction of the power settings. Angles of descent are stable as recommended by ICAO and aircraft certication [START_REF] Erkelens | Research into new noise abatement procedure for the 21st century[END_REF][START_REF] Caep | Generic presentation on the Bamanced Approach[END_REF][START_REF] Icao | Review of noise abatement procedure. Research and development and implementation results[END_REF] in favor for JHB method because the continuous descent approach with one constant stage showing the eciency and performance of the aircraft approach. On the one hand, noise level decrease is conrmed. 6.5 to 9.3 dB reduction is obtained in favor of JHB method. On the other hand, when we compare the measured noise J 0 at a distance of 2 km under the ight path for a standard trajectory approach with the level J obtained with the optimal trajectory given by JHB, change varies with the altitude of approach between 4% to 11% of J 0 -J J 0 . This is because optimization model, in particular the cost function, does not integrate all non-propulsive noise sources and because of optimization model makes noise reduction possible. The ight rate descent is varying between 896 and 1165 ft/mn which is close to the one recommended by ICAO and practiced by the airline companies (1000 ft/mn). The obtained JHB trajectory could be accepted into the airline community for a number of reasons. The soft JHB one-segment approach puts the aircraft in an appropriate envelope with margins for wind uncertainties and errors. There is no question of vortex separation and problems of intercepting a false glide-slope, given that it must be intercepted from above. With autopilot or ight director coupling, this approach would be acceptable for use in regular air carrier service.

ε k+1 = ε k /a, a > 1
Comparison between the described methods applied to our optimal control problem, conrms that the feasible errors are between e -07 and e -09 (table 2). Indeed, in terms of analysis of aircraft noise reduction because of the trajectory is optimal, the problem is more in favor of the case of several observers and the JHB method.

Feasible error (Direct method) 7.39e -07 Feasible error(Indirect method) 5.6e -012to1.50e -07 HJB computed errors 2.3e -09

Table 2 : Comparison of feasible errors

Although computing power has increased substantially making complex problems more practical for large projects, JHB optimization method oers a substantial advantage in detail over DIM methods with much less computer time and less discretization complexity. Optimization model is expected to replace empirical models for well-established applications such as predicting noise contours around airports and fuel saving. It is, now, practical for a wide range of situations where additional details are necessary. To conclude, the obtained results conrmed the good formulation of this problem of optimization and its eective resolution. They also provided good values, in particular for the ight parameters whose maximum values were maintained free during the processing. Further research is needed to include airframe noise sources, and air-brake systems. The cost function must integrate objectives like reduction of pollutant emissions linked to fuel consumption and air trac constraints.

Conclusion

The objective of this paper, qualifying the best applied numerical method solving commercial aircraft trajectory optimization model taking into account noise sources, fuel consumption, constraints and extreme limits, has been reached. First, we described detailed theoretical considerations and algorithms solving the obtained Hamilton-Jacobi-Bellman HJB equation. Second, we carried out comparisons with direct and indirect methods. The OCP problem was transformed into new equalities-constrained as a viscosity problem constituting an original dynamic system extension. Among the existing methods solving the HJB equation a feedback control method using dynamic programming has been developed. Compared to the direct and indirect methods, we show that HJB dynamic method is characterized by its fast computing times and its eciency. It provides the best optimized ight paths called the Shortest and Fastest Continuous Descent Approach (SF-CDA) which is able to reduce commercial aircraft annoyances and fuel consumption. It is a two-segment approach conrmed as an optimized ight path reducing aircraft environmental impacts. Results show that the HJB method is well appropriated for aircraft trajectory optimization problem and could be implemented. Technically, because of its processing speed and algorithm eciency, it can be better interfaced with the in-ight management system respecting airspace system regulation constraints. SF-CDA approach could benet both airlines and communities. Further research is needed to consider non-propulsive sources and air trac regulation.

Figure 1 :

 1 Figure 1: Obtained solution for 100 discretization points

  The number of nonzeros are respectively 3975 in Jacobian, and 80693 in Hessian. A locally optimal solution of the objective function has been found: the nal objective value is 196.7 with a nal feasibility error (abs / rel) of 3.18 10-3 / 1.8 10-9, and a nal optimality error (abs / rel) equal to 4.[START_REF] Scott | The Method of Characteristics with applications to Conservation Laws[END_REF] 10-6 / 7.50 10-7) through 38 function evaluations, 39 gradients and 38 Hessian evaluations. The total program time is equal to 4.341 sec (4.344 CPU time / standard PC).

Table 1 :

 1 Table below summarizes the obtained results for ε. The calculated feasibility errors versus ε

	ε	Feasible error
	1	5.6e -12
	0.2	4.0e -12
	0.04	7.3e -14
	0.008	8.8e -12
	0.0016	1.1e -11
	0.000032	2.3e -07
	6.4e -05	2.5e -07
	1.28e -05	2.1e -07