Olivier Cots

Joseph Gergaud

Damien Goubinat

Boris Wembe

SINGULAR VERSUS BOUNDARY ARCS FOR AIRCRAFT TRAJECTORY OPTIMIZATION IN CLIMBING PHASE

Keywords: 2020 Mathematics Subject Classification. 49K15, 49M05, 70Q05. 2022 Optimal control with state constraints, singular arcs, geometric control, homotopy method, aircraft trajectory optimization

published or not. The documents may come

Singular versus boundary arcs for aircraft trajectory optimization in climbing phase

Introduction

The climbing phase is a normalized phase where most of the civil aircrafts follow the CAS/Mach procedure. This procedure splits the climbing phase in two parts, the first one is an arc at constant Calibrated AirSpeed (CAS) and the second one is an arc at constant Mach number. The CAS constant is chosen among the ones which maximize the climbing slope, the vertical speed or the climbing rate [START_REF] Hull | Fundamentals of Airplane Flight Mechanics[END_REF][START_REF] Verriere | Mécanique du[END_REF]. The Mach constant depends most of the time on the cruise Mach number but this value is also closed to the one which maximizes the vertical speed. The climbing path resulting from a CAS/Mach couple may be seen as an approximation of an optimal flight, minimizing a cost made by a convex combination of the time-to-climb and the fuel consumption, that we call the cost index. Even though historically a lot of studies was done to optimize the climbing flight path [START_REF] Barman | Fixed-range optimum trajectories for short-haul aircraft[END_REF][START_REF] Miele | General solutions of optimum problems in nonstationnary flight[END_REF][START_REF] Miele | Optimum flight paths of turbojet aircraft[END_REF], the CAS/Mach procedure is still used mainly because of its simplicity. The CAS and the Mach numbers are fairly easy quantities to compute since they are defined thanks to the differential pressure which is directly available from the Pitot tubes.

In this article, we study first the minimum time-to-climb trajectory submitted to state constraints which are related to the CAS and the Mach number. The dynamics of this climbing phase is depicted by the following four-dimensional dynamical system:

dh dt = v sin γ (1)
dv dt = T (h) m - 1 2 ρ(h)Sv 2 m C d (C l ) -g 0 sin γ (2)
dm dt = -C s (v) T (h) (3)
dγ dt = 1 2 ρ(h)Sv m C l - g 0 v cos γ, (4)
where the state variable is composed of the altitude, the true air speed, the mass and the flight path angle of the aircraft. We refer to [START_REF] Goubinat | Contrôle optimal géométrique et méthodes numériques : application à un problème de montée d'un avion[END_REF] for more details about the dynamics. The altitude h is given in meter (m), the true air speed v in meter per second (m.s -1 ), the mass m in kilogram (kg) and the flight path angle γ in radian (rad). In this model, the lift coefficient C l may be considered as the control variable. The BADA [START_REF] Poles | Base of Aircraft DAta (BADA) aircraft performance modelling report[END_REF] model is chosen to represent the aircraft performance model and provides the following expressions:

C s (v) := C s,1 1 + v C s,2 , T (h) := C T,1 1 - h C T,2 + h 2 C T,3 , C d (C l ) := C d,1 + C d,2 C 2 l ,
where the constants C T,i , C s,i and C d,i depend on the flight phase and on the aircraft. The International Standard Atmospheric (ISA) model is used to represent the atmospheric model from the sea level till the end of the tropopause which here is considered at 11 000 meters. This model provides the evolution of the pressure P , the temperature Θ and the air density ρ with respect to the altitude through the following expressions:

P (h) := P 0 Θ(h) Θ 0 g 0 βR
, Θ(h) := Θ 0 -βh, ρ(h) := P (h) RΘ(h) .

The remaining data are positive constants: g 0 is the gravitational constant at the sea level, S the wing area, R the specific constant of air, β the thermal gradient and P 0 , Θ 0 represent the pressure and the temperature at the sea level. See Table 1 for the chosen values of the constant parameters for the numerical experiments.

Taking into account the flight path angle γ in the dynamics introduces numerical instabilities [START_REF] Cots | Direct and indirect methods in optimal control with state constraints and the climbing trajectory of an aircraft[END_REF][START_REF] Cots | Time-optimal aircraft trajectories in climbing phase and singular perturbations[END_REF] which are known as singular perturbations [START_REF] Ardema | Singular perturbations in flight mechanics[END_REF][START_REF] Calise | Singular Perturbations and Time Scales in Guidance and Control of Aerospace Systems: A Survey[END_REF][START_REF] Moissev | Problèmes mathématiques d'analyse des systèmes[END_REF]] and which come from large time constant differences between the state variables. A comparison of the time constants [START_REF] Ardema | Singular perturbations in flight mechanics[END_REF][START_REF] Goubinat | Contrôle optimal géométrique et méthodes numériques : application à un problème de montée d'un avion[END_REF] shows that the dynamics (1)-(4) contains slow (the mass m) and fast (the flight path angle γ) variables. The altitude h and the true air speed v are fast compare to the mass but slow compare to the flight path angle. In this article, we consider h and v as slow variables. Roughly speaking, the singular perturbation phenomenon arise when a system of differential equations contains at least one small parameter multiplying the derivative of one or more state variables. To emphasize the presence of a singular perturbation in the flight path angle dynamics, we introduce a parameter ε > 0 such that eq. ( 4) is replaced by:

ε dγ dt = 1 2 ρ(h)Sv m C l - g 0 v cos γ. ( 4

')

Let (P ε t f ) denote the minimum time-to-climb problem with the additional artificial parameter ε > 0. From the control theory, for a fixed value of ε > 0, the candidates as minimizers are selected among a set of BC-extremals, solution of a Hamiltonian system given by the Pontryagin Maximum Principle (PMP) [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]. The application of the PMP leads to define a Boundary Value Problem denoted (BVP ε ), in terms of state and adjoint variables, which can be solved using indirect multiple shooting methods [START_REF] Bulirsch | Introduction to numerical analysis[END_REF]. It is well known that multiple shooting increases numerical stability and a good alternative would be to use direct multiple shooting [START_REF] Bock | A multiple shooting algorithm for direct solution of optimal control problems[END_REF] to solve the optimal control problem. Difficulties in solving the singularly perturbed boundary value problem (BVP ε ) arise because the solution exhibits narrow regions of very fast variation. In this case, it may be difficult to determine a mesh for the nodes of the multiple shooting method that will give an accurate numerical solution and to give an initial guess that will lead to the convergence of the underlying Newton-like algorithm. In this context, an alternative approach would be to solve the boundary value problem using automatic mesh refinement [START_REF] Cash | The numerical solution of nonlinear two-point boundary value problems using iterated deferred correction -a survey[END_REF][START_REF] Cash | An automatic continuation strategy for the solution of singularly perturbed nonlinear boundary value problems[END_REF], or to use homotopy techniques [START_REF] Cots | Time-optimal aircraft trajectories in climbing phase and singular perturbations[END_REF][START_REF] Cots | Homotopic Approach for Turnpike and Singularly Perturbed Optimal Control Problems[END_REF].

Another point of view to deal with singular perturbations is to approximate the solution using asymptotic expansions. An overview of methods to deal with singular perturbations may be found in [START_REF] Ardema | Singular perturbations in flight mechanics[END_REF][START_REF] Calise | Singular Perturbations and Time Scales in Guidance and Control of Aerospace Systems: A Survey[END_REF][START_REF] Moissev | Problèmes mathématiques d'analyse des systèmes[END_REF]. We consider in this article only the zero-order term of the asymptotic expansion. To uniformly approximate the solution of a singularly perturbed boundary value problem, one has to compute at least two different approximations (the so-called inner and outer solutions) which are accurate only for part of the range of the time variable, and usually they are valid on different length-scales. The Method of Matched Asymptotic Expansion consists in computing the inner and outer solutions and then in matching them to get a uniform approximation which satisfies the boundary conditions. In this article, we consider only the zero-order approximation of the outer solution since it has been shown for the minimum time-to-climb problem with no state constraints that it is a good approximation of the solution [START_REF] Cots | Time-optimal aircraft trajectories in climbing phase and singular perturbations[END_REF]. Note that in the aircraft dynamics the zero-order reduction is equivalent to the quasi-steady approximation of the flight which is commonly performed [START_REF] Espin | Aircraft Trajectory Optimization Using Singular Optimal Control Theory[END_REF][START_REF] Hull | Fundamentals of Airplane Flight Mechanics[END_REF][START_REF] Merrer | Optimisation de trajectoire d'avion pour la prise en compte du bruit dans la gestion du[END_REF][START_REF] Nguyen | Singular arc time-optimal climb trajectory of aircraft in a two-dimensional wind field[END_REF][START_REF] Verriere | Mécanique du[END_REF].

The reduced-order optimal control problem denoted (P t f ) is obtained putting ε to 0 in eq. (4') and considering γ as the new control variable. The problem (P t f ) has one state variable less than (P ε t f ), ε > 0, and can be tackled with the tools from geometric optimal control theory. We refer to [START_REF] Ardema | Solution of the Minimum Time-to-Climb Problem by Matched Asymptotic Expansions[END_REF][START_REF] Ardema | Singular perturbations in flight mechanics[END_REF][START_REF] Goubinat | Contrôle optimal géométrique et méthodes numériques : application à un problème de montée d'un avion[END_REF] for details about the reduction process. In particular, one important result is that the reduced-order boundary value problem obtained from (BVP ε ) putting ε to 0 is equivalent to the boundary value problem obtained after applying the maximum principle to the reduced-order dynamical system obtained from (1)-(4') putting ε to 0, when we consider at the end γ as the new control variable. Besides, the reduced-order dynamical system with γ as the control variable is affine with respect to the control if we consider the small-angle approximation, that is if we replace cos γ by 1 and sin γ by γ assuming γ is small. In this case, the associated pseudo-Hamiltonian is of the form H(x, p, u) = H 0 (x, p) + u H 1 (x, p) where x := (h, v, m) is the reduced state, p is the adjoint variable, u the control (i.e. the flight path angle γ), and where H 0 and H 1 are two Hamiltonian lifts. In addition, to complete the definition of the optimal control problem, the control has to satisfy a constraint of the form u ∈ [u min , u max ] and the state has to satisfy two constraints denoted φ(x) ≤ φ max and ψ(x) ≤ ψ max , where φ represents the CAS, ψ represents the Mach, and where φ max and ψ max may be seen as given speed limitations.

Main contributions and results. The minimum time-to-climb problem with no state constraints, or from another point of view, with state constraints but with φ max and ψ max big enough, is analyzed in [START_REF] Cots | On the minimum time optimal control problem of an aircraft in its climbing phase[END_REF][START_REF] Cots | Time-optimal aircraft trajectories in climbing phase and singular perturbations[END_REF]. In [START_REF] Cots | On the minimum time optimal control problem of an aircraft in its climbing phase[END_REF], the influence of the initial mass m 0 and the final true air speed v f is studied and it is shown in particular that the trajectories are of the form σ ± σ s σ ± , where σ -represents a bang arc where u(t) = u min along the arc, where σ + represents a bang arc with u(t) = u max and σ s represents a singular arc where u(t) ∈ (u min , u max ). The first result of this article is that these trajectories from the state unconstrained case are compliant with the so-called Maximum Operating Speed (VMO) and Maximum Operating Mach (MMO), that is φ(x(t)) < VMO and ψ(x(t)) < MMO along the trajectories. Thus, with the constraints VMO and MMO, the trajectories are of the form Bang-Singular-Bang. Another contribution of this article is to analyze the influence of the bounds φ max and ψ max on the structure of the trajectories for the minimum time-to-climb problem. We propose a methodology introduced in [START_REF] Cots | Geometric and numerical methods for a state constrained minimum time control problem of an electric vehicle[END_REF], based on differential homotopy [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF] and geometry techniques [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF] to classify the different structures with respect to φ max and ψ max . Remarkably, the CAS/MACH procedure appears in our classification. In a final part, we present a comparison between trajectories that follow the simple CAS/Mach procedure and minimum cost index solutions (with VMO and MMO bounds), this criterion being commonly used in operations. Our numerical analysis shows that the CAS/MACH trajectories are close to the optimal solutions up to 0.3%, on average, regarding to the cost index criterion (see Table 2 for details). These optimal solutions being of the form Bang-Singular-Bang.

The paper is organized as follows. The physical model with the statements of the optimal control problems are given in Section 2.1 while the necessary conditions of optimality for both the minimum time-to-climb and the fuel consumption problems are given in Sections. 2.2-2.5. A small time analysis which describes the behavior of hyperbolic trajectories with state constraints of order one is given in Sections 3.1 and 3.2. We illustrate in Section 3.3, how this analysis reveals the CAS/Mach procedure. This analysis is then applied to the minimum time-to-climb problem in Section 4 combined with numerical methods (indirect multiple shooting and homotopy with monitoring) to classify the optimal structures with respect to the bounds φ max and ψ max . A comparison between minimal cost index solutions (with singular arcs) and the CAS/Mach procedure is then proposed in Section 5. Finally, Section 6 concludes the article.

Preliminary general results and theoretical background

Statement of the optimal control problems

We restrict the dynamics to the vertical motion of the aircraft. A complete description of the motion can be found in [START_REF] Goubinat | Contrôle optimal géométrique et méthodes numériques : application à un problème de montée d'un avion[END_REF]. The aircraft is subjected to four forces, the Drag #-D, the Lift #-L, the Thrust #-T and its own weight #-P . A nonlinear point mass representation is used and we consider that all the forces apply on the center of gravity of the aircraft. We assume that the thrust is collinear to the velocity vector #-V , which means that the angle of attack is neglected here, and that the aircraft evolves in an horizontal constant wind field. The application of the first dynamics principle assuming that the earth is Galilean, provides the four eqs. (1)- [START_REF] Barman | Fixed-range optimum trajectories for short-haul aircraft[END_REF].

During the flight, the lift coefficient C l depends on the variation of the angle of attack. Since this quantity is not taken into account here, it is quite natural to consider the lift coefficient C l as the control variable. However, according to the work presented in [START_REF] Ardema | Singular perturbations in flight mechanics[END_REF][START_REF] Calise | Singular Perturbations and Time Scales in Guidance and Control of Aerospace Systems: A Survey[END_REF] this four-dimensional dynamics contains slow (m) and fast (γ) variables. The time scale separation between the slow and fast variables is handled by a singular perturbation analysis which consists here in computing a zero-order approximation of the solution. The reduction process has two main features. In a first step, we substitute eq. ( 4) by the quasi-steady approximation:

0 = 1 2 ρSv m C l - g 0 v cos γ ⇒ C l = 2mg 0 ρSv 2 cos γ.
In a second time, the lift coefficient is replaced by the previous expression in eqs. ( 1)-( 3) and the flight path angle γ is taken as the new control variable. We also consider that the flight path angle remains small (γ varies from u min := 0 to u max := 0.262 rad) during the climbing so we set cos γ ≈ 1 and sin γ ≈ γ. These considerations lead to the new affine control system:

ẋ(t) = F 0 (x(t)) + u(t) F 1 (x(t)), (5)
where

x := (h, v, m), u := γ, F 0 (x) :=      0 T (h) m - 1 2 ρ(h)Sv 2 m C d,1 -2 mg 2 0 ρ(h)Sv 2 C d,2 -C s (v)T (h)      and F 1 (x) :=    v -g 0 0    .
According to [START_REF] Cots | Time-optimal aircraft trajectories in climbing phase and singular perturbations[END_REF], this reduced dynamics is a sufficient approximation of the initial dynamics considering the state unconstrained minimum time-to-climb problem.

The true air speed of an aircraft is quite difficult to measure during the flight, that is why historically, the concept of CAS was introduced. This speed can be computed using the differential ratio of pressure and can be expressed, through Bernoulli's equations, as a function of the true air speed. The Mach speed, for its part, is defined by the ratio between the true air speed and the speed of the sound. Here, the CAS is denoted by φ and the Mach by ψ and the corresponding expressions are given by:

φ(x) := 2 R Θ 0 κ P (h) P 0 κ v 2 2R Θ(h) + 1 κ -1 + 1 κ -1 , ψ(x) := v γ air R Θ(h) ,
where γ air := 1.4 J.K -1 is the heat capacity of the air and κ := γ air /(1 -γ air ). From these expressions, we define the state constraints c 1 (x) ≤ 0 and c 2 (x) ≤ 0 with:

c 1 (x) := φ(x) -φ max , (6)
c 2 (x) := ψ(x) -ψ max , (7)
where φ max (resp. ψ max ) can be fixed to VMO (resp. MMO) or to smaller specified values.

In this article we are interested in the minimization of different criteria: the final time t f , the fuel consumption ∆ m := m 0 -m f , where m 0 is the given initial mass and m f the final mass, and a convex combination of the final time and the fuel consumption: g α (t f , m f ) := α t f + (1 -α)(m 0 -m f ). The criterion g α , for α ∈ [0 , 1], is called the cost index. The initial state will be fixed to the realistic values x 0 := (h 0 , v 0 , m 0 ) := (3480, 128.6, 69000). We define the set U t of admissible controls, that is the subset of {u | u : [0 , t] → [u min , u max ] measurable} such that the corresponding trajectory solution of eq. ( 5) with x(0) = x 0 is well defined over [0 , t]. The minimum cost index problem can be summarized by:

(P α )                      min (t f ,u)∈D g α (t f , m(t f )), ẋ(t) = F 0 (x(t)) + u(t) F 1 (x(t)), u(t) ∈ U, t ∈ [0 , t f ] a.e., x(0) = x 0 , c 1 (x(t)) ≤ 0, t ∈ [0 , t f ], c 2 (x(t)) ≤ 0, t ∈ [0 , t f ], b(x(t f )) = 0, where U := [u min , u max ], D := (t f , u) t f ≥ 0, u ∈ U t f and b(x) := (h -h f , v -v f )
with h f := 9144 and v f := 191. Thus, the final altitude and true air speed are fixed while the final mass is free. The minimum time-to-climb problem, denoted (P t f ), is defined as (P α | α=1 ) while the minimum fuel consumption problem, denoted (P ∆m ), is defined as (P α | α=0 ).

Notation. In this paper, we follow the presentation of [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF], which exhibits the role of the Lie and Poisson brackets that we define below. Let F 0 , F 1 be two smooth vector fields and c a scalar smooth function on R n . The Lie derivative of c along F 0 denoted F 0 • c is simply the directional derivative of c at x along F 0 (x), given by (

F 0 • c)(x) := c (x) F 0 (x). The Lie bracket between F 0 and F 1 is given by [F 0 , F 1 ] := F 0 • F 1 -F 1 • F 0 , with (F 0 • F 1 )(x) := F 1 (x) F 0 (x)
. Denoting p the adjoint variable and denoting H 0 (x, p) := p, F 0 (x) and H 1 (x, p) := p, F 1 (x) the Hamiltonian lifts of F 0 and F 1 , then the Poisson bracket of H 0 and H 1 is given by

{H 0 , H 1 } := # - H 0 • H 1 , where # - H 0 := ( ∂H0 ∂x , - ∂H0
∂p ) is the Hamiltonian system associated to H 0 . We also use the notation H 01 (resp. F 01 ) to write the bracket {H 0 , H 1 } (resp. [F 0 , F 1 ]) and so forth. Besides, since H 0 and H 1 are two Hamiltonian lifts, we have Chapter 2] for the expressions of F 01 , F 001 and F 101 associated to eq. ( 5).

{H 0 , H 1 } = p, [F 0 , F 1 ] . Remark 2.1. See [24,

Boundary arcs and assumptions

We consider in the following of Section 2 the problem (P α ) but with only one state constraint, since, according to the numerical results from the next sections, the constraints c 1 and c 2 are active at the same time only at isolated times. In other words, we never have c 1 (x(t)) = c 2 (x(t)) = 0 along a time interval of non-empty interior. Due to these considerations, we present the necessary conditions of optimality considering only one scalar state constraint denoted c. The optimal control problem may be written in the form:

min (t f ,u)∈D g α (t f , m(t f , x 0 , u)) subject to b(x(t f , x 0 , u)) = 0, c(x(t, x 0 , u)) ≤ 0, ∀ t ∈ [0 , t f ],
where x(•, x 0 , u) is the solution of eq. ( 5) with the initial condition x(0, x 0 , u) = x 0 .

We call a boundary arc, labeled σ c , an arc defined on an interval I := [a , b] (not reduced to a singleton), such that c(σ c (t)) = 0, for every t ∈ I. The times a and b are called the entry-and exit-time of the boundary arc; a and b are also termed junction times. An arc σ is said to have a contact point with the boundary at t ∈ [0 , t f ] if c(σ( t )) = 0 and c(σ(t)) < 0 for t = t in a neighborhood of t . A subarc σ on which c(σ(t)) < 0 is called an interior arc.

The generic order of the constraint c is the integer m such that

F 1 • c = F 1 • (F 0 • c) = • • • = F 1 • (F m-2 0 • c) = 0 and F 1 • (F m-1 0 • c) = 0.
If the order of a boundary arc σ c is m, then its associated feedback control can be generically computed by differentiating m times the mapping t → c(σ c (t)) and solving with respect to u the linear equation:

c (m) = F m 0 • c + u F 1 • (F m-1 0 • c) = 0.
The boundary feedback control denoted u c is given by

u c := - F m 0 • c F 1 • (F m-1 0 • c) . Let t → σ c (t), t ∈ [t 1 , t 2 ] ⊂ [0 , t f ],
be a boundary arc associated to u c (•). We introduce the assumptions:

(A 1 ) (F 1 • (F m-1 0 • c))(σ c (t)) = 0 for every t ∈ [t 1 , t 2 ]
, with m the order of the constraint.

(A 2 ) u c (t) ∈ [u min , u max ] for t ∈ [t 1 , t 2 ], i.e. the boundary control is admissible. (A 3 ) u c (t) ∈ (u min , u max ) for t ∈ (t 1 , t 2 ), i.e. u c is not saturating on (t 1 , t 2 ).
Remark 2.2. These assumptions are numerically checked a posteriori.

The maximum principle

We recall hereinafter the necessary conditions due to [START_REF] Jacobson | New necessary conditions of optimality for control problems with state-variable inequality constraints[END_REF][START_REF] Maurer | On optimal control problems with bounded state variables and control appearing linearly[END_REF]. Firstly, we define for (x, p, u, η)

∈ H : R n × (R n ) * × R × R the pseudo-Hamiltonian: H(x, p, u, η) := p, F 0 (x) + u F 1 (x) + η c(x)
where η is the Lagrange multiplier of the constraint c and n := 3 is the dimension of the state. Now, let us consider ( tf , ū(•)) ∈ D an optimal solution with associated trajectory x(•). Let us assume that the set of contact and junction times with the boundary, denoted T , is finite. Assume also that the optimal control is piecewise smooth and that along each boundary arc, assumptions (A 1 ) and (A 2 ) are satisfied. Then, we have the following necessary optimality conditions:

(1) There exists a function η(•) ≤ 0, a real number p 0 ≤ 0 and a function of bounded variation p(•) ∈ BV ([0 , tf ], (R n ) * ) such that for t ∈ [0 , tf ] a.e.:

ẋ(t) = ∂H ∂p (x(t), p(t), ū(t), η(t)), ṗ(t) = - ∂H ∂x (x(t), p(t), ū(t), η(t)).
(2) The maximization condition holds for t ∈ [0 , tf ] a.e.:

H(x(t), p(t), ū(t), η(t)) = max u∈U H(x(t), p(t), u, η(t)). (8)
(3) The boundary conditions b(x( tf )) = 0 are satisfied and we have the transversality conditions pm ( tf ) = -p 0 (1-α) and since

t f is free, if ū(•) is continuous at time tf , then, H(x( tf ), p( tf ), ū( tf ), η( tf )) = -p 0 α. ( 4

) The function η(•) is continuous on the interior of the boundary arcs and we have the complementarity condition

η(t) c(x(t)) = 0, ∀t ∈ [0 , tf ]. (5) For any τ ∈ T we have H[τ + ] = H[τ -]
, where [τ ] stands for (x(τ ), p(τ ), ū(τ ), η(τ )), and p(τ

+ ) = p(τ -) -ν τ c (x(τ ))
, where ν τ ≤ 0 is called a jump.

Remark 2.3. Either p 0 = 0 (abnormal case), or p 0 can be set to -1 by homogeneity (normal case). We restrict our study to the normal case. Definition 2.4.

• We call an extremal a quadruple (x(•), p(•), u(•), η(•)) satisfying items 1, 2, 4 and 5. It is called a BC-extremal if it satisfies also item 3.

• Let H 0 (x, p) := p, F 0 (x) and H 1 (x, p) := p, F 1 (x) denote the Hamiltonian lifts of F 0 and F 1 . Then, along any extremal, we call Φ(t) := H 1 (x(t), p(t)) the switching function.

• It follows from ( 8) that along any extremal, we have

u(t) = u min if Φ(t) < 0 and u(t) = u max if Φ(t) > 0.
We say that a trajectory x(•) restricted to a sub-interval I ⊂ [0 , t f ], not reduced to a singleton, is a bang arc if u(•) is constant on I, taking values in {u min , u max }. The trajectory is called bang-bang if it is the concatenation of a finite number of bang arcs. • We say that a trajectory x(•) restricted to a sub-interval I ⊂ [0 , t f ], not reduced to a singleton, is a singular arc if it is an interior arc and if the associated extremal lift satisfies Φ(t) = 0, for every t ∈ I. • We say that an extremal is a bang, singular, boundary or interior extremal if the associated trajectory is respectively a bang, singular, boundary or interior arc.

Definition 2.5. A bang arc such that u(•) ≡ u min (resp. u max ) is labeled σ -(resp. σ + ). A singular arc is labeled σ s while a boundary arc associated to the constraint c is labeled σ c . Besides, we denote by σ 1 σ 2 an arc σ 1 followed by an arc σ 2 .

Remark 2.6. Along a boundary arc, the maximization condition [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF] with assumption (A 3 ) imply Φ = 0 on the interior of the boundary arc. Besides, the adjoint vector may be discontinuous at τ ∈ T .

Parameterization of the singular extremals

Relaxing the control bounds, singular trajectories are parameterized by the constrained Hamiltonian system:

ẋ = ∂ p H, ṗ = -∂ x H, 0 = ∂ u H = H 1 , with H 1 (x, p) = p, F 1 (x) the Hamiltonian lift of F 1 .
The constraint H 1 = 0 has to be differentiated at least twice along a singular extremal to compute the control. This gives:

H 1 = H 01 = H 001 + u H 101 = 0,
along any singular extremal. A singular extremal along which H 101 = 0 is called of minimal order and the corresponding singular control is given by:

u s (z) := - H 001 (z) H 101 (z) ,
with z := (x, p). In this case, we have the following additional necessary condition of optimality deduced from the high-order maximum principle [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF]. If the singular control is not saturating along the singular extremal, then the generalized Legendre-Clebsch condition must hold along the singular extremal, that is:

∂ ∂u ∂ 2 ∂t 2 ∂H ∂u = H 101 ≥ 0. (9)
(z) := H 0 (z) + u s (z) H 1 (z).
Then, the singular extremals of minimal order are the solutions of

ż(t) = # - H s (z(t)), starting from the set {z | H 1 (z) = H 01 (z) = 0}.
For three-dimensional systems, the singular control may be written in feedback form only with respect to the state x. Indeed,setting D ξ (x) := det(F 1 (x), F 01 (x), F ξ (x)), the singular control is given in feedback form by

u s (x) = - D 001 (x) D 101 (x) ,
whenever D 101 (x) = 0 and p = 0, since, along a singular extremal

p, F 1 (x) = p, F 01 (x) = p, F 001 (x) + u F 101 (x) = 0.

Parameterization of the boundary extremals

We may find in [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF][START_REF] Maurer | On optimal control problems with bounded state variables and control appearing linearly[END_REF] the determination of the multiplier η and the jump ν τ together with the analysis of the junction conditions, which is based on the concept of order and related to the classification of extremals. We give next some results only for the case m = 1 since the constraints c 1 and c 2 are of order 1, according to numerical experiments that we have realized as follows.

First, the differentiation of c 1 (x(t)) = φ(x(t)) -φ max with respect to the time t leads to

ċ1 (x(t)) = ∂φ ∂h (x(t)) ḣ(t) + ∂φ ∂v (x(t)) v(t) and ∂ ċ1 ∂u (x) = v ∂φ ∂h (x) -g 0 ∂φ ∂v (x) = (F 1 • c 1 )(x).
Let us define V * := 140 m.s -1 as the minimal authorized speed in our case study 1 . Then, we check thanks to the software Sage [START_REF] Casamayou | Calcul mathématique avec Sage[END_REF], that (

F 1 • c 1 )(x) does not vanish for any (h, v) ∈ [3400 , 11000] × [V * , V M O]
, so, for our case study, c 1 is a state constraint of order 1. The same reasoning holds for c 2 .

Remark 2.9. We refer to [START_REF] Goubinat | Contrôle optimal géométrique et méthodes numériques : application à un problème de montée d'un avion[END_REF]Chapter 3] for the exact and quite long expressions of F 0 • c, F 1 • c and H 01 for both constraints c 1 and c 2 .

For a first-order constraint denoted c, assuming (A 1 ) and (A 3 ), we have the following result from [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF].

Proposition 2.10. Let m = 1. Then:

(1) along the boundary, the control and the multiplier are given by

u c (x) = - (F 0 • c)(x) (F 1 • c)(x) and η c (z) = H 01 (z) (F 1 • c)(x)
, z = (x, p).

(2) if the control is discontinuous at a contact or a junction time τ between a bang arc and the boundary then the jump ν τ = 0. (3) we have

ν τ = Φ(τ -) (F 1 • c)(x(τ ))
at an entry point and

ν τ = - Φ(τ + ) (F 1 • c)(x(τ ))
at an exit point.

Likewise the singular case, in the state constrained case, we have the following result excerpted from [START_REF] Cots | Geometric and numerical methods for a state constrained minimum time control problem of an electric vehicle[END_REF]Prop. 4.5] which is useful to define the numerical shooting method, see Section 4.1.

Proposition 2.11. Let c(x) ≤ 0 be a smooth scalar state constraint of order 1 and assume the open subset

Ω := {x | (F 1 • c)(x) = 0} is not empty. On Ω we define the Hamiltonian H c (z) := H 0 (z) + u c (x) H 1 (z) + η c (z) c(x), z = (x, p),
where u c and η c are defined in Prop. 2.10. Then, the boundary extremals of order 1 contained in H 1 = 0 are the solutions of the system ż

(t) = # - H c (z(t)), starting from {z = (x, p) | c(x) = H 1 (z) = 0}.

Small time analysis of the minimum time-to-climb problem

The minimum time-to-climb problem with no state constraints is analyzed in [START_REF] Cots | On the minimum time optimal control problem of an aircraft in its climbing phase[END_REF][START_REF] Cots | Time-optimal aircraft trajectories in climbing phase and singular perturbations[END_REF]. In [START_REF] Cots | On the minimum time optimal control problem of an aircraft in its climbing phase[END_REF], the influence of the initial mass m 0 and the final true air speed v f is studied and it is shown in particular that the trajectories are of the form σ ± σ s σ ± , where σ s is an hyperbolic singular extremal. We thus restrict the small time analysis, presented in Section 3.2, to hyperbolic singular extremals, that we introduce in Section 3.1.

Generic classification of extremals near the switching surface

Along a singular extremal we have ∂ u H ≡ 0. In our particular case of a single-input affine control system this condition becomes Φ(t) = H 1 (z(t)) = 0. Let us define the switching surface

Σ := {z := (x, p) ∈ R n × (R n ) * | H 1 (z) = 0} ,
and the set

Σ s := {z := (x, p) ∈ R n × (R n ) * | H 1 (z) = H 01 (z)
= 0} containing all the singular extremals. Let us introduce the notation Φ + (resp. Φ -) if the control along a bang extremal is u max (resp. u min ). The first and second derivatives of Φ ± are given by:

Φ± (t) = H 01 (z(t)), Φ+ (t) = H 001 (z(t)) + u max H 101 (z(t)), Φ-(t) = H 001 (z(t)) + u min H 101 (z(t)).
A crucial point to analyze the minimum time-to-climb problem, is then to apply the results from [START_REF] Kupka | Geometric theory of extremals in optimal control problems. i. the fold and maxwell case[END_REF] (see also [START_REF] Bonnard | Singular Trajectories and their Role in Control Theory[END_REF]) to classify the extremal curves near the switching surface Σ. We have the following:

(1) Ordinary switching time. It is a time t such that two bang arcs switch with Φ( t) = 0 and Φ( t) = H 01 (z( t)) = 0. According to the maximum principle, near Σ, the extremal is of the form

σ -σ + if Φ( t) > 0 and σ + σ -if Φ( t) < 0.
(2) Fold point. It is a time where a bang arc has a contact of order 2 with Σ. We have three cases (if Φ± = 0) depending on Φ± at the switching time:

• Hyperbolic case: Φ+ > 0 and Φ-< 0. A connection with a singular extremal is possible at Σ s and locally each extremal is of the form σ ± σ s σ ± (by convention each arc of the sequence can be empty). • Parabolic case: Φ+ Φ-> 0. The singular extremal at the switching point is not admissible and every extremal curve is locally bang-bang with at most two switchings, i.e. σ + σ -σ + or σ -σ + σ -. • Elliptic case: Φ+ < 0 and Φ-> 0. A connection with a singular arc is not possible and locally each extremal is bang-bang but with no uniform bound on the number of switchings. When dealing with three-dimensional systems we can give conditions depending only on the state x to classify the three cases around fold points. Assuming D 0 (x) = 0, then the family (F 0 (x), F 1 (x), F 01 (x)) forms a basis of R 3 and there exists (α 0 , α 1 , α 01 ) ∈ R 3 and (β 0 , β 1 , β 01 ) ∈ R 3 such that

F 001 (x) + u min F 101 (x) = α 0 F 0 (x) + α 1 F 1 (x) + α 01 F 01 (x), F 001 (x) + u max F 101 (x) = β 0 F 0 (x) + β 1 F 1 (x) + β 01 F 01 (x).
By linearity of the determinant, we have

D 001 (x) + u min D 101 (x) = α 0 D 0 (x), D 001 (x) + u max D 101 (x) = β 0 D 0 (x),
and since D 0 (x) = 0 then we can compute α 0 and β 0 along the trajectory. Besides, along any singular extremal, we have:

Φ-(t) = p(t) , F 001 (x(t)) + u min F 101 (x(t)) = α 0 (t) H 0 (z(t)), Φ+ (t) = p(t) , F 001 (x(t)) + u max F 101 (x(t)) = β 0 (t) H 0 (z(t)),
with H 0 (z(t)) > 0 in the normal case with the convention of the maximum principle. Denoting t the time when the extremal has a contact of order two with the switching surface Σ, then the singular extremal is hyperbolic if α 0 ( t ) < 0 and β 0 ( t ) > 0, is elliptic if α 0 ( t ) > 0 and β 0 ( t ) < 0, and parabolic if α 0 ( t ) β 0 ( t ) > 0.

Remark 3.1. On [24, Fig. 2.11, p. 59], the numerical analysis is illustrated and it shows a posteriori that α 0 (t) > 0 and β 0 (t) < 0 along the singular arc. Thus, we are in the hyperbolic case and the assumption that D 0 does not vanished is fulfilled.

Analysis of the hyperbolic case with a one-order state constraint

The time-optimal trajectories are concatenations of bang, singular and boundary arcs. One of the main difficulties is to determine the number and the sequence of these arcs. We can see these concatenations as a sequence of patterns (or patches), each pattern being made with only few concatenations. These patterns depend on the classification presented in Section 3.1 but also on the order of the state constraint. Since we are dealing with order one state constraints and hyperbolic singular extremals, the number of possible patterns is drastically reduced. Besides, a (long) time-optimal trajectory may be seen as a concatenation of small timeoptimal trajectories since it is necessarily optimal on every sub-interval of time. Finally, the pattern we are looking for are thus given by small time-optimal trajectories that we analyze in this section.

Remark 3.2. In this section, the main result is given by Prop. 3.4 which explains the structures of the timeoptimal trajectories. In the following, we make a brief and quite technical summary of some results excerpted from [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF][START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF], in order to introduce the notations of Prop. 3.4, and which can be skipped in a first read.

Remark 3.3. In [START_REF] Cots | On the minimum time optimal control problem of an aircraft in its climbing phase[END_REF], it is shown that the state-unconstrained time-optimal trajectories are of the form σ ± σ s σ ± , i.e. the optimal sequences are made with only one single pattern.

The aircraft dynamics is a multi-scale system and since the time constant of the mass is a hundred, resp. a thousand, times greater than the time constant of the altitude, resp. the speed, we can neglect the evolution of the mass when considering short time intervals and thus, in this case, we can reduce our three dimensional dynamics into a planar dynamics. We present hereinafter some results from [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF][START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF] about small time analysis considering only a planar single-input control system, in the hyperbolic case and with an order one scalar state constraint c(q(t)) ≤ 0, of the form:

q(t) = F 0 (q(t)) + u(t) F 1 (q(t)),
where |u(t)| ≤ 1 and q(t) := (x(t), y(t)) ∈ R 2 . Let us take q 0 such that c(q 0 ) = 0 and let us identify q 0 to the origin. Let us assume that F 0 (q 0 ) and F 1 (q 0 ) are linearly independent, and that the constraint c is of order one, i.e. (F 1 • c)(q 0 ) = 0. Then replacing, if necessary, u by -u, we can find a local diffeomorphism preserving q 0 = 0 and transforming the constrained system into the following generic model:

ẋ(t) = 1 + y(t) a(q(t)), ẏ(t) = b(q(t)) + u(t), y(t) ≤ 0.
To describe the hyperbolic case, we consider also that det(F 1 (q 0 ), F 01 (q 0 )) = a(q 0 ) = 0 and we assume that the set S := q ∈ R 2 det(F 1 (q), F 01 (q)) = 0 containing the singular trajectories is a simple curve (this assumption is not restrictive since our approach is local) that we approximate by a straight line in our small time model. The equations of the system become:

ẋ(t) = 1 + y(t) (a y(t) + b x(t)), ẏ(t) = c + u(t), y(t) ≤ 0,
where S is identified to q ∈ R 2 2 a y + b x = 0 and we assume that a = 0 and b = 0. Note that the generalized Legendre-Clebsch condition is related to the sign of a. The singular control at q 0 is given by u s (q 0 ) = -c -b/2a

q 0 q S σ - σ + σ c • • b > 0 Figure 1.
Small time representation of different trajectories in the neighborhood of a point q 0 belonging to the singular set S and such that c(q 0 ) = 0, in the case b > 0. Two trajectories joining q 0 and q are represented: one boundary arc σ c and a bang-bang trajectory σ -σ + . Using the clock form and the Stokes theorem, we can show that the trajectory σ c is time-optimal.

and the boundary control is simply u c (q) = -c. Taking the constraint on the control into account, we see that the condition of admissibility implies |c + b/2a| ≤ 1 and |c| ≤ 1. Finally, the hyperbolic case corresponds to:

a < 0, |c + b/2a| < 1, |c| < 1 and b = 0,
that is the strict generalized Legendre-Clebsch condition is satisfied, the singular and boundary controls are strictly admissible, and the singular set S is not identified with the boundary set y = 0. For the state unconstrained problem, the singular arc is optimal, each optimal trajectory has at most two switchings and the local optimal synthesis is of the form σ ± σ s σ ± . For the state constrained problem, the socalled clock form and Stokes theorem (see [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF][START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF] for details) can be used to conclude on the optimatility of the boundary arc when x ≥ 0 and its non-optimality when x < 0, for the case b > 0. For b < 0, the boundary arc is optimal when x ≤ 0 and non-optimal when x > 0. See Fig. 1 for a representation of the bang, singular and boundary arcs in the hyperbolic case with b > 0, and for a comparison between boundary and bang-bang arcs.

When b > 0, the small time optimal synthesis joining two points q 1 and q 2 on the boundary y = 0, with q 1 and q 2 on each side of S and q 1 on the left, is then of the form σ -σ s σ c . In this case, each optimal curve in a neighborhood of q 0 has at most three switchings and the local optimal synthesis is of the form σ ± σ s σ c σ ± . For b < 0, the local optimal synthesis is of the form σ ± σ c σ s σ ± . We summarize this analysis in the following proposition excerpted from [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF]. Proposition 3.4. Under our assumptions, in the hyperbolic case each small time optimal trajectory has at most three switchings. Moreover,

(1) For b < 0, a boundary arc is optimal iff x ≤ 0 and each optimal arc has the form σ ± σ c σ s σ ± .

(2) For b > 0, a boundary arc is optimal iff x ≥ 0 and each optimal arc has the form σ ± σ s σ c σ ± .

Remark 3.5. Prop. 3.4 gives two different patterns σ ± σ s σ c σ ± and σ ± σ c σ s σ ± , depending on the local model, for small time optimal trajectories around the boundary. These two patterns generalize the pattern σ ± σ s σ ± , presented in Section 3.1, in the hyperbolic state unconstrained case.

A first insight on the standard CAS/Mach procedure

In this section, we explain on an example why the CAS/Mach procedure makes sense. Let us take from [18, Fig. 1] the state unconstrained trajectory of the form σ -σ s σ + , for the minimum time-to-climb problem, with data given in Sections 1 and 2.1. Note that in [START_REF] Cots | On the minimum time optimal control problem of an aircraft in its climbing phase[END_REF], u min = -0.262 instead of 0 but this has a small influence and it does not change the structure σ -σ s σ + . From this excerpted trajectory, denoted x(•), we compute and display on Fig. 2, the state constraints c 1 (x(t)) = φ(x(t)) -φ max and c 2 (x(t)) = ψ(x(t)) -ψ max , for φ max = 160 and ψ max = 0.7. These realistic values of φ max and ψ max are chosen to emphasize the following comment: the boundary c 1 = 0 is reached during the first bang arc, before the singular arc, while the boundary c 2 = 0 is reached during the singular arc, after reaching c 1 = 0 and before the last bang arc. From Fig. 2, one can also notice that φ(x(•)) and ψ(x(•)) reach their maxima respectively at the switching times t 1 and t 2 . Let us consider now other values for φ max and ψ max . Assume we fix φ max = φ(x(t 1 )) -ε, ε > 0 small, and ψ max > ψ(x(t 2 )). Considering these values, the state unconstrained trajectory is not admissible anymore for the state constrained case since in a small neighborhood of t 1 , we have c 1 (x(t)) > 0. The trajectory is thus slightly modified and the updated optimal sequence becomes σ -σ c1 σ s σ + corresponding to the first case of Prop. 3.4. On the contrary, if φ max > φ(x(t 1 )) and ψ max = ψ(x(t 2 )) -ε, ε > 0 small, then the updated optimal trajectory is of the form σ -σ s σ c2 σ + corresponding to the second case of Prop. 3.4. Finally, if we set φ max = φ(x(t 1 )) -ε 1 and ψ max = ψ(x(t 2 )) -ε 2 , both ε 1 > 0 and ε 2 > 0 sufficiently small, then the state constrained trajectory is of the form σ -σ c1 σ s σ c2 σ + , that is the boundary arc σ c1 appears before the boundary arc σ c2 and this sequence is made of two patterns. In other words, the trajectory contains one arc at constant CAS which is before the arc at constant Mach, and this is consistent with the actual CAS/Mach procedure. Let us recall that this procedure splits the climbing phase in two parts, the first one is an arc at constant CAS and the second one is an arc at constant Mach number. So, comparing to the CAS/Mach procedure there is an additional singular arc joining the boundary arcs, but, according to the results presented in Section 4, for smaller values of φ max and ψ max , the singular arc vanishes and we get the CAS/Mach procedure. Remark 3.6. In Section 4, we analyze the influence of φ max and ψ max on the structure of the trajectories by deforming the BC-extremal associated to the state unconstrained trajectory, taken here as example. Just note that for all the numerical case studies presented in this article, whatever the values of φ max and ψ max , the boundary arc σ c1 is always before σ c2 .

The minimum time-to-climb problem: numerical methods and results

Introduction to multiple shooting methods

We present in this section the indirect multiple shooting method [START_REF] Bulirsch | Introduction to numerical analysis[END_REF] that we use to solve the necessary conditions of optimality given by the maximum principle presented in Section 2. We describe the method on only one example and we refer to [START_REF] Goubinat | Contrôle optimal géométrique et méthodes numériques : application à un problème de montée d'un avion[END_REF] for more details. Let us consider the same example as the one presented in the previous section and let us fix φ max > φ(x(t 1 )) and ψ max = 0.7. In this case, the resulting trajectory is of the form σ -σ s σ c2 σ + .

The unknowns of the shooting method are the junction times t 1 , t 2 and t 3 , the final time t f , the jumps ν t2 and ν t3 at the junction times t 2 and t 3 and the initial costate denoted p 0 . With p 0 , t 1 , t 2 , t 3 , ν t2 and ν t3 , since we know the structure of the trajectory, we can retrieve the states-costates at the junction times t 1 , t 2 and t 3 simply by integration, starting from x(0) = x 0 and applying the controls u min , u s , u c2 and u max , respectively on [0 , t 1 ], [t 1 , t 2 ], [t 2 , t 3 ] and [t 3 , t f ]. However, to improve numerical stability, we add the states-costates z 1 , z 2 and z 3 at the junction times to the unknowns of the shooting method. The first thing to notice now is that the jumps are zero. Indeed, Φ(t - 2 ) = 0 and so ν t2 = 0 (from Prop. 2.10) since the switching function Φ(•) = 0 along the singular arc σ s . Now, combining Prop. 2.10 and item 5 of the maximum principle, we have:

ν t3 = - Φ(t + 3 ) (F 1 • c 2 )(x(t 3 )) ≤ 0.
But, Φ(t + 3 ) ≥ 0 since the last bang arc is positive. So, assuming (F 1 • c 2 )(x(t 3 )) < 0, then ν t3 = 0. Note that we can check a posteriori that this condition is verified. Finally, grouping all together, we write

y := (p 0 , t 1 , t 2 , t 3 , t f , z 1 , z 2 , z 3 ) ∈ R 25
the unknown of the shooting method.

Remark 4.1. For the jump ν t3 to be not zero, it is necessary that u c2 (x(t 3 )) = u max , which may happen in only few particular cases that we do not encounter in the numerical experiments all through this article. Now, we need to describe the shooting equations and define the shooting function. First of all, we define the following Hamiltonians:

H + (z) := H 0 (z) + u max H 1 (z), H -(z) := H 0 (z) + u min H 1 (z), H s (z) := H 0 (z) + u s (x) H 1 (z), H c2 (z) := H 0 (z) + u c2 (x) H 1 (z) + η c2 (z) c 2 (x)
, where z = (x, p) and where u s , u c2 and η c2 are given in Section 2. Note that we can replace u s (x) = -D 001 (x)/D 101 (x) by u s (z) = -H 001 (z)/H 101 (z). Then, we define the exponential mapping e tϕ (q 0 ) as the solution at time t of the Cauchy problem q(s) = ϕ(q(s)), q(0) = q 0 , where ϕ is any dynamical system and where the state here is denoted q. We define also the canonical projections π x (z) = x and π pm (z) = p m , recalling that z = (x, p), x = (h, v, m) and writing p := (p h , p v , p m ). The multiple shooting function, denoted S c2 , associated to the structure σ -σ s σ c2 σ + is then defined by:

S c2 (y) :=               H 1 (e t1 # - H-(x 0 , p 0 )) H 01 (e t1 # - H-(x 0 , p 0 )) c 2 (π x (e (t2-t1) # - Hs (z 1 )))
Ψ(e (t f -t3) # - H+ (z 3 ))

z 1 -e t1 # - H-(x 0 , p 0 ) z 2 -e (t2-t1) # - Hs (z 1 ) z 3 -e (t3-t2) # - Hc 2 (z 2 )               ∈ R 25 ,
where Ψ(z) := b(π x (z)), π pm (z), H + (z) + p 0 , p 0 = -1, and the multiple shooting method consists in finding a zero of the multiple shooting function S c2 , i.e. in solving S c2 (y) = 0. Remark 4.2. We can replace the first three equations by

H 1 (z 1 ) = H 01 (z 1 ) = c 2 (π x (z 2 )) = 0. Let ȳ := (p 0 , t 1 , t 2 , t 3 , t f , z1 , z2 , z3 ) be a zero of the shooting function S c2 . Then, to ȳ is associated a unique BC-extremal (x(•), p(•), ū(•), η(•)) if 0 ≤ t 1 ≤ t 2 ≤ t 3 ≤ t f and if H 1 (z(t)) < 0 for almost every t ∈ [0 , t 1 ] and H 1 (z(t)) > 0 for almost every t ∈ [ t 3 , t f ]
. These two conditions are not included in the shooting equations but they are easily checked a posteriori and should be satisfied if the structure is the right one. Let us explain now why these shooting equations define a BC-extremal of the form σ -σ s σ c2 σ + . The first two equations in S c2 (y) = 0 impose the extremal to be singular on [ t 1 , t 2 ] according to Prop. 2.7, the third equation impose a boundary arc on [ t 2 , t 3 ] according to Prop. 2.11, noticing that H 1 (e ( t 2 -t 1)

# -

Hs (z 1 )) = 0, while the fourth equation Ψ = 0 contains the limit and transversality conditions. The last equations are the so-called matching conditions.

We use the HamPath software [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF][START_REF] Cots | Geometric and numerical methods for a state constrained minimum time control problem of an electric vehicle[END_REF], to obtain a (numerical) zero, denoted ȳ, of the shooting function with high accuracy: S c2 (ȳ) ≈ 1.511 × 10 -11 . The switching times are t 1 ≈ 88.61s, t 2 ≈ 455.7s and t 3 ≈ 651.46s. The final time is t f ≈ 661.37s and the initial adjoint vector is p0 ≈ (3.180 × 10 -2 , 6.594 × 10 -1 , -2.300 × 10 -1 ). Note that the HamPath software is based on a Newton-like algorithm which is initialized thanks to the state unconstrained trajectory. Fig. 3 depicts the control and the constraint c 2 along the resulting state constrained trajectory.

Introduction to differential path following methods with monitoring

The shooting method is used to solve a single optimal control problem. To solve a one-parameter family of optimal control problems, e.g. for different values of ψ max , we use differential path following techniques [START_REF] Allgower | Introduction to numerical continuation methods[END_REF][START_REF] Gergaud | Sur la résolution numérique de problèmes de contrôle optimal à solution bang-bang via les méthodes homotopiques[END_REF] with arc length parameterization (also called homotopy method). Let h : R N × R → R N , h(y, λ), denote an homotopic function with λ the homotopic parameter. For example, one can consider the homotopic function defined by S c2 , with y := (p 0 , t 1 , t 2 , t 3 , t f , z 1 , z 2 , z 3 ), N = 25 and λ := ψ max . That is to say, in this example h(y, λ) = S c2 (y) with λ = ψ max which appears in the constraint function c 2 . We are interested in solving h = 0 for λ ∈ [λ 0 , λ 1 ], (λ 0 , λ 1 ) ∈ R 2 given.

The classical difficulties about homotopic methods consist in assuring that a curve in h -1 (0) exists, is sufficiently smooth and will intersect the fixed target λ 1 in finite length starting from λ = λ 0 . As a first result we have the following. Suppose h is continuously differentiable and that we know y 0 such that h(y 0 , λ 0 ) = 0 and assume rank ∂h ∂y (y 0 , λ 0 ) = N.

Suppose also that 0 is a regular value of h. Then, a continuously differentiable curve crossing (y 0 , λ 0 ) and transverse to λ = λ 0 exists and is either diffeomorphic to a circle or the real line. The curves in h -1 (0) are disjoints, and we call each branch of h -1 (0) a path of zeros. To compute numerically these paths of zeros, we use the homotopic method from the HamPath software, which is based on a Predictor-Corrector algorithm with a high order and variable step size Runge-Kutta scheme for the prediction and with a classical simplified Newton method for the correction. We present in this section how we combine homotopy and monitoring (that we define hereinafter) to obtain, for a fixed climbing scenario, a cartography of possible time-optimal structures with respect to the values of the bounds of the state constraints c 1 and c 2 , i.e. on the values of φ max and ψ max . We do not give all the details and we refer to [START_REF] Cots | Geometric and numerical methods for a state constrained minimum time control problem of an electric vehicle[END_REF] for a more detailed description of the methodology to obtain such a cartography. We start from the state unconstrained trajectory, denoted x(•), presented in Section 3.3 and of the form σ -σ s σ + . The maximal values of φ and ψ along the trajectory are reached respectively at the switching times t 1 and t 2 . One can notice that φ(x(t 1 )) < V M O and ψ(x(t 2 )) < M M O, so for this trajectory, the state constraints c 1 and c 2 are not violated if we fix φ max = V M O and ψ max = M M O. The idea is then to deform by homotopy, making φ max or ψ max varying, the BC-extremal associated to this trajectory and then to detect by monitoring when a change in the structure occurs. The monitoring consists in checking some conditions after each step of prediction and correction. Technically, we only have to code the conditions and then the HamPath software automatically check them after each prediction-correction step and stop the homotopy process if at least one condition is violated. Here are the three different monitorings we use:

M1 check if c 1 is violated at the entry point of σ s . If yes, then we add a σ c1 arc; M2 check if c 2 is violated at the exit point of σ s . If yes, then we add a σ c2 arc; M3 check if σ s has a positive length. If not, then the arc is removed. When a change is detected, we update the homotopic function accordingly to the new structure and repeat the process. We thus have to limit the range of values for φ max and ψ max . We choose to build the cartography for

(φ max , ψ max ) ∈ [φ 0 , V M O] × [ψ 0 , M M O]
, where V M O ≈ 180, M M O ≈ 0.82, φ 0 := max(φ(x(0)), φ(x(t f ))) ≈ 107 and ψ 0 := max(ψ(x(0)), ψ(x(t f ))) ≈ 0.63, with t f the final time. With these chosen minimal values ψ 0 and φ 0 , whatever φ max and ψ max , the trajectory will always start with an arc σ -and end with an arc σ + because the state constraints depend only on h and v which are fixed at the initial and final times to some given values which do not depend on φ max and ψ max . This explains why we say that we study the deformation of the BC-extremal for a fixed climbing scenario.

Classification of BC-extremals

Let simply recall that we start from the BC-extremal associated to the state unconstrained trajectory of the form σ -σ s σ + , denoted x(•), with t 1 and t 2 the switching times and t f the final time.

Step 1. Let φ c1 := φ(x(t 1 )), ψ c2 := ψ(x(t 2 )) and let us label δ c1 , δ c2 and δ s , a touch point respectively with the state constraint c 1 = 0, c 2 = 0 and the singular set Σ s . Then, we have the following straightforward first result:

φ max ψ max structure φ c1 ψ c2 σ -δ c1 σ s δ c2 σ + φ c1 (ψ c2 , M M O] σ -δ c1 σ s σ + (φ c1 , V M O] ψ c2 σ -σ s δ c2 σ + (φ c1 , V M O] (ψ c2 , M M O] σ -σ s σ +
Step 2. Let fix first ψ max = M M O. For φ max = φ c1 -ε, ε > 0 small enough, the BC-extremal has a structure of the form σ -σ c1 σ s σ + . Then, we perform an homotopy on λ = φ max from λ 0 = φ c1 to λ 1 = φ 0 which stops because of monitoring M3 around λ ≈ 129.8 =: φ s . Hence, at λ = φ s the structure is σ -σ c1 δ s σ + . The second homotopy from λ 0 = φ s to λ 1 = φ 0 with a structure σ -σ c1 σ + is not stopped by the monitoring, thus, no change in the structure is detected anymore. Let fix now (φ max , ψ max ) = (φ s , M M O). The structure is of the form σ -σ c1 δ s σ + . Then, we define ψ s := max ψ(•) ≈ 0.734 along the corresponding trajectory, and we have the following:

φ max ψ max structure [φ 0 , φ s ) MMO σ -σ c1 σ + φ s MMO σ -σ c1 δ s σ + (φ s , φ c1 ) MMO σ -σ c1 σ s σ + φ s ψ s σ -σ c1 δ s δ c2 σ + φ s (ψ s , M M O] σ -σ c1 δ s σ +
Step 3. The same methodology is used each time the homotopy is stopped by a monitoring. The new structure is then guessed thanks to the results from Section 3. The final classification is given on the left sub-graph of Fig. 4.

Discussion on the classification

The left sub-graph of Fig. 4 presents the resulting cartography of the BC-extremal structures for (φ max , ψ max ) ∈

[φ 0 , V M O] × [ψ 0 , M M O],
for the minimum time-to-climb problem. The right sub-graph of Fig. 4 gives the time-to-climb related to the trajectories from left sub-graph. Of course, the state unconstrained trajectory of the form σ -σ s σ + minimizes the time-to-climb and the minimal value is t f,min := 658 s. Note that as a remarkable fact, we retrieve the CAS/Mach procedure for some values of φ max and ψ max . Up to this point, we can say that the CAS/Mach procedure is not the optimal structure for the minimum time-to-climb problem with the Maximum Operating Speed (VMO) and the Maximum Operating Mach (MMO) speeds as constraints. However, when decreasing these constraints, the CAS/Mach procedure may be optimal since the classification from Fig. 4 gives in this case BC-extremals of the form σ -σ c1 σ c2 σ + . Remark 4.3. In the previous paragraph, we say "may be optimal" since we only check necessary conditions of optimality given by the maximum principle, that is we only compute BC-extremals. The question about local or global optimality is difficult in this context and go beyond the scope of this article. Remark 4.4. According to Section 5.3, for the solutions from Fig. 4, the fuel consumption is minimal for (φ max , ψ max ) ≈ (128.9, 0.661) where the structure follows the CAS/Mach procedure.

5.

The CAS/Mach versus the Singular Arc procedure This procedure usually splits the climbing phase in two parts, the first one is an arc σ c1 at constant CAS and the second one is an arc σ c2 at constant Mach number. In our case, we have seen that we need a first negative bang arc σ -and we need a final positive bang arc σ + to reach the limit conditions (b(x(t f )) = 0) which define the cruise phase. We consider the cost index g α (see Section 2.1) as the objective function and we want to optimize the CAS and Mach numbers, that is φ max and ψ max , followed respectively along the arcs σ c1 and σ c2 . The corresponding finite dimensional optimization problem may be summarized this way:

(CM α ) min X∈R 15 f CM α (X) subject to g CM α (X) = 0,
where the decision variable is X = (t 1 , t 2 , t 3 , t f , x 1 , x 2 , x 3 , φ max , ψ max ), where the cost function is given by

f CM α (X) := α t f + (1 -α) m 0 -π m e (t f -t3)F+ (x 3 ) = g α t f , π m e (t f -t3)F+ (x 3 ) ,
with π m (x) = m, recalling that x = (h, v, m), and where m 0 is the initial mass, which is fixed. The equality constraints are defined by:

g CM α (X) :=          
b(e (t f -t3)F+ (x 3 ))

x 1 -e t1F-(x 0 ) x 2 -e (t2-t1)Fc 1 (x 1 )

x 3 -e (t3-t2)Fc 2 (x 2 ) φ max -φ(x 1 ) ψ max -ψ(x 2 )           ∈ R 13 ,
recalling that x 0 = (h 0 , v 0 , m 0 ) is fixed (see Section 2.1 for the numerical values).

Remark 5.1. One can check a posteriori that 0 ≤ t 1 ≤ t 2 ≤ t 3 ≤ t f or include these inequality constraints into the optimization problem formulation.

The Singular Arc (SA) procedure

The Singular Arc (SA) procedure, which is, according to the next section, better than the CM procedure, is defined by trajectories of the form σ -σ s σ + . The cost function is still the cost index and the equality constraints are simply given by the limit and matching conditions. This gives the following NLP problem:

(SA α ) min X∈R 9 f SA α (X) subject to g SA α (X) = 0,
where the decision variable is here X = (t 1 , t 2 , t f , x 1 , x 2 ), where the cost function is given by and where the equality constraints are defined by:

f SA α (X) := α t f + (1 -α) m 0 -π m e (t f -t2)F+ (x 2 ) = g α t f , π m e (t f -t2)F+ (x 2 ) ,
g SA α (X) :=    b(e (t f -t2)F+ (x 2 )) x 1 -e t1F-(x 0 ) x 2 -e (t2-t1)Fs (x 1 )    ∈ R 8 .
Remark 5.2. One can also check a posteriori that 0 ≤ t 1 ≤ t 2 ≤ t f or include these inequality constraints into the optimization problem formulation.

Remark 5.3. Let us denote by (SA t f ) the problem (SA α ) with α = 1 and (CM t f ) the problem (CM α ) with α = 1. According to the results from Section 4 (see Fig. 4 in particular), we have that (SA t f ) is better than (CM t f ).

The minimal fuel consumption problem

We have already compared (SA α ) and (CM α ) when α = 1. Before comparing them for all possible values of α, we present a brief analysis when α = 0, that is for the minimum fuel consumption problem. We denote by (SA ∆m ) the problem (SA α ) with α = 0 and (CM ∆m ) the problem (CM α ) with α = 0.

Fuel consumptions for (SA t f ) versus (CM ∆m ). For the minimum time-to-climb problem, the classification of the different structures is given on the left sub-graph of Fig. 4. Associated to these BC-extremals, we have the final time t f presented on the right sub-graph of Fig. 4. One can notice that the minimal value (t f,min = 658s) is of course for the structure following the SA procedure, and the associated fuel consumption is ∆ m ≈ 882kg. For the BC-extremals from Fig. 4, one can also compute the fuel consumption, which is given on Fig. 5. Even if it is not clear from the figure, the minimal value of fuel consumption is obtained when (φ max , ψ max ) ≈ (128.9, 0.661), that is for a structure following the CM procedure. To obtain these optimal values we solve the problem (CM ∆m ) with the fmincon function from the Matlab software and we use the ode45 function to compute the numerical integrations. The numerical solution we obtain satisfies the first-order optimality conditions with a tolerance of 10 -7 . The switching times are t 1 ≈ 32s, t 2 ≈ 450s and t 3 ≈ 676s while we have a fuel consumption ∆ m ≈ 863kg, a final time t f ≈ 677s and the CAS and Mach optimal values are (φ max , ψ max ) ≈ (128.9, 0.661). Hence, there is a gap of 19s on the final time and 19kg on the fuel consumption between the solution of (SA t f ) and (CM ∆m ). That is, the decrease of the bounds φ max and ψ max of the state constraints c 1 and c 2 allows to save up 19 kg of fuel with a 19 seconds longer flight.

Optimal fuel consumption. In the previous paragraph, we have compared (SA t f ) with (CM ∆m ). As expected, solving (SA t f ) gives a better final time and solving (CM ∆m ) gives a better fuel consumption. But we can do better. Let us come back to the original problem of minimizing the fuel consumption, that is to the optimal control problem (P ∆m ) defined in Section 2.1. To obtain the best possible fuel consumption, we set φ max = V M O and ψ max = M M O. In this case, thanks to the direct collocation method from the Bocop software [START_REF] Bonnans | Bocop -a collection of examples[END_REF], we determine that the structure σ -σ s σ + is still relevant, that is the solution we obtain follows the SA procedure. This means that the problem (SA ∆m ) gives a better solution than the problem (CM ∆m ). To get a BC-extremal solution of (P ∆m ) with φ max = V M O and ψ max = M M O, we use the solution from the Bocop software to initialize the shooting method. This technique is classical, see [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF][START_REF] Goubinat | Contrôle optimal géométrique et méthodes numériques : application à un problème de montée d'un avion[END_REF] for details and examples. We then obtain a trajectory for which the switching times are t 1 ≈ 47s and t 2 ≈ 668s, the final time is t f ≈ 675s and the fuel consumption is ∆ m ≈ 860kg. There is a gap of 3kg on the fuel consumption with the solution from (CM ∆m ) and a difference about 2s on the final time. Hence, even if (SA ∆m ) is better than (CM ∆m ), the difference is very small. Remark 5.4. The Bocop software transforms an infinite dimensional optimal control problem into a finite dimensional optimization problem, applying a full time discretization of the state variables, the control variables, the dynamics and the constraints.

Let us recap in the following table, the fuel consumptions and the final times for the different problems we have encountered so far:

(SA t f ) (CM t f ) (SA ∆m ) (CM ∆m ) t f 658s 660s 675s 677s ∆ m 882kg 884kg 860kg 863kg
From this table, it is clear that the SA procedure is slightly better than the CM procedure in terms of time-toclimb and fuel consumption. In the next section we analyze the case when α ∈ [0 , 1]. This will complete the comparison between (CM α ) and (SA α ).

The minimal cost index problem

We are interested now in the so-called cost index criterion which is a convex combination of the time-to-climb and the fuel consumption. Instead of comparing directly (CM α ) and (SA α ) for α ∈ [0 , 1], we solve (P α ) with φ max = V M O, ψ max = M M O, and we show that the trajectories we get follow the SA procedure. Doing this, we conclude at the end of this section, that not only the SA procedure is better than the CM procedure but it is also the best solution we have. The problem (CM α ) is solved with the fmincon function from the Matlab software and we use the ode45 function to compute the numerical integrations. The optimal control problem (P α ) is solved with the Bocop software.

We present in Table 2, for different values of α, the optimized CAS/Mach couple (that is φ max and ψ max ), the fuel consumption and the time-to-climb, associated to the solutions of problems (CM α ) and (P α ). Since all the trajectories from problem (P α ) with φ max = V M O, ψ max = M M O follow the SA procedure, it is clear from Table 2 that the SA procedure is better than the CM procedure. Indeed, for any α, the final time and the fuel consumption are better, so is the cost index. It is interesting to notice that not only the objective function is better (that is the cost index) but also each part of the criterion (that is t f and ∆ m ). This fact is visible on Fig. 6, where the left sub-graph depicts the evolution of the fuel consumption with respect to α while the right sub-graph presents the time-to-climb.

To conclude, this numerical investigation indicates that the SA procedure of the form σ -σ s σ + is the best policy but the actual CAS/Mach procedure of the form σ -σ c1 σ c2 σ + provides good sub-optimal trajectories regarding the cost index criterion g α . The two curves are fairly similar, the relative gap between the plain and the dotted curves is about 0.3%.

Conclusion

This article is about aircraft trajectory optimization during the climbing phase considering the cost index criterion. The cost index criterion is a convex combination of the time-to-climb and the fuel consumption, parameterized by α. This parameter α is tuned by each airline. The problem is first modeled as an optimal control problem in Mayer form with a single-input affine control system and with pure state constraints. The two state constraints give CAS and Mach speeds limitations while the control variable is the flight path angle of the aircraft. In Section 4, for the minimum time-to-climb problem (α = 0), we classify the BC-extremals structures with respect to the bounds of the state constraints, that is φ max and ψ max , using the small time analysis developed in Section 3, combined with indirect multiple shooting and homotopy methods with monitoring. This classification emphasizes the role of the SA and CM procedures that we define and compare in Section 5. The CM procedure is the actual CAS/Mach procedure of the form σ -σ c1 σ c2 σ + while the SA procedure has the simple form σ -σ s σ + . Fixing the structure, that is considering the SA or CM procedure, the optimal control problem can be reduced to a simpler optimization problem in finite and small dimension, where the adjoint vector p is not needed. This is possible in particular because of the parameterization of the control in the feedback form u s (x), u c1 (x) and u c2 (x) respectively along the arcs σ s , σ c1 and σ c2 . The numerical investigation shows that the SA procedure is better than the actual CM procedure, in terms of time-to-climb and fuel consumption when considering the cost index criterion, but the difference is small. Besides, thanks to the numerical results obtained with the direct collocation method, considering the optimal control problem (P α ) with φ max = V M O and ψ max = M M O, we have also that the SA procedure is possibly the best policy.

It is worth to mention that in this article we consider only the flight path angle as control variable. This means that the thrust is assumed to be constant. It would be interesting to consider the thrust as a control variable and in this case, we would have a bi-input control system still affine in the control and more complex structures may appear. A preliminary study has been done in [START_REF] Cots | Direct and indirect methods in optimal control with state constraints and the climbing trajectory of an aircraft[END_REF] but the singular perturbation phenomenon has not been taken into account, which makes the numerical resolution much more difficult.

Remark 2 . 8 .

 28 Assuming D 0 (x) = 0, the generalized Legendre-Clebsch condition (9) becomes D 0 D 101 ≥ 0 (when α = 0) and on the set H 1 = H 01 = 0 we have D 101 = 0 ⇒ H 101 = 0 and (H 101 = 0 and p = 0) ⇒ D 101 = 0.

Figure 2 .

 2 Figure 2. The plain lines depict the evolution of the state constraints c 1 and c 2 along the state unconstrained hyperbolic trajectory of the form σ -σ s σ + , with φ max = 160 and ψ max = 0.7. The time t 1 (resp. t 2 ) represents the switching time between the negative bang (resp. the singular) arc and the singular (resp. the positive bang) arc.

Figure 3 .

 3 Figure 3. Evolution of the control and the constraint c 2 along the trajectory associated to the zero ȳ of the shooting function S c2 .

Figure 4 .

 4 Figure 4. (Left) Classification of the BC-extremal structures with respect to φ max and ψ max .To the separating blue lines is associated structures with touch points. These blue lines are computed also by homotopy and monitoring. (Right) The final time with respect to φ max and ψ max .

5. 1 .

 1 The CAS/Mach (CM) procedure

Figure 5 .

 5 Figure 5. The fuel consumption with respect to φ max and ψ max for the minimum time-toclimb problem. See the classification on Fig. 4.

  Fuel consumption and time-to-climb for different values of α associated to the solutions of problems (CM α ) and (P α ). For problem (P α ), we set φ max = V M O and ψ max = M M O.

Figure 6 .

 6 Figure 6. Evolution of the fuel consumption (on the left) and of the time-to-climb (on the right) with respect to α. The plain lines are for problem (CM α ) while the dotted lines are for problem (P α ) with φ max = V M O, ψ max = M M O. The two curves are fairly similar, the relative gap between the plain and the dotted curves is about 0.3%.

Table 1 .

 1 Medium-haul aircraft parameters.

	Data Value	Unit	Data Value	Unit
	S	122.6	m 2	C s,1 1.055e -5 kg.s -1 .N -1
	g	9.81	m.s -2 C s,2 441.54	m.s -1
	C T,1 141040	N	R	287.058	J.kg -1 .K -1
	C T,2 14909.9	m	Θ 0	288.15	K
	C T,3 6.997e -10 m -2	β	0.0065	K.m -1
	C d,1 0.0242		P 0	101325	Pa
	C d,2 0.0469

Table 2 .

 2 .4 0.6631 862.8 675.7 860.1 674.0 0.105 129.9 0.6651 862.9 674.3 860.2 672.6 0.158 130.4 0.6670 863.1 673.0 860.4 671.2 0.210 131.0 0.6690 863.4 671.7 860.7 669.9 0.263 131.6 0.6709 863.8 670.5 861.1 668.6 0.316 132.2 0.6729 864.2 669.3 861.6 667.4 0.368 132.8 0.6748 864.8 668.2 862.2 666.3 0.421 133.5 0.6767 865.5 667.1 862.9 665.2 0.474 134.2 0.6787 866.4 666.1 863.7 664.1 0.526 134.9 0.6808 867.3 665.1 864.7 663.2 0.579 135.6 0.6825 868.4 664.2 865.8 662.3 0.631 136.4 0.6844 869.7 663.4 867.1 661.4 0.684 137.2 0.6863 871.1 662.7 868.5 660.7 0.737 138.1 0.6882 872.7 662.0 870.1 660.0 0.790 139.0 0.6901 874.5 661.5 872.0 659.5 0.842 139.9 0.6920 876.6 661.0 874.0 659.0 0.895 140.9 0.6939 878.8 660.7 876.3 658.7 0.947 142.0 0.6959 881.4 660.4 878.8 658.5 1.000 143.2 0.6978 884.3 660.4 881.6 658.4

		(CM α )		(P α )
	α	CAS Mach	∆ m	t f	∆ m	t f
	0	128.9 0.6611 862.7 677.1 860.0 675.4
	0.056 129

Toulouse Univ., INP-ENSEEIHT-IRIT, UMR CNRS 5505,

rue Camichel,

Toulouse, France; e-mail: olivier.cots@irit.fr & joseph.gergaud@irit.fr & boris.wembe@irit.fr 2 Thales Avionics SA, 105 av du General Eisenhower, B.P. 1147, 31047 Toulouse Cedex, France; e-mail: damien.goubinat@fr.thalesgroup.com

V * is closed to the maximal speed below 3480 m: 148.5 m.s -1