Coupled attitude and trajectory optimization for a launcher tilting maneuver

# Coupled attitude and trajectory optimization for a launcher tilting maneuver

Jiamin Zhu***Sorbonne Universités, UPMC Univ Paris 06, CNRS UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France (zhu@ann.jussieu.fr).    Emmanuel TrélatSorbonne Universités, UPMC Univ Paris 06, CNRS UMR 7598, Laboratoire Jacques-Louis Lions, Institut Universitaire de France, F-75005, Paris, France (emmanuel.trelat@upmc.fr).    Max CerfAirbus Defence and Space, Flight Control Unit, 66 route de Verneuil, BP 3002, 78133 Les Mureaux Cedex, France (max.cerf@astrium.eads.net).
###### Abstract

We study the minimum time control problem of the launchers. The optimal trajectories of the problem may contain singular arcs, and thus chattering arcs. The motion of the launcher is described by its attitude kinematics and dynamics and also by its trajectory dynamics. Based on the nature of the system, we implement an efficient indirect numerical method, combined with numerical continuation, to compute numerically the optimal solutions of the problem. Numerical experiments show the efficiency and the robustness of the proposed method.

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

## 1 Introduction

The usual way to approach the control of the guidance of satellites consists of considering separately the attitude and the trajectory dynamics. However, for the launchers like Ariane or Pegasus, the trajectories are controlled by their attitude angles, and it is desirable to determine the optimal control subject to the coupled dynamical system. Thus, we consider the time minimum control of the attitude reorientation of a rocket taking in consideration both the attitude movement and the trajectory dynamics. We denote this problem as . Our objective is to design an efficient numerical strategy to solve the problem .

There are two main types of numerical approaches to solve an optimal control problem. On the one part, the direct methods (see, e.g., [2]) consist of discretizing the state and the control and thus of reducing the problem to a nonlinear optimization problem (nonlinear programming) with constraints. On the other part, the indirect methods consist of numerically solving a boundary value problem obtained by applying Pontryagin maximum principle (PMP, see [10]), by means of a shooting method. Both direct and indirect methods are not easy to initialize successfully, it is required to combine them with other theoretical or numerical approaches (see the survey [11]). The numerical continuation is a powerful tool to be combined with the indirect shooting method based on the Pontryagin Maximum Principle. The idea of continuation is to solve a problem step by step from a simpler one by parameter deformation (see e.g. [1]). For example, in [3, 6, 8], the continuation method is used to solve difficult orbit transfer problems. Here, we will use numerical continuation to numerically solve the problem.

Indeed, the coupling of the attitude movement (fast) with the trajectory dynamics (slow) in the problem generates difficulties in numerical approaches, especially in the indirect methods, where the Newton-like methods are used to solve the boundary value problem. Moreover, this property makes the numerical approaches extremely difficult to be initialized. Another difficulty in numerically solving the problem is the chattering phenomenon. 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., [5, 7, 13]). We proved in [14, 15] that there exists such a chattering phenomenon in the problem , and the chattering phenomenon is a bad news when applying numerical approaches.

In view of these difficulties, we propose an efficient numerical continuation procedure to compute optimal and sub-optimal trajectories for the problem . 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, our indirect approach generates sub-optimal solutions, by stopping the continuation procedure before its failure due to chattering.

Our objective is to design and implement a numerical method having the following qualities: general (vehicle features, terminal conditions), robust, fast and automatic, so that it could be used efficiently by an inexperienced user on a large scope of applications of the real world. Numerical experiments indicates that this approach happens to meet to all the expected features.

The paper is organized as follows. In Section 2, we establish the model of the rocket movement and we formulate the problem . In Section 3, the PMP is applied to the problem and some theoretical results on the extremals are introduced. In Section 4, two auxiliary problems are introduced and the continuation procedure is stated. Some numerical tips improving the proposed numerical strategy are also presented. Finally in Section 5, numerical results are given.

## 2 Problem Statement

We will consider two types of launchers: one is a classical launcher like Ariane rocket, the other one is an airborne launcher like Pegasus rocket (being launched horizontally from an airplane). In the case of an Ariane-type launcher, the aerodynamical forces are very small compared with its thrust and the gravity, and so we ignore them. However, in case of a Pegasus-type launcher, we will need to add the aerodynamical forces to the model. Throughout the paper, we make the following assumptions:

• The Earth is a sphere and is fixed in the inertial space.

• The position and the mass of the rocket are constant during the maneuver.

• The rocket is an axial symmetric cylinder.

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

### 2.1 Coordinates and Model

All the coordinate systems introduced here are Cartesian coordinate systems.

The launch frame (reference frame) is fixed around the launch point (where the rocket is launched). The axis is pointing radially outwardly (normal to the local tangent plane), and the axis points to the North.

The body frame is defined as follows. The origin of the frame is fixed around the mass center of the rocket, the axis is along the axis-symmetric axis of the rocket, and the axis is in the cross-section. The body frame can be derived by three ordered unit single-axis rotations from the launch frame,

 SRRy(θ)−−−→∘Rx(ψ)−−−→∘Rz(ϕ)−−−→Sb,

where is the pitch angle, is the yaw angle, is the roll angle and means to rotate the frame around the axis by an angle . Therefore, the transfer matrix from to is .

The velocity frame is fixed around the mass center of rocket. The axis is parallel to the velocity of the rocket , and the axis is normal to the velocity, pointing to the direction of the lift force of the rocket. This frame can be derived by two unit single-axis rotations from the launch frame,

 SRRx(ξ)−−−→∘Ry(κ)−−−→Sv,

where is the flight path angle and is the bank angle of the flight. Denote by the module of the velocity vector, i.e., . since , we have

 cosξ=vx/v,tanκ=−vy/vz.

This frame will be used to introduce aerodynamical forces (lift and drag ) into the model.

#### Attitude motion

The dynamical and kinematic equations of the attitude are

 ˙θ=(ωxsinϕ+ωycosϕ)/cosψ,˙ψ=ωxcosϕ−ωysinϕ,˙ϕ=(ωxsinϕ+ωycosϕ)tanψ,˙ωx=−¯bu2,˙ωy=¯bu1. (1)

where is the control input expressed in frame , is the angular velocity of the launcher expressed in frame (assuming that there is no rotation along the asymptotical axis of the launcher), and is the angular acceleration constant depending on (the maximum angle between the thrust and the axis ), (the inertia around the axis ), (the length of the rocket). Note that . See [15] for details of this model.

#### Trajectory dynamics

The drag and lift are calculated by , , where is the air density, is the reference surface, and are constant coefficients of the drag and the lift. In the velocity frame , the components of the vectors and are expressed as , where and . By using the transfer matrix from to , we get

 (→D)R=cxmv2(−cosξ,−sinξsinκ,sinξcosκ)⊤,(→L)R=czmv2(−sinξ,−cosξsinκ,cosξcosκ)⊤,

and so we get the trajectory dynamical equations in frame as

 ˙vx=asinθcosψ+gx−cxv2cosξ−czv2sinξ,˙vy=−asinψ+gy−cxv2sinξsinκ−czv2cosξsinκ,˙vz=acosθcosψ+gz+cxv2sinξcosκ+czv2cosξcosκ, (2)

where is the velocity of the launcher expressed in , is the gravity vector expressed in , is the thrust acceleration. For the Ariane rocket, we do not consider the aerodynamical forces, i.e., .

### 2.2 Minimum Time Control Problem (MTCP)

Defining the state variable and the control variable as

 x=(vx,vy,vz,θ,ψ,ϕ,ωx,ωy),u=(u1,u2).

The initial conditions are defined by , , , , , , and ,

 vx(0)=vx0,vy(0)=vy0,vz(0)=vz0,θ(0)=θ0,ψ(0)=ψ0,ϕ(0)=ϕ0,ωx(0)=ωx0,ωy(0)=ωy0. (3)

The desired final velocity is required to be parallel to the axis , according to . The constraints on the final conditions are defined by , , , and .

 vzfsinψf+vyfcosθfcosψf=0,vzfsinθf−vxfcosθf=0,θ(tf)=θf,ψ(tf)=ψf,ϕ(tf)=ϕf,ωx(tf)=ωxf,ωy(tf)=ωyf. (4)

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. The angle of attack, when the wind is set to zero, is defined as the angle between the velocity and the rocket body axis.

We set , and we define the target set

 M1={(vx,vy,vz,θ,ψ,ϕ,ωx,ωy)∈R8 ∣ vzsinψf+vycosθfcosψf=0,vzsinψf+vycosθfcosψf=0,θ=θf,ψ=ψf,ϕ=ϕf,ωx=ωxf,ωy=ωyf}.

The minimum time control problem consists of steering the system (1)-(2) from to the final target in minimum time , with controls satisfying the constraint .

For convenience, the velocity will be defined by polar coordinates, i.e., , , , where and are “pitch” and “yaw” angles of the velocity vector.

## 3 Theoretical Results

According to the Pontryagin Maximum Principle (PMP), the Hamiltonian of the optimal control problem is

 H(x(t),p(t),p0,u(t))=⟨p(t),f(x(t),u(t))⟩+p0,

where is the adjoint variable, is the derivative of the state variable , i.e., and is a real number. Here we assume (see [15]). The differential equation of the adjoint variable is given by the partial derivative of the Hamiltonian as

 ˙p(t)=−∂H∂x(x(t),p(t),p0,u(t)). (5)

The maximization condition of the PMP yields, almost everywhere on ,

 u(t)=(h1(t),h2(t))√h1(t)2+h2(t)2=Φ(t)∥Φ(t)∥,

whenever . Here we have and . We call (as well as its components) the switching function.

Moreover, we have the transversality condition , where is the tangent space to at the point , i.e.,

 pvysinψf=pvxsinθfcosψf+pvzcosθfcosψf, (6)

and, the final time being free and the system being autonomous, we have also .

The quadruple is called an extremal lift of . An extremal is said to be normal (resp., abnormal) if (resp., ). We say that an arc (restriction of an extremal to a subinterval ) is regular if along . Otherwise, the arc is said to be singular.

A switching time is a time at which , that is, both and vanish at time . An arc that is a concatenation of an infinite number of regular arcs over a compact time interval 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.

When and , the explicit form of equation (5) is given by

 ˙pvx=0,˙pvy=0,˙pvz=0,˙pθ=−acosψ(pvxcosθ−pvzsinθ),˙pψ=asinψsinθpvx+acosψpvy+acosθsinψpvz−sinψ(ωxsinϕ+ωycosϕ)/cos2ψpθ−(ωxsinϕ+ωycosϕ)/cos2ψpϕ,˙pϕ=−(ωxcosϕ−ωysinϕ)/cosψpθ+(ωxsinϕ+ωycosϕ)pψ−tanψ(ωxcosϕ−ωysinϕ)pϕ,˙pωx=−sinϕ/cosψpθ−cosϕpψ−sinψsinϕ/cosψpϕ,˙pωy=−cosϕ/cosψpθ+sinϕpψ−sinψcosϕ/cosψpϕ. (7)

We have shown in [15] that the singular extremals of the problem are normal ones, i.e., , and that they are of intrinsic order two. In that case the singular junction can only be realized by chattering, i.e., if the singular junction happens at time and the control is regular on and is singular on , then the control has to switch an infinite number of times on , such that .

The concatenation of regular arcs are classified in [15] by their contact with the switching surface (filled by switching points), and we know that if the switching points of the problem are of order one, then the control will turn an angle when passing the switching surface. We will see this phenomenon in the numerical results.

Actually, when and , the derivative of is no longer zero. However, the terms induced by the air draft and lift do not change the Lie bracket configuration of the system, and thus the results obtained in [15] are still valid.

Note again that the chattering causes hard problems when applying numerical methods. In that case, we can no longer obtain an optimal solution numerically. However, a sub-optimal solution can be derived by stopping the continuation procedure proposed in Section 4.2 before chattering occurs.

## 4 Auxiliary Problems and Continuation Procedure

To construct a suitable continuation procedure, we introduce two auxiliary problems. The first one is the problem of order zero, denoted by . The solution of this problem can be easily computed. Then we just need to plug this simple, low-dimensional solution in higher dimension, in order to initialize an indirect method for the more complicated problem . The second one is the regularized problem, which is used to overcome the chattering issues. Since the solution of problem is contained in the singular surface (filled by the singular solutions) for the problem , passing directly from the solution of problem to the problem makes the optimal extremals to contain a singular arc (and thus chattering arcs), and the shooting method will certainly fail due to the numerical integration of discontinuous Hamiltonian system.

### 4.1 Two auxiliary problems

#### Problem of order zero.

We define the problem of order zero, as a “subproblem” of the complete problem , in the sense that we consider only the trajectory dynamics (without aerodynamical forces) and that we assume a perfect control meaning that the attitude angles (Euler angles) take the desired values instantaneously. Thus, the attitude angles are considered as control inputs in that simpler problem. Denoting the rocket axial symmetric axis as 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)=→v0,→v(tf)//→w,∥→w∥=1,mintf,

where is a given vector that refers to the desired target velocity direction. This problem is easy to solve, and the solution derived by applying the PMP is

 →e∗=1a(k→w−→v0tf−→g),tf=−a2+√a22−4a1a32a1,→pv=−p0a+⟨→e∗,→g⟩→e∗.

with , , , and .

Since the vector is expressed in the launch frame as , the Euler angles can be calculated by assuming as

 θ∗=atan2(e∗x,e∗z)=sign(e∗x)⎧⎨⎩arctan(|e∗x/e∗z|),e∗z>0,π/2,e∗z=0,(π−arctan(|e∗x/e∗z|)),e∗z<0, (8)

where is the sign function, and

 ψ∗=arcsin(−e∗y). (9)

The assumption of generally works for the Ariane-type rockets since the maneuvers are generally quasi planar and the change in the yaw angle is small. However, for the Pegasus-type rockets, which could make large yaw angle maneuvers, we have to consider also the case leading to

 θ∗=atan2(−e∗x,−e∗z),ψ∗=−sign(e∗y)(π−arcsin(|−e∗y|)). (10)

Both 10 and (8) and (9) are possible solutions. We will choose the pair minimizing the value .

Indeed, when , , the Euler angles are not well defined, and if and are given near , we need to “change” them to a be near . This is done in section 4.3 below by transforming of the reference frame to a new reference frame . The singularities of Euler angles are also treated in section 4.4 by calculating the limit values of the vector field at these singular points.

Given a real number , the optimal solution of the problem actually corresponds to a singular solution of with the terminal conditions given by

 vx(0)=vx0,vy(0)=vy0,vz(0)=vz0,θ(0)=θ∗,ψ(0)=ψ∗,,ϕ(0)=ϕ∗,ωx(0)=0,ωy(0)=0, (11)
 vz(tf)sinψf+vy(tf)cosθfcosψf=0,vz(tf)sinθf−vx(tf)cosθf=0, (12)
 θ(tf)=θ∗,ψ(tf)=ψ∗,,ϕ(tf)=ϕ∗,ωx(tf)=0,ωy(tf)=0. (13)

It is worth nothing that this solution is lying in the singular surface of the problem meaning that it is on the “highway” between two fixed points in the state space. On this “highway”, the system goes the most rapidly towards to aimed final point or manifold. Indeed, we observe in the numerical simulations that the singular extremals of the problem acts a role similar to that of the stable points in the “turnpike” phenomenon as described in [12]: the optimal trajectories first tend to reach the singular surface (to have a greater speed in transferring its state), then stay in the singular surface for a while until it is sufficiently close to the target submanifold , and finally get off the singular surface to join the target submanifold. Note that the singular arc is not necessary if the state is sufficiently close to the target.

In view of (11)-(13), a natural continuation strategy is to simply vary the terminal conditions step by step until they correspond to . However, as we have mentioned, the chattering phenomenon causes the failure of this simple strategy, and thus we introduce another auxiliary problem, the regularized problem, in which the singular arcs of the problem become regular arcs, and thus the difficulty caused by chattering is bypassed.

#### Regularized problem.

Let be arbitrary. The regularized problem  consists of minimizing the cost functional

 Cγ=tf+γ∫tf0(u21+u22)dt, (14)

for the system (1)-(2), under the control constraints , , with terminal conditions (3)-(4). Note that, here, we replace the constraint with the constraints . The advantage, for this intermediate optimal control problem with the cost (14), is that the extremal controls are continuous.

The Hamiltonian is

 Hγ=⟨p,f(x,u)⟩+p0(1+γu21+γu22), (15)

and according to the PMP, the optimal controls are

 u1(t)=sat(−1,−¯bpωy(t)/(2γp0),1),u2(t)=sat(−1,¯bpωx(t)/(2γp0),1), (16)

where the saturation operator is defined by if ; if ; and if .

Indeed, an extremal of can be embedded into the problem  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 (8) and (9), with terminal conditions given by (11)-(13). Moreover, it is not a singular extremal for the problem  . The extremal equations for  are the same than for , as well as the transversality conditions.

Note that the difficulty caused by the chattering still exists when trying to go from the regularized problem back to the problem if there exists a singular arc in the optimal trajectory. However, in that case, by stopping at some point during the last continuation step (see below), an acceptable sub-optimal trajectory will be found for the problem .

### 4.2 Continuation procedure

The ultimate objective is to compute the optimal solution of the problem , starting from the explicit solution of .

#### Numerical strategy

We proceed as follows:

• First, we embed the solution of into  . For convenience, we still denote by the problem seen in high dimension.

• Then, we pass from to by means of a numerical continuation procedure, involving four continuation parameters: the first three parameters , and/or are used to pass continuously from the optimal solution of to the optimal solution of the regularized problem  , for some fixed , and the last parameter is then used to pass to the optimal solution of (see Figure 1).

The unknowns of the shooting problem are , , , , , , , and . When and , we have that , and are constants, hence by using the transversality condition (6), we can easily reduce the number of the unknowns by one. In view of more general cases, we will not do this reduction and we will use the transversality condition as a part of the shooting function. In the following, we denote the continuation step corresponding to parameter by -continuation.

#### λ1-continuation and λ2-continuation

The parameter 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 , , and , are calculated through equation (8)-(9).

The shooting function for the -continuation is defined by

 Sλ1=(pωx(tf),pωy(tf),pθ(tf),pψ(tf),pϕ(tf),Hγ(tf),vz(tf)sinψf+vy(tf)cosθfcosψf,vz(tf)sinθf−vx(tf)cosθfpvysinψf−(pvxsinθfcosψf+pvzcosθfcosψf)),

where with is calculated from (15) and and are given by (16).

Note that we can use as shooting function owing to  . For the problem , if , then together with and , the final point of the extremal is then lying on the singular surface and this will cause the failure of the shooting. However, for problem  , even when belongs to the singular surface, the shooting problem can still be solved.

Initializing with the solution of , we can solve this shooting problem with , and we get a solution of  with the terminal conditions (11)-(12) (the other states at being free). Then, by continuation, we make vary from to , and in this way we get the solution of  for . With this solution, we can integrate extremal equations (1), (2) and (7) to get the values of the state variable at . We denote , , , and the “natural” conditions obtained at the final time.

In a second step, we use the continuation parameter to act on the final conditions, in order to make them pass from the “natural” values , , , and , to the desired target values , , , and . The shooting function is

 Sλ2=(ωx(tf)−(1−λ2)ωxe−λ2ωxf,ωy(tf)−(1−λ2)ωye−λ2ωyf,θ(tf)−(1−λ2)θe−λ2θf,ψ(tf)−(1−λ2)ψe−λ2ψf,ϕ(tf)−(1−λ2)ϕe−λ2ϕf,vz(tf)sinψf+vy(tf)cosθfcosψf,vz(tf)sinθf−vx(tf)cosθf,pvysinψf−(pvxsinθfcosψf+pvzcosθfcosψf,Hγ(tf)).

Solving this problem by making vary from to , we obtain the solution of  with the terminal conditions (3)-(4).

#### λ3-continuation

The parameter is added to the terms induced by the aerodynamic forces in system (1)-(2), i.e., replace and by and . It is the only parameter that acts on the differential system, and it is used to differ the Ariane-type and Pegasus-type rockets. We let during the -continuation and the -continuation, which means during these first two steps of the continuation procedure, we do not consider the influence of the drag and the lift. Then, if we are dealing with a Pegasus-type rocket, we vary from to to derive a solution of problem  with drag and lift, and the shooting function remains the same as for the -continuation. Otherwise, for an Ariane-type rocket, we skip -continuation, i.e., we keep and go ahead to the -continuation. By doing this, the program is suitable for both types of rockets.

#### λ4-continuation

Finally, in order to compute the solution of , we use the continuation parameter to pass from  to . We add the parameter to the Hamiltonian and to the cost functional (14) as follows:

 Cγ=tf+γ∫tf0(u21+u22)(1−λ4)dt,
 H(tf,λ4)=⟨p,f(x,u)⟩+p0+p0γ(u21+u22)(1−λ4).

Then, according to the PMP, the extremal controls are given by , , where

 u1e=¯bpωy−2p0γ(1−λ4)+¯bλ4√p2ωx+p2ωy,u2e=−¯bpωx−2p0γ(1−λ4)+¯bλ4√p2ωx+p2ωy.

The shooting function is the same as , replacing with . The solution of is then obtained by making vary continuously from to .

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

As we have mentioned, it is possible that an optimal trajectory meets the singularities of the Euler angles, i.e., , for some . A coordinate system transformation can be made to change the set terminal conditions and avoid meeting the singularity. We will use two numerical tricks as follows. The first one is to use a new reference frame , instead of , such that the terminal conditions in the new reference frame are easier to solve. The new frame is set by rotating the initial frame . This amounts to a nonlinear state transformation, and leads to a preconditioner that makes the proposed continuation procedure more robust. The second trick is to overcome the singularities of the Euler angles by redefining the vector fields at the singular points.

### 4.3 Change of frame

We define a new coordinate which is derived by two single-axis rotations from the frame , i.e.,

 SRRy(α)−−−→∘Rx(β)−−−→SR′,

and so the transfer matrix from to is . Denoting the Euler angles of with respect to the new frame as , and , we get the transfer matrix from to as