Jiamin Zhu

Emmanuel Trélat

Max Cerf

Minimum time control of the rocket attitude reorientation associated with orbit dynamics

Keywords: Coupled attitude orbit problem, optimal control, Pontryagin maximum principle, shooting method, continuation, chattering arcs

In this paper, we investigate the minimal time problem for the guidance of a rocket, whose motion is described by its attitude kinematics and dynamics but also by its orbit dynamics. Our approach is based on a refined geometric study of the extremals coming from the application of the Pontryagin maximum principle. Our analysis reveals the existence of singular arcs of higher-order in the optimal synthesis, causing the occurrence of a chattering phenomenon, i.e., of an infinite number of switchings when trying to connect bang arcs with a singular arc.

We establish a general result for bi-input control-affine systems, providing sufficient conditions under which the chattering phenomenon occurs. We show how this result can be applied to the problem of the guidance of the rocket. Based on this preliminary theoretical analysis, we implement efficient direct and indirect numerical methods, combined with numerical continuation, in order to compute numerically the optimal solutions of the problem.

Introduction

The optimal control of orbit transfer (see, e.g., [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] Bryson | Applied optimal control: optimization, estimation and control[END_REF][START_REF] Marec | Optimal Space trajectories[END_REF]) and attitude reorientation (see, e.g., [START_REF] Bilimoria | Time-optimal three-axis reorientation of a rigid spacecraft[END_REF][START_REF] Seywald | Singular control in minimum time spacecraft reorientation[END_REF][START_REF] Shen | Time-optimal control of axisymmetric rigid spacecraft using two controls[END_REF]) for spacecrafts have been extensively studied in the past few decades. The optimal control problem of orbit transfer focuses mostly on how to move the spacecraft from one orbit or point to another orbit or point by using minimum energy, while the optimal control problem of attitude reorientation is mainly devoted to determine how to change the pointing direction of the spacecraft in minimum time. In the existing literature, these two optimal control problems are considered separately in general. From the engineering point of view, for most satellites, it is appropriate to design separately the control laws for the orbit movement and for the attitude movement. However, for the rockets, the trajectory is controlled by its attitude angles: the way to make the rocket follow its nominal trajectory is to change its attitude angles, and therefore it is desirable to be able to determine the optimal control subject to the coupled dynamical system. Though the control of the coupled problem was also studied in many previous works (see, e.g., [START_REF] Lara | Dynamic Stabilization of L2 Periodic Orbits Using Attitude-Orbit Coupling Effects[END_REF][START_REF] Guzzetti | Coupling Attitude and Orbital Motion of Extended Bodies In The Restricted Circular 3-Body Problem: A Novel Study On Effects And Possible Exploitations[END_REF][START_REF] Knutson | Coupled Orbit and Attitude Dynamics for Spacecraft Comprised of Multiple Bodies in Earth-Moon Halo Orbits[END_REF]), it does not seem that the problem has been investigated in the optimal control framework so far.

In this paper, we consider the time minimum control of the attitude reorientation coupled with the orbit dynamics of a rocket, denoted in short (MTCP). The chattering phenomenon that may occur according to the terminal conditions under consideration, makes in particular the problem quite difficult. Chattering means that the control switches an infinite number of times over a compact time interval. Such a phenomenon typically occurs when trying to connect bang arcs with a higher-order singular arc (see, e.g., [START_REF] Fuller | An Optimum Non-Linear Control System[END_REF][START_REF]Chattering arcs and chattering controls[END_REF][START_REF] Zelikin | Theory of Chattering Control: with applications to Astronautics[END_REF][START_REF] Zhu | Planar tilting maneuver of a spacecraft: singular arcs in the minimum time problem and chattering[END_REF]). In [START_REF] Zhu | Planar tilting maneuver of a spacecraft: singular arcs in the minimum time problem and chattering[END_REF], we studied the planar version of (MTCP), where the system consists of a single-input control-affine system, and we established as well the occurrence of a chattering phenomenon and that the chattering extremals are locally optimal in C 0 topology. * A second important difficulty in (MTCP) is due to the coupling of the attitude movement with the orbit dynamics. Indeed the system contains both slow (orbit) and fast (attitude) dynamics. This observation will be particularly important in order to design appropriate numerical approaches.

In order to analyze the extremals of the problem, we use geometric optimal control theory (see [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Schättler | Geometric Optimal Control: Theory, Methods and Examples[END_REF][START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF]). The Pontryagin maximum principle and the geometric optimal control, especially the concept of Lie bracket, will be used in this paper in order to establish an existence result of the chattering phenomenon. More precisely, based on the Goh and generalized Legendre-Clebsch conditions, we prove that there exist optimal chattering arcs when trying to connect a regular arc with a singular arc of order two.

There exist various numerical approaches to solve an optimal control problem. The direct methods (see, e.g., [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF]) consist of discretizing the state and the control and thus of reducing the problem to a nonlinear optimization problem (nonlinear programming) with constraints. The indirect methods consist of numerically solving a boundary value problem obtained by applying Pontryagin maximum principle (PMP, see [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]), by means of a shooting method. There exist also mixed methods that discretize the PMP necessary conditions and use then a large-scale optimization solver (see, e.g., [START_REF] Bérend | An interior-point approach to trajectory optimization[END_REF]). Since these numerical approaches are not easy to initialize successfully, it is required them to combine with other theoretical or numerical approaches (see the survey [START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF]). Here, we will use numerical continuation, which has proved to be very powerful tool to be combined with the PMP. For example, in [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF][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], the continuation method is used to solve difficult orbit transfer problems.

However, due to the chattering phenomenon, numerical continuation combined with shooting cannot give an optimal solution to the problem for certain terminal conditions for which the optimal trajectory contains a singular arc of higher-order. In that case, we propose sub-optimal strategies by using direct methods computing approximate piecewise constant controls. It is noticeable that our indirect approach can also be adapted to generate sub-optimal solutions, by stopping the continuation procedure before its failure due to chattering. This approach happens to be faster than the direct approach, and appears as an interesting alternative for practice.

From the engineer point of view, the theoretical analysis as well as the numerical strategies and the way to design them (in particular, the design of the problem of order zero) are strongly based on the fact that the orbit movement is much slower than the attitude movement.

The paper is organized as follows. In Section 2, we describe the mathematical model of the system consisting of the attitude dynamics, of the attitude kinematics, and of the orbit dynamics. In Section 3, we recall the Pontryagin maximum principle and some higher necessary conditions of optimality (Goh and generalized Legendre-Clebsch conditions) for bi-input control affine systems, where the two-dimensional control takes its values in the Euclidean unit disk. Based on these necessary conditions of optimality, we establish a result on the existence of the optimal chattering extremals. In Section 4, we analyze the regular and singular extremals of the problem (MTCP). For the regular extremals, we classify the switching points and state some useful properties. For the singular extremals (which are of order two), we show that the chattering phenomenon occurs for the problem (MTCP) by using the results given in the previous section. In Section 5, we propose a numerical approach to solve the problem (MTCP) by implementing numerical continuation combined with shooting. Numerical results are given in Section 6.

Model and problem statement

The problem is to control the attitude movement coupled with the orbit dynamics in the launching ascent stage for a rocket. In this paper, we take the system parameters of the rocket Ariane 5. In order to keep the stability of the rocket along the flight, the attitude maneuver should be moderate, i.e., at most ±20 degrees, and then it is possible to use Euler angles to model the attitude of the engine. In this section, we first define the coordinates systems, and then we give the equations of the attitude dynamics, of the attitude kinematics and of the orbit dynamics. The model consists of eight ordinary differential equations: three for the components of the velocity vector, three for the Euler angles and two for the components of the angular velocity vector.

Coordinate systems

Throughout the paper, we make the following assumptions:

• The Earth is a sphere and is fixed in the inertial space, i.e., the angular velocity of the Earth is zero, which means that ω ei = 0.

• The position of the rocket remains the same during the maneuver of the rocket.

• The rocket is an axial symmetric cylinder.

• The aero-dynamical forces are zero.

• The rocket engine cannot be shut off during the flight and the module of the thrust force is constant, taking its maximum value, i.e., T = T max .

The unit single-axis rotation maps R i (σ): R → R 3×3 , for σ ∈ R, i = x, y, z are defined by

R x (σ) =   1 0 0 0 cos σ sin σ 0 -sin σ cos σ   , R y (σ) =   cos σ 0 -sin σ 0 1 0 sin σ 0 cos σ   , R z (σ) =   cos σ sin σ 0 -sin σ cos σ 0 0 0 1   .
For a given vector e ∈ R 3 , taking R i (σ) e means to rotate the vector e with respect to the axis i by an angle of σ. With this definition, we next introduce the coordinate frames that will be used throughout the paper.

The Earth frame S E = (x E , ŷE , ẑE ) is fixed around the center of the Earth O. The axis ẑE points to the North pole, and the axis xE is in the equatorial plan of the Earth pointing to the equinox.

The launch frame S R = (x R , ŷR , ẑR ) is fixed around the launch point O R (where the rocket is launched). The axis xR is normal to the local tangent plane, pointing to the launch direction (here we assume that the rocket is vertically launched, i.e., the launch direction is perpendicular with the local tangent plane), and the axis ẑR points to the North. As shown in Figure 1 (a), the launch frame is derived from the Earth frame by two ordered unit single-axis rotations R z (ℓ R ) and R y (-L R ),

S E Rz(ℓR) ----→ • Ry(-LR) ------→ S R
where ℓ R and L R are the longitude and latitude of the launch point, respectively.

The body frame S b = (x b , ŷb , ẑb ) is defined as follows. The origin of the frame O b is fixed around the mass center of the rocket, the axis ẑb is along the longitudinal axis of the rocket, and the axis xb is in the cross-section (see Figure 2 (a)). The body frame can be derived by three ordered unit single-axis rotations from the launch frame, as shown in Figure 1 where θ is the pitch angle, ψ is the yaw angle and φ is the roll angle. Therefore, the transformation matrix from S R to S b is

(b), S R Ry(θ) ----→ • Rx(ψ) ----→ • Rz (φ) ----→ S b (a) (b)
L bR =R z (φ)R x (ψ)R y (θ) =   cos θ cos φ + sin θ sin ψ sin φ cos ψ sin φ -sin θ cos φ + cos θ sin ψ sin φ -cos θ sin φ + sin θ sin ψ cos φ cos ψ cos φ sin θ sin φ + cos θ sin ψ cos φ sin θ cos ψ -sin ψ cos θ cos ψ   , (1)
and the transformation matrix from

S b to S R is L Rb = L -1 bR = L ⊤ bR .

Attitude dynamic equations

The attitude dynamics are written in vectorial form in the body frame S b as

d dt (I ω) b = -( ω) b ∧ (I ω) b + ( M ) b , (2)
where I is the inertia matrix, ω is the absolute angular velocity vector, i.e., the angular velocity of the rocket with respect to the inertial space, and M is the control torque introduced by the rocket thrust. The index (•) b means that the vectors are expressed in the body frame S b .

Setting (I) b = diag(I x , I y , I z ), ( ω) b = (ω x , ω y , ω z ) ⊤ and ( M ) b = (M x , M y , M z ) ⊤ , (2) gives      I x ωx = (I y -I z )ω y ω z + M x , I y ωy = (I z -I x )ω x ω z + M y , I z ωz = (I x -I y )ω x ω y + M z . (3)
The control torque M is the cross product of the thrust vector T and of its moment arm L. The moment arm is the vector from the center of mass O b to the force acting point O F , given here by ( L) b = (0, 0, -l) ⊤ , where l is the distance between the points O b and O F . Moreover, as shown in Figure 2 (a), the thrust force vector is

T = (-T sin µ cos ζ, -T sin µ sin ζ, T cos µ) ⊤ ,
where

T = T max , µ ∈ [0, µ max ], and ζ ∈ [-π, π]. The control torque is then ( M ) b = ( L) b ∧ ( T ) b = (-T l sin µ sin ζ, T l sin µ cos ζ, 0) ⊤ .
By assumption, the rocket is axial symmetric, and hence According to the parameters of the rocket engine, µ max is less than 10 degrees and thus the error between sin µ and µ is less than 0.5%. Therefore, in the model we make the approximation sin µ ≃ µ and we define u 1 = μ cos ζ and u 2 = μ sin ζ and μ = µ/µ max with b = bµ max . Hence ωx = -bu 2 , ωy = bu 1 .

I x = I y . Assume that ω z (0) = 0,
(4)

Attitude kinematics equations

Since ω is the angular velocity vector of the rocket with respect to the inertial space, it is equal to the sum of the angular velocity ω bR of the rocket with respect to the launch frame, of the angular velocity ω RE of the launch frame with respect to the Earth frame, of the angular velocity ω Ee of the Earth frame with respect to the Earth, and of the angular velocity ω ei of the Earth with respect to the inertial space. According to the assumptions and definitions of the frames, it is easy to see that the last three terms are zero, and thus ω = ω bR . Therefore, based on the definition of the body frame, the relationship between angular velocity and Euler angles are

  ω x ω y ω z   = L bR   0 θ 0   +   ψ cos φ -ψ sin φ 0   +   0 0 φ  ,
where L bR is given by [START_REF] Allgower | Introduction to numerical continuation methods SIAM[END_REF]. Then the equations of the attitude kinematics are θ = (ω x sin φ + ω y cos φ)/ cos ψ, ψ = ω x cos φ -ω y sin φ, φ = (ω x sin φ + ω y cos φ) tan ψ. [START_REF] Bilimoria | Time-optimal three-axis reorientation of a rigid spacecraft[END_REF] Therefore, the two equations of (4) and the three equations of (5) describe the attitude movement. Note that, if ψ = π/2 + kπ, k ∈ N, then the Euler angles defined above are not well defined (usual singularities of the Euler angles). We assume in this paper that the maneuvers are small enough, so that these singularities will not be encountered. †

Orbit dynamics equations

The equation of the orbit dynamics in vectorial form is

d(m V ) R dt = (m g) R + L Rb ( T ) b + m( ω) R ∧ ( V ) R -2m( ω ei ∧ V ) R -m( ω ei ∧ ( ω R ∧ r)) R , (6)
where the notation (•) R means that the vector is expressed in the launch frame S R . The vector V is the velocity of the rocket with respect to the Earth frame S E , and its components in the launch frame are v x , v y and v z . The vector ( g) R = (g x , g y , g z ) ⊤ can be approximated by ( g) R ≈ (-g 0 , 0, 0) ⊤ , where g 0 is a real number representing the standard gravity (g 0 = 9.8). Here, m is the mass of the rocket, and it is assumed to be constant. The third term of the right-hand side of ( 6) is due to the rotation of the reference frame, the fourth term is the Coriolis force and the last term is the centrifugal force. These three terms are equal to zero because the launch frame is fixed in the inertial frame.

Due to the fact that the control angle µ is very small in practice (physical constraints imposed by the rocket engine), we assume that the thrust force is along the body longitudinal axis. According to the previous assumptions, the equation of the orbit dynamics (6) becomes vx = a sin θ cos ψ + g x , vy = -a sin ψ + g y , vz = a cos θ cos ψ + g z ,

where a = T max /m is constant. Note that this additional assumption is made also because the attitude of the rocket is controlled by only a part of the rocket engines, and so the total thrust remains almost parallel to the rocket longitudinal axis.

Minimum time control problem (MTCP)

Model. The system (4)-( 5)-( 7) has two control inputs u 1 and u 2 , and we obtain the system vx = a sin θ cos ψ + g x , vy = -a sin ψ + g y , vz = a cos θ cos ψ + g z , θ = (ω x sin φ + ω y cos φ)/ cos ψ, ψ = ω x cos φ -ω y sin φ, φ = (ω x sin φ + ω y cos φ) tan ψ,

ωx = -bu 2 , ωy = bu 1 . (8)
Defining the state variable x = (v x , v y , v z , θ, ψ, φ, ω x , ω y ), we write the system (8) as the bi-input control-affine system ẋ

= f (x) + u 1 g 1 (x) + u 2 g 2 (x), (9)
where the controls u 1 and u 2 satisfy the constraint u 2 1 + u 2 2 1, and the vector fields f , g 1 and g 2 are defined by

f = (a sin θ cos ψ + g x ) ∂ ∂v x + (-a sin ψ + g y ) ∂ ∂v y + (a cos θ cos ψ + g z ) ∂ ∂v z + (ω x sin φ + ω y cos φ)/ cos ψ ∂ ∂θ + (ω x cos φ -ω y sin φ) ∂ ∂ψ + tan ψ(ω x sin φ + ω y cos φ) ∂ ∂φ , g 1 = b ∂ ∂ω y , g 2 = - b ∂ ∂ω x . (10)
occur. Moreover, even though a singularity would occur, we could handle the situation by choosing a new launch frame S ′ R in which the new yaw angle is close to zero, i.e., ψ ′ (0) ≈ 0, ψ ′ (t f ) ≈ 0.

Terminal conditions and system parameters. Let v x0 , v y0 , v z0 , θ 0 , ψ 0 , φ 0 , ω x0 , ω y0 , θ f , ψ f , φ f , ω x f and ω y f be real numbers. The initial conditions are fixed to

v x (0) = v x0 , v y (0) = v y0 , v z (0) = v z0 , θ(0) = θ 0 , ψ(0) = ψ 0 , φ(0) = φ 0 , ω x (0) = ω x0 , ω y (0) = ω y0 . (11)
The desired final velocity is required to be parallel to the body axis ẑb , according to ( V (t f )) R ∧ (ẑ b (t f )) R = 0, and therefore, the constraints on the final conditions are

v z f sin ψ f + v y f cos θ f cos ψ f = 0, v z f sin θ f -v x f cos θ f = 0, θ(t f ) = θ f , ψ(t f ) = ψ f , φ(t f ) = φ f , ω x (t f ) = ω x f , ω y (t f ) = ω y f . (12)
Note that the parallel condition on the final velocity is due to the fact that most rockets are planned to maintain a zero angle of attack along the flight. This condition refers to the first two equalities of [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF]. The angle of flight, when the air wind is set to zero, is defined as the angle between the velocity and the rocket body axis.

Minimum Time Control Problem (MTCP). We set

x 0 = (v x0 , v y0 , v z0 , θ 0 , ψ 0 , φ 0 , ω x0 , ω y0 ) ∈ R 8
, and we define the target set (submanifold of R 8 )

M 1 ={(v x , v y , v z , θ, ψ, φ, ω x , ω y ) ∈ R 8 | v z sin ψ f + v y cos θ f cos ψ f = 0, v z sin ψ f + v y cos θ f cos ψ f = 0, θ = θ f , ψ = ψ f , φ = φ f , ω x = ω x f , ω y = ω y f }.
The minimum time control problem (MTCP) consists of steering the bi-input control-affine system (9) from x(0) = x 0 to the final target M 1 in minimum time t f , with controls satisfying the constraint

u 2 1 + u 2 2 1.
3 Some general results for bi-input control-affine systems

In this section, we focus on the chattering phenomenon for bi-input control-affine systems with control constraints and with commuting controlled vector fields. The results that we are going to give are general and will be used in the next section to analyze the problem (MTCP).

We consider the following general framework. Let M be a smooth manifold of dimension n, let x 0 ∈ M be arbitrary, and let M 1 be a submanifold of M . We consider on M the minimal time control problem

         min t f , ẋ(t) = f (x(t)) + u 1 (t)g 1 (x(t)) + u 2 (t)g 2 (x(t)), u = (u 1 , u 2 ) u(t) 2 = u 1 (t) 2 + u 2 (t) 2 1, x(0) = x 0 , x(t f ) ∈ M 1 , t f 0 free, (13)
where f , g 1 and g 2 are smooth vector fields on M . According to classical results (see, e.g., [START_REF] Cesari | Optimization -theory and applications[END_REF][START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF]), there exists at least one optimal solution (x(•), u(•)), defined on [0, t f ].

Application of the Pontryagin maximum principle

According to the Pontryagin maximum principle (in short, PMP, see [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]), there must exist an absolutely continuous mapping p(•) defined on [0, t f ] (called adjoint vector), such that p(t) ∈ T * x(t) M (cotangent space) for every t ∈ [0, t f ], and a real number p 0 0, with (p(•), p 0 ) = 0, such that ẋ(t) = ∂H ∂p (x(t), p(t), p 0 , u(t)), ṗ(t) = -∂H ∂x (x(t), p(t), p 0 , u(t)), almost everywhere on [0, t f ], where H(x, p, p 0 , u) = h 0 (x, p) + u 1 h 1 (x, p) + u 2 h 2 (x, p) + p 0 is the Hamiltonian of the optimal control problem [START_REF] Cesari | Optimization -theory and applications[END_REF]. Here, we have set h 0 (x, p) = p, f (x) , h 1 (x, p) = p, g 1 (x) , and h 2 (x, p) = p, g 2 (x) . The maximization condition of the PMP yields, almost everywhere on [0, t f ],

u(t) = (h 1 (t), h 2 (t)) h 1 (t) 2 + h 2 (t) 2 = Φ(t) Φ(t) , (14)
whenever Φ(t) = (h 1 (t), h 2 (t)) = (0, 0). We call Φ (as well as its components) the switching function. Note that Φ is continuous. Here and throughout the paper, we denote by h i (t) = h i (x(t), p(t)), i = 0, 1, 2 with a slight abuse of notation. Moreover, we have the transversality condition p(t f ) ⊥ T x(t f ) M 1 , where T x(t f ) M 1 is the tangent space to M 1 at the point x(t f ), and, the final time t f being free and the system being autonomous, we have also h 0 (x(t), p(t)) + Φ(t) + p 0 = 0, ∀t ∈ [0, t f ].

The quadruple (x(•), p(•), p 0 , u(•)) is called an extremal lift of x(•). An extremal is said to be normal (resp., abnormal) if p 0 < 0 (resp., p 0 = 0).

We say that an arc (restriction of an extremal to a subinterval I) is regular if Φ(t) = 0 along I. Otherwise, the arc is said to be singular. Note that a singular extremal may be both normal or abnormal. We will see in Section 4.2 that the singular extremals of the problem (MTCP) must be normal.

A switching time is a time t at which Φ(t) = (0, 0), that is, both h 1 and h 2 vanish at time t. Note that, at a switching time, the expression ( 14) is not well defined. If only one component of Φ(t) vanishes, then the optimal control is still regular. An arc that is a concatenation of an infinite number of regular arcs is said to be chattering. The chattering arc is associated with a chattering control that switches an infinite number of times, over a compact time interval. A junction between a regular arc and a singular arc is said to be a singular junction.

Computation of singular arcs, and necessary conditions for optimality

We next define the order of a singular control, since it is important to understand and explain the occurence of chattering. This concept is related to the way singular controls are computed, and since it is a bit technical to define, we start with a preliminary quite informal discussion. Here and throughout the paper, we use the notation adf.g = [f, g] (Lie bracket of vector fields) and adh i .h j = {h i , h j } (Poisson bracket of Hamiltonian functions).

Preliminary informal discussion. In order to compute singular controls, the usual method is to differentiate several times the switching function, until the control appears in a nontrivial way. If Φ(t) = 0 for every t ∈ I, then h 1 (t) = h 2 (t) = 0, and, differentiating in t, we get, using the Poisson bracket, ḣ1 = {h 0 , h

1 } + u 2 {h 2 , h 1 } = 0 and ḣ2 = {h 0 , h 2 } + u 1 {h 1 , h 2 } = 0 along I.
According to the Goh condition (see [START_REF] Goh | Necessary conditions for singular extremals involving multiple control variables[END_REF], see also Lemma 1), if the singular arc is optimal and the associated singular control is not saturating, then the Goh condition {h 1 , h 2 } = p, [g 1 , g 2 ](x) = 0 must be satisfied along I. Therefore we get that ḣ1 = {h 0 , h

1 } = p, [f, g 1 ](x) = 0 and ḣ2 = {h 0 , h 2 } = p, [f, g 2 ](x) = 0 along I. d dt {h 2 , {h 0 , h 1 }} = d dt {h 1 , {h 0 , h 2 }} = {h 0 , {h 1 , {h 0 , h 2 }}} = 0.
This is a new constraint along the singular arc. Note that the time derivative of this constraint is equal to zero and hence does not induce any additional constraint. Due to higher-order necessary conditions for optimality (see Lemma 1), an optimal singular control cannot appear in a nontrivial way with an odd number of derivatives, therefore we must have {h 2 , ad 2 h 0 .h 1 } = 0 and {h 1 , ad 2 h 0 .h 2 } = 0 along I. Accordingly, h

(3) i = 0, i = 1, 2, gives three additional constraints along the singular arc {h 0 , ad 2 h 0 .h 1 } = 0, {h 0 , ad 2 h 0 .h 2 } = 0, and

{h 2 , ad 2 h 0 .h 1 } = -{h 1 , ad 2 h 0 .h 2 } = 0. Differentiating the third constraint with respect to t, we have d dt {h 1 , ad 2 h 0 .h 2 } = {h 0 , {h 1 , ad 2 h 0 .h 2 }} = 0, and thus {h 2 , ad 3 h 0 .h 1 } = {ad 2 h 0 .h 1 , adh 0 .h 2 } and {h 1 , ad 3 h 0 .h 2 } = {ad 2 h 0 .h 2 , adh 0 .h 1 }.
Then, by differentiating the first two constraints with respect to t, we get

h (4) 1 = ad 4 h 0 .h 1 + u 1 {h 1 , ad 3 h 0 .h 1 } + u 2 {ad 2 h 0 .h 1 , adh 0 .h 2 } = 0, h (4) 2 = ad 4 h 0 .h 2 + u 1 {ad 2 h 0 .h 2 , adh 0 .h 1 } + u 2 {h 2 , ad 3 h 0 .h 2 } = 0. Assuming that {h i , ad 3 h 0 .h i } < 0, i = 1, 2 (generalized Legendre-Clebsch condition, see below) and that det ∆ 2 = det {h 1 , ad 3 h 0 .h 1 } {ad 2 h 0 .h 1 , adh 0 .h 2 } {ad 2 h 0 .h 2 , adh 0 .h 1 } {h 2 , ad 3 h 0 .h 2 } = 0
along I, the singular control is given by

u 1 = -(ad 4 h 0 .h 1 ){h 2 , ad 3 h 0 .h 2 } + (ad 4 h 0 .h 2 ){h 2 , ad 3 h 0 .h 1 } / det ∆ 2 , u 2 = (ad 4 h 0 .h 1 ){h 1 , ad 3 h 0 .h 2 } -(ad 4 h 0 .h 2 ){h 1 , ad 3 h 0 .h 1 } / det ∆ 2 .
We say, then, that the singular control u = (u 1 , u 2 ) is of intrinsic order two.

Precise definitions. Now, following [START_REF] Gabasov | High order necessary conditions for optimality[END_REF], let us give a precise definition of the order of a singular control.

Definition 1. The singular control u = (u 1 , u 2 ) defined on a subinterval I ⊂ [0, t f ] is said to be of order q if 1. the first (2q -1)-th time derivative of h i , k = 0, 1, do not depend on u and

d k dt k (h i ) = 0, k = 0, 1, • • • , 2q -1,
2. the 2q-th time derivative of h i , k = 0, 1, depends on u linearly and

∂ ∂u i d 2q dt 2q (h i ) = 0, det ∂ ∂u d 2q dt 2q Φ = 0, i = 1, 2, along I.
The control u is said to be of intrinsic order q if, moreover, the vector fields satisfy

[g i , ad k f.g i ] ≡ 0, k = 1, • • • , 2q -2, i = 1, 2.
The condition of a nonzero determinant guarantees that the optimal control can be computed from the 2q-th time derivative of the switching function. Note that, in the definition, it is required that the two components of the control have the same order.

We next recall the Goh and generalized Legendre-Clebsch conditions (see [START_REF] Goh | Necessary conditions for singular extremals involving multiple control variables[END_REF][START_REF] Kelley | Singular extremals[END_REF][START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF]). It is worth noting that in [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF], the following higher-order necessary conditions are given even when the components of the control u have different orders.

Lemma 1. (higher-order necessary conditions) Assume that a singular control u = (u 1 , u 2 ) defined on I is of order q. Suppose that u is optimal and is not saturating, i.e. u < 1. Then the Goh condition

∂ ∂u j d k dt k (h i ) = 0, k = 0, 1, • • • , 2q -1, i, j = 1, 2, i = j,
must be satisfied along I. Moreover, the matrix of which the (i, j)-th component

(-1) q ∂ ∂u j d 2q dt 2q (h i ), i, j = 1, 2,
is symmetric and nonpositive definite along I (generalized Legendre-Clebsch Condition).

In the problem (MTCP), as we will see, it happens that singular controls are of intrinsic order 2, and that [g 1 , g 2 ] = 0, [g 1 , [f, g 2 ]] = 0, and [g 2 , [f, g 1 ]] = 0, so that the conditions given in the above definition yield

[g 1 , [f, g 1 ]] = 0, [g 2 , [f, g 2 ]] = 0, [g 1 , ad 2 f.g 1 ] = 0, [g 2 , ad 2 f.g 2 ] = 0, p, [g 1 , ad 3 f.g 1 ](x) = 0, p, [g 2 , ad 3 f.g 2 ](x) = 0, and p, [g 1 , ad 3 f.g 1 ](x) p, [g 2 , ad 3 f.g 2 ](x) -p, [g 2 , ad 3 f.g 1 ](x) p, [g 1 , ad 3 f.g 2 ](x) = 0,
and we have the following higher-order necessary conditions, that will be used in the study of the problem (MTCP).

Corollary 1. We assume that the optimal trajectory x(•) contains a singular arc, defined on the subinterval I of [0, t f ], associated with a non saturating control u = (u 1 , u 2 ) of intrinsic order 2.

If the vector fields satisfy

[g 1 , g 2 ] = 0, [g i , [f, g j ]] = 0, for i, j = 1, 2, then the Goh condition p(t), [g 1 , adf.g 2 ](x(t)) = 0, p(t), [g 1 , ad 2 f.g 2 ](x(t)) = p(t), [g 2 , ad 2 f.g 1 ](x(t)) = 0,
and the generalized Legendre-Clebsch condition (in short, GLCC )

p(t), [g i , ad 3 f.g i ](x(t)) 0, i = 1, 2, p(t), [g 1 , ad 3 f.g 2 ](x(t)) = p(t), [g 2 , ad 3 f.g 1 ](x(t))
must be satisfied along I. Moreover, we say that the strengthened GLCC is satisfied if we have a strict inequality above, that is, p(t), [g i , ad 3 f.g i ](x(t)) < 0.

Corollary 1 follows from Lemma 1 and from the arguments developed in the previous informal discussion. It will be used in Section 4.2.

We next investigate the singular junctions for the problem [START_REF] Cesari | Optimization -theory and applications[END_REF], and the chattering phenomenon.

Chattering phenomenon

One can find in [START_REF] Mcdanell | Necessary conditions joining optimal singular and nonsingular subarcs[END_REF] some results on the junction between an optimal regular arc and an optimal singular arc, for single-control affine systems, among which a result stating that, if the singular arc is of even order and if the control is discontinuous at the junction, then the junction must be nonanalytical (meaning that the control is not piecewise analytic in any neighborhood of the junction). In [START_REF] Schättler | Geometric Optimal Control: Theory, Methods and Examples[END_REF][START_REF] Zelikin | Theory of Chattering Control: with applications to Astronautics[END_REF], it is proved that such a nonanalytical junction between a regular arc and a singular arc of intrinsic order two causes chattering (see also [START_REF] Zhu | Planar tilting maneuver of a spacecraft: singular arcs in the minimum time problem and chattering[END_REF]). When the control takes values in the unit disk, explicit analytic expressions for some optimal trajectories of linear-quadratic problems were given, e.g., in [START_REF]Chattering arcs and chattering controls[END_REF][START_REF] Zelikin | Theory of Chattering Control: with applications to Astronautics[END_REF]. However these results cannot be applied to (MTCP) because the control system is bi-input and the cost functional is the time; they are anyway a good source of inspiration to establish the results of that section. The following result is valid for general bi-input control-affine systems.

Theorem 1. Consider the optimal control problem (13). Let (x(•), p(•), p 0 , u(•)) be an optimal extremal lift on [0, t f ]. We assume that u is singular of order two along an open interval I ⊂ [0, t f ], and we denote this control by u s = (u 1s , u 2s ). We assume that u s (t) < 1 (i.e., the singular control does not saturate the constraint) and that ∂ ∂u1 d 4 dt 4 h 2 (x(t), p(t)) = 0 along I. Then the optimal control u must switch infinitely many times at the junction with the singular arc. In other words, there is a chattering phenomenon, which is due to the connection of a regular arc with a singular arc of higher-order.

Proof. Since the singular control is of order two, it follows from the definition that

∂ ∂u i d k dt k h i (x(t), p(t)) = 0, k = 0, • • • , 3, i = 1, 2, ∂ ∂u i d 4 dt 4 h i (x(t), p(t)) = 0.
Thus, we get from ∂ ∂ui d 4 dt 4 h i (x(t), p(t)) = 0 and Lemma 1 that

∂ ∂u j d k dt k h i (x(t), p(t)) = 0, k = 0, • • • , 3, i, j = 1, 2, i = j, and ∂ ∂u i d 4 dt 4 h i (x(t), p(t)) < 0, ∂ ∂u 1 d 4 dt 4 h 2 (x(t), p(t)) = ∂ ∂u 2 d 4 dt 4 h 1 (x(t), p(t)),
along the singular arc I. By assumption, we have

∂ ∂u1 d 4
dt 4 h 2 (x(t), p(t) = 0, and hence we can write h

(4) i (x(t), p(t)) = a i0 (x(t), p(t)) + u is a ii (x(t), p(t)) with a ii (x(t), p(t)) = ∂ ∂ui d 4 dt 4 h i (x(t), p(t)) < 0.
Without loss of generality, we consider a concatenation of a singular arc with a regular arc at time τ ∈ I. Assume that for some ε > 0 the control u is singular along (-ε + τ, τ ), and that, along (τ, τ + ε), the control u = (u 1 , u 2 ) is given by

u i = h i / Φ 0, i = 1, 2.
It can be easily seen from the assumption that u s < 1 that there exists at least one component of the singular control that is smaller than the same component of the regular control, i.e., u ks < u k for k = 1 or k = 2. Then, it follows that

h (4) k (τ ) = a k0 (x(τ ), p(τ )) + u k (τ )a kk (x(τ ), p(τ )) a k0 (x(τ ), p(τ )) + u ks (τ -)a kk (x(τ ), p(τ )) = h (4) k (τ -) = 0. (17)
Hence the switching function h k has a local maximum at t = τ and is nonpositive along the interval (τ, τ + ε). It follows from the maximization property of the Hamiltonian that u k 0. Thus, together with the assumption u k = 0. It follows that u j = 1 and u js < u j for j = k. Replacing k in [START_REF] Gergaud | Low thrust minimum fuel orbital transfer: an homotopic approach[END_REF] with j, we get u j 0, which gives a contradiction with u j = 1.

If, instead, we assume u i = h i / Φ 0, i = 1, 2 over (τ, τ + ε), then there must exists a control component u ks that is larger than u ks , i.e., u ks > u k , and then we obtain h

(4) k (τ ) h (4)
k (τ -) = 0, which yields u k 0, and thus u k = 0. As in the previous case, it follows that u j = -1 and u js > u j for j = k, and thus we get from h

(4) k (τ ) h (4)
k (τ -) = 0 that u j 0, which yields a contradiction. Then, if we assume u i = h i / Φ < 0 and u j = h j / Φ > 0, i, j = 1, 2, i = j, we will have either u is < u i which gives a contradiction with the sign of u i , or u is u i and u js < u j which gives a contradiction with the sign of u j . A similar reasoning can be done for regular-singular type concatenations.

Recall that the extremal is said singular if Φ(t) = h 2 1 (t) + h 2 2 (t) = 0, t ∈ I. Thus, the obtained contradiction indicates that the concatenation of a singular arc with a regular arc violates the PMP and thus there exists a chattering arc when trying to connect a regular arc with a singular arc.

Remark 1. Note that, in this result, we have assumed that u s < 1. In the (nongeneric) case where the singular control saturates the constraint, in order to get the same result we need to assume that the strengthened GLCC is satisfied at the junction point, i.e., a ii (x(τ ), p(τ )) < 0, and the control is discontinuous at the singular junction.

In addition, we have assumed that

∂ ∂u1 d 4 dt 4 h 2 (x(t), p(t)) = 0. Actually, if ∂ ∂u1 d 4
dt 4 h 2 (x(t), p(t)) = 0, then singular and regular extremals can be connected without chattering. For example, [START_REF] Gergaud | Low thrust minimum fuel orbital transfer: an homotopic approach[END_REF] gives

h (4) k (τ ) = a k0 (x(τ ), p(τ )) + u k (τ )a kk (x(τ ), p(τ )) + u m (τ )a km (x(τ ), p(τ )) a k0 (x(τ ), p(τ )) + u ks (τ -)a kk (x(τ ), p(τ )) + u m (τ )a km (x(τ ), p(τ )) = h (4) k (τ -) + a km (x(τ ), p(τ ))(u m (τ ) -u ms (τ )) = a km (x(τ ), p(τ ))(u m (τ ) -u ms (τ )).
where

a km (x(τ ), p(τ )) = ∂ ∂um d 4 dt 4 h k (x(t), p(t)), k, m = 1, 2, k = m.
In contrast to the previous reasoning, now the fact that a km (x(τ ), p(τ ))(u m -u ms ) > 0 does not raise any more a contradiction.

In the next section, we analyze the regular, singular and chattering extremals for the problem (MTCP) by using the results presented previously.

Geometric analysis of the extremals of (MTCP)

In this section, we classify the switching points by their contact with the switching surface, and we establish that the optimal singular arcs of the (MTCP), if they exist, cause chattering.

Regular extremals

Normal extremals. Here, we consider normal extremals and we take p 0 = -1. Let us consider the system (9), with the vector fields f , g 1 and g 2 defined by [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF]. Denoting the adjoint vector by p = (p vx , p vy , p vz , p θ , p ψ , p φ , p ωx , p ωy ), the adjoint equations given by the PMP are ṗvx = 0, ṗvy = 0, ṗvz = 0, ṗθ = -a cos ψ(p vx cos θ -p vz sin θ), ṗψ = a sin ψ sin θp vx + a cos ψp vy + a cos θ sin ψp vzsin ψ(ω x sin φ + ω y cos φ)/ cos

with the transversality condition

p vx (t f ) sin θ f cos ψ f -p vy (t f ) sin ψ f + p vz (t f ) cos θ f cos ψ f = 0. ( 19
)
The switching function is Φ(t) = (h 1 (t), h 2 (t)) = ( bp ωy (t), -bp ωx (t)) and is of class C 1 . The switching manifold Γ is the submanifold of R 16 of codimension two defined by

Γ = {z = (x, p) ∈ R 16 | p ωx = p ωy = 0}.
Let us fix an arbitrary reference regular extremal z(•) = (x(•), p(•)) of the problem (MTCP).

If z(•) never meets Γ, then the extremal control is given by ( 14) along the whole extremal.

If it meets Γ then there is a singularity to be analyzed. It is not even clear if the extremal flow is well defined when crossing such a point (we could lose uniqueness). Let us assume that the extremal z(•) meets Γ at some time t 0 , and we set z 0 = z(t 0 ) = (x 0 , p 0 ). Following [START_REF] Bonnard | Geometric optimal control of elliptic Keplerian orbits[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][START_REF] Kupka | Generalized Hamiltonians and optimal control: a geometric study of extremals[END_REF], we classify the regular extremals by their contact with the switching surface, i.e., if Φ (k-1) (t 0 ) = 0 and Φ (k) (t 0 ) = 0 for some k ∈ N * , then the point z 0 is said to be a point of order k. Without loss of generality, we assume that t 0 = 0.

Let us analyze the singularity occurring at points of order 1, 2, 3 and 4 for the problem (MTCP).

Points of order 1.

Lemma 2. We assume that z 0 is of order 1. Then the reference extremal is well defined in a neighborhood of t = 0, in the sense that there exists a unique extremal associated with the control u = (u 1 , u 2 ) passing through the point z 0 . The control turns with an angle π when passing through the switching surface Γ, and is locally given by

u 1 (t) = a 1 a 2 1 + a 2 2 t |t| + o(1), u 2 (t) = a 2 a 2 1 + a 2 2 t |t| + o(1), with a 1 = {h 0 , h 1 }(z 0 ), a 2 = {h 0 , h 2 }(z 0 ).
Proof. In the problem (MTCP), the vector fields g 1 and g 2 (defined by ( 10)) commute, i.e., [g 1 , g 2 ] = 0. By differentiating the switching function, we get ḣ1 = {h 0 , h 1 } and ḣ2 = {h 0 , h 2 } along the extremal, and then, since the point is of order 1, we have, locally, h 1 (t) = a 1 t + o(t) and h 2 (t) = a 2 t + o(t), with a 2 1 + a 2 2 = 0, and we also have

p ωx (t) = -a 2 t/ b + o(t), p ωy (t) = a 1 t/ b + o(t), p ωx (t) 2 + p ωy (t) 2 = a 2 1 + a 2 2 |t|/ b + o(t).
The expression of the optimal control in the lemma follows. At the crossing point, both control components change their sign, i.e., u i (0 -) = -u i (0 + ), i = 1, 2, which means that the control direction turns with an angle π when crossing Γ.

Points of order 2. Lemma 3. We assume that z 0 is of order 2. Then the reference extremal is well defined in a neighborhood of t = 0, in the sense that there exists a unique extremal associated with the control u = (u 1 , u 2 ) passing through the point z 0 . Moreover, the switching function is of class C ∞ in the neighborhood of t = 0, the control is of class C ∞ , and we have

u 1 (t) = α 1 α 2 1 + α 2 2 + o(1), u 2 (t) = α 2 α 2 1 + α 2 2 + o(1),
with α 1 = {h 0 , {h 0 , h 1 }}(z 0 ) and α 2 = {h 0 , {h 0 , h 2 }}(z 0 ).

Proof. The vector fields f , g 1 and g 2 , defined by [START_REF] Bryson | Applied optimal control: optimization, estimation and control[END_REF], are such that [g i , [f, g j ]] = 0, for i, j = 1, 2 (see Lemma 5). Then, according to the calculations done in Section 3.2, we have ḧ1 = {h 0 , {h 0 , h 1 }} and ḧ2 = {h 0 , {h 0 , h 2 }}, and hence the functions t → h i (t), i = 1, 2 are of class C 2 at 0. Locally, we have h

1 (t) = 1 2 α 1 t 2 + o(t 2 ) and h 2 (t) = 1 2 α 2 t 2 + o(t 2
). The expression of the optimal control follows and the control is continuous.

Differentiating again the switching function, we have h

(3) Lemma 5). Since the control is continuous, the switching function is at least of class C 3 at 0. Hence, locally we can write h

1 (t) = ad 3 h 0 .h 1 + u 2 {h 2 , ad 2 h 0 .h 1 } and h (3) 2 (t) = ad 3 h 0 .h 2 + u 1 {h 1 , ad 2 h 0 .h 2 }, because [g i , [f, [f, g i ]]] = 0, for i = 1, 2 (see
1 (t) = 1 2 α 1 t 2 + 1 6 β 1 t 3 + o(t 3 ) and h 2 (t) = 1 2 α 2 t 2 + 1 6 β 2 t 3 + o(t 3
), where

β i = ad 3 h 0 .h i (z 0 ) + α j / α 2 1 + α 2 2 {h j , ad 2 h 0 .h i }(z 0
) for i, j = 1, 2 and i = j. We infer that the control is at least of class C 1 at 0, with u1 (0

) = 1 6 (α 2 2 β 1 -2α 1 α 2 β 2 )/(α 2 1 + α 2 1 ) 3/2 and u2 (0) = 1 6 (α 2 1 β 2 -2α 1 α 2 β 1 )/(α 2 1 + α 2 1 ) 3/2
. We get the smoothness by an immediate induction argument: for k > 2, assuming that the switching function is of class C k and that the control is of class C k-2 , then the (k + 1)-th time derivative of the switching function can be written as h

(k+1) i = u (k-2) i term 1 + term 2 ,
where term i , i = 1, 2 are terms involving time derivatives of the control of order lower than k -2. Hence the switching function is of class C k+1 since u is of class C k-2 , and

h i = k+1 p=2 a i,p t p + o(t k+1
) where all coefficients can be computed explicitly. The (k -1)-th time derivative of the control can be computed using the coefficients a i,p , p = 2, • • • , k + 1, and hence the control is of class C k-1 . The result follows.

which leads to R 0 (α 4 -35α 2 + 24) = {h 1 , b 1 }(z 0 ) and R 0 (-10α 3 + 50α)t/|t| = 0. It follows that, R 0 = -{h 1 , b 1 }(z 0 )/126 with α ∈ {- √ 5, √ 5} and {h 1 , b 1 } < 0.
It is clear that the uniqueness of the extremal when crossing the point z 0 does not hold true anymore. The switching function Φ(t) converges to (0, 0) when t → 0, while the control switches infinitely many times when t → 0. Indeed, we will see in the following section that this situation is related to the chattering phenomenon.

Abnormal extremals. Abnormal extremals correspond to p 0 = 0 in the PMP. We suspect the existence of optimal abnormal extremals in the problem (MTCP) for certain (nongeneric) terminal conditions. In the planar version of the problem (MTCP) studied in [START_REF] Zhu | Planar tilting maneuver of a spacecraft: singular arcs in the minimum time problem and chattering[END_REF], if the optimal control switches at least two times then there is no abnormal minimizer. We expect that the same property is still true here. We are able to prove that the singular extremals of the problem (MTCP) are normal (see section 4.2), however, we are not able to establish a clear relationship between the number of switchings and the existence of abnormal minimizers as in [START_REF] Zhu | Planar tilting maneuver of a spacecraft: singular arcs in the minimum time problem and chattering[END_REF]. Thus, in our numerical simulations further, we will assume that there is at least one normal extremal for problem (MTCP) and compute it. Note moreover that Lemmas 2, 3 and 4 are also valid for abnormal extremals.

Singular and chattering extremals

Let us compute the singular arcs of the problem (MTCP). According to [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF], the singular trajectories are feedback invariants since they correspond to the singularities of the end-point mapping. This concept is related to the feedback group induced by the feedback transformation and the corresponding control systems are said to be feedback equivalent(see [9, Section 4] for details). Recall that, equivalently, a trajectory x(•) associated with a control u is said to be singular if the differential of the end-point mapping is not of full rank. The end-point mapping

E : R n × R × L ∞ (0, +∞; R) → R n is defined by E(x 0 , t f , u) = x(x 0 , t f , u), where t → x(x 0 , t, u)
is the trajectory solution of the control system associated to u such that x(x 0 , 0, u) = x 0 .

Since the singular trajectories are feedback invariants under feedback transformations, we can replace the vector fields g 1 and g 2 with g 1 = ∂ ∂ωy and g 2 = ∂ ∂ωx . It can be easily shown that the singular control computed with g 1 and g 2 coincides with the one computed with g 1 and g 2 . Let us make precise the Lie bracket configuration of the control system associated with the vector fields f , g 1 and g 2 .

Lemma 5. We have

g 1 = ∂ ∂ω y , g 2 = ∂ ∂ω x , [ g 1 , g 2 ] = 0, adf. g 1 = -cos φ/ cos ψ ∂ ∂θ + sin φ ∂ ∂ψ -tan ψ cos φ ∂ ∂φ = 0, adf. g 2 = -sin φ/ cos ψ ∂ ∂θ -cos φ ∂ ∂ψ -tan ψ sin φ ∂ ∂φ = 0, ad 2 f. g 1 = -ω x ∂ ∂φ + a(cos θ cos φ + sin θ sin φ sin ψ) ∂ ∂v x + a sin φ cos ψ ∂ ∂v y -a(cos φ sin θ -sin φ cos θ sin ψ) ∂ ∂v z , ad 2 f. g 2 = ω y ∂ ∂φ + a(cos θ sin φ -sin θ cos φ sin ψ) ∂ ∂v x -a cos φ cos ψ ∂ ∂v y -a(sin φ sin θ + cos φ cos θ sin ψ) ∂ ∂v z , [ g 1 , adf. g 1 ] = [ g 1 , adf. g 2 ] = [ g 2 , adf. g 1 ] = [ g 2 , adf. g 2 ] = 0, ad 3 f. g 1 = ω x Ω 1 / cos ψ ∂ ∂θ -ω x Ω 2 ∂ ∂ψ + ω x tan ψΩ 1 ∂ ∂φ -aω y cos ψ sin θ ∂ ∂v x + aω y sin ψ ∂ ∂v y -aω y cos ψ cos θ ∂ ∂v z , ad 3 f. g 2 = -ω y Ω 1 / cos ψ ∂ ∂θ + ω x Ω 2 ∂ ∂ψ -ω x tan ψΩ 1 ∂ ∂φ -aω x cos ψ sin θ ∂ ∂v x + aω x sin ψ ∂ ∂v y -aω x cos ψ cos θ ∂ ∂v z , [ g 1 , ad 2 f. g 1 ] = [ g 2 , ad 2 f. g 2 ] = 0, [ g 1 , ad 2 f. g 2 ] = -[ g 2 , ad 2 f. g 1 ] = ∂ ∂φ , ad 4 f. g 1 = -a(ω 2 x + ω 2 y )(cos φ cos θ + sin φ sin ψ sin θ) ∂ ∂v x -a cos ψ sin φ(ω 2 x + ω 2 y ) ∂ ∂v y + a(ω 2 x + ω 2 y )(cos φ sin θ -cos θ sin φ sin ψ) ∂ ∂v z + (ω 3 x + ω x ω 2 y ) ∂ ∂φ , ad 4 f. g 2 = -a(ω 2 x + ω 2 y )(cos θ sin φ -cos φ sin ψ sin θ) ∂ ∂v x + a cos φ cos ψ(ω 2 x + ω 2 y ) ∂ ∂v y + a(ω 2 x + ω 2 y )(sin φ sin θ + cos φ cos θ sin ψ) ∂ ∂v z + (-ω 2 x ω y -ω 3 y ) ∂ ∂φ , [ g 1 , ad 3 f. g 1 ] = -a cos ψ sin θ ∂ ∂v x + a sin ψ ∂ ∂v y -a cos ψ cos θ ∂ ∂v z -(ω x sin φ)/ cos ψ ∂ ∂θ -ω x cos φ ∂ ∂ψ -ω x sin φ tan ψ ∂ ∂φ , [ g 1 , ad 3 f. g 2 ] = -2aω y (cos θ sin φ -cos φ sin ψ sin θ) ∂ ∂v x + 2aω y cos φ cos ψ ∂ ∂v y + 2aω y (sin φ sin θ + cos φ cos θ sin ψ) ∂ ∂v z + (-ω 2 x -3ω 2 y ) ∂ ∂φ , [ g 2 , ad 3 f. g 1 ] = -2aω x (cos φ cos θ + sin φ sin ψ sin θ) ∂ ∂v x -2aω x cos ψ sin φ ∂ ∂v y + 2aω x (cos φ sin θ -cos θ sin φ sin ψ) ∂ ∂v z + (3ω 2 x + ω 2 y ) ∂ ∂φ , [ g 2 , ad 3 f. g 2 ] = -a cos ψ sin θ ∂ ∂v x + a sin ψ ∂ ∂v y -a cos ψ cos θ ∂ ∂v z -(ω y cos φ)/ cos ψ ∂ ∂θ + ω y sin φ ∂ ∂ψ -ω y cos φ tan ψ ∂ ∂φ ,
where Ω 1 = ω x cos φ -ω y sin φ and Ω 2 = ω x sin φ + ω y cos φ. Moreover, we have dim Span g 1 , g 2 , adf. g 1 , adf. g 2 , ad 2 f. g 1 , ad 2 f. g 2 = 6.

The proof of Lemma 5 was done with symbolic computations. Lemma 6. In the problem (MTCP), let us assume that (x(•), p(•), p 0 , u(•)) is a singular arc along the subinterval I, which is locally optimal in C 0 topology. Then we have u = (u 1 , u 2 ) = (0, 0) along I, and u is a singular control of intrinsic order two. Moreover, the extremal must be normal, i.e., p 0 = 0, and the GLCC a + g x sin θ cos ψ -g y sin ψ + g z cos θ cos ψ 0,

must hold along I.

Proof. Using Lemma 5, we infer from Φ = 0 and Φ = 0 that p, g 1 (x) = p ωy = 0, p, g 2 (x) = p ωx = 0,

p, adf. g 1 (x) = -p θ cos φ/ cos ψ + p ψ sin φ -p φ tan ψ cos φ = 0, p, adf. g 2 (x) = -p θ sin φ/ cos ψ -p ψ cos φ -p φ tan ψ sin φ = 0, (22)
and from Φ = 0, that p, ad 2 f. g 1 (x) = -ω x p φ + a(cos θ cos φ + sin θ sin φ sin ψ)p vx + a sin φ cos ψp vy -a(cos φ sin θsin φ cos θ sin ψ)p vz = 0, p, ad 2 f. g 2 (x) = ω y p φ + a(cos θ sin φsin θ cos φ sin ψ)p vx -a cos φ cos ψp vy -a(sin φ sin θ + cos φ cos θ sin ψ)p vz = 0,

along the interval I. Since dim Span g 1 , g 2 , adf. g 1 , adf. g 2 , ad 2 f. g 1 , ad 2 f. g 2 = 6, the six equations in ( 22)-( 23) are independent constraints along the singular arc. Therefore, writing Φ (3) = 0, we get from Corollary 1 that p, [ g 1 , ad 2 f. g 2 (x)] =p φ = 0, p, [ g 2 , ad 2 f. g 1 (x)] = -p φ = 0, p, ad 3 f. g 1 (x) =p θ ω x Ω 1 / cos ψ -p ψ ω x Ω 2 + p φ ω x tan ψΩ 1 -p vx aω y cos ψ sin θ + p vy aω y sin ψ -p vz aω y cos ψ cos θ = 0,

p, ad 3 f. g 2 (x) = -p θ ω y Ω 1 / cos ψ + p ψ ω x Ω 2 -p φ ω x tan ψΩ 1 -p vx aω x cos ψ sin θ + p vy aω x sin ψ -p vz aω x cos ψ cos θ = 0, (24)
These four constraints are dependent: they reduce to two functionally independent constraints. Hence, with ( 22)-( 23)-( 24), we have 8 independent constraints along I. Now, we infer from ( 22)-( 23)-( 24) that p ωx = p ωy = 0 and p θ = p ψ = p φ = 0 and

p vx = tan θp vz , p vy = -tan ψ/ cos θp vz , (25)
along I. Differentiating again, we get that θ = ψ = 0. It follows that ω x = ω y = 0. Using that H = 0 along any extremal, we get

p vz = -p 0 cos θ cos ψ a + g x sin θ cos ψ -g y sin ψ + g z cos θ cos ψ , (26)
Substituting ( 25) and ( 26) into p, ad 4 f. g 1 and p, ad 4 f. g 2 , we get p, ad 4 f. g 1 = p, ad 4 f. g 2 = 0, p, [ g 1 , ad 3 f. g 2 ](x) = p, [ g 2 , ad 3 f. g 1 ](x) = 0, p, [ g 1 , ad 3 f. g 1 ](x) = p, [ g 2 , ad 3 f. g 2 ](x) = -ap vz cos ψ cos θ .

According to Definition 1, to prove that u is of intrinsic order two, it suffices to prove that p, [ g i , ad 3 f. g i ](x) = 0 along I. We prove it by contradiction. If p, [ g i , ad 3 f. g i ](x) = 0, then necessarily p vz = 0 and this would lead to p vx = p vy = 0. It follows then from H = 0 that p 0 = 0. We have obtained that (p, p 0 ) = 0, which is a contradiction. The fact that u = (u 1 , u 2 ) = (0, 0) simply follows from the fact that

u 1 = -{h 1 , ad 3 h 0 .h 1 }/ad 4 h 0 .h 1 , u 2 = -{h 2 , ad 3 h 0 .h 2 }/ad 4 h 0 .h 2 .
Besides, if p 0 = 0, then p vz = 0 and p vx = p vy = 0, which leads to (p, p 0 ) = 0 and thus to a contradiction as well. Therefore, p 0 < 0 (i.e., the singular arc is normal), and then [START_REF] Kelley | Singular extremals[END_REF] follows by applying the GLCC of Corollary 1.

We define the singular surface S, which is filled by singular extremals of the problem (MTCP), by

S = (x, p) | ω x = ω y = 0, p θ = p ψ = p φ = p ωx = p ωy = 0, p vx = tan θp vz , p vz = -p 0 cos θ cos ψ a + g x sin θ cos ψ -g y sin ψ + g z cos θ cos ψ , p vy = -tan ψ/ cos θp vz . ( 27
)
We will see, in the next section, that the solutions of the problem of order zero (defined in Section 5.1.1) live in this singular surface S.

The following result, establishing chattering for the problem (MTCP), is a consequence of Theorem 1, Lemma 5 and Lemma 6.

Corollary 2. For the problem (MTCP), any optimal singular arc cannot be connected with a nontrivial bang arc. There is a chattering arc when trying to connect a regular arc with an optimal singular arc. More precisely, let u be an optimal control, solution of (MTCP), and assume that u is singular on the sub-interval (t 1 , t 2 ) ⊂ [0, t f ] and is regular elsewhere. If t 1 > 0 (resp., if t 2 < t f ) then, for every ε > 0, the control u switches an infinite number of times over the time interval

[t 1 -ε, t 1 ] (resp., on [t 2 , t 2 + ε]).
As mentioned in Section 4.1, points of order 4 are related to the chattering phenomenon. Indeed, from [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF], the switching function for the problem (MTCP) has the form

Φ(t) = R 0 (t 1 -t) 4 e iαln(t1-t) + o((t 1 -t) 4 ) if t < t 1 , R 0 (t -t 2 ) 4 e iαln(t-t2) + o((t -t 2 ) 4 ) if t > t 2 ,
in the neighborhood of the singular junction. The extremal converges to the switching surface when t → t 1 , t < t 1 , and escapes from the singular surface when t > t 2 . In-between, the control switches an infinite number of times if

t → t 1 , t > t 1 , or if t → t 2 , t < t 2 .
This result is important for solving the problem (MTCP) in practice. Indeed, when using numerical methods to solve the problem, the chattering control is an obstacle to convergence, especially when using an indirect approach (shooting). The existence of the chattering phenomenon in the problem (MTCP) explains well why the indirect methods may fail for certain terminal conditions.

Note that, in the planar version of the problem (MTCP) studied in [START_REF] Zhu | Planar tilting maneuver of a spacecraft: singular arcs in the minimum time problem and chattering[END_REF]), one can give sufficient conditions on the initial conditions under which the chattering phenomenon does not occur. Unfortunately, we are not able to derive such conditions in the general problem (MTCP).

Numerical approaches

In this section, we design two different numerical strategies for solving the problem (MTCP): one is based on combining indirect methods with numerical continuation, and the other is based on a direct transcription approach. The first one may be successfully implemented when dealing with solutions without chattering arcs, and the second one is more appropriate to compute solutions involving chattering arcs. However, both approaches are difficult to initialize successfully because the problem (MTCP) is of quite high dimension, is highly nonlinear, and moreover, as a main reason, the system consists of fast (Euler angles and angular velocity) and slow (orbit velocity) dynamics at the same time.

The occurrence of chattering arcs is an obstacle to convergence. Especially for indirect methods, the chattering phenomenon raises an important difficulty due to the numerical integration of the discontinuous Hamiltonian system. Direct transcription approaches provide a sub-optimal solution of the problem that has a finite number of switchings based on a (possibly rough) discretization. Actually, in case of chattering, we are also able to provide a sub-optimal solution with our indirect approach, by stopping the continuation before it would fail due to chattering. Though the suboptimal solutions provided in this way may be "less optimal" compared with those given by a direct approach, in practice they can be computed in a much faster way and also much more accurately.

Indirect method and numerical continuation

The idea of this continuation procedure is to use the (easily computable) solution of a simpler problem, that we call hereafter the problem of order zero, in order then to initialize an indirect method for the more complicated problem (MTCP). Then we are going to plug this simple, low-dimensional problem in higher dimension, and then come back to the initial problem by using appropriate continuations.

This method actually gives an optimal solution with high accuracy. The problem of order zero defined below is used as the starting problem because the orbit movement is much slower compared with the attitude movement and it is easy to solve explicitly. As well, it is worth noting that the solution of the problem of order zero is contained in the singular surface S filled by the singular solutions for the problem (MTCP), defined by [START_REF] Martinon | Using switching detection and variational equations for the shooting method[END_REF].

Two auxiliary problems

Problem of order zero. We define the problem of order zero, denoted by (OCP0), as a "subproblem" of the complete problem (MTCP), in the sense that we consider only the orbit dynamics and that we assume that the attitude angles (Euler angles) can be driven to the target values instantaneously. Thus, the attitude angles are considered as control inputs in that simpler problem. Denoting the rocket longitudinal axis as e and considering it as the control vector (which is consistent with the attitude angles θ, ψ), we formulate the problem as follows:

˙ V = a e + g, V (0) = V 0 , V (t f )// w, w = 1, min t f ,
where w is a given vector that refers to the desired target velocity direction. This problem is easy to solve, and the solution is the following.

Lemma 7. The optimal solution of (OCP0) is given by

e * = 1 a k w -V 0 t f -g , t f = -a 2 + a 2 2 -4a 1 a 3 2a 1 , p v = -p 0 a + e * , g e * .
with k = V 0 , w + g, w t f , a 1 = a 2 -g, w w -g 2 , a 2 = 2( V 0 , w g, w -V 0 , g ), and

a 3 = -V 0 , w w -V 0 2 .
Proof. The Hamiltonian is H = p 0 + p v (a e+ g), and we have ˙ p v = 0, with p v = (p vx , p vy , p vz ) ⊤ , and H = 0 along any extremal. It follows that p v = 0 (indeed otherwise we would get also p 0 = 0, and thus a contradiction). Hence there are no singular controls for this problem. The maximization condition of the PMP yields e * = p v / p v , and hence the optimal control is a constant vector. Moreover, according to the final condition V (t f )// w, the transversality condition is p v ⊥ w, hence e * , w = 0, and using V (t f ) = V 0 + (a e + g)t f = k w we get that e * = 1 a ( k w-V0 t f -g). It follows from the transversality condition that k = V 0 , w + g, w t f . The expression of t f follows, using that e * 2 = 1. Using that H = 0, we get p v = -p 0 a+ e * , g e * .

Since the vector e is expressed in the launch frame as ( e) R = (sin θ cos ψ,sin ψ, cos θ sin ψ) ⊤ , the Euler angles θ * ∈ (-π, π) and ψ * ∈ (-π/2, π/2) are given by

θ * =      arctan(e * 1 /e * 3 ), e * 3 > 0, π/2 sign(e * 1 ), e * 3 = 0 (π -arctan(|e * 1 |/|e * 3 |))sign(e * 1 ), e * 3 < 0, , ψ * = -arcsin(e * 2 ), (28)
where e * i is the i-th component of e * , for i = 1, 2, 3. Using the definition [START_REF] Martinon | Using switching detection and variational equations for the shooting method[END_REF] of the singular surface S, we check that the optimal solution of (OCP0) is contained in S with θ = θ * , ψ = ψ * and φ = φ * (φ * is any real number). Therefore, the relationship between (OCP0) and (MTCP) is the following. Lemma 8. The optimal solution of the problem (OCP0) actually corresponds to a singular solution of (MTCP) with the terminal conditions given by

v x (0) = v x0 , v y (0) = v y0 , v z (0) = v z0 , θ(0) = θ * , ψ(0) = ψ * , , φ(0) = φ * , ω x (0) = 0, ω y (0) = 0, (29)
v z (t f ) sin ψ f + v y (t f ) cos θ f cos ψ f = 0, v z (t f ) sin θ f -v x (t f ) cos θ f = 0, (30)
θ(t f ) = θ * , ψ(t f ) = ψ * , , φ(t f ) = φ * , ω x (t f ) = 0, ω y (t f ) = 0. (31)
Due to this result, a natural idea of numerical continuation strategy consists of deforming continuously (step by step) the terminal conditions given in Lemma 8, to the terminal conditions (11)-( 12) of the problem (MTCP).

However, because of the chattering phenomenon, we cannot make converge the shooting method in such a strategy. More precisely, when the terminal conditions are in the neighborhood of the singular surface S, the optimal extremals are likely to contain a singular arc (and thus chattering arcs). In that case, the shooting method will certainly fail due to the difficulty of numerical integration of discontinuous Hamiltonian system. Hence, we introduce hereafter an additional numerical trick and we define the following regularized problem, in which we modify the cost functional with a parameter γ, so as to overcome the problem caused by chattering.

Regularized problem. Let γ > 0 be arbitrary. The regularized problem (OCPR) γ consists of minimizing the cost functional

C γ = t f + γ t f 0 (u 2 1 + u 2 2 ) dt, (32)
for the bi-input control-affine system (9), under the control constraints -1 u i 1, i = 1, 2, with terminal conditions (11)- [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF]. Here, we replace the constraint u 2 1 + u 2 2 1 (i.e., u takes its values in the unit Euclidean disk) with the constraint that u takes its values in the unit Euclidean square. Note that we use the Euclidean square (and not the disk) because we observed that our numerical simulations worked better in this case. The advantage, for this intermediate optimal control problem with the cost [START_REF] Schättler | Geometric Optimal Control: Theory, Methods and Examples[END_REF], is that the extremal controls are then continuous.

The Hamiltonian is

H γ = p, f (x) + u 1 p, g 1 (x) + u 2 p, g 2 (x) + p 0 (1 + γu 2 1 + γu 2 2 ), (33)
and according to the PMP, the optimal controls are

u 1 (t) = sat(-1, -bp ωy (t)/(2γp 0 ), 1), u 2 (t) = sat(-1, bp ωx (t)/(2γp 0 ), 1), (34)
where the saturation operator sat is defined by sat(-1,

f (t), 1) = -1 if f (t) -1; 1 if f (t) 1; and f (t) if -1 f (t) 1.
As mentioned previously, one of the motivations for considering the intermediate problem (OCPR) γ is that the solution of (OCP0) is a singular trajectory of the full problem (MTCP), and hence, passing directly from (OCP0) to (MTCP) causes difficulties due to chattering (see Corollary 2). The following result shows that when we embed the solutions of (OCP0) into the problem (OCPR) γ , they are not singular. Lemma 9. An extremal of (OCP0) can be embedded into the problem (OCPR) γ by setting

u(t) = (0, 0), θ(t) = θ * , ψ(t) = ψ * , φ(t) = φ * , ω x (t) = 0, ω y (t) = 0, p θ (t) = 0, p ψ (t) = 0, p φ (t) = 0, p ωx (t) = 0, p ωy (t) = 0,
where θ * and ψ * are given by [START_REF] Marec | Optimal Space trajectories[END_REF], with terminal conditions given by ( 29) and ( 30)- [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]. Moreover, it is not a singular extremal for the problem (OCPR) γ . The extremal equations for (OCPR) γ are the same than for (MTCP), as well as the transversality conditions.

Proof. It is easy to verify that the embedded extremal is an extremal of the problem (OCPR) γ and that the transversality conditions are the same. The control is computed from (34) which maximizes the Hamiltonian H γ , and we have H γ = 0 with p 0 = -1. It follows from the PMP that the extremal equations are the same than for (MTCP). Then, for the problem (OCPR) γ , we have

∂ 2 Hγ ∂u 2 i = γp 0 . Note that, in this case, the control u i , i = 1, 2 is singular if ∂ 2 Hγ ∂u 2 i = 0.
Hence there is no normal singular extremal for the problem (OCPR) γ . From Lemma 7, it is easy to see that p 0 = 0 and thus the extremals of (OCP0) are not singular extremals of (OCPR) γ .

Strategy for solving (MTCP)

Continuation procedure. The ultimate objective is to compute the optimal solution of the problem (MTCP), starting from the explicit, simple to compute, solution of (OCP0). We proceed as follows:

• First, according to Lemma 9, we embed the solution of (OCP0) into (OCPR) γ . For convenience, we still denote by (OCP0) the problem (OCP0) seen in high dimension.

• Then, we pass from (OCP0) to (MTCP) by means of a numerical continuation procedure, involving three continuation parameters: the first two parameters λ 1 and λ 2 are used to pass continuously from the optimal solution of (OCP0) to the optimal solution of the regularized problem (OCPR) γ , for some fixed γ > 0, and the third parameter λ 3 is then used to pass to the optimal solution of (MTCP) (see Figure 3). The parameter λ 1 is used to act, by continuation, on the initial conditions, according to

θ(0) = θ * (1 -λ 1 ) + θ 0 λ 1 , ψ(0) = ψ * (1 -λ 1 ) + ψ 0 λ 1 , φ(0) = φ * (1 -λ 1 ) + φ 0 λ 1 , ω x (0) = ω * x (1 -λ 1 ) + ω x0 λ 1 , ω y (0) = ω * y (1 -λ 1 ) + ω y0 λ 1 ,
where ω * x = ω * y = 0, φ * = 0, and θ * , ψ * are calculated through equation ( 28). Using the transversality condition [START_REF] Guzzetti | Coupling Attitude and Orbital Motion of Extended Bodies In The Restricted Circular 3-Body Problem: A Novel Study On Effects And Possible Exploitations[END_REF] and the extremal equations ṗvx = 0, ṗvy = 0 and ṗvz = 0, the unknown p vy can be expressed in terms of p vx and p vz as p vy = (p vx sin θ f cos ψ f + p vz cos θ f cos ψ f )/ sin ψ f , and hence the unknowns of the shooting problem are reduced to p vx , p vz , p θ (0), p ψ (0), p φ (0), p ωx (0), p ωy (0) and t f . The shooting function S λ1 for the λ 1 -continuation is defined by

S λ1 = p ωx (t f ), p ωy (t f ), p θ (t f ), p ψ (t f ), p φ (t f ), H γ (t f ), v z (t f ) sin ψ f + v y (t f ) cos θ f cos ψ f , v z (t f ) sin θ f -v x (t f ) cos θ f ,
where H γ (t f ) with p 0 = -1 is calculated from [START_REF] Seywald | Singular control in minimum time spacecraft reorientation[END_REF] and u 1 and u 2 are given by [START_REF] Shen | Time-optimal control of axisymmetric rigid spacecraft using two controls[END_REF]. In fact, from Lemma 6, we know that a singular extremal of problem (MTCP) must be normal, and since we are starting to solve the problem from a singular extremal, here we assume that p 0 = -1.

Note that we can use S λ1 as shooting function thanks for (OCPR) γ . For problem (MTCP), if S λ1 = 0, then together with ω x (t f ) = 0 and ω y (t f ) = 0, the final point (x(t f ), p(t f )) of the extremal is then lying on the singular surface S defined by [START_REF] Martinon | Using switching detection and variational equations for the shooting method[END_REF] and this will cause the fail of the shooting. However, for problem (OCPR) γ , even when x(t f ) ∈ S, the shooting problem can still be solved.

Initializing with the solution of (OCP0), we can solve this shooting problem with λ 1 = 0, and we get a solution of (OCPR) γ with the terminal conditions ( 29)- [START_REF] Moré | The MINPACK project[END_REF] (the other states at t f being free). Then, by continuation, we make λ 1 vary from 0 to 1, and in this way we get the solution of (OCPR) γ for λ 1 = 1. With this solution, we can integrate extremal equations ( 8) and [START_REF] Goh | Necessary conditions for singular extremals involving multiple control variables[END_REF] to get the values of the state variable at t f . Then denote θ e := θ(t f ), ψ e := ψ(t f ), φ e := φ(t f ), ω xe := ω x (t f ) and ω ye := ω y (t f ).

In a second step, we use the continuation parameter λ 2 to act on the final conditions, in order to make them pass from the values θ e , ψ e , φ e , ω xe and ω ye , to the desired target values θ f , ψ f , φ f , ω xf and ω yf . The shooting function is

S λ2 = ω x (t f ) -(1 -λ 2 )ω xe -λ 2 ω x f , ω y (t f ) -(1 -λ 2 )ω ye -λ 2 ω y f , θ(t f ) -(1 -λ 2 )θ e -λ 2 θ f , ψ(t f ) -(1 -λ 2 )ψ e -λ 2 ψ f , φ(t f ) -(1 -λ 2 )φ e -λ 2 φ f , v z (t f ) sin ψ f + v y (t f ) cos θ f cos ψ f , v z (t f ) sin θ f -v x (t f ) cos θ f , H γ (t f ) .
Solving this problem by making vary λ 2 from 0 to 1, we obtain the solution of (OCPR) γ with the terminal conditions (11)- [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF].

Finally, in order to compute the solution of (MTCP), we use the continuation parameter λ 3 to pass from (OCPR) γ to (MTCP). We add the parameter λ 3 to the Hamiltonian H γ and to the cost functional [START_REF] Schättler | Geometric Optimal Control: Theory, Methods and Examples[END_REF] as follows:

C γ = t f + γ t f 0 (u 2 1 + u 2 2 )(1 -λ 3 ) dt, H(t f , λ 3 ) = p, f + p, g 1 u 1 + p, g 2 u 2 + p 0 + p 0 γ(u 2 1 + u 2 2
)(1 -λ 3 ). Then, according to the PMP, the extremal controls are given by u i = sat(-1, u ie , 1), i = 1, 2, where

u 1e = bp ωy -2p 0 γ(1 -λ 3 ) + bλ 3 p 2 ωx + p 2 ωy , u 2e = -bp ωx -2p 0 γ(1 -λ 3 ) + bλ 3 p 2 ωx + p 2 ωy .
The shooting function S λ3 is defined as S λ2 , replacing H γ (t f ) with H γ (t f , λ 3 ). The solution of (MTCP) is then obtained by making vary λ 3 continuously from 0 to 1.

Remark 2. Note that the above continuation procedure fails in case of chattering (see Corollary 2), and thus cannot be successful for any possible choice of terminal conditions. In particular, if chattering occurs then the λ 3 -continuation is expected to fail for some value λ 3 = λ * 3 < 1. But in that case, with this value of λ 3 , we have generated a sub-optimal solution of the problem (MTCP), which appears to be acceptable and very interesting for practice. Moreover, the overall procedure is very fast and accurate. Note that the resulting sub-optimal control is continuous.

Direct method

We now propose a direct approach for solving the problem (MTCP), where the control is approximated by a piecewise constant control over a given time subdivision. The solutions derived from such a method are therefore sub-optimal, in particular when the control is chattering (and in such a case the number of switchings is limited by the time step). Note that this approach is much more computationally demanding than the indirect one.

Since the initialization of a direct method may also raise some difficulties, we propose the following strategy. The idea is to start from the solution of the problem (MTCP) with less terminal requirements, which is easy to obtain with a direct method, and then we introduce step by step the final conditions [START_REF] Cerf | Continuation from a flat to a round Earth model in the coplanar orbit transfer problem[END_REF] of the problem (MTCP). We implement this direct approach with the software BOCOP and its batch optimization option (see [START_REF] Bonnans | Bocop-A collection of examples[END_REF]).

• Step 1: we solve the problem (MTCP) with initial conditions [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF] and final conditions

ω y (t f ) = 0, θ(t f ) = θ f , v z (t f ) sin θ f -v x (t f ) cos θ f = 0.
These final conditions are the ones of the planar version of (MTCP) in which the motion of the spacecraft is 2D (see [START_REF] Zhu | Planar tilting maneuver of a spacecraft: singular arcs in the minimum time problem and chattering[END_REF] for details). Numerical simulations show that, with such terminal conditions, the problem (MTCP) is easy and fast to solve by means of a direct method (a constant initial guess for the discretized variables suffices to ensure convergence).

• Then, in Steps 2, 3, 4 and 5, we add successively (and step by step) the final conditions v z (t f ) sin At the end of this process, we have obtained the solution of the full problem (MTCP). Note again that this direct approach is much slower than the indirect one, and that the resulting control has many numerical oscillations (see numerical results in Section 6.2).

ψ f + v y (t f ) cos θ f cos ψ f = 0, ψ(t f ) = ψ f , φ(t f ) = φ f

Numerical results

The structure of the rocket is presented in Figure 2 (b). We assume that the thrust I is flexible, i.e., it can turn ±6 • in all directions, and its thrust is around T att = 1400 kN. The other thrusts are fixed with a total thrust T tot = 1 × 10 5 kN. The rocket mass is 800 t, the length of the rocket is l r = 50 m and its radius is r r = 2.5 m. Considering the rocket as a cylinder, we have I x = I y = m(3r 2 r + l 2 r )/12 and I z = mr 2 r /2. The parameters a and b in ( 7) and ( 4) are therefore a = T tot /m ≈ 12 and b = Tattlr 2Ix µ max ≈ 0.02. During the atmospheric ascent phase, the velocity of the rocket remains between several hundreds m/s and around 1000 m/s. Let v = v 2

x + v 2 y + v 2 z be the modulus of the velocity, and let ψ v and φ v be the flight path angles that we use to calculate the components of the velocity in S R frame, i.e., v x = v sin θ v cos ψ v , v y = -v sin ψ v and v z = v cos θ v cos ψ v . In this section, the initial values of the angles θ v and ψ v are chosen equal to the initial values of the angles θ and ψ. This means that, before the maneuver, the rocket is on a trajectory with angle of attack equal to zero.

In the numerical simulations, we set v x0 = v 0 sin θ 0 cos ψ 0 , v y0 = -v 0 sin ψ 0 , v z0 = v 0 cos θ 0 cos ψ 0 . The terminal values are given in Table 1.

(TC1): ω x0 = ω y0 = 0, θ 0 = 75 • , ψ 0 = 0.5 • , φ 0 = 0 • ω x f = ω y f = 0, θ f = 85 • , ψ f = 5 • , φ f = 0 • (TC2): ω x0 = ω y0 = 0, θ 0 = 70 • , ψ 0 = 0.5 • , φ 0 = 0 • ω x f = ω y f = 0, θ f = 85 • , ψ f = 5 • , φ f = 0 • (TC3): ω x0 = ω y0 = 0, θ 0 = 85 • , ψ 0 = 0.5 • , φ 0 = 0 • ω x f = ω y f = 0, θ f = 75 • , ψ f = 5 • , φ f = 0 • Table 1: Terminal conditions
Note that v 0 is the module of velocity at time 0. In the next two subsections, we will choose different values of v 0 , and so here we do not assign to it a specific value. Moreover, we set γ = 50 as the weight of the L 2 -norm control term in the cost functional of problem (OCPR) γ .

Numerical results without chattering

In this section, we solve the problem (MTCP) using the indirect method combined with the numerical continuation described in Section 5.1. This method is implemented using a predictorcorrector (PC) continuation method (see [START_REF] Allgower | Introduction to numerical continuation methods SIAM[END_REF]). Note that the PC continuation method is discrete, in contrast with differential methods used in [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF], for which the Jacobian of the homotopy method must be computed. The Fortran routines hybrd.f (see [START_REF] Moré | The MINPACK project[END_REF]) and dop853.f (see [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF]) are used, respectively, for solving the shooting problem (Newton method) and for integrating the ordinary differential equations (with prediction).

The Euler angle θ is usually called the pitch angle, and a pitching up maneuver designates a maneuver with terminal condition θ f > θ 0 , while a pitching down maneuver designates a maneuver with terminal condition θ f < θ 0 .

Pitching up maneuvers. We set v 0 = 1000 m/s and we use the numerical values denoted by (TC1) in Table 1. The components of the state variable are reported on Figure 4. The optimal control, the adjoint variables p ωx (t) and p ωy (t) and the modulus of the switching function Φ(t) = b(p ωy , -p ωx ) are reported on Figure 5. We observe that the optimal control switches twice, at times 8.8 s and 25.8 s. These two switching points are of order 1 (i.e., Φ(t) = 0 and Φ(t) = 0). Accordingly with Lemma 2, the control turns with an angle π at those points. Let us give another numerical example, taking the same terminal conditions as previously except for v 0 , and we take v 0 = 1500 m/s. The time history of the state, of the optimal control and of the switching function are reported on Figures 6 and7. One can see on Figure 7 that the optimal control turns two more times with an angle π due to two switching points of order one.

Pitching down maneuvers. We set v 0 = 1500 m/s and we use the numerical values denoted by TC3 in Table 1. The optimal solution is drawn on Figures 8 and9.

The shorter maneuver time t f indicates that it is easier to turn clockwise the axis of the velocity vector than to turn it anti-clockwise. This corresponds to the intuition. The reason is that the total force induced by the gravity force tends to reduce v x , i.e., it helps the velocity to turn clockwise, and so together with the rocket thrust force, the maneuver time is less than that of the anti-clockwise case.

Note that the derived time history of the adjoint variable (for both pitching up and pitching down maneuvers) do not have the same order of magnitude, i.e., p ωx and p ωy are ten times larger than p θ and p ψ , and are thousand times larger than p vx , p vy and p vz . This indicates again that the shooting method is difficult to initialize successfully. We note that the indirect strategy proposed in Section 5.1.2 is efficient also because the smallest adjoint variables p vx , p vy and p vz are already quite accurately estimated thanks to the problem of order zero.

Numerical results with chattering arcs

Sub-optimal solution by the indirect approach. On Figures 10 and11 is given a sub-optimal solution of (MTCP) with the terminal conditions (TC2) of Table 1 and v 0 = 2000 m/s. Due to chattering, the continuation parameter λ 3 stops at value λ * 3 = 0.98 (see Remark 2). Observing from Figure 11, the switching function passes four times the switching surface Γ and remains small between times 26.5 and 40.3. The control, instead of bang-bang or singular, is continuous. The cost of this trajectory is 69.3 and the final time t f = 66.0 s. Sub-optimal solution by the direct approach. With the same terminal conditions as above, we now use the direct method described in Section 5.2. Numerical simulations show that the initialization step for the direct method procedure is quite robust (a constant initial guess is enough). The results are reported on Figures 12,13 and 14. We observe that, when t ∈ [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF]43], the control oscillates much with a modulus less than 1: this indicates that there is a singular arc in the "true" optimal trajectory, and therefore chattering according to Corollary 2.

Note that, along the singular arc, the variables ω x , ω y , p ωx , p ωy , p θ , p ψ and p φ are almost equal to 0, and we check that this arc indeed lives on the singular surface S defined by [START_REF] Martinon | Using switching detection and variational equations for the shooting method[END_REF]. Therefore, it turns out that there is a singular arc in the optimal trajectory, causing chattering at the junction with regular arcs.

The maneuver time is t f = 65.4 s. Compared with that of the sub-optimal solution derived from the indirect strategy, only 0.6 s are gained with the direct method. The direct approach is hundreds of times slower than the indirect approach and the obtained control presents many oscillations, which is not much appropriate for a practical use. towards the singular surface and then to get back to the target submanifold. We suspect that this is due to a turnpike phenomenon as described in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF], at least when the required transfer time is quite large.

Conclusion

We have studied the time optimal control of the rocket attitude motion combined with the orbit dynamics. The problem (MTCP) is of interest because of the coupling of guidance and navigation systems. However, this problem is difficult to solve because of the occurrence of the chattering phenomenon for certain terminal conditions.

Using geometric control, we have established a chattering result for bi-input control-affine systems. We have also classified the switching points for the extremals of the problem (MTCP), according to the order of vanishing of the switching function, showing the behavior of the control at the singularities.

In order to compute numerically the solutions of problem (MTCP), we have implemented two approaches. The indirect approach, combining shooting and numerical continuation, is timeefficient when the solution does not contain any singular arcs. For certain terminal conditions, the optimal solution of (MTCP) involves a singular arc that is of order two, and the connection with regular arcs can only be done by means of chattering. The occurrence of chattering causes the failure of the indirect approach. For such cases, we have proposed two possible numerical alternatives. Since our indirect approach involves three continuations, one of them being concerned with a continuation on the cost function (and thus on the Hamiltonian and the control), we have proposed, as a first alternative, to stop this last continuation before its failure: in such a way, we obtain a sub-optimal solution, which seems to be very acceptable for a practical use. The second alternative is based on a direct approach, and then we obtain as well a sub-optimal solution having a finite number of switchings, this finite number being limited by the chosen step of the subdivision in the discretization scheme. In any case, the direct strategy is much more time consuming than the indirect approach. Note that, in both cases, it is not required to know a priori the structure of the optimal solution (in particular, the number of switchings).

As an open issue, one may consider to add atmospheric forces in the model. Since the magnitude of the aero-forces is low (at least, it should be much smaller than the rocket thrust), we expect this extension to be doable, for instance by means of an additional continuation.

Figure 1 :

 1 Figure 1: Coordinate systems and relations.

  and let b = T max l/I x . Then (3) gives ωx = -b sin µ sin ζ, ωy = b sin µ cos ζ, with ω z ≡ 0. 1 flexible thrust I 4 thrust fix Configuration of the engines (b) structure of the rocket (a) thrust in the body frame

Figure 2 :

 2 Figure 2: Thrust in the body Frame

2

 2 ψp θ -(ω x sin φ + ω y cos φ)/ cos 2 ψp φ , ṗφ = -(ω x cos φ -ω y sin φ)/ cos ψp θ + (ω x sin φ + ω y cos φ)p ψ tan ψ(ω x cos φ -ω y sin φ)p φ , ṗωx =sin φ/ cos ψp θcos φp ψsin ψ sin φ/ cos ψp φ , ṗωy =cos φ/ cos ψp θ + sin φp ψsin ψ cos φ/ cos ψp φ ,

Figure 3 :

 3 Figure 3: Continuation procedure.

  , and ω x (t f ) = ω xf , and for each new step we use the solution of the previous one as an initial guess.

Figure 4 :Figure 5 :

 45 Figure 4: State variable and the optimal control with (TC1) and v 0 = 1000.

Figure 6 :

 6 Figure 6: Time histories of state with (TC1) and v 0 = 1500 m/s.

Figure 7 :

 7 Figure 7: Time histories of control and switching function with (TC1) and v 0 = 1500 m/s.

Figure 8 :

 8 Figure 8: Time histories of state with TC3 and v 0 = 1500.

Figure 9 :

 9 Figure 9: Time histories of control and switching function with TC3 and v 0 = 1500.

Figure 4 , 6 , 8 ,

 468 10 and 12, we note that the attitude angles first tend to reach the values θ * deg, ψ * (i.e., θ * = 176.9 deg and ψ * = 18.5 deg for Figures 4 and 6; θ * = -17.4 deg and ψ * = 25.2 deg for Figure 8; θ * = 176.1 deg and ψ * = 11.2 deg for Figures 10 and 12

Figure 10 :

 10 Figure 10: Time histories of state with (TC2) and v 0 = 2000.

Figure 11 :

 11 Figure 11: Time histories of control and switching function with (TC2) and v 0 = 2000.

Figure 12 :

 12 Figure 12: State variable x(t) (v 0 = 2000 m/s).

Figure 13 :

 13 Figure 13: Optimal control and Φ(t) (v 0 = 2000 m/s).

Figure 14 :

 14 Figure 14: Adjoint state p(t) (v 0 = 2000 m/s).

† In fact, for a launcher, the attitude maneuvers are generally planar, i.e., ψ(0) ≈ 0, ψ(t f ) ≈ 0, and no singularities

Acknowledgment. The second author acknowledges the support by FA9550-14-1-0214 of the EOARD-AFOSR.

Let us now assume that the vector fields g 1 and g 2 commute, i.e., [g 1 , g 2 ] = 0. By differentiating again, we get ḧ1 = {h 0 , {h 0 , h 1 }} + u 1 {h 1 , {h 0 , h 1 }} + u 2 {h 2 , {h 0 , h 1 }} = 0, ḧ2 = {h 0 , {h 0 , h 2 }} + u 1 {h 1 , {h 0 , h 2 }} + u 2 {h 2 , {h 0 , h 2 }} = 0. If det ∆ 1 = det {h 1 , {h 0 , h 1 }} {h 2 , {h 0 , h 1 }} {h 1 , {h 0 , h 2 }} {h 2 , {h 0 , h 2 }} = 0 along I, then

and we say that the control u = (u 1 , u 2 ) is of order 1 (also called minimal order in [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Chitour | Singular trajectories of control-affine systems[END_REF]). Note that u 1 and u 2 must moreover satisfy the constraint

Now, if {h 1 , {h 0 , h 1 }} = 0 and {h 2 , {h 0 , h 2 }} = 0 along I, then the singular control is of higher order, and we must have {h i , {h 0 , h j }} = 0, i, j = 1, 2, i = j according to the Goh condition (see [START_REF] Goh | Necessary conditions for singular extremals involving multiple control variables[END_REF][START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF], see also Definition 1 and Lemma 1), and hence we go on differentiating. Assuming that

and we get

Points of order 3. Lemma 4. We assume that z 0 is of order 3. Then the reference extremal is well defined in a neighborhood of t = 0, in the sense that there exists a unique extremal associated with the control u = (u 1 , u 2 ) passing through the point z 0 . If b i = ad 3 h 0 .h i (z 0 ) = 0 and c = {h 2 , ad 2 h 0 .h 1 }(z 0 ) = 0, then the switching function is of class C 3 in the neighborhood of t = 0 and the control turns with an angle π when passing through the switching surface Γ locally, and we have

where,

Proof. Using ( 10) and ( 16), we have h

(3)

= ad 3 h 0 .h 2u 1 {h 2 , ad 2 h 0 .h 1 } (see also Lemma 5). Assuming that we have locally h i (t) = 1 6

. By substituting the control into the expression of h

i , we get

|t| . The result follows by solving these two equations for t > 0 and t < 0.

Points of order 4. We assume that z 0 is of order 4. If {h 2 , ad 2 h 0 .h 1 }(z 0 ) = 0, then 0 = h

2 can be a regular control according to the value of If

then the extremal is well defined in a neighborhood of t = 0 and the control is continuous when passing the switching surface Γ.

see the proof of Lemma 6 further), and we have h