E Trélat

B Bonnard

U Boscain

J.-B Caillau

J Gergaud

F Giannessi

D G Hull

A Rapaport

M Sigalotti

Optimal control and applications to aerospace: some results and challenges

Keywords: optimal control, Pontryagin Maximum Principle, second-order conditions, conjugate point, numerical methods, shooting method, orbit transfer, atmospheric re-entry, geometric optimal control, optimal synthesis, continuation / homotopy method, dynamical systems, mission design 49J15, 93B40, 93B27, 93B50, 65H20, 90C31, 37N05, 37N35

This article surveys the classical techniques of nonlinear optimal control such as the Pontryagin Maximum Principle and the conjugate point theory, and how they can be implemented numerically, with a special focus on applications to aerospace problems. In practice the knowledge resulting from the maximum principle is often insufficient for solving the problem in particular because of the well-known problem of initializing adequately the shooting method. In this survey article it is explained how the classical tools of optimal control can be combined with other mathematical techniques to improve significantly their performances and widen their domain of application. The focus is put onto three important issues. The first is geometric optimal control, which is a theory that has emerged in the 80's and is combining optimal control with various concepts of differential geometry, the ultimate objective being to derive optimal synthesis results for general classes of control systems. Its applicability and relevance is demonstrated on the problem of atmospheric re-entry of a space shuttle. The second is the powerful continuation or homotopy method, consisting of deforming continuously a problem towards a simpler one, and then of solving a series of parametrized problems to end up with the solution of the initial problem. After having recalled its mathematical foundations, it is shown how to combine successfully this method with the shooting method on several aerospace problems such as the orbit transfer problem. The third one consists of concepts of dynamical system theory, providing evidence of nice properties of the celestial dynamics that are of great interest for future mission design such as low cost interplanetary space missions. The article ends with open problems and perspectives.

Introduction: optimal control problems in aerospace

The purpose of this article is to provide a survey of the main issues of optimal control theory and of some geometric results of modern geometric nonlinear optimal control, with a specific focus on applications to aerospace problems. The goal is here not only to report on some classical 1 techniques of optimal control theory (in particular, Pontryagin Maximum Principle, conjugate point theory, associated numerical methods) but also to show how this theory and these methods can be significantly improved by combining them with powerful modern techniques of geometric optimal control, of the theory of numerical continuation, or of dynamical system theory. I will illustrate the different approaches under consideration with different classical but nontrivial aerospace problems: the minimal time or minimal consumption orbit transfer problem with strong or low thrust, the minimal total thermal flux atmospheric re-entry problem of a space shuttle, and space mission design using the dynamics around Lagrange points. On these examples, I will attempt to put in evidence the limits of the classical techniques of optimal control, which are in general not sufficient to provide an adequate solution to the problem, and I will show how these techniques can be considerably enhanced by combining them with some mathematical considerations that are sometimes quite deep but are an efficient (and most of the times, superior) alternative to numerical refinement procedures or other computational efforts. In particular, I will focus on three approaches that have been successfully combined with classical optimal control and that are in my opinion of primary importance in aerospace applications. The first one is geometric optimal control, which started to be developed in the early 80's and has widely proved its superiority over the classical theory of the 60's. The main objective of geometric optimal control is to develop general techniques for general classes of nonlinear optimal control problems, using in particular the concept of Lie bracket to analyze the controllability properties of nonlinear control systems and the regularity properties of optimal trajectories, and to provide optimal synthesis results. I will show how recent results of geometric optimal control can be used to provide a deep geometric insight into the atmospheric re-entry problem and lead to an efficient solution. The second technique focused on is the numerical continuation procedure, which is far from new but has been quite neglected until recently in optimal control probably because of its difficulty to be implemented efficiently and quite systematically. In the last ten years however much progress has been done that permits to understand better how this powerful procedure can be successfully applied and I will show its particular relevance on the orbit transfer problem. The third technique mentioned, that I believe to be of particular interest for future aerospace applications, is the combination with dynamical system theory. Deep mathematical results from this theory permit to put in evidence nice properties of the celestial dynamics due to Lagrange points and gravitational effects, which are particularly interesting in view of designing low cost space missions for future interplanetary exploration. This article is addressed not only to mathematicians wanting to know more about these geometric or mathematical issues associated with concrete applications, but also to engineers already acquainted with classical techniques of optimal control, wishing to get more familiar with the more modern approaches of geometric control and other mathematical notions that have demonstrated significant enhancements in classical aerospace problems, or to experts of nonlinear control wishing to learn about applications of this discipline to nontrivial examples in aerospace.

The mathematical notions whose combination to optimal control has proved its relevance are mainly from (elementary) differential geometry; they are introduced and explained step by step, although they are mostly known by many readers, and can be skipped at the first reading.

The article is structured as follows. Section 2 surveys the most well-known theoretical and numerical aspects of optimal control. It is recalled how first-order necessary conditions lead to the famous Pontryagin Maximum Principle and how it can be used in practice. Second-order optimality conditions, leading to the conjugate point theory, are then briefly surveyed. The possible numerical approaches are then described and discussed, and their limits are underlined. Section 3 is devoted to show how techniques of geometric optimal control can be used in order to provide an efficient solution to the atmospheric re-entry problem. In Section 4, the continuation method is first described, and a theoretical foundation is provided in terms of sensitivity analysis. It is then shown how it can be combined with a shooting method in order to solve different problems, such as the orbit transfer problem. Finally, in Section 5 the focus is made on the properties of the dynamics around Lagrange points and on potential applications to mission design. Throughout the article, and in the conclusion, open problems and challenges are described.

Before going further, let us start with the motivating example of the orbit transfer problem with low thrust, where the system under consideration consists of the controlled Kepler equations q(t) = -q(t) µ r(t) 3 +

T (t) m(t) , ṁ(t) = -β T (t) , (1)
where q(t) ∈ IR 3 is the position of the engine at time t, r(t) = q(t) , T (t) is the thrust at time t, and m(t) is the mass, with β = 1/I sp g 0 . Here g 0 is the usual gravitational constant and I sp is the specific impulsion of the engine. The thrust is submitted to the constraint

T (t) T max , (2)
where the typical value of the maximal thrust T max is around 0.1 Newton, for low-thrust engines.

The well-known orbit transfer problem consists of steering the engine from a given initial orbit (e.g. an initial excentric inclinated orbit) to a final one (e.g., the geostationary orbit). Controllability properties, ensuring the feasibility of the problem, have been studied in [START_REF] Bonnard | Geometric optimal control of elliptic Keplerian orbits[END_REF][START_REF] Bonnard | Mécanique céleste et contrôle de systèmes spatiaux[END_REF], based on a careful analysis of the Lie algebra generated by the vector fields of the system [START_REF] Bonnard | Geometric optimal control of elliptic Keplerian orbits[END_REF].

For this control problem, one is moreover interested in realizing this transfer for instance in minimal time, or minimizing the fuel consumption. Then the problem turns into an optimal control problem, of the form settled hereafter.

A short survey on optimal control: theory and numerics

Let n and m be nonzero integers. Consider on IR n the control system ẋ(t) = f (t, x(t), u(t)),

where f : IR n ×IR m -→ IR n is C 1 , and where the controls are bounded and measurable functions, defined on intervals [0, T (u)] of IR + , and taking their values in a subset U of IR m . Let M 0 and M 1 be two subsets of IR n . Denote by U the set of admissible controls, so that the corresponding trajectories steer the system from an initial point of M 0 to a final point in M 1 . For such a control u, the cost of the corresponding trajectory x u (•) is defined by

C(t f , u) = t f 0 f 0 (t, x u (t), u(t))dt + g(t f , x u (t f )), (4)
where f 0 : IR × IR n × IR m -→ IR and g : IR × IR n → IR are C 1 . We investigate the optimal control problem of determining a trajectory x u (•) solution of (3), associated with a control u on [0, t f ],

such that x u (0) ∈ M 0 , x u (t f ) ∈ M 1 , and minimizing the cost C. The final time t f can be fixed or not.

When the optimal control problem has a solution, we say that the corresponding control (or the corresponding trajectory) is minimizing or optimal.

The orbit transfer problem mentioned previously enters into this category of problems. The state is then x = (q, q, m) ∈ IR 7 , the control is the thrust, the set U of constraints on the control is the closed ball of IR 3 , centered at the origin, with radius T max . If one considers the minimal time problem, then one can choose f 0 = 1 and g = 0, and if one considers the minimal fuel consumption problem, then one can choose f 0 (t, x, u) = β u and g = 0.

Existence results

This is not our aim to provide elaborate existence results for optimal controls. It should just be noted that usual existence results require some convexity on the dynamics since their proof usually relies on weak compactness properties. We refer to [START_REF] Cesari | Optimization -theory and applications[END_REF] for a survey of existence results in optimal control. The result given below, whose early version was obtained by Filippov in [START_REF] Filippov | On certain questions in the theory of optimal control[END_REF], is standard.

Theorem 1. Assume that U is compact, that there may be state constraint c 1 (x) 0, . . . , c r (x) 0, where c 1 , . . . c r are continuous functions on M , and that M 0 and M 1 are compact subsets of M such that M 1 is accessible from M 0 . Assume that there exists b > 0 such that every trajectory steering the system from M 0 to M 1 is bounded by b on [0, t(u)] in C 0 norm, and that the set

Ṽ (t, x) = {(f 0 (t, x, u) + γ, f (t, x, u)) | u ∈ U, γ 0}
is convex, for every t ∈ IR and every x ∈ IR n . Then there exists an optimal control u defined on [0, t(u)] such that the corresponding trajectory steers the system from M 0 to M 1 in time t(u) and minimizing the cost.

Even though existence results would certainly deserve many further interesting discussions, this is not the objective of this paper to report on that subject. It can however be noticed that, if the set U is unbounded then in general existence results lead to optimal controls that are not necessarily in L ∞ ([0, t(u)], U ), leading then to a potential gap with the usual necessary conditions reported hereafter, which assume that the optimal control is essentially bounded. This gap may cause the so-called Lavrentiev phenomenon and raises the question of studying the regularity of the optimal control (see [START_REF] Sarychev | First-and second-order integral functionals of the calculus of variations which exhibit the Lavrentiev phenomenon[END_REF][START_REF] Sarychev | Lipschitzian regularity of minimizers for optimal control problems with control-affine dynamics[END_REF][START_REF] Torres | Lipschitzian regularity of the minimizing trajectories for nonlinear optimal control problems[END_REF][START_REF] Frankowska | Lipschitzianity of optimal trajectories for the Bolza optimal control problem[END_REF] where such issues are investigated).

First-order optimality conditions

The set of admissible controls on [0, t f ] is denoted U t f ,IR m , and the set of admissible controls on [0, t f ] taking their values in U is denoted U t f ,U .

Definition 1. The end-point mapping E : IR n × IR + × U → IR n of the system is defined by E(x 0 , T, u) = x(x 0 , T, u), where t → x(x 0 , t, u) is the trajectory solution of (3), corresponding to the control u, such that x(x 0 , 0, u) = x 0 .

The set U t f ,IR m , endowed with the standard topology of L ∞ ([0, t f ], IR m ), is open, and the end-point mapping is C 1 on U t f ,IR m (and C k whenever the dynamics are C k ).

Note that, in terms of the end-point mapping, the optimal control problem under consideration can be written as the infinite dimensional minimization problem min{C(t f , u) | x 0 ∈ M 0 , E(x 0 , t f , u) ∈ M 1 , u ∈ L ∞ (0, t f ; U )}.

(5)

Definition 2. Assume that M 0 = {x 0 }. A control u defined on [0, t f ] is said singular if and only if the differential ∂E ∂u (x 0 , t f , u) is not of full rank. Singular controls are of high importance in optimal control theory. At this step their potential influence can be stressed by noting that, in the above constrained minimization problem, the set of constraints is a local manifold around a given control u provided u is nonsingular. It is wellknown in constrained optimization theory that, in order to derive the usual necessary conditions exposed hereafter, it is in general needed that the set of constraints be (at least locally) a manifold. Hence at this step it can be easily understood that the existence of minimizing singular controls is a potential source of problems.

Lagrange multipliers

Assume for one moment that we are in the simplified situation where M 0 = {x 0 }, M 1 = {x 1 }, T is fixed, and U = IR m . That is, we consider the optimal control problem of steering the system (3) from the initial point x 0 to the final point x 1 in time T and minimizing the cost (4) among controls u ∈ L ∞ ([0, T ], IR m ). In that case, the optimization problem (5) reduces to min

E(x 0 ,T,u)=x 1 C(T, u). (6)
According to the well-known Lagrange multipliers rule (and using the C 1 regularity of our data), if u is optimal then there exists (ψ,

ψ 0 ) ∈ IR n × IR \ {0} such that ψ.dE x 0 ,T (u) = -ψ 0 dC T (u). ( 7
)
Note that, if one defines the Lagrangian L T (u, ψ, ψ 0 ) = ψE x 0 ,T (u) + ψ 0 C T (u), then this firstorder necessary condition for optimality is written in the usual form as

∂L T ∂u (u, ψ, ψ 0 ) = 0. ( 8
)
Here the main simplification was U = IR m , that is, we considered the case without control constraints. In the general case the situation is more intricate to deal with control constraints and even more when there are state constraints. When there are some control constraints, one possibility could be to take these constraints into account directly in the Lagrangian, with some additional Lagrange multiplier, as it is done e.g. in [START_REF] Bryson | Applied optimal control[END_REF]. This method leads however to weaker results than the Pontryagin Maximum Principle stated hereafter. Actually the method used by Pontryagin in order to take into account control constraints is stronger and consists of considering needle-like variations (see also Remark 2 further).

When there are some state constraints, it is also still possible to express a first-order necessary condition in terms of Lagrange multipliers as above but this has to be done in distributions spaces and the Lagrange multipliers must be expressed as Radon measures (see e.g. [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF][START_REF] Gamkrelidze | Principles of optimal control theory[END_REF][START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF][START_REF] Ioffe | Theory of extremal problems[END_REF][START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]).

Whatever simplified or general case we may consider, the first-order condition ( 8) is in this form not much tractable for practical purposes. This first-order condition can be in some sense parametrized along the trajectory, and this leads to the famous Pontryagin Maximum Principle.

Pontryagin Maximum Principle

The following statement is the most usual Pontryagin Maximum Principle, valuable for general nonlinear optimal control problems (3)-( 4), with control constraints but without state constraint. Usual proofs rely on a fixed point argument and on the use of Pontryagin cones (see e.g. [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF][START_REF] Lee | Foundations of optimal control theory[END_REF][START_REF] Alekseev | Optimal control[END_REF][START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF]).

Theorem 2. If the trajectory x(•), associated to the optimal control u on [0, t f ], is optimal, then it is the projection of an extremal (x(•), p(•), p 0 , u(•)) (called extremal lift), where p 0 0 and p(•) : [0, t f ] → IR n is an absolutely continuous mapping called adjoint vector, with (p(•), p 0 ) = (0, 0), such that

ẋ(t) = ∂H ∂p (t, x(t), p(t), p 0 , u(t)), ṗ(t) = - ∂H ∂x (t, x(t), p(t), p 0 , u(t)),
almost everywhere on [0, t f ], where H(t, x, p, p 0 , u) = p, f (t, x, u) + p 0 f 0 (t, x, u) is the Hamiltonian, and there holds

H(t, x(t), p(t), p 0 , u(t)) = max v∈U H(t, x(t), p(t), p 0 , v) (9)
almost everywhere on [0, t f ]. If moreover the final time t f to reach the target M 1 is not fixed, then one has the condition at the final time

t f max v∈U H(t f , x(t f ), p(t f ), p 0 , v) = -p 0 ∂g ∂t (t f , x(t f )). (10)
Additionally, if M 0 and M 1 (or just one of them) are submanifolds of IR n locally around x(0) ∈ M 0 and x(t f ) ∈ M 1 , then the adjoint vector can be built in order to satisfy the transversality conditions at both extremities (or just one of them)

p(0) ⊥ T x(0) M 0 , p(t f ) -p 0 ∂g ∂x (t f , x(t f )) ⊥ T x(t f ) M 1 , (11)
where T x M i denotes the tangent space to M i at the point x.

The relation between the Lagrange multipliers of the previous section and (p(•), p 0 ) is that the adjoint vector can be constructed so that (ψ, ψ 0 ) = (p(t f ), p 0 ) up to some multiplicative scalar. In particular, the Lagrange multiplier ψ is unique (up to a multiplicative scalar) if and only if the trajectory x(•) admits a unique extremal lift (up to a multiplicative scalar).

If p 0 < 0 then the extremal is said normal, and in this case, since the Lagrange multiplier is defined up to a multiplicative scalar, it is usual to normalize it so that p 0 = -1. If p 0 = 0 then the extremal is said abnormal.

It can be also noted that, in the absence of control constraints, abnormal extremals project exactly onto singular trajectories (it is evident using (8)).

Remark 1. With respect to this relation, it can be noted that in the normal case the Lagrange multiplier ψ (or the adjoint vector p(t f ) at the final time) coincides up to some multiplicative scalar with the gradient of the value function (solution of the Hamilton-Jacobi equation); see e.g. [START_REF] Clarke | The relationship between the maximum principle and dynamic programming[END_REF] for precise results.

Remark 2. In the case U = IR m (no control constraints), the maximization condition [START_REF] Bryson | Applied optimal control[END_REF] implies in particular ∂H ∂u (t, x(t), p(t), p 0 , u(t)) = 0 almost everywhere on [0, t f ]. In this form, the Pontryagin Maximum Principle is exactly the parametrized version of the first-order necessary condition (8) of the simplified case. Note that, in the absence of control constraints, the proof is quite obvious and can be found e.g. in [START_REF] Bonnard | The role of singular trajectories in control theory[END_REF][START_REF] Trélat | Contrôle optimal : théorie & applications, Vuibert, Collection[END_REF]. Note also that the maximimization condition implies as well that the quadratic form ∂ 2 H ∂u 2 (t, x(t), p(t), p 0 , u(t)) is nonpositive almost everywhere on [0, t f ]. These two conditions remain however weaker than the maximization condition [START_REF] Bryson | Applied optimal control[END_REF]. Indeed, those two conditions are local, whereas the maximization condition ( 9) is global. In the proof of the general version of the Pontryagin Maximum Principle, needle-like variations of the control are the main tool in order to derive the strong condition (9) (note that that a short proof of the Pontryagin Maximum Principle is provided in the general case, with needle-like variations and with a conic implicit function theorem, in [START_REF] Haberkorn | Convergence results for smooth regularizations of hybrid nonlinear optimal control problems[END_REF]).

Remark 3. The scalar p 0 is a Lagrange multiplier associated with the instantaneous cost. It may happen that it is equal to 0 and these cases are called abnormal. Abnormal extremals are not detected with the usual Calculus of Variations approach, because this approach postulates at the very beginning that, in a neighborhood of some given reference trajectory, there are other trajectories having the same terminal points, whose respective costs can be compared (and this leads to Euler-Lagrange equations). But this postulate fails whenever the reference trajectory is isolated: it may indeed happen that there is only one trajectory joining the terminal points under consideration (see [START_REF] Agrachev | On abnormal extremals for Lagrange variational problems[END_REF]). In this case, in some sense there is no optimization to do any more. Indeed since the trajectory joining the desired extremities is unique, then obviously it will be optimal, for every optimization criterion we may consider. These cases may appear to be quite trivial, but actually in practice this issue is far from being obvious because a priori, given some extremities, we are not able to say if the resulting problem can be solved with a normal extremal (that is, with a p 0 = -1). It could happen that it is not: this is the case for instance for certain initial and final conditions in the well-known minimal-time attitude control problem (see [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF], where such abnormals are referred to as exceptional singular trajectories).

Hence, when applying the Pontryagin Maximum Principle, we must distinguish between two extremal flows: the normal flow, with p 0 < 0 (and in general we normalize to p 0 = -1), and the abnormal one, for which p 0 = 0.

In many situations, where some qualification conditions hold, abnormal extremals do not exist in the problem under consideration, but in general it is impossible to say whether, given some initial and final conditions, these qualification conditions hold or not. Furthermore it can be noted that, when they exist, extremities of projections of abnormal extremals do not fill much of the state space (see [START_REF] Rifford | Morse-Sard type results in sub-Riemannian geometry[END_REF][START_REF] Trélat | Some properties of the value function and its level sets for affine control systems with quadratic cost[END_REF] for precise statements). Remark 4. An ultimate remark on the multiplier p 0 is about its sign. According to the convention chosen by Pontryagin, we consider p 0 0. If, instead, we adopt the convention p 0 0, then we have to replace the maximization condition (9) with a minimization condition. This just consists in considering the opposite of the adjoint vector.

If there is no integral term in (4) (that is, f 0 = 0), this does not change anything for the considerations on p 0 . The p 0 still appears in the transversality condition [START_REF] Gamkrelidze | Principles of optimal control theory[END_REF]. Note that, if M 1 = IR n (no condition on the final point), then this condition leads to p(t f ) = p 0 ∂g ∂x (t f , x(t f )), and then necessarily p 0 = 0 (otherwise the adjoint vector (p(t f ), p 0 ) would be zero and this would contradict the assertion of the Pontryagin Maximum Principle) and in that case we can normalize to p 0 = -1.

Generalizations

The Pontryagin Maximum Principle withstands many possible generalizations.

First, it can be expressed for control systems evolving on a manifold, that is, control systems of the form (3) with a dynamics f : M × N -→ T M , where M (resp., N ) is a smooth manifold of dimension n (resp., m). This situation can be of interest if for instance the system is evolving on an energy surface (it is often the case in aerospace) and/or if the controls take their values in the unit sphere (this situation occurs often in aerospace as well, for instance when the control models a thrust of constant modulus). We refer the reader e.g. to [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF] for such a version on manifolds.

The Pontryagin Maximum Principle can be generalized to wider classes of functionals and boundary conditions; for instance, periodic boundary conditions (see [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF]), systems involving some delays (see [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF][START_REF] Leitmann | The calculus of variations and optimal control. An introduction[END_REF]), intermediate conditions (see [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Arutyunov | Necessary optimality conditions for optimal control problems with intermediate constraints[END_REF][START_REF] Dmitruk | Maximum Principle for optimal control problems with intermediate constraints[END_REF]), or, more generally, hybrid systems where the dynamics may change along the trajectory, accordingly to time and/or state conditions (see [START_REF] Branicky | A unified framework for hybrid control: model and optimal control theory[END_REF][START_REF] Garavello | Hybrid necessary principle[END_REF][START_REF] Bressan | Optimal control problems on stratified domains[END_REF][START_REF] Shaikh | On the hybrid optimal control problem: theory and algorithms[END_REF][START_REF] Dmitruk | The hybrid maximum principle is a consequence of Pontryagin maximum principle[END_REF][START_REF] Haberkorn | Convergence results for smooth regularizations of hybrid nonlinear optimal control problems[END_REF]). In particular in this last case, when the system crosses a given boundary then a jump condition must hold on the adjoint vector, which means that the adjoint vector is no more continuous (but is however piecewise absolutely continuous). It can be also generalized to nonsmooth systems (see [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF]).

Probably the most difficult generalization is when there are some state constraints. In that case, we impose to the trajectories to lie in a given part of the space. From the mathematical point of view the situation is more intricate since the adjoint vector becomes a measure, evolving in some distribution space. What can be probably considered as the most general Maximum Principle statement has been derived in [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF], which can be applied to very general (possibly nonsmooth) control systems with state/control constraints. Note that a version of that result for smooth systems (with state and control constraints) has been written in [START_REF] Vinter | Optimal control. Systems & Control: Foundations & Applications[END_REF]. As explained in [START_REF] Bonnard | Mécanique céleste et contrôle de systèmes spatiaux[END_REF], the difficulty in practice is that since the adjoint vector is a vectorial measure it may admit an accumulation of atoms; in other words the measure does not admit necessarily piecewise a density, for instance there may occur an accumulation of touching points with the boundary of the state domain. To overcome this mathematical difficulty the usual argument consists of assuming that, in practice, such an accumulation does not occur, and the trajectory is a "regular" succession of free arcs (i.e., arcs inside the authorized state domain) and of boundary arcs (i.e., arcs that are at the boundary of the authorized state domain), with possibly some isolated touching points (at which the trajectory touches the boundary). This regular structure was already assumed by Pontryagin and his co-authors in their seminal book [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF], however in this book they moreover restricted to so-called first-order constraints. In few words, and roughly speaking, the order of a state constraint is the number of times one must differentiate the equality constraint along the trajectory with respect to time in order to make appear the control. The general case has been treated in [START_REF] Jacobson | New necessary conditions of optimality for control problems with statevariable inequality constraints[END_REF] in the case where the Hessian of the Hamiltonian is nondegenerate, and in [START_REF] Maurer | On optimal control problems with bounded state variables and control appearing linearly[END_REF] for the case where the controls appear linearly in the system (note that the jump conditions have been clarified in [START_REF] Arutyunov | Optimality conditions. Abnormal and degenerate problems[END_REF][START_REF] Bonnans | Well-posedness of the shooting algorithm for state constrained optimal control problems with a single constraint and control[END_REF][START_REF] Bonnans | Stability and sensitivity analysis for optimal control problems with a first-order state constraint and application to continuation methods[END_REF][START_REF] Bonnans | Revisiting the analysis of optimal control problems with several state constraints[END_REF][START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF][START_REF] Bonnard | Une approche géométrique du contrôle optimal de l'arc atmosphérique de la navette spatiale[END_REF]). Note that a nice survey on the Pontryagin Maximum Principle for control systems with general state constraints has been written in [START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF]. In any case, under this regular structure assumption, as in the hybrid situation the adjoint vector is absolutely continuous by parts, with jump conditions at the boundary (the discontinuity of the adjoint vector is along the gradient of the frontier). It can be noted however that for state constraints of order greater than or equal to three the bad chattering phenomenon of accumulation of touching points may occur in a typical way (see [START_REF] Robbins | Junction phenomena for optimal control with state-variable inequality constraints of third order[END_REF]).

Practical use of the Pontryagin Maximum Principle

In practice in order to compute optimal trajectories with the Pontryagin Maximum Principle the first step is to make explicit the maximization condition. An usual assumption to make this step feasible is to assume the so-called strict Legendre assumption, that is, to assume that the Hessian ∂ 2 H ∂u 2 (t, x, p, p 0 , u) is negative definite. Under that assumption, a standard implicit function argument permits to end up, at least locally, with a control u expressed as a function of x and p. This assumption is for instance obviously satisfied for normal extremals if one considers control affine systems with a cost that is quadratic in u. Assume for example that we are in the normal case (p 0 = -1). Then, plugging the resulting expression of the control in the Hamiltonian equations, and defining the reduced (normal) Hamiltonian by H r (t, x, p) = H(t, x, p, -1, u(x, p)), it follows that every normal extremal is solution of

ẋ(t) = ∂H r ∂p (t, x(t), p(t)), ṗ(t) = - ∂H r ∂x (t, x(t), p(t)), (12)
and this leads to define the (normal) exponential mapping.

Definition 3. The exponential mapping is defined by exp x 0 (t, p 0 ) = x(t, x 0 , p 0 ), where the solution of (12) starting from (x 0 , p 0 ) at t = 0 is denoted as (x(t, x 0 , p 0 ), p(t, x 0 , p 0 )).

In other words, the exponential mapping parametrizes the (normal) extremal flow. This mapping generalizes the concept of Riemannian ou sub-Riemannian exponential mapping. Indeed in Riemannian geometry the extremal equations correspond exactly to the cotangent formulation of the geodesics equations. Note as well that the equations [START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF] are the cotangent version of the usual Euler-Lagrange equations of the calculus of variations 1 .

The abnormal extremal flow can be parametrized as well provided that there holds such a kind of Legendre assumption in the abnormal case.

At the opposite, when the Hessian of the Hamiltonian considered above is degenerate, the situation can be far more intricate. A typical example is when one considers the minimal time problem for single-input control affine systems ẋ(t) = f 0 (x(t)) + u(t)f 1 (x(t)) without constraint on controls. In that case, the maximization condition leads to ∂H ∂u = 0, that is, there must hold p(t), f 1 (x(t)) = 0 along the corresponding extremal. By the way, note that, since the optimal control takes here its values in the interior of the domain of constraints, it is necessarily singular. To compute the control, the method consists of differentiating two times this relation with respect to t, which leads at first to p(t), [f 0 , f 1 ](x(t)) = 0 and then at p(t), where [•, •] denotes the Lie bracket of vector fields. This permits as well to express the optimal control u(t) as a function of x(t) and p(t), provided that the quantity p(t), [f 1 , [f 0 , f 1 ]](x(t)) does not vanish along the extremal. The latter condition is called strong generalized Legendre-Clebsch condition. We refer the reader to [START_REF] Bonnard | The role of singular trajectories in control theory[END_REF] for more details on this theory. It can also be shown that this kind of computation is valid in a "generic" situation (see [START_REF] Bonnard | Generic properties of singular trajectories[END_REF][START_REF] Chitour | Propriétés génériques des trajectoires singulières[END_REF][START_REF] Chitour | Genericity results for singular curves[END_REF][START_REF] Chitour | Singular trajectories of control-affine systems[END_REF]).

[f 0 , [f 0 , f 1 ]](x(t)) + u(t) p(t), [f 1 , [f 0 , f 1 ]](x(t)) = 0,
Remark 5. Note the important fact that the normal extremals are distinguished from the abnormal ones by a binary variable, namely, the variable p 0 ∈ {0, -1}. In the case where an abnormal flow is well defined, we then have to deal with two extremal flows. Intuitively, however, is is expected that the abnormal flow fills less space than the normal flow, in the sense that almost every point of the accessible set should be reached by a normal extremal. This kind of statement is however difficult to derive. There exist some results for control-affine systems and for controlaffine systems without drift, that assert that the end-points of projections of abnormal extremals fill only a negligible part of the state space (see [START_REF] Rifford | Morse-Sard type results in sub-Riemannian geometry[END_REF][START_REF] Rifford | On the stabilization problem for nonholonomic distributions[END_REF][START_REF] Trélat | Some properties of the value function and its level sets for affine control systems with quadratic cost[END_REF]). In [START_REF] Rifford | Morse-Sard type results in sub-Riemannian geometry[END_REF] in particular it is proved that for control-affine systems without drift, with quadratic cost (framework of sub-Riemannian geometry), the image normal exponential mapping is everywhere dense in the state space under usual the usual Chow condition. Remark 6. Note that the Pontryagin Maximum Principle is nothing else but a far-reaching version of the Lagrange multipliers necessary condition derived formerly. It is thus only a firstorder necessary condition for optimality, asserting that if a trajectory is optimal then it should be sought among projections of extremals joining the initial set to the final target. Conversely, the projection of a given extremal is not necessarily (locally or globally) optimal. This leads us to next section on second-order optimality conditions.

Second-order optimality conditions

Throughout this section we assume that we are in the simplified situation where M 0 = {x 0 }, M 1 = {x 1 } and U = IR m . Also, in order to consider second-order derivatives, we assume that the dynamics are at least C 2 .

Abstract conditions

In this simplified situation, we have seen that the usual first-order necessary condition for optimality is [START_REF] Frankowska | Lipschitzianity of optimal trajectories for the Bolza optimal control problem[END_REF], that is, the vanishing of the differential of the Lagrangian.

In this simplified situation where there is no constraint on the control, conditions of order two are also standard in terms of the Lagrangian. Defining as usually the intrinsic second order derivative Q T of the Lagrangian by

Q T = ∂ 2 L T ∂ 2 u (u, ψ, ψ 0 ) | ker ∂L T ∂u ,
it is well-known that a second-order necessary condition for optimality is that Q T be nonpositive (recall the agreement ψ 0 0), and a second-order sufficient condition for local optimality is that Q T be negative definite.

In this form, these conditions are not convenient for practical purposes. Fortunately, in the same way that Lagrange multipliers conditions can be parametrized into the Pontryagin Maximum Principle, the above second-order conditions can be parametrized as well along the extremals and this leads to the theory of conjugate points, briefly sketched in next section. LiuSussmann Remark 7. The above quadratic form is the one considered in the simplified situation. Such abstract conditions have been widely generalized in the literature (see e.g. [START_REF] Maurer | First and second order necessary and sufficient optimality conditions for infinite-dimensional programming problems[END_REF]).

Conjugate points

Under the strict Legendre assumption mentioned formerly, it happens that the quadratic form Q T is negative definite whenever T > 0 is small enough. This leads naturally to the following definition.

Definition 4. The first conjugate time t c along x(•) is defined as the infimum of times t > 0 such that Q t has a nontrivial kernel.

Under the strict Legendre assumption, there holds t c > 0, and this first conjugate time characterizes the (local) optimality status of the trajectory (see [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF][START_REF] Bonnard | The role of singular trajectories in control theory[END_REF][START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF][START_REF] Liu | Shortest paths for sub-Riemannian metrics on rank-two distributions[END_REF][START_REF] Sarychev | The index of second variation of a control system[END_REF]).

Theorem 3. The trajectory x(•) is locally optimal (in L ∞ topology) on [0, t] if and only if t < t c .

The following result is crucial for practical computations of conjugate times.

Theorem 4. The time t c is a conjugate time along x(•) if and only if the mapping exp x 0 (t c , •) is not an immersion at p 0 (that is, its differential is not injective).

Its proof can be found in [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF][START_REF] Bonnard | Mécanique céleste et contrôle de systèmes spatiaux[END_REF]. Essentially it states that computing a first conjugate time reduces to compute the vanishing of some determinant along the extremal. Indeed, the fact that the exponential mapping is not an immersion can be translated in terms of so-called vertical Jacobi fields. Note however that the domain of definition of the exponential mapping requires a particular attention in order to define properly these Jacobi fields according to the context: normal or abnormal extremal, final time fixed or not. A more complete exposition can be found in the survey article [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF], which provides also some algorithms to compute first conjugate times in various contexts (however, always in the case where the control can be expressed as a smooth function of x and p) and some practical hints for algorithmic purposes (see also section 2.4.2).

Generalizations, open problems and challenges

A first remark is that the conjugate point theory sketched previously can be generalized in the case where the initial and final sets M 0 and M 1 are not necessarily restricted to a single point. In that case, the notion of conjugate point must be replaced with the notion of focal point. The theory and the resulting algorithms remain however similar (see [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF][START_REF] Bonnard | Mécanique céleste et contrôle de systèmes spatiaux[END_REF]).

It should be also stressed that the above conjugate point theory holds only in the "smooth case", that is, whenever the optimal controls under consideration can be expressed as smooth functions of x and p (thus, essentially, when there is no constraint on the control, although this assumption involves some possible cases where there are some control constraints), and without any state constraint. In this theory, the definition and the computation of conjugate points are based on second-order conditions which do not involve in particular bang-bang situations where the control is discontinuous and consists of successive arcs saturating the constraints.

In this case where the extremal controls are continuous, the literature on first and/or second order sufficient conditions is vast (see also [START_REF] Bliss | Lectures on the Calculus of Variations[END_REF][START_REF] Hestenes | Application of the theory of quadratic forms in Hilbert spaces to the calculus of variations[END_REF][START_REF] Maurer | First and second order sufficient optimality conditions in mathematical programming and optimal control[END_REF][START_REF] Maurer | Second order sufficient conditions for optimal control problems with free final time: The Riccati approach[END_REF][START_REF] Milyutin | Calculus of Variations and Optimal Control[END_REF][START_REF] Zeidan | The Riccati equation for optimal control problems with mixed state-control constraints: Necessity and sufficiency[END_REF] and references therein), and there exist also numerical procedures to test second order sufficient conditions that are based on the Riccati equation (of course these procedures are equivalent to the one described previously in terms of the exponential mapping). We refer also the reader to [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF][START_REF] Agrachev | On abnormal extremals for Lagrange variational problems[END_REF][START_REF] Arutyunov | Optimality conditions. Abnormal and degenerate problems[END_REF][START_REF] Trélat | Asymptotics of accessibility sets along an abnormal trajectory[END_REF]] (and references therein) for extensions of such theories to the abnormal case.

A conjugate time theory has been developed in the bang-bang case. In [START_REF] Sussmann | Envelopes, conjugate points and optimal bang-bang extremals[END_REF], the definition of conjugate points for bang-bang situations is based on the theory of envelopes and is used for the study of the structure of locally optimal bang-bang trajectories for the minimal time control problem for single-input control-affine systems in dimension two and three. These results were generalized in [START_REF] Krener | The structure of small-time reachable sets in low dimensions[END_REF]. Second order necessary and/or sufficient conditions of optimal control problems with nonlinear control systems and discontinuous controls have been developed in [START_REF] Agrachev | Strong optimality for a bang-bang trajectory[END_REF][START_REF] Dunn | Second-order optimality conditions in sets of L ∞ functions with range in a polyhedron[END_REF][START_REF] Osmolovskii | Transformation of quadratic forms to perfect squares for broken extremal[END_REF][START_REF] Osmolovskii | Equivalence of second order optimality conditions for bang-bang control problems. Part 1: main results[END_REF][START_REF] Poggiolini | State-local optimality of a bang-bang trajectory: a Hamiltonian apprach[END_REF][START_REF] Sarychev | First-and second-order sufficient optimality conditions for bang-bang controls[END_REF] and references therein. A different point of view is provided in [START_REF] Agrachev | Symplectic methods for optimization and control[END_REF][START_REF] Noble | Sufficient conditions for relative minima of broken extremals in optimal control theory[END_REF] (see also [START_REF] Agrachev | Strong optimality for a bang-bang trajectory[END_REF][START_REF] Poggiolini | State-local optimality of a bang-bang trajectory: a Hamiltonian apprach[END_REF]) in terms of extremal fields, and consists of embedding the reference trajectory into a local field of broken extremals (corresponding to piecewise continuous controls). The occurrence of a conjugate point is then related with an overlap of the flow near the switching surface. In [START_REF] Maurer | Optimization methods for the verification of second order sufficient conditions for bang-bang controls[END_REF][START_REF] Maurer | Second order sufficient conditions for time-optimal bang-bang control problems[END_REF] optimization methods are given to test second order sufficient optimality conditions for such bang-bang situations. The idea is to reduce the problem to the finitedimensional subproblem consisting of moving the switching times and a second variation is defined as a certain quadratic form associated to this subproblem. Then, finding a conjugate time consists in testing the positivity of that quadratic form (and this can only happen at a switching time). We refer the reader to [START_REF] Silva | Asymptotic approach on conjugate points for minimal time bang-bang controls[END_REF] where a brief survey with a unifying point of view of all these (apparently different) approaches has been written in the introduction.

This idea of parametrizing the trajectories by their switching times also extends to the case where there are state constraints, and has been carried out in [START_REF] Zeidan | The Riccati equation for optimal control problems with mixed state-control constraints: Necessity and sufficiency[END_REF][START_REF] Augustin | Computational sensitivity analysis for state constrained optimal control problems[END_REF][START_REF] Malanowski | Sensitivity analysis for parametric control problems with control-state constraints[END_REF][START_REF] Malanowski | Sensitivity analysis for state constrained optimal control problems[END_REF][START_REF] Maurer | Second-order sufficient conditions for optimal control problems with mixed control-state constraints[END_REF], and in these references the test of positivity of the resulting quadratic form permits to derive sufficient conditions for local optimality. The concept of conjugate point is however not so straightforward to define, since one has to investigate stability properties of the structure. Indeed the loss of positivity could be due to a change in the global structure of the trajectory. Hence the definition of a conjugate point is not clear.

The situation is similar for trajectories that involve both bang and singular arcs. At the moment there is no general conjugate point theory which would involve such cases. Note however that there are some partial results such as in [START_REF] Poggiolini | Sufficient optimality conditions for a bang-singular extremal in the minimum time problem[END_REF] where the authors study the particular case of a bang-singular extremal.

This is an open (and important) problem to derive a complete conjugate point theory that would consist of any possible smooth, bang, singular, or boundary arcs.

From the algorithmic point of view, note that, although the theory of conjugate times in the bang-bang case has been well developed, the computation of conjugate times in the bang-bang case is difficult in practice with the algorithms of the previously mentioned references (see in particular [START_REF] Maurer | Optimization methods for the verification of second order sufficient conditions for bang-bang controls[END_REF][START_REF] Maurer | Second order sufficient conditions for time-optimal bang-bang control problems[END_REF] and references therein). Besides, in the smooth case, as explained in the previous section efficient tools are available (see [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]). In [START_REF] Silva | Smooth regularization of bang-bang optimal control problems[END_REF][START_REF] Silva | Asymptotic approach on conjugate points for minimal time bang-bang controls[END_REF] a regularization procedure is proposed which allows the use of these tools for the computation of the first conjugate time of a bang-bang situation for a single-input control-affine system, by showing the convergence of the conjugate time of the regularized system to the conjugate time of the initial bang-bang system. This result is of interest because it provides an efficient way to compute conjugate points in the bang-bang case. It is an open problem to extend that kind of result to more general systems and more general situations.

Numerical methods in optimal control

It is usual to distinguish between two kinds of numerical approaches in optimal control: direct and indirect methods. Roughly speaking, direct methods consist in discretizing the state and the control and thus reduce the problem to a nonlinear optimization problem with constraints. Indirect methods consist of solving numerically the boundary value problem derived from the application of the Pontryagin Maximum Principle, and lead to the shooting methods.

Direct methods

They are the most evident methods when one addresses practically an optimal control problem. By discretizing both the state and the control, the problem reduces to a nonlinear optimization problem in finite dimension, or nonlinear programming problem, of the form

min Z∈C F (Z), (13)
where Z = (x 1 , . . . , x N , u 1 , . . . , u n ), and

C = {Z | g i (Z) = 0, i ∈ 1, . . . , r, g j (Z) 0, j ∈ r + 1, . . . , m}. (14)
There exist many ways to carry out such discretizations. In any case, one has to choose finite dimensional representations of the control and of the state, and then express in a discrete way the differential equation representing the system. Once all static or dynamic constraints have been transcribed into a problem with a finite number of variables, one is ought to solve the resulting optimization problem with constraints, using some adapted optimization method. Let us first explain hereafter one possible very simple way of such a discretization. Consider

a subdivision 0 = t 0 < t 1 < • • • < t N = t f of the interval [0, t f ].
Controls are discretized in such a way that they are piecewise constant on this subdivision (with values in U ). Moreover, we choose a discretization process of ordinary differential equations, for instance choose here (to simplify) the classical explicit Euler method. Setting h i = t i+1 -t i , we obtain

x i+1 = x i + h i f (t i , x i , u i ).
There exist of course an infinite number of possible variants. On one hand, one may discretize the set of admissible controls by piecewise constant, or piecewise affine controls, or splines, etc. On the other hand, there exist many methods in order to discretize ODEs, such as Euler methods (implicit or explicit), middle point, Heun, Runge-Kutta, Adams-Moulton, etc (see for instance [START_REF] Stoer | Introduction to Numerical Analysis[END_REF]). The choice of the method is guided by the problem under consideration. Here we choose the Euler method for the simplicity of its writing but in practice it should be avoided because it is too much rough.

The previous discretization process leads to the nonlinear programming problem min C(x 0 , . . . , x N , u 0 , . . . , u N ),

x i+1 = x i + h i f (t i , x i , u i ), u i ∈ U, i = 0, . . . , N -1,
which is a problem of the form [START_REF] Ioffe | Theory of extremal problems[END_REF].

Note that this kind of method is easy to implement in the sense that it does not require a precise a priori knowledge of the optimal control problem. Moreover, it is easy to take into account some state constraints or any other kinds of constraints in the optimal control problem. In this sense, this approach is not sensitive to the model.

From a general point of view, this amounts to choosing a discretization of the control, of the state, in some finite dimension spaces, using specific Galerkin bases. Typically, one may choose piecewise polynomial approximations. The ODE and state or control constraints are satisfied only at the discretization points. We thus arrive at a finite dimensional nonlinear optimization problem of the form [START_REF] Ioffe | Theory of extremal problems[END_REF].

The numerical implementation of such a nonlinear programming problem is standard. It can be achieved, for instance, using a penalty method, or more usually using a SQP method (sequential quadratic programming). The aim of these methods is to reduce the problem to subproblems, easier to treat, without constraints, by using penalty functions to handle constraints, or applying Kuhn-Tucker necessary conditions for optimization problems with constraints. For the problem ( 13)-( 14), SQP methods permit to solve iteratively the optimality system derived from the application of Kuhn-Tucker conditions, using Newton or quasi-Newton methods. At each step of the iteration process, one chooses a quasi-Newton method in order to estimate the Hessian of the Lagragian function associated to the nonlinear programming problem, and one solves a quadratic programming subproblem based on a quadratic approximation of the Lagrangian function. For more details, see e.g. [START_REF] Fletcher | Practical Methods of Optimization[END_REF][START_REF] Gill | Practical optimization[END_REF][START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF].

From the point of view of practical implementation, in the last years much progress has been done in the direction of combining automatic differentiation softwares (such as the modeling language AMPL, see [START_REF] Fourer | AMPL: A modeling language for mathematical programming[END_REF]) with expert optimization routines (such as the open-source package IPOPT, see [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF], carrying out an interior point optimization algorithm for large-scale differential algebraic systems, combined with a filter line-search method). With such tools it has become very simple to implement with only few lines of code difficult (nonacademic) optimal control problems, with success and within a reasonable time of computation. Even more, web sites such as NEOS (http://neos-server.org/neos/) propose to launch online such kinds of computation: codes can be written in a modeling language such as [START_REF] Fourer | AMPL: A modeling language for mathematical programming[END_REF] (or others) and can be combined with many optimization routines (specialized either for linear problems, nonlinear, mixed, discrete, etc). Note that there exists a large number (open-source or not) of automatic differentiation softwares and of optimization routines, it is however not our aim to provide a list of them. They are easy to find on the web.

For an excellent survey on direct methods with a special interest to applications in aerospace, we refer the reader to [START_REF] Betts | Survey of numerical methods for trajectory optimization[END_REF][START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF].

There exist many approaches to discretize an optimal control problem. Many of them are described in [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF] (in this survey book sparsity issues, very important in practice, are also discussed). Among these different approaches, we quote the following.

Speaking in a general way, collocation methods consist of choosing specific points or nodes on every subinterval of a given subdivision of the time interval. From the point of view of the discretization spaces, these methods consist in approximating the trajectories (and/or the control functions) by polynomials on each subinterval. Then the collocation conditions state that the derivatives of the approximated state match exactly with the dynamics at the nodes mentioned previously. Note that Runge-Kutta discretizations are a particular case. We refer the reader to [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF][START_REF] Hargraves | Direct trajectory optimization using nonlinear programming and collocation[END_REF] for the use of such methods in aerospace.

In spectral or pseudospectral methods, the above nodes are chosen as the zeros of special polynomials, such as Gauss-Legendre or Gauss-Lobatto polynomials. Equivalently, these polynomials serve as a basis of approximation spaces for the trajectories and the controls. Since they share nice orthogonality properties, the collocation conditions turn into constraints that are easily tractable for numerical purposes. We refer the reader to [START_REF] Elnagar | Pseudospectral Chebyshev optimal control of constrained nonlinear dynamical systems[END_REF][START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF][START_REF] Ross | Legendre pseudospectral approximations of optimal control problems[END_REF] and to the references therein for more details on these approaches.

There exist also some probabilistic approaches, such as the method described in [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMIrelaxations[END_REF] (see also [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]), which consists of expressing the optimal control problem in measure spaces and then of seeking the optimal control as an occupation measure, which is approximated by a finite number of its moments. This method uses algebraic geometry tools and permits as above to reduce the optimal control problem to some finite dimensional optimization problem (more precisely, involving LMI inequalities).

Remark 8. Another approach to optimal control problems, which can be considered (although it can be discussed) as a direct method, consists of solving the Hamilton-Jacobi equation satisfied by the value function, that is, the optimal cost for the optimal control problem of reaching a given point. The Hamilton-Jacobi equation is of the form

∂S ∂t + H 1 (x, ∂S ∂x ) = 0,
where H 1 (x, p) = max u∈U p, f (x, u) -f 0 (x, u) is the reduced normal Hamiltonian obtained from the Pontryagin Maximum Principle (see [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF][START_REF] Cannarsa | Semiconcave functions, Hamilton-Jacobi equations, and optimal control[END_REF][START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF] for the theoretical problem of existence and/or uniqueness of viscosity solutions). Although it can be discussed, we tend to place numerical methods for solving the Hamilton-Jacobi equation among direct methods. Consider a discretization (x ī) of the space, where ī = (i 1 , . . . , i n ) ∈ Z n , assumed to be regular for simplification, and consider a regular subdivision (t j ) of the time interval. Denote by h = (h 1 , . . . , h n ) the space step, and by k = t j+1 -t j the time step. Let S ī,j denote the approximate value of S(t j , x ī). One should approximate ∂S ∂xp (x ī) with a backward (resp. forward) finite difference whenever f p (x ī, u) is positive (resp. negative). For every real number a, set a + = max(a, 0) = a+|a| 2 , a -= min(a, 0) = a-|a| 2 . For every p ∈ {1, . . . , n}, denote by e p = (0, . . . , 1, . . . , 0), the "1" being in p-th position. We get the explicit scheme

0 = S ī,k+1 -S ī,k k + max u∈U n p=1 f p (x, u) + S ī,k -S ī-ep,k h p + f p (x, u) - S ī+ep,k -S ī,k h p -f 0 (x ī, u) .
There exist many discretization processes. The previous explicit finite difference scheme is the most simple, but it is not much efficient. We can choose schemes of larger order. Notice the existence of level set methods (see [START_REF] Sethian | Level set methods and fast marching methods[END_REF]), which consist of computing at each step of the iteration process the level set of the value function S solution of the Hamilton-Jacobi equation. Very efficient in low dimension, these methods consist of making evolve the wave front of the value function, starting from an initial point or set. The algorithmic complexity is a linear function of the number of discretization points. These methods have been implemented in a very efficient way for problems of low dimension (typically 3). The construction of such schemes is however not easy, and, in function of the equation, one should be able to elaborate a stable and consistant scheme, see [START_REF] Sethian | Level set methods and fast marching methods[END_REF] for examples. Notice that, as previously, taking into account some state constraints is not a problem. One should indeed impose to the value function to take an infinite (numerically, very large) value on the forbidden domain.

Once switching curves have been localized, one can refine the discretization in the neighborhood of these curves, in order to get a better accuracy.

Indirect methods

The idea is as follows. Consider the optimal control problem (3)-( 4) and assume, at first, that the final time t f is fixed. The Pontryagin Maximum Principle provides a necessary condition for optimality, and states that every optimal trajectory is the projection of an extremal, as explained before. After having made explicit the maximization condition, we reduce the problem to an extremal system of the form ż(t) = F (t, z(t)), where z(t) = (x(t), p(t)), and initial, final, transversality conditions, are of the form R(z(0), z(t f )) = 0. Finally, we get the boundary value problem (BVP)

ż(t) = F (t, z(t)), R(z(0), z(t f )) = 0. ( 15
)
Denote by z(t, z 0 ) the solution of the Cauchy problem ż(t) = F (t, z(t)), z(0) = z 0 , and set G(z 0 ) = R(z 0 , z(t f , z 0 )). The boundary value problem ( 15) is then equivalent to solving G(z 0 ) = 0, that is, one should determine a zero of the function G. This can be achieved in practice by using a Newton like method.

If the final time t f is free, one can reduce the problem to the previous formulation by considering t f as an auxiliary unknown. We thus augment the dimension of the state by considering the additional equation dt f dt = 0. The same trick can be used if the control is bang-bang, in order to determine switching times. It may be however better to use the transversality condition on the Hamiltonian function, when the final time is free. This method sketched above is called the shooting method. It has many possible refinements, among which the multiple shooting method. The latter method consists of subdividing the time interval [0, t f ] in N intervals [t i , t i+1 ], and of considering as unknowns the values z(t i ). One should then take into account gluing conditions at every time t i (continuity conditions). The aim is actually to improve the stability of the method. A classical reference as to the multiple shooting algorithm is [START_REF] Stoer | Introduction to Numerical Analysis[END_REF].

More precisely, applying the maximum principle reduces the problem to a BVP problem of the form

ż(t) = F (t, z(t)) =            F 0 (t, z(t)) if t 0 t < t 1 F 1 (t, z(t)) if t 1 t < t 2 . . . F s (t, z(t)) if t s t t f (16)
where z = (x, p) ∈ IR 2n (p is the adjoint vector), and t 1 , t 2 , . . . , t s ∈ [t 0 , t f ] are either switching times, or, in case of state constraints, junction times with a boundary arc, or contact times with the boundary. Moreover, there are some continuity conditions on the state and on the adjoint vector at switching times. In case of state constraints, one has jump conditions on the adjoint vector, and conditions on the constraint at junction or contact points, see [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF][START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations by multiple-shooting techniques[END_REF][START_REF] Bryson | Applied optimal control[END_REF][START_REF] Jacobson | New necessary conditions of optimality for control problems with statevariable inequality constraints[END_REF][START_REF] Maurer | On optimal control problems with bounded state variables and control appearing linearly[END_REF][START_REF] Miele | Recent advances in the optimization and guidance of aeroassisted orbital transfers[END_REF][START_REF] Pesch | A practical guide to the solution of real-life optimal control problems[END_REF]). Moreover, there are some limit conditions on the state, on the adjoint vector (transversality conditions), and on the Hamiltonian if the final time is free. Note that the final time is a priori unknown. Moreover, in the multiple shooting method, the number s of switchings has to be fixed; it is determined if possible by a preliminary geometric study of the problem.

The multiple shooting method consists in subdividing the interval [t 0 , t f ] in N subintervals, the value of z(t) at the left side of each subinterval being unknown. More precisely, let

t 0 < σ 1 < • • • < σ k < t f be a fixed subdivision of [t 0 , t f ].
At each point σ j , the function z is continuous. One can consider σ j as a fixed switching point, at which there holds z(σ + j ) = z(σ - j ). Now define the nodes {τ 1 , . . . , τ m } = {t 0 , t f } ∪ {σ 1 , . . . , σ k } ∪ {t 1 , . . . , t s }. Finally, we arrive at the BVP problem

• ż(t) = F (t, z(t)) =            F 1 (t, z(t)) if τ 1 t < τ 2 F 2 (t, z(t)) if τ 2 t < τ 3 . . . F m-1 (t, z(t)) if τ m-1 t τ m • ∀j ∈ {2, . . . , m -1} r j (τ j , z(τ - j ), z(τ + j )) = 0 • r m (τ m , z(τ 1 ), z(τ m )) = 0 (17)
where τ 1 = t 0 is fixed, τ m = t f , and the r j represent interior or limit conditions.

The stability of the method is improved by augmenting the number of nodes. It is the advantage of the multiple shooting method, contrary to the simple shooting one, in which errors may grow exponentially in function of t f -t 0 (see [START_REF] Stoer | Introduction to Numerical Analysis[END_REF]). Of course, there are much more unknowns in the multiple shooting method as in the simple shooting method, but note that the numerical integration of the system ( 16) can be parallelized.

Set z + j = z(τ + j ), and denote by z(t, τ j-1 , z + j-1 ) the solution of the Cauchy problem ż(t) = F (t, z(t)), z(τ j-1 ) = z + j-1 . There holds z(τ - j ) = z(τ - j , τ j-1 , z + j-1 ). Interior and limit conditions are written as

∀j ∈ {2, . . . , m -1} r j (τ j , z(τ - j , τ j-1 , z + j-1 ), z + j ) = 0, r m (τ m , z + 1 , z(τ - m , τ m-1 , z + m-1 )) = 0. ( 18
) Set Z = (z + 1 , τ m , z + 2 , τ 2 , . . . , z + m-1 , τ m-1 ) T ∈ IR (2n+1)(m-1) (where z ∈ IR 2n ). Then, conditions (18) hold if G(Z) =      r m (τ m , z + 1 , z(τ - m , τ m-1 , z + m-1 )) r 2 (τ 2 , z(τ - 2 , τ 1 , z + 1 ), z + 2 ) . . . r m-1 (τ m , z(τ - m-1 , τ m-2 , z + m-2 ), z + m-1 )      = 0. ( 19
)
Hence, this amounts to determining a zero of the function G, which is defined on a vector space whose dimension is proportional to the number of switching and subdivision points. The equation G = 0 is then iteratively solved by using a Newton type method.

From the practical implementation point of view, note on the one hand that there exist many variants of Newton methods, among which the Broyden method or the Powell hybrid method are quite competitive (see [START_REF] Stoer | Introduction to Numerical Analysis[END_REF]). On the other hand, note that, as for direct methods, the shooting methods can be combined with automatic differentiation. Here, the use of automatic differentiation can help to generate the Hamiltonian equations of extremals. This is particularly useful when one works on a problem whose model is not completely fixed. In [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] the authors provide the description for the package COTCOT (Conditions of Order Two and COnjugate Times), available for free on the web (http://www.n7.fr/apo/cotcot/), implementing the following issues:

• automatic generation in Fortran of the equations of the Pontryagin Maximum Principle (automatic differentiation with Adifor );

• automatic creation of mex files for Matlab;

• Fortran codes for the numerical integration, the shooting method, and the computation of conjugate times, interfaced with Matlab.

Remark 9. It must be noted that, when implementing a shooting method, the structure of the trajectory should be known in advance, particularly in the case where the trajectory involves singular arcs (see e.g. [START_REF] Maurer | Numerical solution of singular control problems using multiple shooting techniques[END_REF][START_REF] Bonnans | Singular arcs in the generalized Goddard's problem[END_REF]). This remark shows the importance of being able to determine at least locally the structure of optimal trajectories: this is one of the main issues of geometric optimal control theory, as explained further in this article (see sections 3.2 and 3.3).

Remark 10. Proving that the shooting method is feasible amounts to proving that the Jacobian of the mapping G defined by ( 19) is nonzero. To simplify, if one considers the simplified situation of section 2.3, in this case the (single) shooting method is well-posed at time t, locally around p 0 if and only if the exponential mapping exp x 0 (t, •) is an immersion at p 0 . In other words, according to Theorem 4, the shooting method is feasible (well-posed) if and only if the final time under consideration is not a conjugate time. This argument can be generalized to far more general situations. First of all, if the initial and final sets are not restricted to single points, the above argument still holds except that the notion of focal point has to be used instead of conjugate point (see section 2.3.3). Note that a modification of the shooting method is proposed in [START_REF] Aronna | A well-posed shooting algorithm for optimal control problems with singular arcs[END_REF], which consists in adding unknowns to the method (so that there are more unknowns than equations) to overcome partially the problem of a priori structure determination, and then the Newton method must be adapted with the use of pseudo-inverse. In [START_REF] Bonnans | Well-posedness of the shooting algorithm for state constrained optimal control problems with a single constraint and control[END_REF][START_REF] Bonnans | Second-order analysis for optimal control problems with pure state constraints and mixed control-state constraints[END_REF] it is shown that the shooting method is well-posed also in the presence of control and state constraints, provided that a certain second-order coercivity holds; this second-order condition is not translated in terms of conjugate points but this could be probably done if the corresponding conjugate point theory would exist (see section 2.3.3).

An open problem

If we summarize the main issues of the previous direct and indirect approaches, we realize that direct methods consist of 1. discretizing first the differential system, the cost, in order to reduce the optimal control problem to an usual nonlinear minimization problem with constraints (the dimension being as larger as the discretization is finer), 2. and then dualizing, by applying e.g. a usual Lagrange-Newton method to the nonlinear minimization problem (applying Kuhn-Tucker and then a Newton method to solve the resulting optimality system), whereas indirect methods consist of 1. first dualizing the optimal control problem, by applying the Pontryagin Maximum Principle (or, equivalently, the Lagrange multipliers necessary condition for optimality in infinite dimension), 2. and then discretizing, by applying a shooting method (that is, a Newton method composed with a numerical integration method).

In shorter words, direct methods consist of 1) discretize, 2) dualize, and indirect methods consist of the converse: 1) dualize, 2) discretize. It is natural to wonder whether this diagram is commutative or not, under usual approximation assumptions. It happens that, even under usual assumptions of consistency and stability (Lax scheme), it is not. Although it is very simple to see that, under these classical assumptions, the indirect approach is convergent, the direct method may diverge, whenever it was not conveniently designed. That is, although one chooses a convergent method in order to integrate the system, a convergent method in order to discretize the cost, the consistency and stability properties of the numerical schemes are not sufficient to ensure the convergence of the resulting direct method. Very simple counterexamples are provided in [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF].

It is not obvious to obtain simple conditions on the schemes ensuring the convergence of the resulting direct method, and up to now there exist only few positive results. The results of [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF] assert the convergence for "smooth" problems provided that the underlying discretization method be based on a Runge-Kutta method whose all coefficients are positive. The smoothness assumptions mean that the optimal controls under consideration take their value in the interior of the authorized domain of control (so that the maximization condition of the Pontryagin Maximum Principle reduces to ∂H ∂u = 0) and that coercivity second-order conditions hold, ensuring the smoothness of the optimal controls (as in Section 2.3), in the continuous case as well as in the discrete case. This is for instance the case for linear quadratic problems. We refer also the reader to [START_REF] Bonnans | Computation of order conditions for symplectic partitioned Runge-Kutta schemes with application to optimal control[END_REF] for further comments on this result and for other considerations on symplectic integrators. The class of Legendre pseudospectral methods is up to now the other one for which the commutation issues have been proved (see [START_REF] Elnagar | Pseudospectral Chebyshev optimal control of constrained nonlinear dynamical systems[END_REF][START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF][START_REF] Ross | Legendre pseudospectral approximations of optimal control problems[END_REF][START_REF] Ruths | Optimal control of inhomogeneous ensembles[END_REF] and see [START_REF] Ross | A roadmap for optimal control: the right way to commute[END_REF] for a detailed discussion on the commutation properties).

Apart from those few results, up to our knowledge the situation is still open in the general case and there do not exist any simple criteria or any systematic method to build adapted numerical schemes for discretizing the differential system and the cost in order to ensure the convergence of the resulting direct method. As explained above, since in the smooth case the conditions ensuring the commutation of the diagram rely on second-order conditions, the problem is clearly related to the theory of conjugate points, in the sense that, in order to handle the general case, there is need for a general conjugate point theory involving all possible situations (smooth, bang, singular arcs, state constraints). The, numerical schemes should be designed in such a way to ensure coercivity properties in the discretized second-order conditions under consideration.

It can be noted that this discrepancy in the dualization-discretization diagram arises as well in the infinite dimensional setting e.g. when one is interested to carry out practically the so-called HUM method (Hilbert Uniqueness Method), which is roughly speaking the optimal control problem of steering an infinite dimensional linear control system from a given point to a final point by minimizing the L 2 norm of the control (linear quadratic problem in a Hilbert space). In the case of the wave equation a phenomenon of interference of highfrequencies with the mesh has been put in evidence, that causes the divergence of the method (see [START_REF] Zuazua | Propagation, observation, control and numerical approximation of waves approximated by finite difference method[END_REF] and references therein for more details and more possible remedies, see also [START_REF] Labbé | Uniform controllability of semidiscrete approximations of parabolic control systems[END_REF] for a general result of convergence in the parabolic case). The literature is quite abundant for this commutation problem in the infinite dimensional framework, however the situation is still not well understood, in particular for hyperbolic equations where the question is raised as well of deriving a systematic way to build adapted shemes so that discretization and dualization commute.

Comparison between methods

We can sketch a brief comparison between both direct and indirect approaches, although such comments are a bit of caricatural. Anyway, it can be said that direct methods have the following advantages on indirect methods: they do not require any a priori theoretical study, in particular, one does not have to know a priori the structure of switchings; they are more robust, the model can be easily modified, and they are less sensitive to the choice of the initial condition. Moreover it is easy to take into account some constraints of any possible kind. However, it is difficult to reach with direct methods the precision provided by indirect methods. The direct discretization of an optimal control problem often causes several local minima. Direct methods require a large amount of memory, and thus may become inefficient if the dimension of the space is too large or if the problem cannot be easily parallelized or does not have an evident sparse structure.

The advantages of indirect methods are their extremely good numerical accuracy. Indeed since they rely on a Newton method they inherit of the very quick convergence properties of the Newton method. Moreover the shooting methods can, by construction, be parallelized, and their implementation can thus be achieved on a cluster of parallel computers. They however suffer from the following drawbacks: the optimal controls are computed in an open-loop form; they are based on the maximum principle, which gives a necessary condition for optimality only, and thus one should be able to check, a posteriori, the optimal status of the computed trajectory (with conjugate point theory); the method is not soft, in the sense that for instance the structure of switchings has to be known a priori. Furthermore it is not easy to introduce state constraints, because, on one hand, this requires to apply a maximum principle with state constraints, and on the other hand, the presence of state constraints may imply a very intricate structure of the optimal trajectory, in particular the structure of switchings. The main drawback of the shooting methods is that they are difficult to make converge. Indeed since they are based on a Newton method, they suffer from the usual drawback of the Newton method, that is, they may be very difficult to initialize properly. In other words, to make converge a shooting method one should be able to guess good initial conditions for the adjoint vector. Indeed the domain of convergence of the Newton method may happen to be very small, depending on the optimal control problem.

There exist many solutions to overcome the different flaws of both approaches. There is however no universal answer and the choice of the method should be guided by the practical problem under consideration and by the experience (note again the excellent surveys [START_REF] Betts | Survey of numerical methods for trajectory optimization[END_REF][START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF][START_REF] Pesch | A practical guide to the solution of real-life optimal control problems[END_REF][START_REF] Stryk | Direct and indirect methods for trajectory optimization[END_REF]). Speaking however in a general way, a first idea for a reasonable solution consists of combining both direct and indirect approaches, thus obtaining a so-called hybrid method. When one addresses an optimal control problem, one could indeed try at first to implement a direct method. In such a way, one can hope to get a first (maybe rough) approximation of the optimal trajectory and a good idea of the structure of switchings and of the associated adjoint vector. If one wishes more numerical accuracy, one can then carry out an indirect method, hoping that the result provided by the direct method gives a sufficient approximation, thus providing an initial point hopingly belonging to the domain of convergence of the shooting method. Combining in such a way both direct and indirect methods, one can take benefit of the extremely good accuracy provided by the shooting method, reducing considerably the drawback due to the smallness of the domain of convergence. Applying first a direct method, one can obtain an approximation of the adjoint vector. Indeed, the total discretization method consists of solving a nonlinear programming problem with constraints. The Lagrange multipliers associated to this problem give an approximation of the adjoint vector (see [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Stryk | Direct and indirect methods for trajectory optimization[END_REF][START_REF] Grimm | Adjoint estimation from a direct multiple shooting method[END_REF]).

By the way, among the many variants of direct and indirect approaches, we mention here the possibility of designing hybrid methods, neither direct or indirect, consisting of solving the boundary value problem resulting from the application of the PMP, not by a Newton method, but by an optimization method, in which the unknowns may for instance only consist of the initial adjoint vector, and the minimization functional is the cost, seen as a function of the initial adjoint vector (there are many possible various formulations for such problems). Furthermore we should quote the so-called direct multiple shooting method (see [START_REF] Bock | A multiple shooting algorithm for direct solution of optimal control problems[END_REF][START_REF] Gerdts | Direct shooting method for the numerical solution of higher-index DAE optimal control problems[END_REF]), based on constrained nonlinear programming, where the optimization variables are, similarly to the multiple shooting method, the states at some nodes, and where the controls are parametrized over the intervals between the nodes by well chosen functions. The advantage of such an approach is that it can be efficiently parallelized and it has nice sparsity features (see [START_REF] Diehl | Real-time optimization and nonlinear model predictive control of processes governed by differential-algebraic equations[END_REF][START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF] for variants).

In the present article it is our aim to focus on applications of optimal control to aerospace, and in such problems indirect methods are often priviledged because, although they are difficult to make converge, they offer a very good numerical accuracy. Hence in the sequel of that article we will describe several optimal control problems in aerospace, providing some methods in order to make converge the shooting method:

• a geometric insight (geometric optimal control tools) for the problem of atmospheric reentry of a space shuttle (section 3),

• the continuation method for orbit transfer problems (section 4),

• dynamical systems theory for interplanetary mission design (section 5).

3 Geometric optimal control and applications to the atmospheric re-entry problem

In this section we focus on the problem of the atmospheric re-entry of a space shuttle controlled by its bank angle, and where the cost to minimize is the total thermal flux. The engine is moreover submitted to state constraints on the thermal flux, the normal acceleration and the dynamic pressure. It is our aim here to show how results of geometric optimal control can help to make converge the shooting method.

The atmospheric re-entry problem

The precise problem under consideration is the following. We call atmospheric phase the period of time in which the altitude of the engine is between around 20 and 120 kilometers. It is indeed in this range that, in the absence of any motor thrust, the aerodynamic forces (friction with the atmosphere) can be employed to adequately control the space shuttle to as to steer it to a desired final point and meanwhile satisfying the state constraints in particular on the thermal flux. Thus, during this phase the shuttle can be considered as a glider, only submitted to the gravity force and the aerodynamic forces. The control is the bank angle, and the minimization criterion under consideration is the total thermal flux. The model of the control system is standard (see e.g. [START_REF] Miele | Recent advances in the optimization and guidance of aeroassisted orbital transfers[END_REF][START_REF] Cnes | Mécanique spatiale[END_REF]) and is written as

dr dt = v sin γ dv dt = -g sin γ - 1 2 ρ SC D m v 2 + Ω 2 r cos L(sin γ cos L -cos γ sin L cos χ) dγ dt = cos γ - g v + v r + 1 2 ρ SC L m v cos µ + 2Ω cos L sin χ + Ω 2 r v cos L(cos γ cos L + sin γ sin L cos χ) dL dt = v r cos γ cos χ dl dt = v r cos γ sin χ cos L dχ dt = 1 2 ρ SC L m v cos γ sin µ + v r cos γ tan L sin χ + 2Ω(sin L -tan γ cos L cos χ) + Ω 2 r v sin L cos L sin χ cos γ . (20)
Here, r denotes the distance of the center of gravity of the shuttle to the center of the Earth, v is the modulus of its relative velocity, γ is the flight angle (or path inclination, that is, the angle of the velocity vector with respect to an horizontal plane), L is the latitude, l is the longitude, and χ is the azimuth (angle between the projection of the velocity vector onto the local horizontal plane measured with respect to the axis South-North of the planet).

The gravitational force appears with a usual model g(r) = µ 0 r 2 , where µ 0 is the gravitational constant. The aerodynamic forces consist of the drag force, whose modulus is 1 2 ρSC D v 2 , which is opposite to the velocity vector, and of the lift force, whose modulus is 1 2 ρSC L v 2 , which is perpendicular to the velocity vector. Here, ρ = ρ(r) = ρ 0 e -βr is the air density, S is some positive coefficient (reference area) featuring the engine, and C D and C L are the drag and the lift coefficients; they depend on the angle of attack and on the Mach number of the shuttle. Notice that more specific models can be used and that in general the gravity, the air density and the aerodynamic coefficients are tabulated (we refer to [START_REF] Bonnard | Une approche géométrique du contrôle optimal de l'arc atmosphérique de la navette spatiale[END_REF] for precise tabulations used in the study).

The control is the bank angle µ; it acts on the orientation of the lift force and thus its effect may be to make the shuttle turn left or right but also to act on the altitude. It is a scalar control that is assumed to take values in [0, π]. Note that the mass m of the engine is constant along this atmospheric phase since it is assumed that there is no thrust.

Finally, Ω is the angular rotation speed of the planet. In the above model the terms linear in Ω represent the Coriolis force, and the terms proportional to Ω 2 are due to the centripetal force.

The optimal control problem under consideration is to steer the vehicle from initial conditions to final conditions stated in Table 1, in free final time, and moreover the system is submitted to three state constraints:

• a constraint on the (instantaneous) thermal flux: ϕ = C q √ ρv 3 ϕ max ,

• a constraint on the normal acceleration: • a constraint on the dynamic pressure: 1 2 ρv 2 P max , where C q , ϕ max , γ n 0 , γ max

γ n = γ n 0 ρv 2 γ
J(µ) = t f 0 C q √ ρv 3 dt. ( 21
)
Figure 1: Constraints, and Harpold-Graves strategy Note that, if we approximate v -d, then

J(µ) = K v f v 0 v 2 √
d dv (with K > 0), and hence for this approximated criterion the optimal strategy is to maximize the drag d all along the flight. This strategy, described in [START_REF] Harpold | Shuttle entry guidance[END_REF] and usually employed, reduces the problem to the problem of finding a trajectory tracking the boundary of the authorized domain in the following order: thermal flux -normal acceleration -dynamic pressure, as drawn on figure 1. The advantage of this method is that along the boundary arcs the control can be easily expressed in closed-loop (feedback), which is greatly convenient for stabilization issues and for real-time embarked implementation. Anyway this strategy is not optimal for the minimization criterion [START_REF] Haberkorn | Convergence results for smooth regularizations of hybrid nonlinear optimal control problems[END_REF], and it was the aim of [START_REF] Bonnard | Une approche géométrique du contrôle optimal de l'arc atmosphérique de la navette spatiale[END_REF][START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF][START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations by multiple-shooting techniques[END_REF][START_REF] Trélat | Optimal control of a space shuttle and numerical simulations[END_REF] to solve this optimal control problem with geometric considerations.

A version of the Pontryagin Maximum Principle can be applied to that problem but it is then difficult to make converge the resulting shooting method, due to the fact that the domain of convergence is very small and getting a good initial condition of the adjoint vector is a real challenge. Of course, many numerical refinements can be proposed to overcome this initialization problem, and similar optimal control problems have been considered in a number of articles (see e.g. [START_REF] Betts | Survey of numerical methods for trajectory optimization[END_REF][START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF][START_REF] Pesch | A practical guide to the solution of real-life optimal control problems[END_REF][START_REF] Stengel | Optimal guidance for the space shuttle transition[END_REF][START_REF] Naidu | Three-dimensional atmospheric entry problem using method of matched asymptotic expansions[END_REF][START_REF] Kreim | Minimizing the maximum heating of a re-entering space shuttle: an optimal control problem with multiple control constraints[END_REF] with various approaches (direct or indirect). This is indeed a classical problem, but we insist on the fact that our objective is here to show how a result of geometric optimal control can be of some help in order to guess a good initial condition to make converge the shooting method (rather than making it converge through numerical refinements). Note that, without the aid of such a tool, solving this problem with a shooting method is nearly intractable.

Geometric optimal control results and application to the problem

In this section, instead of providing a solution with computational or numerical refinements, our goal is to provide a rough analysis of the control system and show how geometric control can be of some help in order to provide a better understanding of the structure of the system and finally lead to a precise description of the optimal trajectories, then reducing the application of the shooting method to an easy exercise.

Before that, let us first explain what is geometric control. As explained in the introduction of this article, modern optimal control theory combines classical techniques developed in the 60's, typically the Pontryagin Maximum Principle, with other powerful mathematical techniques in order to provide results on the structure of optimal trajectories for general classes of nonlinear control systems. Literally, geometric optimal control is the combination of classical optimal control with geometric methods in system theory. More precisely, it can be described as the combination of the knowledge inferred from the Pontryagin Maximum Principle with geometric considerations such as the use of Lie brackets, of subanalytic sets, of differential geometry on manifolds, of symplectic geometry and Hamiltonian systems, with the ultimate objective of deriving optimal synthesis results, permitting to describe in a precise way the structure of optimal trajectories. In other words, the objective is to derive results saying that, according to the class of control systems we are considering, the optimal trajectories have a precise structure and are of such and such kind. Geometric optimal control has furnished a modern and uniform framework to realize this objective.

The foundations of geometric control can be dated back, first, to the important Chow's Theorem (see [START_REF] Chow | Uber Systeme von linearen partiellen Differentialgleichungen erster Ordnung (German)[END_REF]) on reachable sets of integral curves of families of vector fields, which was not part of the classical calculus of variations theory, and second, to the articles [START_REF] Brunovsky | Every normal linear system has a regular synthesis[END_REF][START_REF] Brunovsky | Existence of regular synthesis for general problems[END_REF], where Brunovsky discovered that it was possible to derive regular synthesis results using geometric considerations for a large class of control systems, yielding a precise description of the structure of optimal trajectories. Since then, many different tools from differential geometry have been introduced in optimal control, forming gradually a package of techniques and knowledge now identified as geometric optimal control :

• the use of differentiable manifolds, extending the field of applications of the Pontryagin Maximum Principle to optimal control problems naturally posed on a manifold or on a Lie group (see e.g. [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Jurdjevic | Geometric control theory[END_REF]), very much encountered in mechanics, robotics, aerospace, quantum control theory, etc;

• Lie brackets and Lie algebras, used to derive important results on accessibility and controllability properties, to derive higher order optimality conditions (see [START_REF] Sussmann | Controllability of nonlinear systems[END_REF][START_REF] Brockett | Lie algebras and Lie groups in control theory[END_REF][START_REF] Hermes | Lie algebras of vector fields and local approximation of attainable sets[END_REF][START_REF] Krener | A generalization of Chow's Theorem and the bang-bang theorem to nonlinear control problems[END_REF][START_REF] Lobry | Contrôlabilité des systèmes non-linéaires[END_REF]) permitting as mentioned above to restrict the set of optimal candidates, and to derive local regularity and optimal synthesis results (see [START_REF]Schättler The local structure of time optimal trajectories in dimension 3[END_REF][START_REF] Krener | The structure of small-time reachable sets in low dimensions[END_REF][START_REF] Agrachev | On the local structure of optimal trajectories in IR 3[END_REF][START_REF] Sigalotti | Local regularity of optimal trajectories for control problems with general boundary conditions[END_REF]);

• stratification and real analyticity theory considerations (see [START_REF] Sussmann | Lie brackets, real analyticity and geometric control, Differential geometric control theory[END_REF][START_REF] Sussmann | Regular synthesis and sufficient conditions for optimality[END_REF][START_REF] Agrachev | On subanalyticity of Carnot-Carathéodory distances[END_REF]), used as well for regularity and optimal synthesis issues;

• singularity theory, providing a starting point to the classification of extremals (see [START_REF] Ekeland | Discontinuités de champs hamiltoniens et existence des solutions optimales en calcul des variations[END_REF][START_REF] Kupka | Geometric theory of extremals in optimal control problems. I. The fold and Maxwell case[END_REF][START_REF] Kupka | The ubiquity of Fuller's phenomenon[END_REF][START_REF] Bonnard | The role of singular trajectories in control theory[END_REF][START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF][START_REF] Caillau | Minimum time control of the restricted three-body problem[END_REF]) et permitting to study how trajectories may lose optimality (see [START_REF] Agrachev | On sub-Riemannian caustics and wave fronts for contact distributions in the three-space[END_REF][START_REF] Boscain | Morse properties for the minimum time function on 2D manifolds[END_REF][START_REF] Agrachev | Two-dimensional almost-Riemannian structures with tangency points[END_REF]);

• the use of the Miele clock form (see [START_REF] Miele | Extremization of linear integrals by Green's theorem[END_REF][START_REF] Hermes | Functional analysis and time optimal control[END_REF]), widely generalized with the framework of symplectic geometry and methods, the latter being used to provide sufficient conditions for optimality either in terms of conjugate points, Maslov index, envelope theory, extremal field theory, or of optimal synthesis (see [START_REF] Sussmann | Envelopes, conjugate points and optimal bang-bang extremals[END_REF][START_REF] Agrachev | Symplectic methods for optimization and control[END_REF][START_REF] Bonnard | The role of singular trajectories in control theory[END_REF][START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Noble | Sufficient conditions for relative minima of broken extremals in optimal control theory[END_REF][START_REF] Poggiolini | State-local optimality of a bang-bang trajectory: a Hamiltonian apprach[END_REF]);

• fine considerations from differential geometry, e.g. the concepts of conjugate or cut locus, of Jacobi curves or of curvature, used to provide global optimality results (see e.g. [START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF][START_REF] Agrachev | Geometry of Jacobi curves[END_REF][START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF]);

• sub-Riemannian metrics (see e.g. [START_REF] Brockett | Control theory and singular Riemannian geometry[END_REF][START_REF] Bellaïche | Sub-Riemannian geometry[END_REF][START_REF] Montgomery | A tour of subriemannian geometries, their geodesics and applications[END_REF]), much used for applications to robotics and more recently to aerospace problems;

• and many other notions and mathematical concepts, borrowed from differential geometry and related areas.

Typically, one should keep in mind the following idea. The aim of using these geometric tools is to provide a complement to the Pontryagin Maximum Principle whenever its application alone happens to be insufficient to adequately solve an optimal control problem, due to a lack of information. As explained in details in section 2, the PMP is a first-order condition for optimality and its aim is to select a set of trajectories that are candidate to be optimal. Apart from the ultimate goal of providing a complete optimal synthesis, one of the objectives of geometric control is to derive higher order optimality conditions in order to restrict more the set of candidate optimal trajectories. Second-order conditions have been briefly reviewed in section 2.3 and their connection to conjugate point theory has been put in evidence. More generally the objective of higher order conditions, in terms of conjugate points, Maslov index techniques, envelope theory, Miele clock form or more generally symplectic tools, is to select among the extremals derived from the PMP, those who are candidates to be indeed optimal. When this selection is achieved in a so nice way that there exists only one possible way to steer the system from the initial to the final prescribed conditions, one speaks of an optimal synthesis, although this wording underlies some more regularity properties, in particular regularity properties of the selected extremals ensuring that they form an extremal field (see e.g. [START_REF] Brunovsky | Every normal linear system has a regular synthesis[END_REF][START_REF] Brunovsky | Existence of regular synthesis for general problems[END_REF][START_REF] Sussmann | Regular synthesis and sufficient conditions for optimality[END_REF]). We refer the reader interested in a deeper insight on geometric control issues to the survey [START_REF] Sussmann | Shortest paths for the Reeds-Shepp car: a worked out example of the use of geometric techniques in nonlinear optimal control[END_REF] and to the textbooks [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Bonnard | The role of singular trajectories in control theory[END_REF][START_REF] Jurdjevic | Geometric control theory[END_REF][START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF][START_REF] Bullo | Geometric control of mechanical systems. Modeling, analysis, and design for simple mechanical control systems[END_REF][START_REF] Bressan | Introduction to the mathematical theory of control[END_REF][START_REF] Schättler | Geometric optimal control, theory, methods, examples[END_REF] and the many references therein. Note that we do not mention here the many geometric issues related with stabilization, that are outside of the scope of this article.

Let us now show how, starting from a simple remark on the structure of the control system [START_REF] Trélat | Contrôle optimal : théorie & applications, Vuibert, Collection[END_REF], results from geometric control theory can be applied and then help to guess the precise structure of optimal trajectories, and ultimately make the application of a shooting method much easier. In a first step, let us assume that the rotation of the planet can be neglected, that is, Ω = 0. Note that, at the end, this is not the case and the effect of the Coriolis force happens to be necessary to reach the desired final conditions. Anyway, assuming that Ω = 0, the control system (20) gets simpler and in particular the three first differential equations become

dr dt = v sin γ dv dt = -g sin γ - 1 2 ρ SC D m v 2 dγ dt = cos γ - g v + v r + 1 2 ρ SC L m v cos µ. (22)
Setting u = cos(µ), x=(r, v, γ), and defining the vector fields in IR 3

f 0 (x) = v sin γ ∂ ∂r -(g sin γ + 1 2 ρ SC D m v 2 ) ∂ ∂v + cos γ - g v + v r ∂ ∂γ , f 1 (x) = 1 2 ρ SC L m v ∂ ∂γ ,
the system ( 22) can be rewritten as the single-input control-affine system in IR 3

ẋ(t) = f 0 (x(t)) + u(t)f 1 (x(t)), |u| 1. (23)
Ignoring the coordinates (L, l, χ), the induced optimal control problem is to steer the above three-dimensional control system, from a given x(0) to a final target (γ(t f ) is free but the two first coordinates are fixed), with a control satisfying the constraint |u| 1, and moreover, under the three state constraints on the thermal flux, normal acceleration and dynamic pressure (which depend only on x), by minimizing the cost [START_REF] Haberkorn | Convergence results for smooth regularizations of hybrid nonlinear optimal control problems[END_REF].

Reparametrizing by the instantaneous cost (that is, dividing the equations ( 22) by ϕ = C q √ ρv 3 , and setting as a new time s = ϕt), we end up with the minimal time problem for a single-input control-affine system, with the constraint |u| 1 on the control, and with pure state constraints of the form c i (x) 0, i = 1, 2, 3.

Besides, there exist results coming from geometric optimal control theory, providing a qualitative description of minimal time trajectories for control systems of the form [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF], within small time, in small dimension (two and three), and under generic assumptions. We refer the reader to [START_REF] Agrachev | On the local structure of optimal trajectories in IR 3[END_REF][START_REF] Boscain | Optimal syntheses for control systems on 2-D manifolds[END_REF][START_REF] Krener | The structure of small-time reachable sets in low dimensions[END_REF][START_REF] Piccoli | Classification of generic singularities for the planar time-optimal synthesis[END_REF][START_REF]Schättler The local structure of time optimal trajectories in dimension 3[END_REF][START_REF] Sussmann | The structure of time-optimal trajectories for single-input systems in the plane: the C ∞ nonsingular case[END_REF] for precise results. For instance, in dimension three, in the absence of state constraint, it has been derived in [START_REF] Krener | The structure of small-time reachable sets in low dimensions[END_REF] that, if the vector fields of the system are such that f 0 , f 1 and their Lie bracket [f 0 , f 1 ] are linearly independent at x 0 , then minimal time trajectories starting from x 0 are locally bang-bang with at most two switchings. Moreover, denoting by x + (resp., x -) an arc corresponding to the control u = 1 (resp., u = -1), the small time accessible set is homeomorphic to a cone whose boundary consists of all trajectories of the form x + x -and x -x + (i.e., concatenations of two bang arcs). Furthermore, using the Miele clock form (see [START_REF] Miele | Extremization of linear integrals by Green's theorem[END_REF][START_REF] Hermes | Functional analysis and time optimal control[END_REF][START_REF] Moreno | Optimal time control of bioreactors for the wastewater treatment[END_REF][START_REF] Bonnard | Une approche géométrique du contrôle optimal de l'arc atmosphérique de la navette spatiale[END_REF] for the use of this form in the plane, and see [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Bonnard | The role of singular trajectories in control theory[END_REF] for a generalization in a symplectic context), according to the sign of some coefficient only depending on the Lie structure of the vector fields at x 0 , it can be shown that, locally around x 0 , trajectories x + x -x + (starting from x 0 ) are minimal time whereas trajectories x -x + x -are maximal time (or conversely according to that sign).

Motivated by the atmospheric re-entry problem, this kind of result has been extended in [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF] to the case where there are state constraints. Actually in this reference local optimal syntheses are derived for systems in dimension two and three, with one or several state constraints. This classification involves many cases, depending on the order of the state constraints under consideration and cannot be sketched in few words. In the case of the atmospheric re-entry problem, all state constraints are of order two, since one has to differentiate two times with respect to t the relations characterizing a boundary arc to make appear the control. The results of [START_REF] Bonnard | Optimal control with state constraints and the space shuttle re-entry problem[END_REF], combined with numerical simulations, then lead to the following result.

Proposition 1. The optimal trajectory for the simplified three-dimensional model ( 22) is of the kind x -x + x f lux x + x acc x + , where x + (resp., x -) is an arc corresponding to the control u = 1 (resp., u = -1), and x f lux (resp., x acc ) denotes a boundary arc saturating the constraint on the thermal flux (resp. on the normal acceleration).

Note that, since the abovementioned geometric results are of local nature, to make them global they must be combined with numerical simulations and possibly with conjugate point arguments, as it was done for the atmospheric re-entry problem in [START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations by multiple-shooting techniques[END_REF]. In this latter reference, it was also shown, using perturbation arguments, how this result in dimension three could be used in order to provide an approximation of the optimal trajectory for the true problem in dimension six. In this perturbation argument, the parameter Ω is in some sense viewed as a small parameter, but to justify properly the argument it must also be observed that the simplified three-dimensional system is almost a projection onto IR 3 of the complete system in dimension six (we refer to [START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations by multiple-shooting techniques[END_REF] for details). Anyway, what is important is that the strategies announced in Proposition 1 provide a good approximation of the optimal trajectories of the complete problem in dimension six. Now, the point is that it is very easy to make converge a shooting method for the simplified problem in dimension three. Indeed, since one knows precisely the structure of the optimal trajectory, the trajectory to be determined can be parametrized only with its switching times, and hence the shooting problem reduces to a problem with only five unknowns (which are the switching times). The resulting optimal trajectory can then serve as a good initial guess for seeking the optimal trajectory of the complete problem in dimension six. Moreover, it is possible to derive as well a good approximation of the initial adjoint vector in dimension six, by completing the Lagrange multiplier of the optimal solution in dimension three with zeros (it is shown in [START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations by multiple-shooting techniques[END_REF][START_REF] Trélat | Optimal control of a space shuttle and numerical simulations[END_REF] that it is indeed a good approximation because for Ω = 0 the optimal trajectory of the three-dimensional problem can be viewed as a singular trajectory of the sixdimensional problem with corank three).

Open challenges

It has been shown previously how a result of geometric optimal control theory on local optimal syntheses can help to make converge a shooting method, or at least can simplify its implementation by describing precisely the structure of the optimal trajectory (e.g. as a succession of bang, singular or boundary arcs, in a precise order). As briefly surveyed previously, these results exist only for control-affine systems in small dimension (essentially, two and three). Note that, in dimension three, more general situations have been considered in [START_REF] Agrachev | On the local structure of optimal trajectories in IR 3[END_REF] for single-input control-affine systems, providing a precise local structure of optimal trajectories having a finite number of switching times and involving possible singular arcs. These results have been generalized in the deep article [START_REF] Sigalotti | Local regularity of optimal trajectories for control problems with general boundary conditions[END_REF], in which the author studies the local structure of minimal time trajectories for single-input control-affine systems with a starting point in a submanifold S. It is shown that, if the codimension of S is less than or equal to four, then generic minimal time trajectories starting from S are concatenations of at most seven between bang and singular arcs. This result can be applied to four-dimensional systems when S is a point.

For larger dimensions the situation is far more intricate, not only due to the possible occurence of singular trajectories, but also to the generic occurence of Fuller phenomena (see [START_REF] Kupka | The ubiquity of Fuller's phenomenon[END_REF][START_REF] Zelikin | Theory of chattering control[END_REF][START_REF] Robbins | Junction phenomena for optimal control with state-variable inequality constraints of third order[END_REF]) in which case an infinite number of switchings may occur in a compact time interval.

It would be however useful to derive such local optimal synthesis results for systems in larger dimensions, however necessarily under strong assumptions in particular to avoid the Fuller phenomenon, for instance in view of providing alternative ways of making converge the shooting method for the orbit transfer problem. Concerning the latter problem note that, in order to generate the accessible set for the orbit transfer problems, the authors of [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF] have used tools of Riemannian geometry to determine the cut and conjugate loci on a complete two-surface of revolution in order to infer the global structure of the extremals of the problem.

Note that the results of geometric optimal control mentioned in the previous section essentially rely on a careful analysis of the extremal flow using second-order conditions or a Hamiltonian approach and hence are strongly related to the concept of conjugate time (surveyed previously in this paper). Hence, the methods e.g. of [START_REF] Agrachev | Symplectic methods for optimization and control[END_REF][START_REF] Agrachev | Strong optimality for a bang-bang trajectory[END_REF][START_REF] Noble | Sufficient conditions for relative minima of broken extremals in optimal control theory[END_REF][START_REF] Osmolovskii | Equivalence of second order optimality conditions for bang-bang control problems. Part 1: main results[END_REF][START_REF] Poggiolini | State-local optimality of a bang-bang trajectory: a Hamiltonian apprach[END_REF][START_REF] Poggiolini | Sufficient optimality conditions for a bang-singular extremal in the minimum time problem[END_REF][START_REF] Sarychev | First-and second-order sufficient optimality conditions for bang-bang controls[END_REF] that are applicable to bang-bang situations in any dimension should permit to derive local optimal syntheses in larger dimension under additional assumptions, and as well for control-affine systems with more than one control (although it can be expected that the situation is much more complicated). Note however that, according to the results of [START_REF] Chitour | Propriétés génériques des trajectoires singulières[END_REF][START_REF] Chitour | Genericity results for singular curves[END_REF][START_REF] Chitour | Singular trajectories of control-affine systems[END_REF], generic (in the Whitney sense) controlaffine systems do not admit any minimizing singular trajectory whenever the number of controls is more than two (more precisely it is shown in these references that such generic control-affine systems do not admit any trajectories satisfying the Goh necessary condition derived in [START_REF] Goh | Necessary conditions for singular extremals involving multiple control variables[END_REF]). For a first result concerning the classification of extremals for control-affine systems with two controls, we quote the recent article [START_REF] Caillau | Minimum time control of the restricted three-body problem[END_REF], with an application to the minimum time control of the restricted three-body problem. [START_REF] Filippov | On certain questions in the theory of optimal control[END_REF] The continuation method and applications

The continuation method

The objective of continuation or homotopy methods is to solve a problem step by step from a simpler one by parameter deformation. There exists a well-developed theory and many algorithms and numerical methods implementing these ideas, and the field of applications encompasses Brouwer fixed point problems, polynomial and nonlinear systems of equations, boundary value problems in many diverse forms, etc. We refer the reader to the survey articles [START_REF] Allgower | Simplicial and continuation methods for approximating fixed points and solutions to systems of equations[END_REF][START_REF] Allgower | Continuation and path following[END_REF][START_REF] Rheinboldt | Numerical continuation methods: a perspective, Numerical analysis[END_REF][START_REF] Watson | Probability-one homotopies in computational science[END_REF] or to the textbook [START_REF] Allgower | Numerical continuation methods. An introduction[END_REF] for a complete report on these theories and methods.

Here we will use the continuation or homotopy approach in order to solve the shooting problem resulting from the application of the Pontryagin Maximum Principle to an optimal control problem. More precisely, the method consists of deforming the problem into a simpler one that we are able to solve (without any careful initialization of the shooting method), and then of solving a series of shooting problems, step by step, to come back to the original problem. In practice the choice of an adapted parameter deformation of the problem is done according to an intuition or an heuristics with respect to the physical meaning of the different parameters entering into the problem, and thus may require considerable physical insight into the problem. The homotopic parameter λ can be a physical parameter (or several) of the problem, or an artificial one. Some examples are provided in the sequel. The deformation should also be chosen in such a way to enjoy sufficient regularity conditions, making the homotopy method converge. Notice that not only the simpler problem should be chosen according to an heuristics, but also the path between the simpler problem and the original problem. When the homotopic parameter λ is a real number and when the path is linear in λ (meaning that, in some coordinates, the path consists of a convex combination of the simpler and of the original problem, with λ ∈ [0, 1]), the homotopy method is rather called a continuation method in the literature. The continuation method consists then of tracking a set of zeros, as the parameter λ is increased monotonically from 0 to 1 (starting from the simpler known solution). Numerical continuation is well-known in numerical analysis and has been applied to a wide field of various problems. It can fail whenever the path of zeros which is tracked has bifurcation points or more generally singularities, or whenever this path fails to exist globally and does not reach λ = 1. Homotopy methods generalize continuation methods, in the sense that the parameter λ is not necessarily increased monotonically from 0 to 1, dealing with the possible occurence of bifurcations or singularities, and in the sense that the parameter λ is not necessarily a real number but can be considered in more general spaces (it can be a real number, or a vectorial number, or even a parameter evolving in some general Banach space); furthermore, in general homotopy methods the path can be nonlinear and considered in various spaces.

For the moment, for the sake of simplicity we focus on the continuation method and consider a real parameter λ ∈ [0, 1] (we comment further on homotopy methods). Let us provide shortly the basic arguments showing the feasibility of the continuation method. From the theoretical point of view, regularity properties require at least that the optimal solution be continuous, or differentiable, with respect to the parameter λ that is expected to increase monotonically in [0, 1]. This kind of property is usually derived using an implicit function argument, which is encountered in the literature as sensitivity analysis. Let us explain what is the general reasoning of sensitivity analysis, in the simplified framework of section 2.2.1, that is, assuming that M 0 = {x 0 }, M 1 = {x 1 } and U = IR m . We are faced with a family of optimal control problems, parametrized by λ, that can be as in [START_REF] Sarychev | Lipschitzian regularity of minimizers for optimal control problems with control-affine dynamics[END_REF] written in the form of min

E x 0 ,T,λ (u λ )=x 1 C T,λ (u). ( 24
)
According to the Lagrange multipliers rule, if u λ is optimal then there exists (ψ λ , ψ 0 λ ) ∈ IR n ×IR\ {0} such that ψ λ dE x 0 ,T,λ (u λ ) + ψ 0 λ dC T,λ (u) = 0. Assume that there are no minimizing abnormal extremals in the problem. Under this assumption, since the Lagrange multiplier (ψ λ , ψ 0 λ ) is defined up to a multiplicative scalar we can definitely assume that ψ 0 λ = -1. Then, we are seeking (u λ , ψ λ ) such that F (λ, u λ , ψ λ ) = 0, where the function F is defined by

F (λ, u, ψ) = ψdE x 0 ,T,λ (u) -dC T,λ (u) E x 0 ,T,λ (u) -x 1 = ∂L T,λ ∂u (u, ψ) E x 0 ,T,λ (u) -x 1 ,
where L T,λ (u, ψ) = ψE x 0 ,T,λ (u) -C T,λ (u) is the Lagrangian, as defined in section 2.2.1. Let ( λ, uλ, ψλ) be a zero of F . Assume that F is of class C 1 . If the Jacobian of F with respect to (u, ψ), taken at the point ( λ, uλ, ψλ), is invertible, then according to a usual implicit function argument one can solve the equation F (λ, u λ , ψ λ ) = 0, and the solution (u λ , ψ λ ) depends in a C 1 way on the parameter λ. Note that this standard argument from sensitivity analysis is at the base of the well-known Lagrange-Newton method in optimization.

Let us now analyze the invertibility condition of the Jacobian of F with respect to (u, ψ). This Jacobian matrix is

Q T,λ dE x 0 ,T,λ (u) * dE x 0 ,T,λ (u) 0 , (25)
where Q T,λ is the quadratic form considered in section 2.3.1, defined by

Q T,λ = ∂ 2 L T,λ ∂ 2 u (u, ψ, ψ 0 ) | ker ∂L T,λ ∂u
, and dE x 0 ,T,λ (u) * is the transpose of dE x 0 ,T,λ (u). The matrix [START_REF] Trélat | Some properties of the value function and its level sets for affine control systems with quadratic cost[END_REF] (which is a matrix of operators) is called sensitivity matrix in sensitivity analysis. It is an easy exercise to prove that this sensitivity matrix is invertible if and only if the linear mapping dE x 0 ,T,λ (u) is surjective and the quadratic form Q T is nondegenerate. Having in mind the definitions given previously in this article, the meaning of these assumptions is the following. The surjectivity of dE x 0 ,T,λ (u) exactly means that the control u is not singular (see Definition 2). The nondegeneracy of Q T,λ is exactly related with the concept of conjugate point (see Definition 4). Note that, as long as we do not encounter any conjugate time along the continuation path, the extremals that are computed are locally optimal. It follows that, to ensure the surjectivity of dE x 0 ,T,λ (u) along the continuation process, it suffices to assume the absence of singular minimizing trajectory. Note that, in the simplified problem that we considered, where the controls are unconstrained, singular trajectories are exactly the projections of abnormal extremals. Therefore, we conclude that, as long as we do not encounter any minimizing singular control nor conjugate point along the continuation procedure, then the continuation method works locally and the extremal solution (u λ , ψ λ ) which is locally computed as above is of class C 1 with respect to the parameter λ. These two assumptions are the basic ones ensuring the existence of a local solution in the continuation procedure. In other words, under these two sufficient conditions, the continuation method is locally feasible.

Before going to global considerations, let us make an ultimate comment on these two assumptions. The absence of conjugate point can be tested numerically: as explained in section 2.3.2, it suffices to test the vanishing of some determinant along the extremal flow (see [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] for details, and see section 2.4.2 where the package COTCOT is cited). As long as this test does not detect any conjugate point along the continuation process, this means that the extremals that are computed are locally optimal. The assumption of the absence of minimizing singular trajectories is of a much more geometric nature. Such results exist for some classes of controlaffine systems under some strong Lie bracket assumptions (see [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Bonnard | The role of singular trajectories in control theory[END_REF][START_REF] Trélat | Global subanalytic solutions of Hamilton-Jacobi type equations[END_REF][START_REF] Trélat | Singular trajectories and subanalyticity in optimal control and Hamilton-Jacobi theory[END_REF]). Moreover, as mentioned in section 3.3, it is proved in [START_REF] Chitour | Propriétés génériques des trajectoires singulières[END_REF][START_REF] Chitour | Genericity results for singular curves[END_REF][START_REF] Chitour | Singular trajectories of control-affine systems[END_REF] that for generic (in the Whitney sense) control-affine systems with more than two controls, there is no minimizing singular trajectory; hence for such kinds of systems the assumption of the absence of minimizing singular trajectory is automatically satisfied.

Remark 11. The implicit function argument given above is on the control, but the continuation method is usually implemented on the exponential mapping (defined in Definition (3)) and consists of tracking a path of initial adjoint vectors doing the job. More precisely, instead of [START_REF] Rifford | Morse-Sard type results in sub-Riemannian geometry[END_REF], one has to solve exp x 0 ,λ (T, p 0,λ

) = x 1 , (26)
where the exponential mapping is parametrized by λ. Note that this is the shooting method in the simplified case. Of course for the method to be locally feasible as above we have the same sufficient conditions.

The previous implicit function arguments permit to ensure the local feasibility of the continuation procedure, locally around a given solution (that is, locally around a given parameter λ). Now to make it global over [0, 1], we ought to ensure that the path of zeros λ → p 0,λ is globally defined on [0, 1] and joins p 0,0 to p 0,1 . It could indeed happen that the path is not globally defined and either reaches some singularity or wanders off to infinity before reaching λ = 1. To eliminate the first possibility, since a limit of optimal controls is optimal as well (see e.g. [START_REF] Trélat | Some properties of the value function and its level sets for affine control systems with quadratic cost[END_REF][START_REF] Haberkorn | Convergence results for smooth regularizations of hybrid nonlinear optimal control problems[END_REF]), we can make the assumption of the absence of minimizing singular trajectory and of conjugate point over all the domain under consideration (not only along the continuation path), and for every λ ∈ [0, 1]. As said before, the absence of singular minimizing trajectory over the whole space is generic for large classes of systems, hence this is a reasonable assumption; however the global absence of conjugate point is a strong assumption. There exist however some other possibilities to tackle singularities2 . To eliminate the second possibility, we ought to provide sufficient conditions ensuring that the tracked paths remain bounded. In other words, considering [START_REF] Leitmann | The calculus of variations and optimal control. An introduction[END_REF], we have to ensure that the initial adjoint vectors p 0,λ that are computed along the continuation procedure remain bounded, uniformly with respect to the homotopic parameter λ. This means that we have to ensure that the exponential mapping is proper, uniformly with respect to λ. The properness of the exponential mapping has been studied in [START_REF] Trélat | Some properties of the value function and its level sets for affine control systems with quadratic cost[END_REF], where it has been proved that, if the exponential mapping is not proper, then there must exist an abnormal minimizer (see also [START_REF] Bonnard | On the role of abnormal minimizers in sub-Riemannian geometry[END_REF], and see [START_REF] Haberkorn | Convergence results for smooth regularizations of hybrid nonlinear optimal control problems[END_REF]Lemma 2.16] for a more general statement). By contraposition, if one assumes the absence of minimizing abnormal extremals, then the required boundedness follows.

Note again that, in the simplified problem that we considered, where the controls are unconstrained, singular trajectories are exactly the projections of abnormal extremals. Hence we have obtained the following result. Proposition 2. In the simplified case where M 0 = {x 0 }, M 1 = {x 1 } and U = IR m , if there is no minimizing singular trajectory nor conjugate points over all the domain, for every λ ∈ [0, 1], then the continuation procedure (26) is globally feasible on [0, 1].

We stress again that, as explained before, the assumption of the absence of minimizing singular trajectory made in this proposition holds generically for large classes of control systems. Therefore this is a reasonable assumption, contrarily to the assumption of the absence of conjugate point over all the domain, which is a too strong assumption. It is however not easy to weaken this assumption. One possibility is to consider a smaller domain, covering the continuation paths under consideration; in this case however one ought to ensure that the tracked continuation path remains in the domain, that is, remains far from its boundary. There does not seem to exist any simple and tractable condition ensuring this fact in general. Note that in [START_REF] Bonnard | The smooth continuation method in optimal control with an application to quantum systems ESAIM Control Optim[END_REF] the authors use the concept of injectivity radius in order to provide estimations of domains in which the continuation method is globally feasible, on an example which is however specific in Riemannian geometry. This simple result stated in Proposition 2 withstands many possible generalizations. For more general optimal control problems Proposition 2 can be extended quite easily, by adapting the above arguments, and in particular the implicit function argument (although this may be a bit technical, for instance whenever there are some state constraints, see [START_REF] Bonnans | Stability and sensitivity analysis for optimal control problems with a first-order state constraint and application to continuation methods[END_REF]).

In any case, this kind of result provides the mathematical foundations ensuring the global feasibility of the continuation method in optimal control. It can be noted that the feasibility of the continuation method has been much studied for other less specific issues in numerical analysis (see e.g. [START_REF] Ficken | The continuation method for functional equations[END_REF][START_REF] Avila | The feasibility of continuation methods for nonlinear equations[END_REF][START_REF] Allgower | Simplicial and continuation methods for approximating fixed points and solutions to systems of equations[END_REF][START_REF] Allgower | Continuation and path following[END_REF][START_REF] Rheinboldt | Numerical continuation methods: a perspective, Numerical analysis[END_REF] and references therein).

In the more general case of homotopies, the parameter λ is not necessarily increasing monotonically from 0 to 1 and we can encounter turning points (see [START_REF] Watson | Probability-one homotopies in computational science[END_REF]). One of the methods, known as differential homotopy (or differential pathfollowing), consists of tracking a path of zeros s → (λ(s), p 0,λ(s) ) satisfying [START_REF] Leitmann | The calculus of variations and optimal control. An introduction[END_REF] for every s. It is then usual to assume that the mapping F has maximal rank (more precisely, that 0 is a regular value of F ) so that the path of zeros evolves on a submanifold (see e.g. [START_REF] Allgower | Numerical continuation methods. An introduction[END_REF] for the details): this kind of implicit function argument permits to establish, as before, the local feasibility of the method; but now the difference is that turning points3 are allowed: the zero p 0,λ is not necessarily a local function of λ. The global feasibility issues require topological considerations such as connectedness features. Note that, if one does not make this assumption that the mapping F has maximal rank, then one is faced with the possible occurence of singularities. As explained previously for the continuation method, assuming the absence of singularities is a too strong assumption in general. In the existing literature there are essentially two approaches to tackle this difficulty. The first approach, of local type, consists of detecting the possible singularities or bifurcations along the zero path. There is a huge literature on this problem and we refer to [START_REF] Allgower | Numerical continuation methods. An introduction[END_REF] for a survey on these methods applied to homotopy procedures. The second approach, of global type, consists of considering a global perturbation of the homotopy function, more precisely, of the simpler problem under consideration, in order to ensure that, with probability one, 0 is a regular value of F . This variant of the method that can be proved to be globally convergent is known as globally convergent probability-one homotopy method. It is based on nontrivial transversality arguments, combined in [START_REF] Chow | Finding zeros of maps: homotopy methods that are constructive with probability one[END_REF] with Sard's Theorem4 and yielding to homotopy methods with a guarantee of success of probability one with respect to the choice of the simpler problem (see [START_REF] Watson | Probability-one homotopies in computational science[END_REF] for a nice survey discussion and the statement of a general result of global convergence). The "almost everywhere" statement of such a result is used to avoid the possible singularities of the curves to be tracked in the homotopy procedure. The last crucial requirement to ensure global feasibility is as before that the tracked paths remain bounded, in order to ensure that the zero paths are globally well-defined and do not wander off to infinity. This properness can be handled as before by assuming the absence of abnormal minimizers (see arguments of the previous discussion). Having in mind these issues, it is then possible to derive results similar to Proposition 2, according to the specific homotopy method under consideration.

Note that the mathematical foundations of the differential homotopy method applied to optimal control are studied in [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF], and more deeply in [START_REF] Caillau | Minimum time control of the restricted three-body problem[END_REF] where the relation between turning points of the path and conjugate points is clearly elucidated. Note that the authors of [START_REF] Caillau | Minimum time control of the restricted three-body problem[END_REF], studying by homotopy a three-body problem, recommend to stop following a path in case a conjugate point (resulting into a turning point) appears, and provide some hints to jump to another path (these hints are however specific to their problem).

From the numerical point of view there exist many methods and strategies in order to implement continuation or homotopy methods (see [START_REF] Deuflhard | A stepsize control for continuation methods and its special application to multiple shooting techniques[END_REF]), and one has to distinguish between differential pathfollowing (see e.g. [START_REF] Haberkorn | Transfert orbital à poussée faible avec minimisation de la consommation: résolution par homotopie différentielle[END_REF] for applications to orbit transfer problems), simplicial methods (see e.g. [START_REF] Martinon | Numerical resolution of optimal control problems by a piecewise linear continuation method[END_REF] for similar applications), predictor-corrector methods, piecewise-linear methods, etc. Extensive documentation about path following methods with theoretical and algorithmic issues can be found in [START_REF] Allgower | Numerical continuation methods. An introduction[END_REF]. Also, many codes can be found on the web, such as the well-known Hompack90 (see [START_REF] Watson | Algorithm777: HOMPACK90: a suite of Fortran 90 codes for globally convergent homotopy algorithms[END_REF]) or the recent Hampath (see [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF]), just to name a few. Many others are listed in [START_REF] Allgower | Numerical continuation methods. An introduction[END_REF] or can be easily found on internet libraries.

Application to the orbit transfer problem with low thrust

Let us turn back to the orbit transfer problem (1), [START_REF] Bonnard | Mécanique céleste et contrôle de systèmes spatiaux[END_REF], mentioned in the introduction as a motivating example. We show here how the minimal time problem of steering this control system from any initial position to some final orbit can be solved by combining a shooting method with a continuation procedure. On this problem one immediately realizes that the main difficulty is the fact that the maximal authorized modulus of thrust is very low. It is then not surprising to observe numerically that the lower is the maximal thrust, the smallest is the domain of convergence of the Newton method in the shooting problem. In these conditions it is natural to carry out a continuation on the value of the maximal thrust, starting with larger values of the maximal thrust (for which the problem is no more realistic, but for which the shooting method is by far easier to make converge), and then decreasing step by step the value of the maximal thrust, down to low realistic values.

This strategy was implemented in [START_REF] Caillau | 3D geosynchronous transfer of a satellite: continuation on the thrust[END_REF] in order to realize the minimal time 3D transfer of a satellite from a low and eccentric inclinated initial orbit towards the geostationary orbit, for an engine of around 1500 kg. Their continuation procedure starts with the orbit transfer problem with the value T max = 60 N, for which the domain of convergence of the shooting function is large enough so that the shooting method can be initialized easily. Then they decrease the value of T max step by step in order to reach down the value T max = 0.14 N. Along this continuation procedure, the authors observe (and prove) that the minimal time t f is right-continuous with respect to the maximal thrust T max , hence, in theory, it could be expected that the minimal time t f obtained at the step k of the continuation procedure is a good initial guess for the step k + 1. However, they note that this strategy is not so much efficient numerically for low thrusts, in the sense that, for low values of T max , the value of T max has to be decreased with very small steps to ensure convergence. The authors then use the remarkable heuristics t f T max Cst, which permits to significantly improve the efficiency of their continuation procedure and to reach down the low value of T max = 0.14 (for which the resulting time of transfer is more than six months).

The minimal fuel consumption orbit transfer problem has also been solved in [START_REF] Gergaud | Homotopy method for minimum consumption orbit transfer problem[END_REF][START_REF] Gergaud | Low thrust minimum fuel orbital transfer: an homotopic approach[END_REF]. It consists of minimizing the cost t f 0 T (t) dt, and the problem is more difficult than the minimal time problem, because the optimal control derived from the Pontryagin Maximum Principle is no more continuous. This lack of continuity implies difficulties to apply the shooting method. To overcome this problem, the authors propose to implement a continuation on the cost functional, parametrized by λ ∈ [0, 1]. More precisely, they propose to minimize the cost

t f 0 (1 -λ) T (t) 2 + λ T (t) dt.
The case λ = 0 corresponds to the minimization of the energy, while λ = 1 corresponds to the original problem (minimization of the consumption). For every λ < 1, the application of the Pontryagin Maximum Principle leads to smooth controls, for which the shooting method can be applied successfully. Also, for λ = 0 the shooting problem is easier to intialize. The authors prove that it is possible to follow a path of solutions starting from λ = 0 and reaching a value of λ very close to 1, which permits then to initialize successfully the shooting method with λ = 1.

It can be noted that the heuristics t f T max Cst has been understood and clearly explained in the papers [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF][START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF]. In these references, based on preliminary results of [START_REF] Epenoy | Optimal low-thrust transfers with constraints: generalization of averaging techniques[END_REF] where the optimal trajectories of the energy minimization problem are approximated using averaging techniques, the averaged Hamiltonian system is explicitly computed and is shown to be a Riemannian problem. The geodesics and their integrability properties are investigated and deeply analyzed, as well as the Riemannian metrics of the averaged system. Since the averaged system is Riemannian, this means, roughly speaking, that optimal trajectories are straight lines up to a change of coordinates. Since the averaged system can serve as a good approximation of the initial system for low values of the maximal thrust (this fact is proved in these references), the heuristics follows. This is one more example where a geometric insight provides a good understanding of the problem, leading to an efficient numerical solving.

Remark 12. In [START_REF] Haberkorn | Convergence results for smooth regularizations of hybrid nonlinear optimal control problems[END_REF] it is shown how one can moreover take into account a shadow cone (eclipse) constraint in the orbit transfer problem. The approach is based on an hybridization of the problem, considering that the controlled vector fields are zero when crossing the shadow cone. A regularization procedure consisting of smoothing the system, combined with a continuation, is also implemented (it is actually the objective of the article to derive convergence properties of smoothing procedures) and compared with other penalization methods as in [START_REF] Ch | Optimal control for engines with electro-ionic propulsion under constraint of eclipse[END_REF].

A continuation approach to the strong thrust orbit transfer problem by flattening the Earth

In this section we sketch shortly an alternative approach to the strong thrust minimal consumption orbit transfer planification problem developed in [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF], consisting of considering at first the problem for a flat model of the Earth with constant gravity, and then of introducing step by step (by continuation) the variable gravity and the curvature of the Earth, in order to end up with the true model. Of course, the fuel efficient orbit transfer of a satellite has been widely studied (see [START_REF] Marec | Optimal space trajectories[END_REF][START_REF] Bonnard | Mécanique céleste et contrôle de systèmes spatiaux[END_REF]), with many possible approaches. For instance in [START_REF] Naidu | Fuel-optimal trajectories for aeroassisted coplanar orbital transfer problem[END_REF][START_REF] Neustadt | A general theory of minimum-fuel space trajectories[END_REF][START_REF] Vinh | Optimal time-free nodal transfers between elliptical orbits[END_REF] it is solved using impulsive orbit transfers in which it is considered that the engine makes instantaneous changes of velocity. In [START_REF] Leitmann | On a class of variational problems in rocket flight[END_REF] a continuous dynamical approach is considered, taking into account thrust limitations, and then one must distinguish between low-thrust (as in the previous section) and high-thrust (as here) problems. Many possible direct methods are described in [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF] to solve this problem. Indirect methods have also been used in spite of the drawback of the initialization difficulty explained previously. For instance in [START_REF] Augros | Computation of optimal coplanar orbit transfers[END_REF] the impulse transfer solution is used to provide a good initial guess to the shooting method for nearly circular initial and final orbits, based on the intuition that the constrained thrust trajectory should be close to impulsive for strong thrusts (see [START_REF] Neustadt | A general theory of minimum-fuel space trajectories[END_REF][START_REF] Gergaud | Orbital transfer: some links between the low-thrust and the impulse cases[END_REF]). A multiple shooting method is used in [START_REF] Oberle | Existence and multiple solutions of the minimum-fuel orbit transfer problem[END_REF], where the number of thrust arcs evolves iteratively. To initialize properly the shooting method, the authors of [START_REF] Gergaud | Low thrust minimum fuel orbital transfer: an homotopic approach[END_REF][START_REF] Martinon | Using switching detection and variational equations for the shooting method[END_REF] implement differential or simplicial homotopies for passing continuously from the easier problem of minimizing the L 2 -norm of the control to the initial fuel consumption problem, however for the low-thrust problem (the approach does not work for strong thrusts). In [START_REF] Dargent | An integrated tool for low thrust optimal control orbit transfers in interplanetary trajectories[END_REF] an impressive and efficient software for solving the orbit transfer problem is described, combining and smoothing processes and taking advantage of numerical interpolation formulas derived in [START_REF] Thorne | Minimum-time continuous thrust orbit transfers[END_REF] whose goal were to initialize successfully the shooting method within certain ranges of orbit values.

Here we described shortly an alternative approach, based on the remark that the problem is extremely easy to solve whenever the Earth is flat with a constant gravity. Then we pass continuously to the initial model. We restrict to the two-dimensional case and consider the coplanar orbit transfer problem with a spherical Earth and a central gravitational field g(r) = µ r 2 . Written in cylindrical coordinates, the control system under consideration is

ṙ(t) = v(t) sin γ(t), φ(t) = v(t) r(t) cos γ(t), v(t) = -g(r(t)) sin γ(t) + Tmax m(t) u 1 (t), γ(t) = v(t) r(t) -g(r(t)) v(t) cos γ(t) + Tmax m(t)v(t) u 2 (t), ṁ(t) = -βT max u 1 (t) 2 + u 2 (t) 2
(27) where the thrust is T (t) = u(t)T max (here in the application T max is large since we consider a strong thrust), and the control is u(t) = (u 1 (t), u 2 (t)) satisfying the constraint u 1 (t) 2 +u 2 (t) 2 1. The optimal control problem under consideration consists of steering the above system from a given initial configuration r(0

) = r 0 , ϕ(0) = ϕ 0 , v(0) = v 0 , γ(0) = γ 0 , m(0) = m 0 , to some point of a specified final orbit r(t f ) = r f , v(t f ) = v f , γ(t f ) = γ f (
other kinds of final conditions have been considered in [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF]), by maximizing the final mass m(t f ). Note that the final time t f must be fixed in this problem, otherwise the optimal control problem would not have any solution (see [START_REF] Gergaud | Homotopy method for minimum consumption orbit transfer problem[END_REF][START_REF] Oberle | Existence and multiple solutions of the minimum-fuel orbit transfer problem[END_REF]) since it is always better in terms of consumption to let the engine turn more around the planet with shorter thrust arcs. The application of the Pontryagin Maximum Principle to this problem leads to a shooting problem with discontinuous controls (consisting of thrust and balistic arcs) that is not easy to solve directly because it is difficult to initialize adequately. In contrast, consider the very simple flat Earth model

ẋ(t) = v x (t), ḣ(t) = v h (t), vx (t) = Tmax m(t) u x (t), vh (t) = Tmax m(t) u h (t) -g 0 , ṁ(t) = -βT max u x (t) 2 + u h (t) 2 , ( 28
)
where x denotes the horizontal variable, h is the altitude, and v x and v h are the corresponding components of the velocity. The control (u x (•), u h (•)) must satisfy the constraint u x (•) 2 +u h (•) 2 1. It happens that the problem of maximizing the final mass max m(t f ) (here, it makes sense to consider a free final time), with initial conditions x(0) = x 0 , h(0

) = h 0 , v x (0) = v x0 , v h (0) = v h0 , m(0) = m 0 , and final conditions h(t f ) = h f , v x (t f ) = v xf , v h (t f ) = 0,
is extremely simple to solve. It can even be solved explicitly, analytically, and the shooting method can be simplified in order to converge automatically and instantaneously, without any careful initialization (see [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF] for details). In view of that it is tempting to try to pass continuously from this simple model to the initial one by acting on the gravity and on the curvature of the planet. Note that, since the coordinates used in [START_REF] Dmitruk | Maximum Principle for optimal control problems with intermediate constraints[END_REF] are Cartesian whereas those in [START_REF] Arutyunov | Necessary optimality conditions for optimal control problems with intermediate constraints[END_REF] are polar, at the end of the continuation procedure a change of coordinates will be required. Evidently, this change of coordinates is x = rϕ, h = r -r T (where r T is the radius of the Earth), v x = v cos γ, v h = v sin γ, and for the control, u x = u 1 cos γ -u 2 sin γ, u h = u 1 sin γ + u 2 cos γ. When passing from polar to Cartesian coordinates, note however that we not take into account an obvious physical feature: in the absence of control (u = 0), in the flat Earth model [START_REF] Dmitruk | Maximum Principle for optimal control problems with intermediate constraints[END_REF] there do not exist any horizontal trajectories (for which h(t) is constant), whereas the round Earth model [START_REF] Arutyunov | Necessary optimality conditions for optimal control problems with intermediate constraints[END_REF] does admit round (Keplerian) orbits (for which r(t) is constant). This still holds even though we transform the flat Earth model with a variable gravity. This is of course due to the model that is too much simplist, and we are going to modify this model accordingly, by introducing some new terms into the dynamics of the flat Earth model, so that there may exist such horizontal trajectories with null thrust.

First of all, let us apply the above change of coordinates to the control system [START_REF] Arutyunov | Necessary optimality conditions for optimal control problems with intermediate constraints[END_REF]. This leads to

ẋ(t) = v x (t) + v h (t) x(t) r T +h(t) , ḣ(t) = v h (t), vx (t) = Tmax m(t) u x (t) -vx(t)v h (t) r T +h(t) , vh (t) = Tmax m(t) u h (t) - µ (r T +h(t)) 2 + vx(t) 2 r T +h(t) , ṁ(t) = -βT max u x (t) 2 + u h (t) 2 . ( 29
)
This control system is exactly the system [START_REF] Arutyunov | Necessary optimality conditions for optimal control problems with intermediate constraints[END_REF], expressed in cylindrical coordinates. With respect to the flat Earth model [START_REF] Dmitruk | Maximum Principle for optimal control problems with intermediate constraints[END_REF], except the fact that the gravity term is variable, we observe the presence of additional terms in the dynamics of x, v x and v h , which can be viewed for the flat Earth model as kinds of correcting terms that permit the possible occurence of horizontal trajectories. In view of that, in order to pass continuously from the flat Earth model [START_REF] Dmitruk | Maximum Principle for optimal control problems with intermediate constraints[END_REF] to the (actually round Earth) model ( 29), we introduce two parameters λ 1 and λ 2 , the first of which is acting on the gravity, and the second of which permits to introduce the correcting terms. Finally, we consider the family of control systems

ẋ(t) = v x (t) + λ 2 v h (t) x(t) r T +h(t) , ḣ(t) = v h (t), vx (t) = Tmax m(t) u x (t) -λ 2 vx(t)v h (t) r T +h(t) , vh (t) = Tmax m(t) u h (t) - µ (r T +λ 1 h(t)) 2 + λ 2 vx(t) 2 r T +h(t) , ṁ(t) = -βT max u x (t) 2 + u h (t) 2 , (30)
parameterized by 0 λ 1 1 and 0 λ 2 1. Then, we implement the following continuation procedure on the resulting family of optimal control problems: implementing a continuation on λ 1 , keeping λ 2 = 0, we first pass from the simplified flat Earth model [START_REF] Dmitruk | Maximum Principle for optimal control problems with intermediate constraints[END_REF] with constant gravity (for λ 1 = λ 2 = 0) to the intermediate flat Earth model with variable gravity (for λ 1 = 1 and λ 2 = 0). Along this first continuation it makes sense to consider free final times. Then, we implement a second continuation on the parameter λ 2 , keeping λ 1 = 1, to pass continuously to the initial model (for λ 1 = λ 2 = 1). Along this second continuation, we fix the final time for the optimal control problems under consideration to the value obtained at the end of the first continuation.

The details of the procedure as well as numerical simulations are provided in [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF], and comparisons are led with classical direct methods.

To end this section, it remains to explain how the change of coordinates acts onto the adjoint vector, in order to come back to the initial cylindrical coordinates after the continuation procedure. Denoting by F the change of variables from Cartesian to cylindrical coordinates, one passes from the adjoint vector in Cartesian coordinates to cylindrical coordinates by applying the transpose of the inverse of the differential of F . This is indeed a general geometric result whose proof is provided in the appendix of [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF].

Solving the atmospheric re-entry problem by continuation

Another approach to solve the atmospheric re-entry problem of Section 3.1 by a shooting method, implemented in [START_REF] Hermant | Optimal control of the atmospheric reentry of a space shuttle by an homotopy method[END_REF], consists of carrying out a continuation on the maximal value of the state constraint on the thermal flux in order to introduce this constraint step by step. The procedure automatically determines the structure of the optimal trajectory, and permits to start from the easier problem without state constraint and to introduce the constraints progressively (see also [START_REF] Berkmann | Abort landing in windshear: optimal control problem with third-order state constraint and varied switching structure[END_REF]). The theoretical foundations which permit to take into account the change of the structure of the trajectory (and hence the number of unknowns in the shooting method) along the continuation procedure were derived in [START_REF] Bonnans | Stability and sensitivity analysis for optimal control problems with a first-order state constraint and application to continuation methods[END_REF] for first-order state constraints and in [START_REF] Hermant | Homotopy algorithm for optimal control problems with a second-order state constraint[END_REF] for second-order state constraints, and permit to prove that, under some appropriate assumptions, the change in the structure of the trajectory is regular is the sense that, when a constraint becomes active along the continuation, only one boundary arc appears. Note indeed that it could happen that infinitely many boundary arcs appear; see for instance [START_REF] Robbins | Junction phenomena for optimal control with state-variable inequality constraints of third order[END_REF] where this phenomenon is shown to be typical for constraints of order more than or equal to three. Here however in the problem under consideration the state constraints are of order two. To take into account this change of structure along the continuation, the usual continuation procedure must be modified accordingly. For the atmospheric re-entry problem with a constraint on the thermal flux, this procedure is described in details in [START_REF] Hermant | Optimal control of the atmospheric reentry of a space shuttle by an homotopy method[END_REF], and permits to recover in a nice way the results of [START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations by multiple-shooting techniques[END_REF].

General Goddard's problem and singular trajectories

Variants of Goddard's problems (see [START_REF] Goddard | A method of reaching extreme altitudes[END_REF]) are investigated in [START_REF] Bonnans | Singular arcs in the generalized Goddard's problem[END_REF][START_REF] Bonnans | Numerical study of optimal trajectories with singular arcs for an Ariane 5 launcher[END_REF] for nonvertical trajectories. The control is the thrust force, and the objective is to maximize a certain final cost, typically, the final mass. Performing an analysis based on the Pontryagin Maximum Principle, it is proved that optimal trajectories may involve singular arcs (along which the norm of the thrust is neither zero nor maximal), that are computed and characterized. Numerical simulations are carried out, both with direct and indirect methods, demonstrating the relevance of taking into account singular arcs in the control strategy. The indirect method combines a shooting method with a continuation method. The continuation approach leads to a quadratic regularization of the problem similar to the one presented in section 4.2 and is a way to tackle with the problem of nonsmoothness of the optimal control. Note that this quadratic regularization has also been used in [START_REF] Gergaud | Homotopy method for minimum consumption orbit transfer problem[END_REF][START_REF] Martinon | Using switching detection and variational equations for the shooting method[END_REF]. To tackle the lack of continuity of the optimal control u, which makes difficult the application of a shooting method, the authors consider a family of optimal control problems indexed by a continuation parameter λ ∈ [0, 1] with minimization criterion

t f 0 ((1 -λ) u(t) 2 + λ u(t)
)dt, so that the case λ = 0 corresponds to the minimization of the energy and λ = 1 to the original problem (minimization of the consumption) under consideration in their articles.

Dynamical systems theory and applications to mission design 5.1 Dynamics around Lagrange points

Consider the so-called circular restricted three-body problem, in which a body with negligible mass evolves in the gravitational field of two masses m 1 and m 2 called primaries and assumed to have circular coplanar orbits with the same period around their center of mass. The gravitational forces exerted by any other planet or body are neglected. In the solar system this problem provides a good approximation for studying a large class of problems. In a rotating frame the equations are of the form

ẍ -2 ẏ = ∂Φ ∂x , ÿ + 2 ẋ = ∂Φ ∂y , z = ∂Φ ∂z with Φ(x, y, z) = x 2 +y 2 2 + (1 -µ)((x + µ) 2 + y 2 + z 2 ) -1/2 + µ((x -1 + µ) 2 + y 2 + z 2 ) -1/2 + µ(1-µ)
2 . These equations have the first integral (called Jacobi first integral) J = 2Φ-( ẋ2 + ẏ2 + ż2 ), hence the solutions evolve on a five-dimensional energy manifold, the topology of which determines the so-called Hill's region of possible motions (see e.g. [START_REF] Koon | Dynamical Systems, the three-body problem and space mission design[END_REF]).

It is well-known that the above dynamics admit five equilibrium points called Lagrange points, the three first of which, denoted L 1 , L 2 and L 3 , being collinear points on the axis joining the centers of the two primaries, and the two last of which, denoted L 4 and L 5 , located in an equilateral way with respect to the primaries (see [START_REF] Szebehely | Theory of orbits: the restricted problem of three bodies[END_REF]). It must be noted that the linearized system around these equilibrium points admits eigenvalues with zero real part, hence the study of their stability is not obvious. It follows from a generalization of a theorem of Lyapunov (due to Moser [START_REF] Moser | On the generalization of a theorem of A. Lyapunov[END_REF]) that, for a value of the Jacobi integral a bit less than the one of the Lagrange points, the solutions have the same qualitative behavior as the solutions of the linearized system around the Lagrange points. It was then established in [START_REF] Meyer | Introduction to Hamiltonian dynamical systems and the N-body problem Applied Math[END_REF] that the three collinear Lagrange points are always unstable, whereas L 4 and L 5 are stable under some conditions (that are satisfied in the solar system for instance for the Earth-Moon system, or for the system formed by the Sun and any other planet).

The dynamics around these Lagrange points have particularly interesting features for space mission design. Using Lyapunov-Poincaré's Theorem, it is shown that there exists a twoparameter family of periodic trajectories around every Lagrange point (see [START_REF] Meyer | Introduction to Hamiltonian dynamical systems and the N-body problem Applied Math[END_REF], see also [START_REF] Bonnard | Mécanique céleste et contrôle de systèmes spatiaux[END_REF]), among which the well-known halo orbits are periodic orbits that are diffeomorphic to circles (see [START_REF] Breakwell | The halo family of 3-dimensional of periodic orbits in the Earth-Moon restricted 3-body problem[END_REF]) whose interest for mission design was put in evidence by Farquhar (see [START_REF] Farquhar | Station-keeping in the vicinity of collinear libration points with an application to a Lunar communications problem[END_REF][START_REF] Farquhar | A halo-orbit lunar station[END_REF]). There exist many other families of periodic orbits (called Lissajous orbits) and quasi-periodic orbits around Lagrange points (see [START_REF] Gómez | Lissajous orbits around halo orbits[END_REF][START_REF] Gómez | Quasihalo orbits associated with libration points[END_REF]). The invariant (stable and unstable) manifolds of these periodic orbits, consisting of all trajectories converging to the orbit (as the time tends to ±∞), are four-dimensional tubes, topologically equivalent to S 3 × IR, in the five-dimensional energy manifold (see [START_REF] Gómez | Connecting orbits and invariant manifolds in the spatial three-body problem[END_REF]). Hence they play the role of separatrices. Therefore they can be used for mission design and space exploration, since a trajectory starting inside such a tube (called transit orbit) stays inside this tube. It can be noted however that the invariant manifolds of halo orbits (which can be really seen as tubes) are chaotic in large time: they do not keep their aspect of tube and behave in a chaotic way, far from the halo orbit (see [START_REF] Koon | Dynamical Systems, the three-body problem and space mission design[END_REF]). In contrast, the invariant manifolds of eight-shaped Lissajous orbits 5 (which are eight-shaped tubes) are numerically shown in [START_REF] Archambeau | Eight Lissajous orbits in the Earth-Moon system[END_REF] to keep their regular structure globally in time. Moreover, in the Earth-Moon system, it is shown that they permit to fly over almost all the surface of the Moon, passing very close to the surface (between 1500 and 5000 kilometers). These features are of particular interest in view of designing low-cost space missions to the Moon. Indeed in the future space exploration the Moon could serve as an intermediate point (with a lunar space station) for farther space missions.

Applications to mission design and challenges

The idea of using the specific properties of the dynamics around Lagrange points in order to explore lunar regions is far from new but has recently received a renewal of interest (see e.g. [START_REF] Farquhar | Utilization of libration points for human exploration in the Sun-Earth-Moon system and beyond[END_REF][START_REF] Koon | Dynamical Systems, the three-body problem and space mission design[END_REF][START_REF] Renk | Exploration missions in the Sun-Earth-Moon system: a detailed view on selected transfer problems[END_REF][START_REF] Topputo | Low energy interplanetary transfers exploiting invariant manifolds of the restricted three-body problem[END_REF]). In [START_REF] Koon | Dynamical Systems, the three-body problem and space mission design[END_REF][START_REF] Mingotti | Combined optimal low-thrust and stable-manifold trajectories to the Earth-Moon halo orbits[END_REF][START_REF] Mingotti | Low-energy, low-thrust transfers to the Moon[END_REF][START_REF] Mingotti | Invariant manifold, low-thrust transfers to low Mars orbits[END_REF], the authors combine the use of low-thrust propulsion with the use of the nice properties of invariant manifolds of periodic orbits around Lagrange points in order to design low-cost trajectories for space exploration. Their techniques consist of stating an optimal control problem that is numerically solved using either a direct or an indirect transcription, carefully initialized with the trajectories of the previously studied system (with no thrust). In such a way they are able to realize a reasonable compromise between fuel consumption and time of transfer, and design trajectories requiring moderate propellant mass and reaching the target within reasonable time.

Hence, in these studies the previously studied circular restricted three-body problem approximation is used to provide an appropriate first guess for carefully initializing an optimal control method (for instance, a shooting method) applied to a more precise model. In view of that, and having in mind the previous methodology based on continuation, it is natural to develop an optimal planification method based on the combination of the dynamics of the three-body problem with a continuation on the value of the maximal authorized thrust. This idea is used in the recent article [START_REF] Caillau | Minimum time control of the restricted three-body problem[END_REF] where a homotopy procedure is carried out on the maximal value of the thrust, starting from a zero value (natural dynamics), and ending with a low value of the thrust. The authors are then able to design minimal time trajectories with low thrust passing from a geostationary orbit around the Earth to a circular lunar one. This idea opens new directions for future investigations and is a promising method for designing efficiently fuel low-consumption space missions. Although the properties of the dynamics around Lagrange points have been widely used for developing planification strategies, up to now, and up to our knowledge they have not been combined with continuation procedures that would permit to introduce, for instance, the gravitational effects of other bodies, or values of the maximal thrust that are low or mild, or other more complex models. This is a challenge for potential future studies.

Note that, in [START_REF] Picot | Shooting and numerical continuation methods for computing time-minimal and energy-minimal trajectories in the Earth-Moon system using low propulsion[END_REF], the author implements a numerical continuation procedure to compute minimal-energy trajectories with low thrust steering the engine from the Earth to the Lagrange point L 1 in the Earth-Moon system, by making a continuation on the gravitational constant of the Moon. The continuation procedure is initialized with the usual Kepler transfer, in which the Moon coincides with the point L 1 , and ends up with a trajectory reaching the point L 1 with a realistic gravitational effect of the Moon.

Another challenge, which is imperative to be solved within next years, is the problem of debris cleaning. Indeed, recently it was observed a drastic growth of space debris in the space around the Earth, in particular near the SSO orbit and polar orbits with altitude between 600 and 1200 km (indeed these orbits are intensively used for Earth observation). These debris are due to former satellites that were abandoned, and now cause high collision risks for future space flights.

It has become an urgent challenge to clean the space at least from its biggest debris in order to stabilize the debris population, otherwise it will soon become impossible to launch additional satellites. At present, all space agencies in the world are aware of that problem and are currently working to provide efficient solutions for designing space debris collecting missions. One of them, currently led at EADS (see [START_REF] Cerf | Multiple space debris collecting mission: debris selection and trajectory optimization[END_REF]), consists of deorbiting five heavy debris per year, selected in a list of debris (in the LEO region) so that the required fuel consumption for the mission is minimized. The problem to be solved turns into a global optimization problem consisting of several continuous transfer problems and of a combinatorial path problem (selection of the debris and of the collecting order). It is not obvious to solve since it must combine continuous optimal control methods with combinatorial optimization, and other considerations that are specific to the problem. The results of [START_REF] Cerf | Multiple space debris collecting mission: debris selection and trajectory optimization[END_REF] (which are valuable for high-thrust engines) provide first solutions in this direction, and open new problems for further investigation. For instance it is an open problem to design efficient space cleaning missions for low-thrust engines, taking benefit of the gravitational effects due to Lagrange points and to invariant manifolds associated with their periodic orbits. Such studies can probably be carried out with appropriate continuation procedures, carefully initialized with trajectories computed from the natural dynamics of the three-body problem.

Conclusion and final remarks 6.1 Optimal control and trajectory optimization

Although the techniques of optimal control surveyed in this article provide a nice way to design efficient trajectories in particular in aerospace problems, their applications require a reasonably simple model. In practice many problems remain difficult due to the complexity of real-life models. For instance in the problem of low-thrust orbit transfer, many problems remain such as the one of taking into account the gravitational perturbations due to the Earth or the Moon, the atmospheric drag, the constraint of crossing the Van Allen barrier as quickly as possible, cone constraints on the control, eclipse constraints, taking into account the launching phase, and the insertion of the problem in a more global one, using multidisciplinary optimization. The eclipse constraint in particular may be viewed as a state constraint and can be handled by modelizing the system as a hybrid system. This problem is called the shadow cone problem. The objective is to develop necessary optimality conditions leading to efficient computation algorithms. Classical approaches are based on penalization methods and there is a challenging problem to use rather shooting methods, based on a Pontryagin approach, which are potentially more efficient from the point of view of the convergence. Of course all the constraints mentioned above are of different nature. Some of them can probaly be treated for instance by using some continuation procedures, but some others are not so well adapted to the use of indirect methods and then, according to the problem, one should then rather use direct methods (see discussions in [START_REF] Betts | Survey of numerical methods for trajectory optimization[END_REF][START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF]). There is a compromise to be found between the complexity of the model under consideration, the robustness of that model, and choice of an adapted numerical method to treat the problem.

The presence of numerous constraints on the controls or on the state makes the trajectography problems difficult, and in certain situations as for instance in the problem of atmospheric re-entry, a preliminary step to optimization is an estimation of the accessible set. A challenging question is then to combine the tools of numerical optimal control and of stabilization with fine geometric techniques developed recently in nonlinear control in order to design an estimation tool of the accessible set.

In many situations the system under consideration is modeled by infinite dimensional systems (PDE's). From the mathematical point of view many difficulties arise, and from the numerical point of view one has to use sophisticated techniques of numerical analysis (on meshes in particular). Applicative issues concern for instance motor outflows, fuel optimal management, the minimization of electromagnetic interferences or other perturbations (acoustic, thermic, etc). The modelization of these problems uses optimal control of PDE's. The numerical implementation of necessary optimality conditions (of Pontryagin Maximum Principle type) causes numerous difficulties, as put in evidence e.g. in [START_REF] Zuazua | Propagation, observation, control and numerical approximation of waves approximated by finite difference method[END_REF]. For hyperbolic equations interference phenomena appear between the mesh and highfrequency modes. Some remedies do exist, such as for instance highfrequencies filtering or the use of multigrid methods, and the objective is to adapt and apply them to the complex systems stemming from aerospace. Note that in certain cases the system in consideration is made of a finite dimensional system coupled with a "quite simple" partial differential equation, for instance the problem of optimizing the trajectories of planes in order to minimize noise pollution. The model should take into account sound propagation, shock waves, and thus, wave-like PDE's. Nowadays the control of such coupled systems is a challenging problem in mathematics and recently some first results have been published, which show how nonlinear couplings may help to recover controllability properties for the system (see [START_REF] Coron | Control and nonlinearity[END_REF]). Many other problems require complex models based on nonlinear PDE's: propulsion problems, thermics, aerodynamics, etc. Due to the complexity of these problems, every optimization procedure is in general impossible, and the questions that arise are in general on nonlinear control, in view of applications concerning in particular the design and dimensioning of space engines.

Pluridisciplinary optimization

In celestial mechanics many issues are still to be investigated in the very interesting field of the dynamics around Lagrange points; in particular it should be done a precise cartography of all invariant manifolds generated by all possible periodic orbits (not only halo or eight-shaped orbits) around Lagrange points in order to make them useful for mission design. The existence of such invariant manifolds indeed makes possible the design of low-cost interplanetary missions. The design of trajectories taking advantage of these corridors, of gravitational effects of celestial bodies of the solar system, of "swing-by" strategies, is a difficult problem related to techniques of continuous and discrete optimization (multidisciplinary optimization). Is is an open challenge to design a tool combining refined techniques of nonlinear optimal control, continuation procedures, mixed optimization, and global optimization procedures.

Many problems are modeled by hybrid systems, that is, systems whose dynamics may evolve with the time and contain discrete variables. An example is the problem of shadow cone constraint, another one is the global launcher problem in which the dynamics change whenever modules fall down. A certain number of deep theoretical results do exist on Pontryagin Maximum Principle in the hybrid case (see [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF][START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF][START_REF] Vinter | Optimal control. Systems & Control: Foundations & Applications[END_REF]) but the question of an efficient numerical implementation is still open in general (see [START_REF] Haberkorn | Convergence results for smooth regularizations of hybrid nonlinear optimal control problems[END_REF]); indeed when one implements a version of hybrid maximum principle, one is then immediately faced with a combinatorial explosion. It is not clear how to adapt efficiently tools of pluridisciplinary optimization to that difficult problem.

Another optimization problem is to determine the optimal placement of actuators, controls, in order to minimize or maximize a certain criterion: for instance, where the retrorockets should be optimally placed for the attitude control of a satellite, what should be the optimal shape of tailpipes in order to guarantee the best possible fluid outflow, where should the injectors be positioned in a motor to maximize combustion, etc. This kind of problem is part of the thematic of optimal design, in which the unknown is no more a vector but a domain. The problems of optimal locations of sensors or actuators in linear partial differential equations have been widely considered in engineering problems (see e.g. [START_REF] Kubrusly | Sensors and controllers location in distributed systems -a survey[END_REF][START_REF] Van De Wal | A review of methods for input/output selection[END_REF] and references therein). Usual popular approaches consist of recasting the optimal location problem for distributed systems as an optimal control problem with an infinite dimensional Riccati equation and then of computing approximations with optimization techniques. These techniques rely on an exhaustive search over a predefined set of possible candidates and are limited with combinatorial difficulties due to the selection problem. We thus recover the usual flaws of combinatorial optimization methods. Many challenging problems fall into this category.

An important problem in aerospace is the optimal design of launchers. The objective is to optimize both the trajectory and the engine (launcher). In the optimal design of a launcher one may seek to optimize thrust levels, the number of floors, of tanks, or to know what type of propellant should be used. Methods classically used consist in splitting the global problem into sub-problems handled with specific methods like genetic algorithms. The development of pluridisciplinary optimization tools should provide some breakthroughs in this domain. Another very important problem which could be treated efficiently with this kind of approach is the problem of space cleaning mentioned previously. We indeed have at our disposal a precise catalog of fragments, wreckage, scraps, and one of the top priorities in the next years is to clean the space from big fragments (essentially coming from old satellites). The problem is to design optimally a space vehicle able to collect in minimal time a certain number of fragments, themselves being chosen in advance in the catalog in an optimal way. This problem combines techniques of continuous optimal control in order to determine a minimal time trajectory between two successive fragments, and techniques of discrete optimization for the best possible choice of the fragments to be collected.

Inverse problems

From a certain number of measures one aims at detecting flaws in the structure of an engine, which are due to shocks, thermic or electromagnetic problems. This is an inverse problem which requires a mesh adapted to the engine, adequately placed sensors providing the measures, and this problem is thus related to the aforementioned ones.

The objective may be also to reconstruct the electromagnetic, thermic or acoustic environment of a launcher (after take-off or along the flight) from measures, in order to protect efficiently the fragile, delicate components of the launcher like computers for instance. This is a difficult inverse problem modeled with nonlinear PDE's. The challenge is to develop a tool permitting to design efficiently a launcher in order to make it more robust, less sensitive to environment perturbations. From the point of view of numerical analysis this requires the development of mesh methods ir spectral methods that are adapted to this specific problem. It can be noted that since fractional derivatives appear naturally in fluid mechanics problems (acoustic in particular), for instance when computing a heat flux getting out from the side of a fluid outflow in function of the time evolution of the internal source, or in the modeling of viscoelastic materials, it is important to develop efficient numerical approximation schemes of fractional derivatives (schemes of Grunwald-Letnikov, or Gear, etc) as in [START_REF] Galucio | An adaptation of the Gear scheme for fractional derivatives[END_REF]. It is a challenge to improve fractional methods with optimal inverse reconstruction procedures and with optimal design issues.

Finally, a last inverse problem is the one of optimal design of sensors. The problem is to determine where the sensors should be placed in order to ensure an optimal observation of the system, for instance in view of ensuring the success of online guidance processes. The applications are numerous in aerospace, and this problem is connected to the previous ones, the measures serving also for instance to reconstruct the electromagnetic or thermic environment of an engine, or to detect flaws. This problem enters into the category of shape optimization problems. In a general way measures are taken to reconstruct an environment. A difficult question is to determine which measures are required in order to optimize the reconstruction and to be able to approximate in the best possible way an inverse problem.

Combined with guidance objectives, these inverse problems may probably be recast in terms of pluridisciplinary optimization, as discussed previously. This results into difficult, complex problems and raises very interesting challenges for the future.

n

  and P max are positive constants. They are drawn on figure 1 in the flight domain, in terms of the drag d = 1 2 SC D m ρv 2 and of v. The minimization criterion is the total thermal flux along the flight

Table 1 :

 1 Boundary conditions

	max n	,

The usual framework of calculus of variations consists of solving the problem of minimizing the action R 1 0 L(t, x(t), ẋ(t))dt among all possible curves x(•). With respect to optimal control, this corresponds to the trivial control system ẋ(t) = u(t), where all moving directions are authorized.

Singularities due to conjugate points may be either detected and then handled with specific methods, or can be removed generically by Sard arguments (see comments further, on homotopy methods).

A turning point is a point of the path of zeros at which λ(s) has a local extremum.

See also[START_REF] Smale | A convergent process of price adjustment and global Newton methods[END_REF] where such arguments are applied to derive the global Newton method by homotopy. The authors of[START_REF] Chow | Finding zeros of maps: homotopy methods that are constructive with probability one[END_REF] make a good historical survey on the use of Sard arguments in homotopy theories.

Eight-shaped Lissajous orbits are the Lissajous orbits of the second kind, in the sense that they are diffeomorphic to a curve having the shape of an eight. They are chiefly investigated in[START_REF] Archambeau | Eight Lissajous orbits in the Earth-Moon system[END_REF].