Thierry Miquel
email: thierry.miquel@enac.fr

SUBOPTIMAL LONGITUDINAL REFERENCE TRAJECTORY COMPUTATION FOR TIME BASED CONTINUOUS DESCENT OPERATIONS

come

Suboptimal longitudinal reference trajectory Computation for time based continuous descent operations

Introduction

Nowadays, environmental impact and efficiency have become the two very important aspects in aviation industry after safety. New operations, such as Continuous Descent Operation (CDO), can significantly reduce the noise impact of landing aircraft by keeping them longer at higher altitudes and by avoiding steps during descent. Such new operations need a collaborative work between airlines and Air Navigation Service Providers (ANSPs) to be defined and operated. CDO is defined as "an aircraft operating technique aided by appropriate airspace and procedure design and appropriate ATC clearances enabling the execution of a flight profile optimized to the operating capability of the aircraft, with low engine thrust settings and, where possible, a low drag configuration, thereby reducing fuel burn and emissions during descent. The optimum vertical profile takes the form of a continuously descending path, with a minimum of level flight segments only as needed to decelerate and configure the aircraft or to establish on a landing guidance system (e.g. ILS)" [START_REF] Miquel | Continuous Descent Operations (CDO) Manual[END_REF].

Since 2008 European stakeholders have initiated an European CDO implementation program. Up to now, Basic CDO is already in operation at UK Heathrow Airport [3]. In Sweden the European project NUP2 has enabled SAS and LFV to operate 4D flight paths, or green approaches, at Stockholm-Arlanda airport with beneficial impacts on environment. Thanks to the collaboration between Air France and French ANSP0, CDOs are operated routinely at Marseille airport in France [START_REF] Guibaud | Evaluation des descentes continues à Marseille-Provence, Séminaires de la DSNA[END_REF].

In United States a program named Partnership for AIR Transportation Noise and Emission Reduction (PARTNER) conducted tests at Louisville

International Airport in 2002 and at Los Angeles International Airport in 2007 [START_REF] Cao | Evaluation of Continuous Descent Approach as a Standard Terminal Airspace Operation[END_REF].

NASA and FAA have been involved in extensive efforts to develop advanced concepts, technologies and procedures for the Next Generation Air Transportation System (NextGen) later on. One aim of NextGen is to develop ground-side automation systems to assist controllers in strategic planning operations. The En Route Descent Advisor (EDA) is one of the Center TRACON Automation System (CTAS) decision support tools under development at the NASA Ames Research Center. EDA generates maneuver advisories for arrival aircraft to meet scheduled arrival times at the arrival meter fix, sometimes 20 -25 minutes ahead of the aircraft's scheduled meter fix arrival time [START_REF] Mueller | Effect of Uncertainty on En Route Descent Advisor (EDA) Predictions[END_REF]. A research [START_REF] Vaddi | 4D Green Trajectory Design for Terminal Area Operations Using Nonlinear Optimization Techniques[END_REF] has also been done to develop ground-side automation to enable 4D-Trajectory-Based Operations (4DTBO) in terminal airspace. This research illustrates a computational framework for the design of 4D-Trajectories (4DTs) based on fundamental flight mechanics and nonlinear trajectory optimization techniques with sample scenarios and open-source models. AIRE Project (Atlantic Interoperability Initiative to Reduce Emissions) is a joint initiative by the European Commission and the FAA to improve energy efficiency and aircraft noise. AURORA project was a project implemented by Airbus, Scandinavian Airlines International (SAS), Swedish ANSP LFV and Stockholm Arlanda Airport. They conducted Continuous Descent Approaches for the first time on SAS transatlantic flights using an Airbus A330.

Another new ATC technique called Point Merge System (PMS) [START_REF] Boursier | Merging Arrival Flows Without Heading Instructions, 7th USA/Europe Air[END_REF] aims to facilitate the merging of traffic from a number of Area Navigation (RNAV) arrival routes. The technique is based upon aircraft flying a quasi-arc, up to 30NM long, with a radius of more than 20NM from the designated merging point. Each arc has a published height that the aircraft must have reached before establishing on the arc and a predefined speed to fly it. In general the arc nearest to the merging point has the highest height while the other has the lowest height so that the external sequencing leg is free from traffic from the internal sequencing leg during descent of the traffic. On April 2011 Oslo became the world's first airport to implement a PMS in their airspace. Other countries and airports will follow.

In this paper, the task of merging an aircraft over a specified meter fix is addressed through a novel clearance in which air traffic control clears an aircraft to track an ad-hoc computed reference trajectory. This enables the aircraft to merge at a specified meter fix at a given time. This kind of application may be envisioned as an enhancement of the use of the ground based Arrival MANager (AMAN), which is a tactical controller assistance system enabling the computation of rendez-vous time at meter fix to meet the runway capacity and absorb the traffic [START_REF] Bourdon | Maestro: Sequenced Traffic Flow Tool[END_REF]. The time computed by the AMAN is envisioned to be used by the air traffic controller to request the aircrew to overfly the meter fix at the desired time. In comparison with current operations, the change is that the controller would communicate decisions on traffic flow organization at a higher level to the pilot rather than requiring the controller to calculate and communicate lower-level guidance instructions.

While keeping the controller responsible for making the traffic flow decisions, the envisioned application involves new avionics capabilities for merging operations including the tracking of height and speed. Thus air traffic control improvement is achieved through a greater involvement of pilots in cooperation with air traffic controllers. This type of application is clearly in the scope of the 4D trajectory concept promoted by the European programme SESAR [START_REF]The ATM Target Concept -D3[END_REF] and by the US programme NextGen [START_REF] Merkle | 3D Path Arrival Management (PAM) Project, High Density Airport Time-based RNAV/RNP NextGen / JPDO Demonstrations, REACT Workshop[END_REF].

The main benefit expected from this application [START_REF] Miquel | 4D Green Approach Trajectory: Illustrative Assessment On Toulouse Approach[END_REF] is to improve flight efficiency by more precise maneuvering resulting from onboard capabilities, and also noise abatement and fuel saving. More precise maneuvers are also expected to increase sector capacity. Indeed moving from radar vectoring to monitoring pre-computed trajectories would contribute to decrease controller's workload, and therefore to increase sector capacity. This paper presents a new methodology to compute a reference trajectory for time based continuous descent operations. The aircraft is considered as a point mass model and focuses on aircraft longitudinal motion including known wind. As far as time constrained operations are assumed, final time as well as final altitude and along track distance to be flown are imposed. We propose a new methodology to compute a reference calibrated airspeed (CAS) and a reference vertical speed to achieve imposed final position and altitude at a prescribed time. The optimal control problem consists in minimizing fuel consumption while ensuring that maximum longitudinal and normal accelerations remain lower than the acceptable level for civil flights. The computed trajectory is a time parametrized trajectory which will be used as a reference trajectory by some envisioned tracking controller installed on board the aircraft. Nevertheless the design of the tracking controller is out of the scope of this paper. Numerical simulations using Bada 3.11 [1] are provided to illustrate the suboptimal trajectory generation method and achieved results.

The structure of the paper is the following: first the equations of motion which will be used are presented as well as the constraints to be satisfied. Then the notion of differentially flat system is introduced and the dimensionless parameters which will be used are defined. The core of the paper consists in describing the proposed methodology to compute the reference trajectory, firstly without and secondly with optimization. Then a scenario is proposed and simulation results are provided.

Equation of motion

In the following, the aircraft is modeled as a point mass and the flight trajectory is strictly confined in a vertical plane on a non-rotating, flat earth. In addition, the aircraft is assumed to fly in an time varying atmospheric wind field comprising of both horizontal and vertical components. The notations are depicted on Figure 1:

Figure 1. aircraft forces and angles

The notations which will be used are the following:

• y is the horizontal position of the aircraft • h is the vertical position (altitude) of the aircraft • V is the airspeed (TAS: true airspeed) of the aircraft

• θ is the pitch angle • q is the pitch rate =q ( 1
)
• α is the angle of attack

• γ is the flight path angle

=- (2)
• g is the gravitational acceleration

• m is the aircraft mass

• L is the lift force L= 1 2  h S V 2 C L (3)
• ρ(h) is the air density which depends on altitude h

• S is the reference wing surface area

• C L is the lift coefficient. • D is the drag force D= 1 2 hS V 2 C D (4)
• C D is the drag coefficient. We will assume a parabolic drag coefficient which is specified as a function of the lift coefficient

C L C D =C D0 C D2 C L 2 (5)
• Following Bada [1] the lift coefficient C L is determined assuming that the lift force L compensates exactly weight mg:

L≈m g ⇔ C L ≈  h S V 2 2 m g (6)
• F is the thrust, which is assumed to be in the opposite direction of the drag

D y V F γ θ α D mg L h 0 w y w h
• w y (t) is the horizontal time varying component of the wind field

• w h (t) is the vertical time varying component of the wind field

Aircraft control is denoted u and is assumed to be the vector gathering thrust F and flight path angle γ:

u=[ F ] T (7)
The state vector will be denoted x; it gathers horizontal position y, vertical position h and airspeed V:

x =[ y h V ] T (8)
The equations of motion which will be used are the following:

{ ẏ = V cosw y t  ḣ = V sin w h t  V = F -D m -g sin   -ẇ y cos -ẇh sin  (9)
Constraints to be satisfied

The following constraints are imposed on the vertical axis:

{ h0= h 0 hT = h f ḣ0= 0 ḣT = 0 (10)
Basically the preceding constraints indicate that the aircraft starts the maneuver at altitude h 0 and finishes it at altitude h f once the maneuver duration is accomplished at time T. Moreover the aircraft is assumed to be level at the beginning and at the end of the maneuver.

As far as the horizontal axis is concerned, the following constraints are imposed:

{ y 0= 0 y T = y f ẏ 0= V 0 cos  asin  w h 0 V 0  w y 0  ẏ T = V f cos  asin  w h T  V f  w y T  (11)
Basically the preceding constraints indicate that the aircraft shall fly a distance equal to y f once the maneuver is finished at time T. Moreover the airspeed is equal to V 0 at the beginning of the maneuver and to V f at the end of the maneuver.

It can be seen that vertical component of the wind w h appears in the two last constraints of [START_REF] Merkle | 3D Path Arrival Management (PAM) Project, High Density Airport Time-based RNAV/RNP NextGen / JPDO Demonstrations, REACT Workshop[END_REF]. Indeed this comes from the fact that as soon as the vertical component w h of the wind is not null, the flight path angle is not null because we impose that the time derivative of the vertical position h is zero at the beginning and at the end of the maneuver, that is ḣ=0 ; inserting this value in the second equation of ( 9) leads to the value of the flight path angle: t =-asin  w h t  V t   . Inserting this value in the first equation of ( 9) and specializing the result to the initial and final conditions leads to the last constraints in [START_REF] Merkle | 3D Path Arrival Management (PAM) Project, High Density Airport Time-based RNAV/RNP NextGen / JPDO Demonstrations, REACT Workshop[END_REF].

Constraints [START_REF]The ATM Target Concept -D3[END_REF] and [START_REF] Merkle | 3D Path Arrival Management (PAM) Project, High Density Airport Time-based RNAV/RNP NextGen / JPDO Demonstrations, REACT Workshop[END_REF] associated with the set of coupled ordinary differential equations ( 9) constitute a Two-point Boundary Value Problem (TPBVP).

Differentially flat system

Flatness was originally introduced 1995 [START_REF] Fliess | Flatness and Defect of Non-linear Systems: Introductory Theory And Examples[END_REF]. Roughly speaking, a flat system is a square input/output system (i.e. a system having the same number of inputs and outputs) for which there exists an output vector such that all states and inputs can be expressed in terms of this output vector and its derivatives. More precisely, consider the following nonlinear system:

ẋ = f  x ,u  (12)
This system is said differentially flat if one can find an output z of the form: z= x , u ,... , u    [START_REF] Fliess | Flatness and Defect of Non-linear Systems: Introductory Theory And Examples[END_REF] where u (s) denotes the s-th order derivative of control u with respect to time, and such that:

{ x= 0  z ,... , z s  u= 1 z ,... , z  s , z  s1  ( 14
)
Output z is called flat output. In addition, system (13) is said Lie-Bäcklund equivalent to the following system (called trivial system) where vector v is the new input, also called pseudo control:

v =z  s1 (15)
Imposing a given arbitrary trajectory to z yields a trajectory for all the system variables x and u, without integrating any differential equation. Remark that the time derivatives involved in the above formulas do not imply to take derivatives of noisy signals since it involves precomputed open-loop time functions.

We will assume that w y (t) and w h (t) are exogenous parameters as well as aircraft mass m and initial altitude h 0 . In the following we will check that the following variables z 1 and z 2 are flat outputs for the system (9):

{ z 1 t = y t -∫ 0 t w y  d  z 2 t =h t -h 0 -∫ 0 t w h  d  (16)
By taking the first time derivative of z 1 and z 2 we get the expression airspeed V and flight path angle γ from (9) as follows:

{ V =  ż1 2  ż2 2 =atan  ż2 ż1  (17)
Using (17) we can express the time derivative of airspeed V:

V = ż1 z1  ż2 z2 V ⇔ V = ż1 z1  ż2 z2  ż1 2  ż2 2 (18)
Finally the use of ( 9) leads to the expressions of thrust F according to the flat outputs z 1 and z 2 and their derivatives:

{ F = m V g sin D m  ẇy cos  ẇh sin   (19)

Dimensionless parameters

We will denote by τ the ratio between the actual time t and the duration T of the time constrained maneuver:

= t T ( 20
)
In addition we will denote by r'(τ) the derivative of the trajectory r(τ) with respect to dimensionless parameter τ and we will identify ṙ t to r ' T  :

d d t r t = ṙ t ≡r ' T  = d d  r T  (21)
As a consequence, the relationships between the integral and the derivative of r'(τ) with the integral and the derivative of r(t) are the following:

= t T ⇒ { ∫ 0 T ṙ t dt =T ∫ 0 1 r ' d  rt = 1 T r ' '  (22)

Reference trajectory computation without optimization

The following expression will be used to build the reference trajectories:

{ r =a 0  a 1 b arctan b  a 2 b arctan -1 barctan  b   r '  = d r d  =a 0  a 1 1b  2  a 2 1b -1 2 r ' ' = d 2 r d  2 =- 2 a 1 b  1b 2  2 - 2 a 2 b -1 1b-1 2  2 (23)
It is clear that whatever parameters a 0 , a 1 and a 2 the reference trajectory r(τ) satisfies:

r 0=0 (24)
As far as parameter b is set, the reference trajectory depends on three parameters, namely a 0 , a 1 and a 2 . They will be set according to the initial and final constraints [START_REF]The ATM Target Concept -D3[END_REF] and [START_REF] Merkle | 3D Path Arrival Management (PAM) Project, High Density Airport Time-based RNAV/RNP NextGen / JPDO Demonstrations, REACT Workshop[END_REF].

More specifically to satisfy constraints [START_REF]The ATM Target Concept -D3[END_REF] on vertical position we set:

z 2 t =r  T  , a 0h , a 1h , a 2h ,b h  (25)
Assuming that parameter b = b h is given we shall have:

z 2 t =ht -h 0 -∫ 0 t w h  d  ⇒ { z 2 T =h f -h 0 -∫ 0 T w h t dt ż2 0= ḣ 0 -w h 0 =-w h 0 ż2 T = ḣ 0-w h T =-w h T  (26)
Constraints [START_REF] Merkle | 3D Path Arrival Management (PAM) Project, High Density Airport Time-based RNAV/RNP NextGen / JPDO Demonstrations, REACT Workshop[END_REF] on horizontal position are treated similarly. Nevertheless it is worth noticing that the derivative of reference trajectory (21) is not used to set the horizontal component of the true airspeed (TAS) but the horizontal component of the calibrated airspeed (CAS), which will be denoted CAS y :

ż1 t ≡CAS y t =r ' T  ,a 0y ,a 1y ,a 2y ,b y  (27)
We have the following equivalence:

TAS  ż1 t , ht ≡V cos  (28)
CAS / TAS conversion formulas are provided within Bada documentation [1].

Assuming that parameter b = b y is given, parameters a 0 , a 1 and a 2 are set as follows:

{ z 1 T = y fcas ż1 0=CAS V 0, h 0  ż1 T =CAS V f , h f  (29)
Where y f'cas is computed such that:

∫ 0 T TAS  ż1 t , ht dt= y f -∫ 0 T w y t dt (30)

Reference trajectory computation with optimization

The preceding computations have been done assuming that parameter b is known, for both the vertical and the horizontal motion. In that section we will assess the influence of parameter b on the reference trajectories and we will use it to minimize fuel consumption while ensuring that the maximum longitudinal and normal accelerations remain lower than the acceptable level for civil flights.

The fuel consumption is computed as follows, where f(t) is the fuel flow:

∫ 0 T f t dt (31)
Following Bada 3.11 [1] model, the fuel flow f is calculated using thrust F and the thrust specific fuel flow η; it is limited to the minimum fuel flow f min if necessary:

f =max   F , f min (32)
For a jet aircraft the thrust specific fuel flow η and the minimum fuel flow f min have the following expression where (C f1 , C f2 , C f3 , C f4 ) are constant parameters which depends on the aircraft type:

{ =C f1  1 V C f2  f min =C f3  1-h C f4  (33)
Moreover, for civil flights maximum longitudinal and normal accelerations are limited in order to take into account passenger comfort. More specifically, and following Bada 3.11 [1], true airspeed (TAS) acceleration (longitudinal) and normal acceleration shall be bellow the following thresholds:

{ V ≤2 ft / s 2 V ≤5 ft / s 2 (34)
Consequently parameter b used for the vertical reference trajectory, which is denoted b h , and parameter b used for the horizontal reference trajectory, which is denoted b y , will be computed in order to satisfy the Two-point Boundary Value Problem (TPBVP) while minimizing fuel consumption and ensuring that the maximum longitudinal and normal accelerations remain lower than the acceptable level for civil flights:

b h ,b y =arg min ∫ 0 T f t dt s.t. { ẏ=V cosw y ḣ=V sin w h h 0=h 0, ḣ 0=0 h T =h f , ḣT =0 y 0=0, V 0 =V 0 y T =y f , V T =V f and { V ≤2 ft / s 2 V ≤5 ft / s 2 (35)
Computations have been done using the routine fminsearch() within the open source platform for numerical computation Scilab [START_REF]Scilab, the open source platform for numerical computation[END_REF]. The cost which has been minimized is the following:

∫ 0 T f t dt  2 10 5  V -2  2 V -5 2  (36)

Scenario

The scenario which has been simulated is the following: the aircraft starts the descent at FL 140 with a conventional airspeed (CAS) of 220 kts; the final approach fix (FAF) is assumed to be 35 NM ahead and shall be overflown at 2500 feet with a conventional airspeed of 170 kts. The maneuver duration is set to T = 540 sec (9 min). The aircraft is assumed to descent in clean configuration.

As far as wind is concerned, we will assume a constant vertical wind w h of +1 kts; the horizontal component w y is assumed to be time dependent: starting at +20 kts this value will decrease to -20 kts after 120 sec of flight; this may for example simulate the change in wind direction when the aircraft turns from downwind to base leg and final approach segment.

w y = { 20 kts ∀t t s 20 cos    t-t s t d  kts ∀t s ≤tt s t d -20 kts ∀ t ≥t s t d where { t s =120 sec t d =30 sec (37)
The simulation have been done using Bada 3.11 [1]. Two aircraft types have been simulated: an Airbus A320 and a Boeing B777-300. For the Airbus A320 model, the initial value of the mass has been set to 46,600 kg whereas for the Boeing B777-300 the initial value of the mass has been set to 187,540 kg (in fact 1-p m min  p m max where m min and m max are the minimum and the maximum mass and p = 20%).

Simulation results

The following figures present for the Airbus A320 model some simulation results assuming no optimization, that is b h = b y = 1 when computing the reference trajectories of the flat outputs z 1 and z 2 .

The following figure represents the reference calibrated airspeed (CAS in kts) with respect to time (in sec): As expected, the conventional airspeed changes from 220 kts to 170 kts. We can see on the preceding figure that the calibrated airspeed starts to increase before decreasing smoothly; the initial increase may be disturbing for the pilot during the approach phase.

The flight level with respect to time (in sec) is depicted hereafter: We can see on the preceding figure the sharp decrease of the ground speed starting at t = 120 sec; this is obviously the effect of the changing wind described in (37). Moreover we have checked that the distance flown is 35 NM, as expected (this is a question of integrating the ground speed).

Finally the aircraft vertical speed (V z in fpm) is shown on the following figure: We can see that the initial and the final values of the vertical speed in not zero. Indeed, we have assumed a constant vertical wind of +1 kts, that is about 101 fpm. In order to have a level flight at the beginning and at the end of the scenario, the vertical speed of the aircraft shall consequently be set at -101 fpm, as indicated by constraints (26). The calibrated airspeed stays close to its initial value at 220 kts before decreasing smoothly towards its final value, that is 170 kts.

As far as flight level with respect to time (in sec) is concerned, we can see on the following figure that compared to the case without optimization of parameters b y and b h the descent is more straight.

Figure 7. flight level for optimized values of bh and by

The following figures depicts the ground speed (GS in kts) with respect to time (in sec): Finally it can be seen on the following figure that the descent from FL 140 to 2500 feet is accomplished with a vertical speed of -1400 fpm: (39)

The shape of the reference calibrated airspeed (CAS) and vertical speed are very similar to those obtained for the Airbus A 320 model.

Conclusion

This paper addresses a specific aspect of air traffic control services, namely the achievement of an orderly and expeditious flow of air traffic under time constrained continuous descent approach. This paper presents a new methodology to compute a reference calibrated airspeed (CAS) and a reference vertical speed for time constrained descent which satisfy length and endpoint constraints.

The proposed approach is based on the shaping of a parametric reference airspeed, which is applied to the controlled calibrated airspeed (CAS) of the aircraft and its vertical speed. An appropriate choice of the available degree of freedom when computing the reference calibrated airspeed (CAS) of the aircraft and its vertical speed (namely parameter b) enables to minimize fuel consumption while ensuring that the maximum longitudinal and normal accelerations remain lower than the acceptable level for civil flights.

Simulation results which include time dependent wind speed illustrate the efficiency of the proposed design with respect to the time and speed constraints.

Future developments include the test of the robustness of the proposed design with respect to unexpected wind. This may be addressed through the periodic update of the computation of reference calibrated airspeed and vertical speed as well as the design of some tracking controller installed on board the aircraft.

The proposed approach can be extended to the case where constraints over multiple fixes are imposed. Indeed, the way to compute the reference speed takes explicitly into account the initial and final values of aircraft's position and speed and can thus be used to accommodate the reference speed to each segment of flight.

Figure 2 .

 2 Figure 2. calibrated airspeed for bh = by = 1

Figure 3 .

 3 Figure 3. flight level for bh = by = 1 As expected, the aircraft descents smoothly from FL 140 with a conventional airspeed (CAS) of 220 kts to 2500 feet. The next figure depicts the ground speed (GS in kts) with respect to time (in sec):

Figure 4 .

 4 Figure 4. ground speed for bh = by = 1

Figure 5 .

 5 Figure 5. aircraft vertical speed for bh = by = 1

For

  this simulation without optimization of parameters b y and b h , the maximum longitudinal acceleration is 0.53 ft/sec 2 whereas the maximum normal acceleration is 0.28 ft/sec 2 . The estimated estimated burnt fuel is 96.7 kg. Now the same figures are shown but the optimization process on parameters b y and b h is used when computing the reference trajectories of the flat outputs z 1 and z 2 . This leads to the following values for b h and b y :{ b y =335.1 b h =36903.6(38)As expected the maximum longitudinal acceleration is 2 ft/sec 2 whereas the maximum normal acceleration is 5 ft/sec 2 for this simulation. The estimated estimated burnt fuel is 74.8 kg, that is a gain of more than 20 kg compared to the non optimized coefficients!The following figures present for the Airbus A320 model the reference calibrated airspeed (CAS in kts) assuming optimized values for parameters b h and b y :

Figure 6 .

 6 Figure 6. calibrated airspeed for optimized values of bh and by

Figure 8 .

 8 Figure 8. ground speed for optimized values of bh and by

Figure 9 .

 9 Figure 9. aircraft vertical speed for optimized values of bh and by For the Boeing B777-300 model, the estimated fuel consumption is 270.9 kg without optimization whereas is decreases to 222.6 kg with optimized coefficients. For this aircraft and the associated scenario, the optimized values for b h and b y are the following: { b y =335.4 b h =36918.9(39)

Acknowledgements

The author wishes to thank P. Louyot for his comments and suggestions and C. Bertrand for her support.

34th Digital Avionics Systems Conference September [START_REF] Fliess | Flatness and Defect of Non-linear Systems: Introductory Theory And Examples[END_REF][START_REF]Scilab, the open source platform for numerical computation[END_REF][15][16][17]2015