F Bonnans

P Martinon
email: martinon@cmap.polytechnique.fr

E Trélat

Singular arcs in the generalized Goddard's Problem *

Keywords: Optimal control, Goddard's problem, singular trajectories, shooting method, homotopy, direct methods

We investigate variants of Goddard's problems for nonvertical trajectories. The control is the thrust force, and the objective is to maximize a certain final cost, typically, the final mass. In this article, performing an analysis based on the Pontryagin Maximum Principle, we prove 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 we use is based on our previous theoretical analysis and consists in combining a shooting method with an homotopic method. The

Introduction

The classical Goddard's problem (see [START_REF] Goddard | A Method of Reaching Extreme Altitudes[END_REF][START_REF] Maurer | Numerical solution of singular control problems using multiple shooting techniques[END_REF][START_REF] Miele | Generalized variational approach to the optimum thrust programming for a vertical flight of a rocket[END_REF][START_REF] Seywald | Goddard problem in presence of a dynamic pressure limit[END_REF][START_REF] Tsiotras | Drag-law effects in the Goddard problem[END_REF]) consists in maximizing the final altitude of a rocket with vertical trajectory, the controls being the norm and direction of the thrust force. Due to nonlinear effects of aerodynamic forces, the optimal strategy may involve subarcs along which the thrust is neither zero nor equal to its maximal value, namely, since the control variable enters linearly in the dynamics and the cost function is over the final cost, singular arcs. A natural extension of this model for nonvertical trajectories is the control system

ṙ = v, v = - D(r, v) m v v -g(r) + C u m , ṁ = -b u , (1)
where the state variables are r(t) ∈ IR 3 (position of the spacecraft), v(t) ∈ IR 3 (velocity vector) and m(t) (mass of the engine). Also, D(r, v) > 0 is the drag component, g(r) ∈ IR 3 is the usual gravity force, and b is a positive real number depending on the engine. The thrust force is Cu(t), where C > 0 is the maximal thrust, and the control is the normalized thrust u(t) ∈ IR 3 , submitted to the constraint

u(t) ≤ 1. (2)
The real number b > 0 is such that the speed of ejection is C/b. Here, and throughout the paper, denotes the usual Euclidean norm in IR 3 .

We consider the optimal control problem of steering the system from a given initial point

r(0) = r 0 , v(0) = v 0 , m(0) = m 0 , (3)
to a certain target M 1 ⊂ IR 7 , in time t f that may be fixed or not, while maximizing a final cost. For the moment, there is no need to be more specific with final conditions and the cost. In real applications, the problem is typically to reach a given orbit, either in minimal time with a constraint on the final mass, or by maximizing the final mass, or a compromise between the final mass and time to reach the orbit. In our numerical experiments we will study the problem of maximizing the final mass (i.e., minimizing the fuel consumption) subject to a fixed final position r(t f ) = r f , the final velocity vector and final time being free.

Depending on the features of the problem (initial and final conditions, mass/thrust ratio, etc), it is known that control strategies that consist in choosing the control so that u(t) is piecewise constant all along the flight, either equal to 0 or to the maximal authorized value 1, may not be optimal, as a consequence of the high values of the drag for high speed. Optimal trajectories may indeed involve singular arcs, and it is precisely the aim to this article to perform such an analysis and prove that the use of singular arcs is relevant in the problem of launchers.

The article is structured as follows. In Section 2, we recall the Pontryagin Maximum Principle, and the concept of singular trajectories. A precise analysis of the optimal control problem is performed in Section 3, where extremals are derived, and singular trajectories are computed. Theorem 1 makes precise the structure of the optimal trajectories. Section 4 is devoted to numerical simulations. The problem is first implemented with indirect methods, based on our theoretical analysis with the maximum principle, and, numerically, our method uses a shooting method combined with an homotopic approach. The homotopic method, leading to a quadratic regularization, permits to tackle with the problem of nonsmoothness of the optimal control. Experiments are also made using direct methods, i.e., by discretizing control variables and solving the resulting nonlinear optimization problem. Less precise than the indirect one, this method permits however to validate our approach by checking that results are consistent with the previously computed solution.

Our results show, as expected, that taking into account singular arcs in the control strategy permits to improve slightly the optimization criterion.

The numerical simulations presented in this paper, using a simplified and more academic model and set of parameters, constitute the first step in the study of a realistic launcher problem.

Preliminaries

In this section we recall a general version of the Pontryagin Maximum Principle (see [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF], and for instance [START_REF] Bryson | Applied optimal control[END_REF] for its practical application), and a definition and characterizations of singular arcs.

Consider the autonomous control system in IR n

ẋ(t) = f (x(t), u(t)), (4)
where Define the cost of a control u on [0, t] by

f : IR × IR n × IR m -→ IR n is of class C 1 ,
C(t, u) = t 0 f 0 (x(s), u(s))ds + g 0 (t, x(t)),
where f 0 : IR n × IR m -→ IR and g 0 : IR × IR n → IR are of class C 1 , and x(•)

is the trajectory solution of (1) associated to the control u.

Consider the optimal control problem of finding a trajectory joining M 0 to M 1 and minimizing the cost. The final time may be free or not.

Pontryagin Maximum Principle

According to the Pontryagin Maximum Principle (see [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]), if the control u ∈ U associated to the trajectory x(•) is optimal on [0, T ], then there exists an absolutely continuous mapping p(•) : [0, T ] -→ IR n called adjoint vector, and a real number p 0 ≤ 0, such that the couple (p(•), p 0 ) is nontrivial, and such that, for almost every t ∈ [0, T ],

ẋ(t) = ∂H ∂p (x(t), p(t), p 0 , u(t)), ṗ(t) = - ∂H ∂x (x(t), p(t), p 0 , u(t)), (5)
where H(x, p, p 0 , u) = p, f (x, u) + p 0 f 0 (x, u) is the Hamiltonian of the optimal control problem. Moreover, the function

t -→ max v∈Ω H(x(t), p(t), p 0 , v)
is constant on [0, T ], and the maximization condition

H(x(t), p(t), p 0 , u(t)) = max v∈Ω H(x(t), p(t), p 0 , v) (6)
holds almost everywhere on [0, T ].

Moreover, if the final time T to join the target set M 1 is free, then

max v∈Ω H(x(t), p(t), p 0 , v) = -p 0 ∂g 0 ∂t (T, x(T )). ( 7
)
for every t ∈ [0, T ].

Furthermore, if M 0 and M 1 (or just one of them) are submanifolds of IR n having tangent spaces in x(0) ∈ M 0 and x(T ) ∈ M 1 , then the adjoint vector can be chosen so as to satisfy the transversality conditions at both extremities (or just one of them)

p(0) ⊥ T x(0) M 0 (8)
and

p(T ) -p 0 ∂g 0 ∂x (T, x(T )) ⊥ T x(T ) M 1 . (9)
An extremal of the optimal control problem is a fourth-tuple (x(•), p(•), p 0 , u(•))

solution of ( 5) and ( 6). If p 0 = 0, then the extremal is said to be abnormal , and if p 0 = 0 then the extremal is said to be normal .

Singular arcs

Given x 0 ∈ IR n and two real numbers t 0 , t 1 , with t 0 < t 1 , denote by U x 0 ,t 0 ,t 1 the set of controls u ∈ L ∞ ([t 0 , t 1 ], Ω1), with Ω1 an open subset of Ω, such that the trajectory t → x(t, x 0 , t 0 , u), solution of (1), associated with the control u on [t 0 , t 1 ], and such that x(t 0 ) = x 0 , is well defined on [t 0 , t 1 ].

Define the end-point mapping E x 0 ,t 0 ,t 1 by E x 0 ,t 0 ,t 1 (u) := x(t 1 , x 0 , t 0 , u), for every u ∈ U x 0 ,t 0 ,t 1 . It is classical that E x 0 ,t 0 ,t 1 : U x 0 ,t 0 ,t 1 → IR n is a smooth map.

A control u ∈ U x 0 ,t 0 ,t 1 is said to be singular if u is a critical point of the end-point mapping E x 0 ,t 0 ,t 1 , i.e., its differential dE x 0 ,t 0 ,t 1 (u) at u is not surjective. In this case, the trajectory x(•, x 0 , t 0 , u) is said to be singular on

[t 0 , t 1 ].
Recall the two following standard characterizations of singular controls (see [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]). A control u ∈ U x 0 ,t 0 ,t 1 is singular if and only if the linearized system along the trajectory x(•, x 0 , t 0 , u) on [t 0 , t 1 ] is not controllable. This is also equivalent to the existence of an absolutely continuous mapping p 1 :

[t 0 , t 1 ] -→ IR n \ {0} such that, for almost every t ∈ [t 0 , t 1 ], ẋ(t) = ∂H 1 ∂p (x(t), p 1 (t), u(t)), ṗ(t) = - ∂H 1 ∂x (t, x(t), p 1 (t), u(t)), ∂H 1 ∂u (x(t), p 1 (t), u(t)) = 0,
where

H 1 (x, p 1 , u) = p 1 , f (x, u)
is the Hamiltonian of the system.

Note that singular trajectories coincide with projections of abnormal extremals for which the maximization condition (6) reduces to ∂H ∂u = 0.

For a given trajectory x(•) of the system (1) on [0, T ], associated to a control u ∈ U x(0),0,T , we say that x(•) involves a singular arc, defined on the subinterval [t 0 , t 1 ] ⊂ [0, T ], whenever the control u |[t 0 ,t 1 ] for the control system restricted to [t 0 , t 1 ] is singular.

In the case when the dynamics f and the instantaneous cost f 0 are linear in the control u, a singular arc corresponds to an arc along which one is unable to compute the control directly from the maximization condition of the Pontryagin maximum principle (at the contrary of the bang-bang situation). Indeed, in this case, the above condition ∂H 1 ∂u = 0 along the arc means that some function (called switching function) vanishes identically along the arc. Then, it is well known that, in order to derive an expression of the control along such an arc, one has to differentiate this relation until the control appears explicitly. It is also well known that such singular arcs, whenever they occur, may be optimal. Their optimal status may be proved using generalized Legendre-Clebsch type conditions or the theory of conjugate points (see [START_REF] Robbins | A generalized Legendre-Clebsch condition for the singular case of optimal control[END_REF][START_REF] Goh | Necessary conditions for singular extremals involving multiple control variables[END_REF], or 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] for a complete second-order optimality theory of singular arcs).

Analysis of the optimal control problem

With respect to the notations used in the previous section, we set

x =       r v m       ∈ IR 3 × IR 3 × IR, f (x, u) =       v -D(r,v) m v v -g(r) + C u m -b u      
, and f 0 = 0. Here, the set Ω of constraints on the control is the closed unit ball of IR 3 , centered at 0. Consider the optimal control problem of minimizing some final cost g 0 (t f , x(t f )), for the control system (1), with initial conditions (3) and final conditions x(t f ) ∈ M 1 in time t f which may be free or not.

We make the following assumption.

Assumption (H). The function g 0 is such that:

• either the final mass m(t f ) is free, and ∂g 0 ∂m = 0,

• or the final time t f is free, and ∂g 0 ∂t = 0.

In the first situation, the target set M 1 ⊂ IR 7 can be written as

M 1 = N 1 × IR,
where N 1 is a subset of IR 6 . A typical example is the problem of maximizing the final mass, for which

g 0 (t, x) = -m. If the final condition is r(t f ) = r 1 and v(t f ) = a, then M 1 = {r 1 } × S(0, a) × IR, where S(0, a)
is the sphere of IR 3 , centered at 0, with radius a.

In the second situation, a typical example is the minimal time problem to reach some target. In this case, g 0 (t, x) = t.

Computation of extremals

According to Section 2.1, the Hamiltonian of the optimal control problem under consideration is

H = p r , v + p v , - D(r, v) m v v -g(r) + C u m -p m b u , (10)
where , denotes the usual scalar product in IR 3 . Here, the adjoint vector is denoted by

p(t) =       p r (t) p v (t) p m (t)       ∈ IR 3 × IR 3 × IR.
In what follows, we assume the mappings D and g to be of class C 1 . Applying Pontryagin's Maximum Principle leads to the adjoint equations

ṗr = 1 m p v , v v ∂D ∂r + p v , ∂g ∂r , ṗv = -p r + 1 m p v , v v ∂D ∂v + D m p v v - D m p v , v v v 3 , ṗm = 1 m p v , - D(r, v) m v v + C u m . (11)
Moreover, if u is an optimal control on [0, t f ], then, for almost every t ∈ [0, t f ], u(t) maximizes the function

Φ t (w) := C m(t) p v (t), w -bp m (t) w , among all possible w ∈ IR 3 such that w ≤ 1.
The next technical lemma is the first step in the analysis of extremals.

Lemma 3.1. If there exists t 0 ∈ [0, t f ] such that p r (t 0 ) = p v (t 0 ) = 0, then p r (t) = p v (t) = 0, and p m (t) = p m (t f ), for every t ∈ [0, t f ]. Moreover, p m (t f ) = 0, and if p m (t f ) > 0 then u(t) = 0 on [0, t f ], otherwise u(t) = 1 on [0, t f ].
Proof. The first statement follows immediately from a uniqueness argument applied to the system [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF]. It follows from the expression of the Hamiltonian function that, if p m (t) > 0, then u(t) = 0, and if p m (t) < 0, then u(t) = 1.

In the first case of Assumption (H), the transversality condition (9) yields in particular

p m (t f ) = p 0 ∂g 0 ∂m (t f , x(t f )).
Therefore, p m (t) cannot be equal to zero (otherwise the adjoint vector (p, p 0 ) would be zero, contradicting the maximum principle). In the second case of Assumption (H), it follows from ( 7) and ( 10) that

p m (t)b u(t) = p 0 ∂g 0 ∂t (t f , x(t f )).
Therefore, similarly, p m (t) cannot be equal to zero. The conclusion follows.

An extremal satisfying the conditions of Lemma 3.1 (ie p r (t) = p v (t) = 0 for every t ∈ [0, t f ]) is called degenerate. For such extremals, the control is either identically equal to zero, or of maximal norm, along the whole trajectory. Such kind of trajectories can be excluded for practical applications and are thus discarded in the sequel.

Lemma 3.2. Consider a nondegenerate extremal. Then:

1. The set T := {t ∈ [0, t f ] | p v (t) = 0} has a finite cardinal.
2. There exists a measurable function α on [0,

t f ], with values in [0, 1],
such that

u(t) = α(t) p v (t) p v (t) , a.e. on [0, t f ]. ( 12
) 3. Set Ψ(t) := C m(t) p v (t) -bp m (t). Then, α(t) =      0 if Ψ(t) < 0, 1 if Ψ(t) > 0.
Proof. If t ∈ T , then by the costate equation [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF], ṗv (t) = -p r (t) is not zero (since the extremal is not degenerate). Therefore T has only isolated points, and hence, has a finite cardinal.

Writing w = αd, with α = w and d of unit norm, we get Φ t (w) =

α C m(t) p v (t), d -bp m (t) . Since p v (t) = 0 a.e.
, points 2 and 3 of the lemma follow immediately from the maximization condition.

The continuous function Ψ defined in Lemma 3.2 is called switching function. In the conditions of the lemma, the extremal control is either equal to 0, or saturating the constraint and of direction p v (t). The remaining case, not treated in this lemma and analyzed next, is the case where the function Ψ vanishes on a (closed) subset I ⊂ [0, t f ] of positive measure.

Remark 3.1. Let [t 0 , t 1 ] be a subinterval of I on which α(t) > 0. Then, the

control u |[t 0 ,t 1 ] is singular.
Indeed, it suffices to notice that, using [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF],

∂Φ t ∂w (u(t)) = C m(t) p v (t) -bp m (t) p v (t) p v (t) = Ψ(t) p v (t) p v (t) ,
and to use the Hamiltonian characterization of singular controls recalled in Section 2.2.

Singular arcs may thus occur in our problem whenever Ψ vanishes, and we next provide an analysis of that case, and show how to derive an expression of such singular controls.

Analysis of singular arcs

Throughout this section, we assume that Using this lemma, and extremal equations [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF], one gets, for a.e. t ∈ I,

Ψ(t) = C m(t) p v (t) -bp m (t) = 0 (13)
Ψ(t) = bC m(t) 2 ( p v (t) u(t) -p v (t), u(t) ) + Ξ(r(t), v(t), m(t), p(t)) = 0, ( 14
)
where the function

Ξ(r, v, m, p) = Db m 2 v p v , v + C m p v p v , p r + p v , v m v ∂D ∂v , p v + ∂D ∂m p v 2 m v - D m p v , v 2 v 3
does not depend on u. From Lemma 3.2, the relation ( 12) holds almost everywhere, and hence the first term of ( 14) vanishes. Therefore,

Ψ(t) = Ξ(r(t), v(t), m(t), p(t)) = 0, (15)
for almost every t ∈ I (actually over every subinterval of positive measure, since the above expression is continuous).

Relations ( 13) and ( 14) are two constraint equations, necessary for the existence of a singular arc. Derivating once more, using Lemma 3.3, leads to Ψ(t) = 0, a.e. on I.

The control u is expected to appear explicitly in this latter relation. However, since calculations are too lengthy to be reported here, we next explain how ( 16) permits to derive an expression for α(t), and hence, from [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF], an expression for u(t). When derivating [START_REF] Oberle | Numerical computation of singular control functions in trajectory optimizatio[END_REF], the terms where the control u appears are the terms containing v, ṗm , and ṁ. Recall that ṁ = -b u ,

that ṗm = 1 m p v , -D(r,v) m v v + C u m
, and that v is affine in u. Hence, since α(t) ≥ 0, it is not difficult to see that this derivation leads to an equation of the form

A(r, v, m, p r , p v , p m )α = B(r, v, m, p r , p v , p m ), (17)
almost everywhere on I. This relation should be "generically" nontrivial, that is, the coefficient A should not be equal to zero. This fact proves to hold true on numerical simulations. We explain below rigorously why this is true generically at least in the case of a scalar control (recall that we deal here with a three-dimensional control). For a scalar control, the control system (1) is of the form q = f 0 (q) + uf 1 (q), [START_REF] Allgower | Numerical Continuation Methods[END_REF] where f 0 and f 1 are smooth vector fields, and q is the state. In this case, it is well known (see e.g. [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]) that, if u is a singular control on I, then there must exist an adjoint vector p such that p, f 1 (q) = 0 on I,

p,

[f 0 , f 1 (q)] = 0 on I, (20)
p,

[f 0 , [f 0 , f 1 (q)] + u p, [f 1 , [f 0 , f 1 (q)] = 0 a.e. on I. (21)
The situation encountered here for 3D Goddard's problem is similar to that case: Equations ( 19), ( 20), [START_REF] Hiltmann | Numerische Lösung von Mehrpunkt-Randwertproblemen und Aufgaben der optimalen Steuerung mit Steuerfunktionen über endlichdimensionalen Räumen[END_REF], are respectively similar to Equations ( 13), ( 15), [START_REF] Martinon | Numerical resolution of optimal control problems by a Piecewise Linear continuation method[END_REF]; Equations ( 19), (20) (similarly, Equations ( 13), ( 15)) are constraint equations, and Equation ( 21) (similarly, Equation ( 16)) permits in general to derive an expression for the control u. The term "generic" employed above can now be made more precise: it is proved in [START_REF] Bonnard | Generic properties of singular trajectories[END_REF] that there exists an open and dense (in the sense of Whitney) subset G of the set of couples of smooth vector fields such that, for every control system [START_REF] Allgower | Numerical Continuation Methods[END_REF] with

(f 0 , f 1 ) ∈ G, the set where p, [f 1 , [f 0 , f 1 (q)
] vanishes has measure zero, and hence Equation ( 21) always permits to derive u. Additionaly, we can notice that the classical one-dimensional Goddard problem can be formulated as a particular case of the general 3D problem described here. In this case, it is well known that the second derivative of the switching function provides the expression of the singular control, so we can safely assume that 17 is nontrivial for the restriction to the 1D problem. Based on these arguments, we should expect the coefficient A of Equation ( 17) to be non zero in general. This is indeed the case in our numerical simulations presented next. Of course, once α(t) has been determined, one has to check (numerically) that 0 ≤ α(t) ≤ 1, so that the constraint u ≤ 1 is indeed satisfied. Here also, numerical simulations show the existence and admissibility of such singular arcs (see Section 4).

Conclusion

We sum up the previous results in the following theorem.

Theorem 1. Consider the optimal control problem of maximizing a fi-nal cost g 0 (t f , x(t f )), for the control system (1), with initial conditions

(3) and final conditions x(t f ) ∈ M 1 . We assume that Assumption (H) holds. Let u be an optimal control defined on [0, t f ], associated to the trajectory (r(•), v(•), m(•)). Then, there exist absolutely continuous mappings

p r (•) : [0, t f ] → IR 3 , p v (•) : [0, t f ] → IR 3 , p m (•) : [0, t f ] → IR,
Ψ(t) = C m(t) p v (t) -bp m (t).
Then,

• if Ψ(t) < 0 then u(t) = 0; • if Ψ(t) > 0 then u(t) = pv(t) pv(t) ;
• if Ψ(t) = 0 on a subset I ⊂ [0, t f ] of positive Lebesgue measure, then Equation ( 15) must hold on I, and

u(t) = α(t) p v (t) p v (t)
a.e. on I, where α(t) ∈ [0, 1] is determined by [START_REF] Schwartz | Theory and implementation of numerical methods based on Runge-Kutta integration for solving optimal control problems[END_REF].

Remark 3.2. The optimal control is piecewise either equal to zero, or saturating the constraint with the direction of p v (t), or is singular. Notice that, in all cases, it is collinear to p v (t), with the same direction.

Remark 3.3 (Optimality status). The maximum principle is a necessary condition for optimality. Second-order sufficient conditions are usually charac-available on the web1 , developed in [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]. This reference provides algorithms to compute the first conjugate time (where the trajectory ceases to be optimal) along a smooth extremal curve, based on theoretical developments of geometric optimal control using second order optimality conditions. The computations are related to a test of positivity of the intrinsic second order derivative or a test of singularity of the extremal flow.

It can be checked as well that every smooth sub-arc of the trajectory is locally optimal between its extremities. However, the problem of proving that the whole trajectory (i.e., a succession of bang and singular arcs) is

locally optimal is open. Up to now no conjugate point theory exists to handle that type of problem. Of course, one could make vary the times of switchings but this only permits to compare the trajectory with other trajectories having exactly the same structure. A sensitivity analysis is actually required to treat trajectories involving singular subarcs.

Numerical experiments

In this section, we provide numerical simulations showing the relevance of singular arcs in the complete Goddard's Problem. For given boundary conditions, the optimal trajectory is first computed using indirect methods (shooting algorithm) combined with an homotopic approach. Then we use a direct method (based on the discretization of the problem) to check the obtained solution. All numerical experiments were led on a standard computer (Pentium 4, 2.6 GHz).

Numerical values of the parameters of the model

We implement the optimal control problem of maximizing m(t f ) for the system (1), with the constraint (2). The equations of motion can be made dimension free by scaling the equations and choosing the model parameters in terms of r(0), m(0), and g 0 . We follow [START_REF] Oberle | Numerical computation of singular control functions in trajectory optimizatio[END_REF] and set the following parameters.

-g 0 = r 0 = m 0 = 1.

-C = 3.5, b = 7.

-Drag D(r, v) = K D v 2 e -500( r -1) with K D = 310.

-Initial and final conditions r 0 = (0.999949994, 0.0001, 0.01), v 0 = (0, 0, 0), m 0 = 1,

r f = (1.01, 0, 0), v f is free, m f is free. t f is free.
With the scaled equations, r = 1 corresponds to a distance of 6378km (Earth radius), v = 1 to a speed of 7910m.s -1 and t = 1 to a time of 806s.

The prescribed final position is therefore at an altitude of 63.78km, and r 2 (0) being small, the trajectory will be nearly planar (this is consistent with the case of a real launch from a nearly equatorial base to a geostationary transfer orbit).

Numerical simulations with indirect methods

In our simulations presented hereafter, we prefer to express the objective of the optimal control problem in the following form.

20

Maximizing m(t f ) is equivalent to minimizing the cost

t f 0 u(t) dt,
and we assume that there are no minimizing abnormal extremals, therefore the adjoint vector can be normalized so that p 0 = -1. The results of the simulations are consistent with this assumption.

According to Section 2.1, the Hamiltonian of the optimal control problem under consideration is

H = p r , v + p v , - D(r, v) m v v -g(r) + C u m -(1 + bp m ) u ,
The only difference with the Hamiltonian in 2.1 for the M ax m(t f ) objective is the additional "-1" in the u term, which leads to the switching function

ψ(t) = C m(t) p v (t) -(1 + bp m (t)),
• if ψ(t) < 0 then u(t) = 0;

• if ψ(t) > 0 then u(t) = pv(t) pv(t) ;

• if ψ(t) = 0 on I ⊂ [0, t f ], then Equation ( 15) must hold on I, the control u is singular, and

u(t) = α(t) p v (t) p v (t)
a.e. on I, where α(t) ∈ [0, 1] is determined by [START_REF] Schwartz | Theory and implementation of numerical methods based on Runge-Kutta integration for solving optimal control problems[END_REF]. We check numerically that 0 ≤ α(t) ≤ 1.

Furthermore, on a singular subarc, derivating the switching function twice yields the expression of α via a relation of the form A(x, p)α = B(x, p), see [START_REF] Schwartz | Theory and implementation of numerical methods based on Runge-Kutta integration for solving optimal control problems[END_REF]. The computations are actually quite tedious to do by hand, and we used the symbolic calculus tool Maple. The expressions of A and B are quite complicated and are not reported here.

The free final time problem is formulated as a fixed final time one via the usual time transformation t = t f s, with s ∈ [0, 1] and t f an additional component of the state vector, such that ṫ f = 0 and t f (0), t f (1) are free, with the associated costate satisfying ṗ t f = -H. All the graphs in the following will use this normalized time interval [0, 1].

Transversality conditions on the adjoint vector yield p v (1) = (0 0 0), p m (1) = 0, and p t f (0) = p t f (1) = 0. The unknown of the shooting function

S is therefore z = (t f , p r (0), p v (0), p m (0)) ∈ R 8 .

Homotopic approach

In the indirect approach, it is necessary to get some information on the structure of the solutions, namely, to know a priori the number and approximate location of singular arcs. To this aim, we perform a continuation (or homotopic) approach, and regularize the original problem by adding a quadratic ( u 2 ) term to the objective, as done for instance in [START_REF] Martinon | Numerical resolution of optimal control problems by a Piecewise Linear continuation method[END_REF][START_REF] Schwartz | Theory and implementation of numerical methods based on Runge-Kutta integration for solving optimal control problems[END_REF]. The general meaning of continuation is to solve a difficult problem by starting from the known solution of a somewhat related, but easier problem. Here, we regularize the cost function by considering an homotopic connection with an energy,

t f 0 ( u(t) + (1 -λ) u(t) 2 ) dt, (22)
where the parameter of the homotopy is λ ∈ [0, 1]. The resulting perturbed problem (P λ ) has a strongly convex Hamiltonian (with respect to u), with a continuous optimal control, and is much easier to solve than (P ) = (P 1 ).

Assuming we have found a solution of (P 0 ), we want to obtain a solution of (P ) for λ = 1 (or at least sufficient information).

The continuation can be conducted manually, by solving iteratively a suitable sequence of problems P (λ k ) from λ = 0 to λ = 1. However, finding such a sequence can be quite difficult in practice, which is why we chose here to perform a full path-following continuation. Extensive documentation about path following methods can be found in [START_REF] Allgower | Numerical Continuation Methods[END_REF], and we use here a simplicial method, which basically builds a piecewise linear approximation of the homotopy path over a triangulation of the research space. The reason behind the choice of this method over a more classical predictor-corrector continuation is that we expect the problem to be ill-conditioned, due to the presence of singular arcs, which is indeed the case in the numerical experiments.

Preliminary continuation on the atmosphere density

In our case, even solving the regularized problem (P 0 ) is not obvious, due to the aerodynamic forces (drag). For this reason, we introduce a preliminary continuation on the atmosphere density, starting from a problem without atmosphere. Technically, this is done by using an homotopy with the modified parameter

K θ D = θK D , θ ∈ [0, 1],
where K D appears in the model of the drag. The shooting method for the problem without atmosphere at θ = 0 converges immediately with the trivial initialization p r (0) = (0.1, 0.1, 0.1), p v (0) = (0., 1, 0.1, 0.1), p m (0) = 0.1 and t f = 0.1. We would like to emphasize the fact that we have here no difficulties to find a starting point for the shooting method.

We then perform the first continuation from θ = 0 to θ = 1, gradually introducing aerodynamic forces. The solution we obtain is sufficient to initialize the shooting method at the beginning of the main homotopy.

Notice that a direct continuation on the atmosphere with the original non regularized problem (P ) fails. During the continuation, the process abruptly diverges at a certain value for θ, certainly due to the appearance of the singular arc.

Main continuation on the quadratic regularization

We now perform the main continuation on the cost [START_REF] Stoer | Introduction to Numerical Analysis[END_REF]. Figure 1 represents the solutions for λ = 0, 0.5 and 0.8. The shape on the switching function and of the control norm graphs are particularly interesting concerning suspicion of singular arcs. Indeed, we observe that, on a certain time interval (roughly [0.1, 0.4]), the switching function comes closer to zero as λ increases, while the control norm keeps values in (0, 1). These facts strongly suggest the appearance of a singular arc. Trying to go further with the continuation as λ → 1 becomes extremely difficult since we lose the singular structure and encounter trajectories with incorrect bang-bang structures. However the knowledge of the solution for λ = 0.8 happens to be sufficient to solve the problem: it provides a good starting point for which the shooting method applied to the original problem (P ) converges.

Shooting method applied to the original problem (P )

When implementing a shooting method (see for instance [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF][START_REF] Bulirsch | Die mehrzielmethode zur numerischen lösung von nichtlinearen randwertproblemen und aufgaben der optimalen steuerung[END_REF][START_REF] Hiltmann | Numerische Lösung von Mehrpunkt-Randwertproblemen und Aufgaben der optimalen Steuerung mit Steuerfunktionen über endlichdimensionalen Räumen[END_REF][START_REF] Stoer | Introduction to Numerical Analysis[END_REF]), the structure of the trajectory has to be known a priori. The structure of the control must be prescribed here by assigning a fixed number of interior switching times that correspond to junctions between nonsingular and sin-gular arcs. These times (t i ) i=1..n switch are part of the shooting unknowns and must satisfy some switching conditions. Each arc is integrated separately, and matching conditions must be verified at the switching times, as drawn on the diagram below.

Unknown: z IVP unknown at t0 (x 1 , p 1 ) ... (x s , p s ) t1 ... ts Value: S Sing (z) Switch cond (t1) M atch cond (t1) ... Switch cond (ts) M atch cond (ts) TC(t f )

Here, matching conditions reduce to imposing state and costate continuity at the switching times.

A switching condition indicates a change of structure, which corresponds here to an extremity of a singular arc. Along such a singular arc, it is required that ψ = ψ = 0. The control is computed using the relation ψ = 0. Therefore, using this expression of the control, switching conditions consist in imposing either ψ = 0 at the extremities of the singular arc, or ψ = ψ = 0 at the beginning of the arc. In our simulations, we choose the latter solution which happens to provide better and more stable results.

The previous results, obtained with an homotopic approach, provide an indication on the expected structure of the optimal trajectory for the original problem (P ). Inspection of Figure 1 suggests to seek a solution involving a singular arc on an interval [t 1 , t 2 ], with t 0 < t 1 < t 2 < t f . As a starting point of the shooting method, we use the solution previously obtained with the homotopy on the cost at λ = 0.8.

The IVP integration is performed with the radau5 code (see [START_REF] Hairer | Solving ordinary differential equations[END_REF]), with absolute and relative tolerances of, respectively, 10 -6 and 10 -6 . The shooting method converges in 17 seconds, with a shooting function of norm 5 10 -4 . In particular, the switching conditions at t 1 are well satisfied, with (ψ, ψ) = (2.00 10 -5 , 3.67 10 -7 ). The condition number for the shooting function is quite high (about 10 12 ), which was expected. The overall execution time of the whole approach (preliminary atmosphere homotopy, regularization homotopy, final shooting) is about 400 seconds.

The solution obtained has the following times

t 1 t f = 1.054679 10 -1 , t 2 t f = 3.

Numerical simulations with direct methods

In order to validate the solution obtained previously with the shooting algorithm, we next implement a direct method. Although direct methods can be very sophisticated (see for instance [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF][START_REF] Waechter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF]), we here use a very rough formulation, since our aim is just to check if the results are consistent with our solution. We discretize the control using piecewise constant functions, and the state is integrated on [0, t f ] with a basic fixed step Runge-Kutta fourth order formula. The values of the control at the discretization nodes, as well as the final time t f , thus become the unknowns of a nonlinear constrained optimization problem, the constraints being the final conditions for the state. To solve the optimization problem, we use the ipopt solver, which implements an interior point algorithm with a filter line-search method (see [START_REF] Waechter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF] for a complete description). We chose spherical coordinates u i = (r, θ, ψ), i = 1 . . . N for the 3D control instead of cartesian coordinates, as this makes the norm constraint u ≤ 1 a simple bound constraint. With N = 100 integration steps, the unknown size is 901 (discretized 3D control + free final time), and we take as initialization

     u i = (0.5, -0.1, 0), i = 1 . . . N t f = 0.2
We use as main options (see IPOPT manual) an error tolerance of 10 -5 , SR1 update to approximate reduced Hessian, and augmented Lagrangian for the line search. The algorithm converges after 173 iterations (and 185 seconds) to a solution with a final time of 0.2188 and a criterion value of 0.3999. This solution is clearly consistent with the results of the shooting method, as shown on Figure 4, which represents the norm of the control for the shooting method solution and the direct method solution, as well as a

Comparison with a bang-bang solution

Recall that the usual launch strategy consists in implementing piecewise controls either saturating the constraint or equal to zero. To prove the relevance of the use of singular controls in the control strategy, we next modify slightly the formulation above in order to find a bang-bang solution. Our aim is to demonstrate that taking into account singular arcs in the control strategy actually improves (as expected) the optimization criterion.

We set a "on-off " structure, with only one switching time t off . The control is chosen so as to satisfy u(t) = 1 for t 0 < t < t off , and u(t) = 0 for t off < t < t f . Here, the unknowns of the optimization problem are t f , t off and the direction of the control at the discretization nodes before t off .

For the initialization, we keep the same values as before for the control direction and final time, and set t off = 0.05. We obtain a solution with t f = 0.2105, t off = 0.0580, and the value of the criterion is 0.4061, which represents a loss of about 1.6% compared to the solution with a singular arc. On this academic example, the gain of the optimal strategy, involving a singular arc, over a pure bang-bang strategy, is quite small. This simplified problem is a first step in the study of a realistic launcher problem, and permits to illustrate the method.

Extension to more realistic physical models

The target application for this study is the computation of optimal trajectories for multi-stage heavy launchers (such as Ariane). In such problems, the angle of attack is set to zero during the whole flight except for a small time period. Thus, lift force is neglected, and only the drag term appears in the equations. Adding a lift term to the model is possible, but would require a new analysis of the problem. Adaptation of the numerical method, however, should be straightforward. Another significant difficulty of using a more realistic model for aerodynamic forces is the use of tabulated physical data, such as atmospherical pressure, sound speed and aerodynamic coefficient (Cx). These tables make it impossible to obtain the analytical expression of ψ, therefore we need an alternate way of finding the value of the singular control.

As we want to ensure that ψ = 0 over a singular arc, we can try to chose the control that minimizes ψ 2 + ψ2 . We obtain convergence for the shooting method and find a similar solution, with the same singular arc. We observe on Figure 5 that the approximate singular control values are actually quite close to the exact singular control (difference is about 10 -3 ).

The drawbacks of this alternate formulation are a heavy increase in the numerical cost (about one order of magnitude), and some loss of accuracy (the norm of the shooting function is here only about 10 -2 compared to 10 -4 before). However, this approach should be applicable to more complicated and realistic physical models.

Conclusions

In this paper we give an analysis of optimal trajectories with singular subarcs for a nonvertical variant of Goddard's problem. The numerical simulations confirm the existence of an optimal trajectory with a singular arc, and show the relevance of using a continuation approach in order to determine the optimal control structure and obtain a suitable estimate of the initial adjoint variables. The simple experiments we made with direct methods are consistent with these results. We also give an alternate method for the singular control, that should be applicable to more complicated and realistic physical models.

  and where the controls are measurable and bounded functions defined on a subinterval [0, t e (u)[ of IR + with values in Ω ⊂ IR m . Let M 0 and M 1 be subsets of IR n . Denote by U the set of admissible controls u, whose associated trajectories are well defined and join an initial point in M 0 to a final point in M 1 , in time t(u) < t e (u).

Lemma 3 . 3 .

 33 Let a, b be real numbers such that a < b, and f : [a, b] → IR be an absolutely continuous function. Let J be a subset of {t ∈ [a, b] | f (t) = 0} of positive Lebesgue measure. Then f ′ (t) = 0 a.e. on J.

Figure 1 :

 1 Figure 1: Main homotopy -u and ψ for λ = 0, 0.5 and 0.8.

Figure 2 :

 2 Figure 2: Solution with singular arc: altitude, speed and mass.

Figure 3 :

 3 Figure 3: Solution with singular arc: control and switching function.

Figure 4 :

 4 Figure 4: Control norm for the shooting and direct method.

Figure 5 :

 5 Figure 5: Control norm and switching function for the alternate formulation.

  and a real number p 0 ≤ 0, such that (p x (•), p v (•), p m (•), p 0 ) is nontrivial, and such that Equations (11) hold a.e. on [0, t f ]. Define the switching function Ψ on [0, t f ]

	by

http://www.n7.fr/apo/cotcot

* Support from the French space agency CNES (Centre National d'Etudes Spatial) is gratefully acknowledged.

terized in terms of conjugate points (see e.g. [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]. Unfortunately standard theories do not apply here for two reasons: first, the equation in m(t) involves the term u(t) which is not smooth; second, the structure of trajectories stated in the theorem involves both bang arcs and singular arcs, and up to now a theory of conjugate points that would treat this kind of trajectory.

We mention however below a trick, specific to the form of our system, which permits to apply the standard theory of conjugate points on every subinterval J of [0, t f ] on which u is singular and 0 < u(t) < 1. Let J be such a subinterval. Then, ṁ = 0 a.e. on J, and the system can be reparametrized by -m(t). Then, denoting q = (r, v), system (1) yields

and consider as new control the control ũ = (v, θ 1 , θ 2 ). Notice that the controls θ 1 and θ 2 are unconstrained, and that v must satisfy the constraint v ≥ 1. However, along the interval J it is assumed that 0 < u(t) < 1, and thus v does not saturate the constraint. Hence, the standard theory of conjugate points applies and the local optimality status of the trajectory between its extremities on J can be numerically checked, for instance using the code COTCOT (Conditions of Order Two and COnjugate times),


P Martinon
email: pierre.martinon@inria.fr

F Bonnans
email: eric.bonnans@inria.fr

J Laurent-Varin
email: julien.laurent-varin@cnes.fr

E Trélat

Numerical study of optimal trajectories with singular arcs for an Ariane 5 launcher

Introduction

We consider a flight mission to the geostationary transfer orbit (GTO) for an Ariane 5 launcher, while maximizing the payload or, as a variant, minimizing the fuel consumption. We first solve the complete flight sequence up to the final orbit, assuming a maximal thrust for all propulsion systems. Then we focus on the atmospheric ascent phase, which has been studied for instance in [START_REF] Seywald | Goddard problem in presence of a dynamic pressure limit[END_REF][START_REF] Calise | Optimization of launch vehicle ascent trajectories with path constraints and coast arcs[END_REF][START_REF] Lu | Closed-loop endoatmospheric ascent guidance[END_REF]. We are more specifically interested in optimal trajectories with singular arcs (flight phases with a non maximal thrust) for the boosters. Due to the presence of tabulated data in the physical model, the exact expression of the singular control cannot be obtained from the time derivatives of the switching function.

An alternate way to compute the singular control is provided, and numerical experiments are carried out for for several launcher variants.

Problem statement

State and control variables

The state variables include the position and speed of the launcher, in dimension 3, as well as the masses of the three fuel-consuming parts of the launcher: m 1 for the two boosters (treated as one propulsion system), m 2 for the first stage, and m 3 for the second stage (m denoting the total mass). This allows us to treat the separations more easily, avoiding the discontinuities on the state variables. Such discontinuities can be treated properly when applying Pontryagin's Minimum Principle (see for instance [START_REF] Pesch | A practical guide to the solution of real-life optimal control problems[END_REF]), but in view of the disappointing numerical results, we got rid of discontinuities by splitting the total mass. The original problem being with free final time, we use the standard transformation to obtain a formulation with a fixed final time.

The control variables include the throttle α ∈ [0, 1] for the boosters (we assume a maximal thrust for the propulsion systems of the stages), and the flight angles θ and ψ (heading and azimuth) giving the thrust direction. The control vector is then

u = (α, θ, ψ) ∈ [0, 1] × [-π, π] 2 .
(1)

Flight dynamics

The dynamics for the masses correspond to the fuel consumption, which depends on the flight phase. The flight dynamics are

   ṙ = v v = 1 m (T (r, u) -D(r, v)) + g(r) ṁi = -α β (for each fuel-consuming part) (2)
with • T (r, u) = α β Isp g 0 -S P z : rocket thrust. β : mass flow rate; Isp : specific impulse S : area at rocket exit; P z : atmospheric pressure (table )  • D(r, v) = q Sr C x (M ): drag due to earth atmosphere. q : dynamic pressure; S r : reference area

C x : drag coefficient (table); M : speed in Machs • g(r) = -µ r 3 (r + J 2 R 2
T r 2 M J2 r): gravity with J 2 corrector term.

Note: the lift term in the aerodynamic forces is here omitted.

Applying Pontryagin's Minimum Principle shows a decoupling for the control variables. Minimizing the Hamiltonian leads to the thrust direction being opposite to the speed costate, while the throttle follows a bang-bang law, with possible switchings or singular arcs.

Complete flight at full thrust

We start by solving the complete flight with stage separations, by setting a full thrust for the boosters. We want here to maximize the payload, and assume all fuel is burnt. As we consider the mass of each component instead of the global mass, the state variables remain continuous at the separations times. Therefore the separation times t 1 and t 2 for the boosters and first stage only need to satisfy the continuity of the Hamiltonian. We introduce as additional shooting unknowns the value of the state and costate at the separation times, namely, (x 1 , p 1 ) and (x 2 , p 2 ), for a better numerical behaviour. We have the corresponding matching conditions that enforce the continuity of x and p at the separation times.

Numerical results

Solving the complete flight problem directly with the shooting method is difficult, due to the lack of a suitable initial guess. We first solve the flight restricted to the first phase (until separation of the boosters) then add the second and third phase. We obtain a solution with a payload of 14623kg and a final time of 1232.75 seconds. Boosters separation occurs after 119 seconds, and 533.375 seconds for the first stage, which is consistent with the constraints set on the final masses (consumption of all fuel).

Remark: The thrust is actually completely fixed during the first 10 seconds of the flight, vertical and maximal.

The first two graphs in Fig. 1 represent the optimal control with the thrust level for the boosters and two stages, here set to 1, and the heading and azimuth angles for the thrust direction; and the masses of the different parts of the launcher. The separations correspond to the discontinuities observed on the total mass. Remark: It should be kept in mind that we have allowed a free direction for the thrust, which is not the case in the real flight where the thrust direction is strongly constrained by the dynamic pressure limit. However, the trajectory we obtained is close to the reference one.

Study of singular arcs for the boosters

We focus now on the study of singular arcs for the boosters, whose throttle is now free in [0, 1], while keeping a maximal thrust for the two stages. We also restrict the flight to the first atmospheric climbing phase until the separation of the boosters. We reformulate the payload maximization criterion into a fuel consumption minimization, which is more suitable for the study of singular arcs. As we are only considering a variable throttle for the boosters (EAP), we set as the new criterion

M in t f 0 α β EAP . ( 3
)
The payload is now fixed and is no longer part of the state variables, and the final masses m 1 , m 2 , m 3 are free. We also replace the final condition on the orbit by new final conditions on the position and speed. We take the values corresponding to the boosters' separation on the trajectory previously computed for the complete flight with full thrust.

An alternate expression for the singular control

A singular arc is characterized by the fact that the switching function ψ (the derivative of the Hamiltonian with respect to the control) vanishes over a whole time interval. The expression of the singular control is usually obtained by differentiating the equation ψ(x, p) = 0 with respect to time until the control appears explicitely (the required order of differentation being always even, see for instance [START_REF] Robbins | A generalized Legendre-Clebsch condition for the singular case of optimal control[END_REF]). However, due to the presence of tabulated data in the thrust and drag terms, the analytic expression of ψ is already not available. We only have the first derivative ψ that depends on ẋ and ṗ. As the idea behind the formal expression of the singular control is that the switching function and its successive time derivatives vanish over a singular arc, we try to enforce this constraint with the terms at our disposal.

Therefore, for the singular control we choose at each time step the value αsing that minimizes ψ and ψ at the next integration step, i.e., αsing = ArgM in [0,1] (ψ(y i+1 ), ψ(y i+1 ))

t 2 = ψ 2 (y i+1 ) + ψ2 (y i+1 ), (4)
where y i+1 denotes the state-costate pair (x, p) obtained after one integration step from the current time. This minimization is currently performed by a BFGS method [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF], starting from an initial value α = 0.5. The numerical experiments showed that taking 0.1 or 0.9 as initial values gives similar results, which is of course reassuring. Furthermore, using the singular control continuity to initialize this minimization step (after the arc entry) seems to give slighltly better results than using a constant initialization.

We tested this formulation on the three dimensional Goddard problem (cf [START_REF] Bonnans | Singular arcs in the generalized Goddard's problem[END_REF]) and found control values quite close to the exact value α * sing obtained from the equation ψ(x, p) = 0, see Fig. 3. The difference between the exact singular control and the alternate formulation is about 10 -3 , but the computational cost is significantly increased. However, the minimization procedure was not optimized and so might be further improved.

Continuation approach for singular arcs

In order to solve a problem with singular arcs, besides the expression of the singular control, we also need some information about the control structure, namely the number and position of the singular arcs. To obtain such information, we approach the original problem by a sequence of regularized problems with strictly convex Hamiltonians (with respect to the throttle), such as in [START_REF] Gergaud | An application of PL continuation methods to singular arcs problems[END_REF][START_REF] Bonnans | Singular arcs in the generalized Goddard's problem[END_REF], for instance. This is done by adding a quadratic term to the criterion, which becomes, for given λ ∈ [0, 1]:

M in t f 0 α β EAP + (1 -λ) α 2 β EAP dt. ( 5
)
This regularization adds a new term (1λ) α 2 β EAP to the Hamiltonian, which gives the modified command law for λ < 1:

   If ψ(x, p) > 0 then α = 0 If ψ(x, p) < -2(1 -λ)β EAP then α = 1 If -2(1 -λ)β EAP < ψ(x, p) < 0 then α = -ψ(x,p) 2(1-λ)βEAP . (6)
The optimal control is now continuous, without singular arcs or switchings.

We perform a discrete continuation on this problem family, and solve an automated sequence of problems (P ) λ starting from λ = 0. We use a very basic algorithm to generate the sequence of problems:

• Solve (P 0 ) • Iterate: set λ k+1 = λ k + h and solve (P λ k+1 ) using (P λ k ) as initialization, with a linear prediction (h = 1 initially, and is decreased in case of failure) • Stop when reaching λ = 1, or maximum number of iterations / minimal stepsize for h.

Notice that we do not actually need to reach λ = 1, but we still need to obtain sufficient information regarding the control structure, as well as a suitable initialization for the shooting method adapted to the singular structure.

However, even the strongly regularized problem (P 0 ) is not easy to solve directly, so we add a second layer of continuation over atmospheric forces, gradually introducing the atmospheric drag for the problem (P 0 ). Numerical experiments indicate that it is usually easy to solve the regularized problem in vacuum and complete the continuation on the drag, which gives a solution for the regularized problem (P 0 ). We can then begin the regularization homotopy in order to determine the control structure. The results of this method are described below, first for the original Ariane 5 launcher, and then for a slightly modified launcher.

Note: All numerical simulations were run on a standard desktop computer (Pentium IV 2.4GHz), using the gfortran compiler. Using a fixed-step RK4 integration (100 steps) and the HYBRD solver (see [START_REF] Garbow | User Guide for Minpack-1[END_REF]) for the shooting method, the total cpu time for the atmosphere and regularization continuations and the final shooting is about 10 minutes.

Study of the original launcher

We first apply the approach described above to the original launcher problem, with a fixed payload m CU = 12610kg. The regularized problem in vacuum can be solved from a very simple starting point: (t f = 100, p r (0) = (-0.1, -0.1, 0.1), p v (0) = (-10 3 , -10 3 , -10 3 ) and p i (0) = -0.1, i = 1 . . . 3). The preliminary continuation is performed without any difficulties and gives a solution of (P 0 ). During the main continuation, we observe (see Fig. 4) that the interval with a non-maximal thrust becomes smaller when the regularization decreases as λ tends to 1. At λ = 0.999, we have a solution with a full thrust throughout the flight except for the last few seconds, but with no sign of a singular arc. Fig. 4 Regularization homotopy: thrust level for λ = 0, 0.5, 0.9, 0.999

We make some simple experiments with a direct method: piecewise constant control with 100 steps, RK4 integration for the state, IPOPT (see [START_REF] Waechter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF]) solver for the resulting optimization problem. We find a similar solution, with a full thrust and a switching to null thrust just before the end of the flight. Thus both methods suggest that the optimal solution to this problem does not present singular arcs.

Study of a modified launcher

Now we modify the launcher parameters in order to increase the drag effect: we increase the reference area S r and the specific impulse Isp EAP . For S r × 2 and Isp EAP × 1.25, we find an optimal trajectory with a singular arc.

As before, the regularized problem with no drag is easily solved, and the first continuation to introduce the drag poses no difficulties. Then we perform the regularization continuation until λ = 0.95, which shows strong hints of a singular arc. The graphs in Fig. 5 show the thrust level and switching function at the solutions for λ = 0, 0.5, 0.9 and 0.95. Contrary to the previous case, we observe a small time interval (around t = 30s) where the thrust level remains in ]0, 1[. Moreover, we see that the switching function ψ comes closer to 0 at the same time interval. These two facts together strongly suggest the presence of a singular arc. What happens at the end of the flight is less clear, as the thrust level again takes values in ]0, 1[, and decreases near t f . This could indicate a second singular arc, maybe followed by an arc with a null thrust, or just a switching to a null thrust arc. Fig. 5 Regularization homotopy: thrust level for λ = 0, 0.5, 0.9, 0.95

The shooting method fails to converge for structures with two singular arcs, but the formulation with one arc actually gives a solution with a singular arc and a switching at the end. We detail below the shooting formulation for a control structure assuming one interior singular arc. In addition to the usual shooting unknowns (final time and initial costate), we have the entry and exit times for the singular arc. The corresponding additional conditions ensure that we enter the switching surface tangentially at t entry , i.e. ψ(x(t entry ), p(t entry ) = ψ(x(t entry ), p(t entry )) = 0, see for instance [START_REF] Pesch | A practical guide to the solution of real-life optimal control problems[END_REF]. Notice that the switching at the end is handled automatically by a switching detection algorithm checking the sign of the switching function during the integration (see for instance [START_REF] Hairer | Dense output, discontinuities, derivatives[END_REF], pages 195-200).

Shooting formulation for one interior singular arc

• Shooting function unknown z ∈ R 12 : z = (t f , p r (0), p v (0), p m (0), t entry , t exit ).

We use the solution from the end of the regularization continuation (λ = 0.986 after 100 iterations) to initialize z.

• Shooting function value S(z) ∈ R 12 :                r(t f ) -r f F inal position v(t f ) -v f F inal speed p t f (t f ) F inal time T C p mi (t f ) , i = 1..3
F inal masses T C ψ(x(t entry ), p(t entry )) Switching conditions ψ(x(t entry ), p(t entry )) at entry time [START_REF] Bonnans | Singular arcs in the generalized Goddard's problem[END_REF] Once again, we test a basic direct method, and find a similar solution. The graph on Fig. 6 shows these results, with the small singular arc (about 2 seconds) clearly visible on both solutions from the indirect and direct methods, around t = 30s. Both solutions also confirm that we have a simple switching near the end of the flight, and not a second singular arc. As expected, the switching function ψ and its first derivative ψ are both close to 0 on the singular arc, around 10 -3 and 10 -4 respectively. Thus the approach of computing the singular control that minimizes ψ 2 + ψ2 worked rather well on this problem.

The criterion value (i.e., final mass of the EAP) is quite close for both solutions, with m 1 (t f ) = 165163kg for the shooting method and m 1 (t f ) = 164990kg for the direct method. Solving the same problem with the direct method with a fixed bang(1)-bang(0) control structure gives a close solution (without the singular arc of course), with a criterion of m 1 (t f ) = 164442kg. It is comforting to observe that the solution with the singular arc is indeed slightly better than the forced bang-bang one, even if the difference is quite small. It seems reasonable to assume that for a different set of parameters that would give a solution with a longer singular arc, the gain compared to the bangbang solution would be more significant. Indeed, we can perform a continuation directly on the solution with a singular arc and further increase the value of the parameter S r . We observe that the solutions exhibit longer singular arcs when S r increases, with a more significant gain with respect to the corresponding solution (Fig. 7).

S

Conclusion

This study indicates that while the original problem for an Ariane 5 launcher seems to admit a simple bang-bang optimal solution, slightly modifying some parameters such as the reference area and specific impulse of the launcher gives an optimal solution with a singular arc. We also introduced a new way of computing the singular control when its analytic expression from the time derivatives of the switching function is not available (here due to the presence of tabulated data in the physical model). We plan to extend this work to the study of a winged launcher while taking into account the lift force and dynamic pressure constraint.

Fig. 1 Fig. 2

 12 Fig.1Complete flight with full thrust -Controls and masses Fig.2shows the evolution of the altitude and speed of the launcher during the flight. The change of the propulsion system at the beginning of the second and third phases is clearly visible on the speed graph. These graphs are consistent with the typical GTO flight described, for instance, in the Ariane 5 User's Manual from Arianespace.

Fig. 3 3D

 3 Fig.3 3D Goddard problem -Approximate and exact singular control

Fig. 6

 6 Fig.6 Solution with singular arc for λ = 1

Fig. 7

 7 Fig.7 Singular arcs for increasing values of S r

Acknowledgement

This study was supported by Centre National d'Etudes Spatiales, Paris.


M Soulier

T Maho

C Muja

J Lo

Ph Guillot

Non-thermal electron cyclotron resonance (ECR) air plasmas characteristics and preliminary results in tomato seeds germination

Keywords: ECR plasma source, OES, Langmuir probe, seed germination

Non-equilibrium plasmas have been studied to be used in agricultural applications. Synthetic air plasmas were generated at 1 Pa in an 80 litres vacuum chamber by an ECR coaxial microwave source at 2.45 GHz. Plasma characteristics were determined by OES and Langmuir probe measurements. Tomato seeds adding showed no influence on plasma characterizations. A 12% increase in germination rate was measured after 14 days following a plasma treatment of 10 min at 50 W.

Introduction

Non-thermal plasmas are emerging alternative treatments relative to seed conservation and improvement on their germination [1]. Nowadays, chemical products are used by industrial agriculture. They imply water and soil pollution which cause biodiversity and human health troubles. Non-thermal plasmas are advantageous because of their composition: heavy particles (ions and neutrals) induce a low-temperature preserving the seeds integrity; and energetic electrons produce reactive species and UVphotons impacting on germination and seed conservation. Different low-pressure and atmospheric pressure plasma sources have been tested on many seeds in this field [2].

In this work, an investigation into the low-pressure plasma action on seeds germination has been carried out with an electron cyclotron resonance microwave plasma source. The approach used to define the plasma treatment conditions and to better understand the mechanisms involved in these processes consisted in different steps. Firstly, a plasma parametric study without seeds was performed to determine the evolution of plasma temperatures (electron, vibration and rotation temperatures) and electronic densities. Secondly, seeds were placed inside the vacuum chamber to observe their influence on plasma characteristics. Finally, the source efficiency on modification of tomato seeds was highlighted with germination experiments.

Experimental setup

Air plasmas were generated at 1 Pa in an 80 litres vacuum chamber presented in Fig. 1 and Fig. 2 where vacuum pumping, gas monitoring and plasma diagnostics were set up through different vacuum passages.

The electron cyclotron resonance coaxial microwave source (ECRsee Fig. 1) powered by a 2.45 GHz solidstate generator was designed to sustain microwave plasmas over several decades of pressure varying from 10 -2 Pa to a few Pa (Aura-wave -SAIREM SAS). This configuration ensures minimal power-losses and perfectly matched impedance from 0 W to 200 W injected power. The ECR source was positioned at the top centre of the reactor through a KF40 passage as shown in Fig. 1 and Fig. 2. At the bottom of the vacuum chamber, a radiofrequency capacitively coupled plasma electrode (CCP) could be added to the ECR plasma source to homogenize the discharge. This electrode is a 30 cm diameter disk feed by a 13.56 MHz generator. The impedance matching system is tunable with two position parameters. Plasma diagnostics were conducted along the reactor central axis at three different height positions (Y = 6 cm, 14 cm and 22 cm from the top). For each height position, diagnostics were also performed radially where X = 0 cm correspond to the central position.

The experimental conditions for the plasma physical characterization are resumed in Table 1. All measurements were realized at 1 Pa with injected microwave power varied between 25 W to 150 W. Axial (Y): 6, 14, 22 cm from the top Radial (X): 0 to 20 cm from the centre Species identification has been conducted on spectra obtained from Ocean Optics HR2000+ optical emission spectrometer measurements. Acquired spectra were calibrated on relative intensity in order to follow the evolution of band emitted intensities as a function of microwave power. 2 shows the nitrogen and oxygen compounds identified from all spectra conditions, 𝑁 2 emissions are predominant and will be used for vibrational and rotational temperature determination. Table 2. Species identified in the ECR plasma.

Emitted species

Wavelength range (nm) 𝑁 2 : 𝐶 3 Π 𝑢 -𝐵 3 Π 𝑔 (SPS) 296.2 -464.9 𝑁 2 : 𝐵 3 Π 𝑔 -𝐴 3

844.6

In order to estimate the rotational and vibrational temperatures, a higher resolution spectrometer (Horiba Jobin Yvon iHR320) has been used. This spectrometer was equipped with an 1800 g/mm grating blazed at 400 nm and with a 0.3 nm resolution. Spectra from 365 nm to 385 nm were collected to measure the 𝑁 2 (𝐶 -𝐵) emission bands at 370.9, 375.4 and 380.4 nm. Prior to measurements, intensity relative calibration has been done with a Deuterium-Halogen calibration light source (DH3plus Ocean Optics). These spectra were fitted with Specair software to determine the vibration and rotation temperatures (resp. TVIB and TROT). The baseline was adjusted, and the measurement data were normalized to be compared with the simulation theoretical spectra. The vibrational temperature was obtained through the 375.4 nm and 380.4 nm maximum intensities ratio. Once TVIB was known, a focus on the 380.4 nm band slope allowed TROT estimation.

A linearly driven double Langmuir probe Impedans ALP system has been used to calculate electron temperatures and densities. The linear translation system permitted the collect of intensity-voltage characteristics from X = 0 cm to 20 cm every 5 cm (see Fig. 2). Electron temperatures and densities were calculated by IVcharacteristics interpretation. Electron temperature 𝑇 𝑒 was determined by equation 1 and electron density 𝑛 𝑒 by equation 2 as following. The diagnostics on plasma spectra, densities and different temperatures (TVIB, TROT and Te) were also replicated with about 500 tomato seeds on the grid inside the reactor to estimate eventual plasma parameters modification on plasma characteristics.

Evaluation on germination rate enhancement by plasma was performed on tomato (Solanum lycopersicum) seeds. Each sample consisted in about 50 seeds exposed to ECR plasmas at different power input and exposition time related in Table 3. After plasma treatment, seeds were stocked between two imbibed towel papers disposed in a 20°C dark chamber during germination experiments. The percentage of germinated seeds in a sample was estimated at 5 days and 14 days after treatment as indicated in the standard rules for germination assessment [3] and compared with the control sample (without plasma treatment). Table 3. Plasma conditions for seed treatments.

Plasma characterization results

Vibration and rotation temperatures as a function of microwave power are presented in Fig. 3. Measurements were taken at the three axial positions. The spectrometer resolution of 0.3 nm allows best-fitted spectra with standard deviations around ± 100 K for TVIB and around ± 50 K for TROT. The closer the spectrometer is to the source, the higher is the vibration and rotation temperatures. TVIB ranges from 5000 K to 8000 K and globally increase with power and TROT remained constant around 400 K. Triangular interpolation has been performed on Langmuir probe measurements for all positions specified in Table 1 in order to issue cartographies of electron temperature and electronic density in the vacuum chamber. The obtained cartographies for electronic temperature are shown in Fig. 4 whereas Fig. 5 corresponds to electronic density cartographies.

Electron temperature globally increases with power between 0.7 to 1.7 eV (i.e. 8000 to 20000 K). The higher temperature zone is found around the source (see Fig. 4), in the principal ionisation area. Although the structure is slightly different, these values are consistent with the ones found in the literature [4]. Fig. 4. Electron temperature cartographies with X, the radial position (0 to 20 cm from the centre) and Y, the axial one (6, 14, 22 cm from the top).

Electronic densities obtained inside the plasma generation zone are in the range of 10 13 cm -3 exceeding by far the microwave critical density which is defined at 7.6×10 10 cm -3 for 2.45 GHz excitation (see Fig. 5). It decreases however with distance and becomes radially uniform for Y = 15 cm (and further), corresponding to seeds position inside the reactor.

Preliminary results in seed germination

Germination rates evaluation at different treatment times was performed as presented on Fig. 7 for a 50 W ECR plasma treatment. The dark points correspond to control sampleseeds unexposed to plasma. On the 5 th and 14 th days after plasma treatments, control seeds' germination rates are lower than plasma-exposed seeds' germination rates. Furthermore, certain exposure durations are more favourable than others in order to improve germination rates. Compared to control sample, 10-and 20-minutes plasma treatment induce an increase about 12 % in germination rates against 4 % raise in the case of 5 and 30 min of plasma treatment at 50 W. In Fig. 8, 10-minutes treatment time was used as a benchmark in order to evaluate microwave power influence on germination rates. As in the case of Fig. 7, increases in germination rates thanks to plasma treatments should be noted. On the 14 th day after plasma treatment, 12 % increase in seed germination was obtained with 50 W and 150 W against 7 % with 100 W injected power. Moreover, the germination index [5] which takes into account germination relative to time, is higher on the 5 th day after plasma treatment with higher power input.

Conclusion

Air plasmas generated at 1 Pa by a microwave ECR source have been characterized in order to optimize plasma treatments for tomato seeds germination. A study of plasma temperatures as a function of microwave power has been conducted from 25 W to 150 W. Electron temperatures are between 0.7 and 1.7 eV (from 8000 to 20000 K). As expected on non-thermal plasmas, they are higher than the vibration temperatures (from 5000 K to 8000 K) and the rotation ones (around 400 K). Electron density cartographies show a maximum density in the range of 10 13 cm -3 inside the plasma generation area and a decrease far from the source.

Tomato seeds presence inside the reactor does not modify plasma characteristics. First experiments on tomato seeds showed a promising 12% increase in germination rate after 14 days following a 50 W -10 min plasma treatment.

The parametric studies undertaken indicated that this germination rate could be increased by choosing optimal microwave power-treatment time values.

Fig. 1 .

 1 Fig. 1. Aura-wave ECR plasma source (left) and the sustained discharge at 1 Pa in air plasma with 150 W power input (right).

Fig. 2 .

 2 Fig. 2. Schematic diagram of the vacuum chamber and the plasma diagnostics positions.

  𝑘 𝐵 : Boltzmann constant 𝑚 𝑒 and 𝑒: electron mass and charge 𝐴 𝑝 : probe area (9.4 µm²) 𝐼 𝑒 : electron current (probe current subtracted by ion current) 𝑉 𝑝 : plasma potential (probe potential when the IV-first derivative is maximal) 𝑉 𝑓 : floating potential (probe potential when there is no probe current)

  cm from the top (onto the grid) Seeds quantity about 50 seeds (about 0.1 g) Seeds treatment time 0 (control), 5, 10, 20 and 30 min

Fig. 3 .

 3 Fig. 3. TVIB (top) and TROT (bottom) determined by the 𝑁 2 (𝐶 -𝐵) emission bands fitted with Specair software.

Fig. 5 .

 5 Fig. 5. Electron density cartographies with X, the radial position (0 to 20 cm from the centre) and Y, the axial position (6, 14, 22 cm from the top).

Fig. 6 .

 6 Fig. 6. Electron density radial distribution at 50 W without (full) and with (empty) seeds inside the reactor.

Fig. 6

 6 Fig. 6 is an example of the seeds influence on plasma characteristics. It represents electron density radial

Fig. 7 .

 7 Fig. 7. Tomato seeds germination rates at 5 and 14 days after 50 W plasma treatmentstreatment time varied as indicated in legend.

Fig. 8 .

 8 Fig. 8. Tomato seeds germination rates at 5 and 14 days after 10 min exposed to plasmatreatment power varied as indicated in legend.

Table 1 .

 1 Plasma experimental conditions.

	Plasma

source ECR Gaz Synthetic air 4.0 Flow 50 sccm Pressure 1 Pa Injected power 25 W to 150 Wstep: 25 W Diagnostic positions

Acknowledgements

This study was supported by the Occitanie region, France. The authors would like to acknowledge SAIREM SAS and the Réseau des Plasmas Froids for their material support (DH3-plus Ocean Optics calibration light source).