Optimal Trajectory Synthesis and Tracking Control for Spacecraft Large Attitude Manoeuvers

Abstract

The classical approach to the problem of synthesizing an optimal attitude manoeuver trajectory, involves the use of the calculus of variations and the use Lagrange multipliers or co-states. The nonlinear large attitude manoeuver trajectory is controlled by a set of nonlinear evolving co-states. In this paper, following a review of the methodologies available for trajectory synthesis followed by tracking control, the optimal trajectory for a typical optimal attitude manoeuver is synthesized by solving for the states and co-states defined by a two point boundary value problem. Gravity gradient torques are included as a matter of course. Following the synthesis of the optimal attitude-rate trajectory, tracking control laws are synthesized by re-formulating the optimal control as a feedback law. The approximate linear tracking feedback controls are evaluated by relating the optimal state and co-state vector by a linear relation. The control laws are synthesized numerically. The problem of optimal attitude orientation trajectory synthesis is also addressed. The methodologies are applied to typical sample problems and results are presented. Quantitative comparisons of the results of the methods are made to the results obtained by the application of other linear and nonlinear methods, to illustrate the key features of the methodologies.

1. Introduction

The need for designing fast attitude and angular rate acquisition manoeuvers for a spacecraft with restricted or low actuation torqueing capacity arises in many space recent applications. Spacecraft are usually equipped with an attitude control system (ACS), which operates in one of two modes; the first mode involves maneuvering for attitude or angular rate acquisition while the second is to ensure stability. In the first mode, the ACS is responsible for acquiring and tracking an attitude or an angular rate state trajectory which could include a steady state. The requirements for this mode are set by the need to remotely capture an orbiting body, de-tumble a spacecraft, synchronize with another orbiting body or point at a specific direction in space. Although a large number of publications have appeared before the end of the last millennium on the subject of attitude stabilization and feedback control, a few recent papers have focused on the construction of optimal maneuvering trajectories synthesis for attitude or angular rate acquisition. There have been some publications related to the synthesis of optimal maneuvering trajectories for attitude or angular rate acquisition during the last two decades of the preceding millennium. Yet some significant advances have been made in the early part of this century. This includes papers by Lee et al. [1], Yoshida et al. [2], Aghili [3, 4], Yang and Wu [5], Liu et al. [6], and Zhang et al. [7], who have considered the maneuvering for attitude or angular rate acquisition problems using classical methodologies. Sharma and Tewari [8] have addressed the issue of nonlinear tracking of spacecraft attitude manoeuvers while Hegrenas et al. [9] have considered the maneuvering for attitude or angular rate acquisition problem by means of explicit model predictive control, via a nonlinear programming approach.

In dealing with the optimal attitude trajectory synthesis, several real world effects such as gravity gradient torques are generally neglected. Neglecting the gravity gradient torques can have a serious effect on the trajectory synthesis problem as (i) the gravity gradient torques can influence the stability of the spacecraft and (ii) they tend to couple the attitude rate dynamics with the quaternion kinematics. For this reason, it is not always advisable to ignore these torques on grounds of “smallness” as even the smallest of these perturbations can not only trigger instability but also induce the bifurcation of the orbit. It is the gravity that is primarily responsible for the orbital motion and the attitude stability of a spacecraft. The importance of gravity gradient torques has also been underscored by Lobo et al. [10].

The classical methodologies for trajectory synthesis compare well with other nonlinear and deterministic artificial intelligence approaches such as those developed by Sands et al. [11], Nakatani and Sands [12] and Baker et al. [13].

In this paper, the optimal trajectory for a typical attitude rate manoeuver is synthesized by solving for the states and co-states defined by a two point boundary value problem. Gravity gradient torques are included as a matter of course. Following the synthesis of the optimal attitude-rate trajectory, tracking control laws are synthesized by re-formulating the optimal control as a feedback law. The approximate linear tracking feedback controller gains are evaluated by relating the optimal state and co-state vector by a linear relation. The feedback control laws are synthesized numerically. The problem of optimal attitude orientation trajectory synthesis is also addressed. The optimization methodologies are applied to typical sample problems and results are presented. Quantitative comparisons of the results of the optimization method are made to the results obtained by the application of other linear and nonlinear control methods, to illustrate the key features of the methodologies.

2. Spacecraft attitude dynamics and quaternion kinematics

In matrix form, when the inertia matrix is not diagonal the equations of attitude motion of chaser spacecraft are

$$\mathbf{I}\dot{\mathbf{\omega }}+\mathbf{\Omega I\omega }=\mathbf{M}+{\mathbf{M}}_{\mathit{gg}}+{\mathbf{M}}_d,$$ E1

where $\mathbf{I}$ is the moment of inertia matrix which is assumed to be

$$\mathbf{I}=\left[\begin{array}{lll}{I}_{11}& I_{12}& I_{13}\\ I_{12}& I_{22}& I_{23}\\ I_{13}& I_{23}& I_{33}\end{array}\right],\mathbf{\omega }\equiv \left[\begin{array}{l}{\omega }_1\\ {\omega }_2\\ {\omega }_3\end{array}\right],\mathbf{\Omega }=\left[\begin{array}{lll}0& -{\omega }_3& {\omega }_2\\ {\omega }_3& 0& -{\omega }_1\\ -{\omega }_2& {\omega }_1& 0\end{array}\right]$$ E2

${\mathbf{M}}_d$ are the disturbance torques and ${\mathbf{M}}_{\mathit{gg}}$ are the gravity gradient torques.

It is important to emphasize that the targets dynamics are irrelevant to us as there is little or no chance of acquiring the target’s inertia properties. However the target’s angular velocity vector is assumed to be given by ${\mathbf{\omega }}_d$ , its attitude quaternion relative to the chaser’s body frame is assumed to be ${\mathbf{q}}_d$ or its relative attitude quaternion $\mathrm{\Delta }\mathbf{q}$ , relative to the chaser’s body frame, can in principle be measured from within the chaser spacecraft.

Expressions for the gravity gradient moment are obtained assuming that z axis of the spacecraft body is nominally pointing to the Earth. The direction vector the center of gravity of the spacecraft pointing to the Earth is given by the last column of ${\mathbf{T}}_{\mathit{BR}}$ , the transformation from the Earth orbiting frame to the body fixed frame of the spacecraft as

$$\mathbf{c}={\left[c_1{c}_2{c}_3\right]}^T.$$ E3

The corresponding cross product operator ${\mathbf{c}}_{\times }$ is defined as

$${\mathbf{c}}_{\times }=\left[\begin{array}{lll}0& -c_3& c_2\\ c_3& 0& -c_1\\ -c_2& c_1& 0\end{array}\right].$$ E4

Hence the gravity gradient moments acting on the spacecraft and manipulator body are:

$${\mathbf{M}}_{\mathit{gg}}=3n^2{\mathbf{c}}_{\times }\mathbf{Ic}\equiv 3n^2{\mathbf{c}}_{\times }\mathbf{Ic}={\left[\begin{array}{lll}{L}_{\mathit{gg}}& M_{\mathit{gg}}& N_{\mathit{gg}}\end{array}\right]}^T.$$ E5

Thus,

$${\mathbf{M}}_{\mathit{gg}}=\left[\begin{array}{l}{L}_{\mathit{gg}}\\ M_{\mathit{gg}}\\ N_{\mathit{gg}}\end{array}\right]=3n^2\left[\begin{array}{l}{c}_2{c}_1{I}_{31}+c_2^2{I}_{32}+c_2{c}_3\left(I_{33}-I_{22}\right)-c_3{c}_1{I}_{21}-c_3^2{I}_{23}\\ c_3{c}_1\left(I_{11}-I_{33}\right)+c_2{c}_3{I}_{12}+c_3^2{I}_{13}-c_1^2{I}_{31}-c_1{c}_2{I}_{32}\\ c_1^2{I}_{21}+c_1{c}_2\left(I_{22}-I_{11}\right)+c_1{c}_3{I}_{23}-c_2^2{I}_{12}-c_2{c}_3{I}_{13}\end{array}\right].$$ E6

If we express the transformation from the orbiting frame to the body coordinates in terms of an attitude quaternion of the chaser spacecraft with components ${\epsilon }_1$ , ${\epsilon }_2$ , ${\epsilon }_3$ and $\eta$ as

$${\mathbf{T}}_{\mathit{BR}}\left(\mathbf{q}\right)=\left[\begin{array}{lll}{\eta }^2+{\epsilon }_1^2-{\epsilon }_2^2-{\epsilon }_3^2& 2\left({\epsilon }_1{\epsilon }_2+{\epsilon }_3\eta \right)& 2\left({\epsilon }_1{\epsilon }_3-{\epsilon }_2\eta \right)\\ 2\left({\epsilon }_1{\epsilon }_2-{\epsilon }_3\eta \right)& {\eta }^2-{\epsilon }_1^2+{\epsilon }_2^2-{\epsilon }_3^2& 2\left({\epsilon }_2{\epsilon }_3+{\epsilon }_1\eta \right)\\ 2\left({\epsilon }_1{\epsilon }_3+{\epsilon }_2\eta \right)& 2\left({\epsilon }_2{\epsilon }_3-{\epsilon }_1\eta \right)& {\eta }^2-{\epsilon }_1^2-{\epsilon }_2^2+{\epsilon }_3^2\end{array}\right],$$ E7

then from the last column, the Earth pointing direction vector is:

$$\mathbf{c}=\left[\begin{array}{l}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{l}2\left({\epsilon }_1{\epsilon }_3-{\epsilon }_2\eta \right)\\ 2\left({\epsilon }_2{\epsilon }_3+{\epsilon }_1\eta \right)\\ {\eta }^2-{\epsilon }_1^2-{\epsilon }_2^2+{\epsilon }_3^2\end{array}\right].$$ E8

The quaternion kinematics satisfies

$$\frac{d\mathbf{q}}{\mathit{dt}}=\frac{1}{2}{\mathbf{A}}_{\mathbf{\omega }}\left(\mathbf{\omega }\right)\mathbf{q},{\mathbf{q}}^T\cdot \mathbf{q}=\mathbf{1}$$ E9

where the quaternion $\mathbf{q}={\left[\begin{array}{llll}{\epsilon }_1& {\epsilon }_2& {\epsilon }_3& \eta \end{array}\right]}^T$ , consists of a vector part, $\mathit{\epsilon }={\left[\begin{array}{lll}{\epsilon }_1& {\epsilon }_2& {\epsilon }_3\end{array}\right]}^T$ and the scalar $\eta$ so,

$$\mathbf{q}=\left[\begin{array}{l}\mathit{\epsilon }\\ \eta \end{array}\right]\text{and}{\mathbf{A}}_{\mathbf{\omega }}=\left[\begin{array}{ll}-\mathbf{\Omega }\left(\mathbf{\omega }\right)& \mathbf{\omega }\\ -{\mathbf{\omega }}^T& 0\end{array}\right],\mathbf{\Omega }\left(\mathbf{\omega }\right)=\left[\begin{array}{lll}0& -{\omega }_3& {\omega }_2\\ {\omega }_3& 0& -{\omega }_1\\ -{\omega }_2& {\omega }_1& 0\end{array}\right].$$ E10

The quaternion kinematics may also be compactly expressed as

$$\frac{d\mathbf{q}}{\mathit{dt}}=\frac{d}{\mathit{dt}}\left[\begin{array}{l}\mathit{\epsilon }\\ \eta \end{array}\right]=\frac{1}{2}\mathbf{\Gamma }\left(\mathbf{q}\right)\mathbf{\omega },$$ E11

$$\mathbf{\Gamma }\left(\mathbf{q}\right)=\left[\begin{array}{l}\eta {\mathbf{I}}_{3\times 3}+\mathbf{S}\left(\mathit{\epsilon }\right)\\ -{\mathit{\epsilon }}^T\end{array}\right],\mathbf{S}\left(\mathit{\epsilon }\right)=\left[\begin{array}{lll}0& -{\epsilon }_3& {\epsilon }_2\\ {\epsilon }_3& 0& -{\epsilon }_1\\ -{\epsilon }_2& {\epsilon }_1& 0\end{array}\right],$$ E12

where ${\mathbf{I}}_{3\times 3}$ is the $3\times 3$ unit matrix. These relations may be inverted as

$$\mathbf{\omega }=2\left[\begin{array}{ll}\eta {\mathbf{I}}_{3\times 3}+{\mathbf{S}}^T\left(\mathit{\epsilon }\right)& -\mathit{\epsilon }\end{array}\right]{\left[\begin{array}{ll}\dot{\mathit{\epsilon }}& \dot{\eta }\end{array}\right]}^T=2{\mathbf{\Gamma }}^{-1}\left(\mathbf{q}\right){\left[\begin{array}{ll}\dot{\mathit{\epsilon }}& \dot{\eta }\end{array}\right]}^T.$$ E13

The desired attitude quaternion relative to the chaser’s body frame which is assumed to be ${\mathbf{q}}_d$ and the relative attitude quaternion $\mathrm{\Delta }\mathbf{q}$ , relative to the chaser’s body frame are related to the chasers attitude by

$${\mathbf{q}}_d=\mathrm{\Delta }\mathbf{q}\otimes \mathbf{q}.$$ E14

Given two quaternions, $\mathrm{\Delta }\mathbf{q}={\left[\begin{array}{llll}\mathrm{\Delta }{q}_1& \mathrm{\Delta }{q}_2& \mathrm{\Delta }{q}_3& \mathrm{\Delta }{q}_4\end{array}\right]}^T$ , $\mathbf{q}={\left[\begin{array}{llll}{\epsilon }_1& {\epsilon }_2& {\epsilon }_3& \eta \end{array}\right]}^T$ , the quaternion product ${\mathbf{q}}_d=\mathrm{\Delta }\mathbf{q}\otimes \mathbf{q}$ is defined as

$${\mathbf{q}}_d=\left[\begin{array}{l}{q}_{1d}\\ q_{2d}\\ q_{3d}\\ q_{4d}\end{array}\right]=\left[\begin{array}{llll}\eta & {\epsilon }_3& -{\epsilon }_2& {\epsilon }_1\\ -{\epsilon }_3& \eta & {\epsilon }_1& {\epsilon }_2\\ {\epsilon }_2& -{\epsilon }_1& \eta & {\epsilon }_3\\ -{\epsilon }_1& -{\epsilon }_2& -{\epsilon }_3& \eta \end{array}\right]\left[\begin{array}{l}\mathrm{\Delta }{q}_1\\ \mathrm{\Delta }{q}_2\\ \mathrm{\Delta }{q}_3\\ \mathrm{\Delta }{q}_4\end{array}\right]=\left[\begin{array}{llll}\mathrm{\Delta }{q}_4& -\mathrm{\Delta }{q}_3& \mathrm{\Delta }{q}_2& \mathrm{\Delta }{q}_1\\ \mathrm{\Delta }{q}_3& \mathrm{\Delta }{q}_4& -\mathrm{\Delta }{q}_1& \mathrm{\Delta }{q}_2\\ -\mathrm{\Delta }{q}_2& \mathrm{\Delta }{q}_1& \mathrm{\Delta }{q}_4& \mathrm{\Delta }{q}_3\\ -\mathrm{\Delta }{q}_1& -\mathrm{\Delta }{q}_2& -\mathrm{\Delta }{q}_3& \mathrm{\Delta }{q}_4\end{array}\right]\left[\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\\ \eta \end{array}\right].$$ E15

Hence it is expressed in matrix form as

$${\mathbf{q}}_d={\mathbf{C}}_0\mathrm{\Delta }\mathbf{q},{\mathbf{C}}_0=\left[\begin{array}{llll}\eta & {\epsilon }_3& -{\epsilon }_2& {\epsilon }_1\\ -{\epsilon }_3& \eta & {\epsilon }_1& {\epsilon }_2\\ {\epsilon }_2& -{\epsilon }_1& \eta & {\epsilon }_3\\ -{\epsilon }_1& -{\epsilon }_2& -{\epsilon }_3& \eta \end{array}\right],{\mathbf{C}}_0=\left[\begin{array}{ll}\eta {\mathbf{I}}_{3\times 3}+{\mathbf{S}}^T\left(\mathit{\epsilon }\right)& {\mathbf{q}}_{1:3}\\ -{\mathbf{q}}_{1:3}& \eta \end{array}\right].$$ E16

Similarly,

$${\mathbf{q}}_d=\left[\begin{array}{l}{q}_{1d}\\ q_{2d}\\ q_{3d}\\ q_{4d}\end{array}\right]=\left[\begin{array}{llll}\mathrm{\Delta }{q}_4& -\mathrm{\Delta }{q}_3& \mathrm{\Delta }{q}_2& \mathrm{\Delta }{q}_1\\ \mathrm{\Delta }{q}_3& \mathrm{\Delta }{q}_4& -\mathrm{\Delta }{q}_1& \mathrm{\Delta }{q}_2\\ -\mathrm{\Delta }{q}_2& \mathrm{\Delta }{q}_1& \mathrm{\Delta }{q}_4& \mathrm{\Delta }{q}_3\\ -\mathrm{\Delta }{q}_1& -\mathrm{\Delta }{q}_2& -\mathrm{\Delta }{q}_3& \mathrm{\Delta }{q}_4\end{array}\right]\left[\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\\ \eta \end{array}\right]=\left[\begin{array}{ll}\mathrm{\Delta }{q}_4{\mathbf{I}}_{3\times 3}+\mathbf{S}\left(\mathrm{\Delta }{\mathbf{q}}_{1:3}\right)& \mathrm{\Delta }{\mathbf{q}}_{1:3}\\ -\mathrm{\Delta }{\mathbf{q}}_{1:3}& \mathrm{\Delta }{q}_4\end{array}\right]\mathbf{q}.$$ E17

The traditional method of defining the attitude of a spacecraft is by the use of Euler angle sequences. The conversion of Euler angles defined as Euler angle sequences, may be converted to an equivalent quaternion set, using well-known conversion formulae, such as, those given by Smeresky et al. [14].

3. Formulation of the optimal angular rate trajectory synthesis problem

The first task is to formulate the optimal control problem, so it can be solved numerically. This is briefly reviewed. The attitude equations of the spacecraft may be expressed in state space form as

$$\frac{d\mathbf{q}}{\mathit{dt}}=\frac{1}{2}{\mathbf{A}}_{\mathbf{\omega }}\left(\mathbf{\omega }\right)\mathbf{q},{\mathbf{q}}^T\cdot \mathbf{q}=\mathbf{1},$$ E18

$$\mathbf{I}\dot{\mathbf{\omega }}+\mathbf{\Omega I\omega }={\mathbf{M}}_c+{\mathbf{M}}_{\mathit{gg}}+{\mathbf{M}}_d,$$ E19

where ${\mathbf{M}}_c$ is the control torque vector acting on the spacecraft. For our purposes it will be assumed that it can be expressed as, ${\mathbf{M}}_c=\mathbf{I}T\mathbf{W}\mathit{u}$ , where T is the scalar magnitude of the specific torque or torque per unit inertia, $\mathbf{W}$ is a symmetric, positive-definite, torque transformation weighting matrix and $\mathit{u}$ defines the direction of the torque vector. It is the attitude steering control input to the spacecraft. Thus Eq. (19) is

$$\dot{\mathbf{\omega }}={\mathbf{I}}^{-1}{\mathbf{M}}_{\mathit{gg}}\left(\mathbf{q}\right)-{\mathbf{I}}^{-1}\mathbf{\Omega I\omega }\mathbf{+}T\mathbf{W}\mathit{u},$$ E20

with

$$\mathit{u}={\left[\begin{array}{lll}sin\alpha cos\beta & cos\alpha cos\beta & sin\beta \end{array}\right]}^T.$$ E21

When one is interested in the problem of finding the steering control

$$\mathit{u}=\mathit{u}\left(t\right),t_0\le t\le t_f,$$ E22

the torque direction time history is sought, such that it minimizes the cost functional:

$$J=0.5{\left.{\left(\mathbf{\omega }\left(t\right)-{\mathbf{\omega }}_d\right)}^T{\mathbf{Q}}_f\left(\mathbf{\omega }\left(t\right)-{\mathbf{\omega }}_d\right)\right|}_{t=t_f}={\left.\mathrm{\Phi }\left\{\mathbf{\omega }\left(t\right)\right\}\right|}_{t=t_f},$$ E23

subject to, Eqs. (18), (20) and (21).

Introducing the single state vector, $\mathbf{x}={\left[\begin{array}{ll}{\mathbf{q}}^T& {\mathbf{\omega }}^T\end{array}\right]}^T$ , so the Eqs. (18), (20) and (22) are expressed as

$$\frac{d}{\mathit{dt}}{\mathbf{x}}^T=\frac{d}{\mathit{dt}}\left[\begin{array}{ll}{\mathbf{q}}^T& {\mathbf{\omega }}^T\end{array}\right]={\mathit{f}}^T,{\mathbf{q}}^T\cdot \mathbf{q}-\mathbf{1}\equiv g=0$$ E24

To solve the optimization problem, seven Lagrangian multipliers or co-states are introduced given by the two vectors ${\mathit{\lambda }}_{\mathbf{q}}\left(t\right)$ and ${\mathit{\lambda }}_{\mathbf{\omega }}\left(t\right)$ , and a scalar ${\lambda }_c$ denoted by a single column vector, $\mathit{\lambda }\left(t\right)$ . Following Bryson and Ho [15], a Hamiltonian function is defined as

$$H={\mathit{\lambda }}^T\left(t\right){\left[\begin{array}{ll}{\mathit{f}}^T\left(t\right)& g\end{array}\right]}^T.$$ E25

Hence,

$$H={\mathit{\lambda }}_{\mathbf{q}}^T\left(\frac{1}{2}{\mathbf{A}}_{\mathbf{\omega }}\left(\mathbf{\omega }\right)\mathbf{q}\right)+{\mathit{\lambda }}_{\mathbf{\omega }}^T\left({\mathbf{I}}^{-1}{\mathbf{M}}_{\mathit{gg}}\left(\mathbf{q}\right)-{\mathbf{I}}^{-1}\mathbf{\Omega I\omega }\right)+{\mathit{\lambda }}_{\mathbf{\omega }}^{T}T\mathbf{W}\mathit{u}\mathit{+}{\lambda }_{c}g.$$ E26

The necessary conditions (Bryson and Ho [15], Conway [16]) for the first variation of $J$ to be zero include the co-state differential equations

$$\frac{d}{\mathit{dt}}{\mathit{\lambda }}^T\left(t\right)=-\frac{\partial H}{\partial \mathbf{x}}.$$ E27

Explicitly the co-state equations are

$${\dot{\mathit{\lambda }}}_{q_i}\left(t\right)=-\frac{\partial H}{\partial q_i}=-\left\{{\mathit{\lambda }}_{\mathbf{q}}^T\left(\frac{1}{2}{\mathbf{A}}_{\mathbf{\omega }}\left(\mathbf{\omega }\right)\frac{\partial }{\partial q_i}\mathbf{q}\right)+{\mathit{\lambda }}_{\mathbf{\omega }}^T\left({\mathbf{I}}^{-1}\frac{\partial }{\partial q_i}{\mathbf{M}}_{\mathit{gg}}\left(\mathbf{q}\right)\right)\right\},$$ E28

$${\dot{\mathit{\lambda }}}_{{\omega }_i}\left(t\right)=-\frac{\partial H}{\partial {\omega }_i}=-\left\{{\mathit{\lambda }}_{\mathbf{q}}^T\left(\frac{1}{2}\frac{\partial }{\partial {\omega }_i}{\mathbf{A}}_{\mathbf{\omega }}\left(\mathbf{\omega }\right)\mathbf{q}\right)\mathit{-}{\mathit{\lambda }}_{\mathbf{\omega }}^T{\mathbf{I}}^{-1}\frac{\partial }{\partial {\omega }_i}\left(\mathbf{\Omega I\omega }\right)\right\},$$ E29

$${\dot{\lambda }}_c=0.$$ E30

The optimality conditions are

$$\frac{\partial H}{\partial \alpha }={\mathit{\lambda }}_{\mathbf{\omega }}^T\cdot T\mathbf{W}{\left[\begin{array}{lll}cos\alpha cos\beta & -sin\alpha cos\beta & 0\end{array}\right]}^T=0,$$ E31

and

$$\frac{\partial H}{\partial \beta }={\mathit{\lambda }}_{\mathbf{\omega }}^T\cdot T\mathbf{W}{\left[\begin{array}{lll}-sin\alpha sin\beta & -cos\alpha sin\beta & cos\beta \end{array}\right]}^T=0.$$ E32

Hence,

$$\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}=\left|\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}\right|{\left[\begin{array}{lll}sin\alpha cos\beta & cos\alpha cos\beta & sin\beta \end{array}\right]}^T=-\left|\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}\right|\mathit{u}.$$ E33

Thus the two-parameter control vector $\mathit{u}$ , can be expressed as

$$\mathit{u}=-\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}/\left|\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}\right|.$$ E34

The choice of the sign in Eq. (34) will depend on the direction of the desired torque, forward or reverse torque. Thus the closed-loop equations of motion are

$$\frac{d}{\mathit{dt}}\mathbf{x}=\mathit{f},{\mathbf{q}}^T\cdot \mathbf{q}-\mathbf{1}\equiv g=0,\mathit{u}=-\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}/\left|\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}\right|.$$ E35

To complete the definition of the optimal solution, the boundary conditions at $t=t_f$ for the co-state system are found by applying the transversality conditions.

The transversality conditions ensure that the initial and final states are selected optimally within the feasible regions of the states. For the transversality conditions, one may write

$${\mathit{\lambda }}_{\mathbf{q}}\left(t_f\right)={\left.\frac{\mathrm{\partial \Phi }\left\{\mathbf{\omega }\left(t\right)\right\}}{\partial \mathbf{q}}\right|}_{t=t_f}=0,$$ E36

$${\mathit{\lambda }}_{\mathbf{\omega }}\left(t_f\right)={\left.\frac{\mathrm{\partial \Phi }\left\{\mathbf{\omega }\left(t\right)\right\}}{\partial \mathbf{\omega }}\right|}_{t=t_f}={\mathbf{Q}}_f\left(\mathbf{\omega }\left(t_f\right)-{\mathbf{\omega }}_d\right).$$ E37

$${\lambda }_c\left(t_f\right)=0.$$ E38

Thus, ${\lambda }_c\left(t\right)=0$ is a feasible solution. If it can be ensured that the constraint on the quaternion, ${\mathbf{q}}^T\cdot \mathbf{q}-\mathbf{1}\equiv g=0$ is satisfied at each and every integration time step, ${\lambda }_c\left(t\right)$ could be set to zero for all time.

The solution to the state and co-state equations, subject to the optimal control defined by Eq. (34) and the boundary conditions defined by Eqs. (36) and (37), may be found by solving a two point boundary value problem (TPBVP). This can be done using MATLAB’s function, bvp4c.m.

Re-considering Eq. (34), the control input vector may be expressed as

$$\mathit{u}=-\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}/\left|\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}\right|=-r^{-1}{\mathbf{B}}^T\mathit{\lambda },$$ E39

where ${\mathbf{B}}^T$ is a projection matrix relating $\mathit{\lambda }$ to $\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}$ and $r=\left|\mathbf{W}{\mathit{\lambda }}_{\mathbf{\omega }}\right|$ . It is important to recognize that the weighting matrix $\mathbf{W}$ also needs to be chosen. On one hand it provides a set of free parameters so one can construct an optimal solution, while on the other it makes the solution much harder to obtain. Its choice is discussed further in Section 8.

4. Feedback implementation of optimal co-states

Of interest at this stage is to be able to implement the controller, obtained in the last section, as a feedback control law. Thus, inspired by linear optimal control, it can be assumed that, locally, the co-state vector $\mathit{\lambda }$ is proportional to the state vector $\mathbf{x}$ . Hence one may express

$$\mathit{\lambda }=\mathbf{Px}.$$ E40

Furthermore it is assumed that the matrix P is slowly varying and hence does not change appreciably as the time t changes from the current time t to $t\pm \mathrm{\Delta }t$ and to $t\pm n\mathrm{\Delta }t$ , $n=2,3,4$ . Thus the matrix P may be obtained by evaluating $\mathit{\lambda }$ and $\mathbf{x}$ in Eq. (40) at the times $t\pm n\mathrm{\Delta }t$ , $n=0,1,\dots ,4$ , where $\mathrm{\Delta }t$ is a reasonably small time step, over which both $\mathit{\lambda }$ and $\mathbf{x}$ change appreciably. The solution for the matrix P is obtained by differencing the data and solving an over determined system of linear equations by a least squares approach over a moving time window. The matrix P is also constrained to be a symmetric non-negative definite matrix. Moreover $\mathbf{x}$ is expressed as

$$\mathbf{x}=\left(\mathbf{x}-{\mathbf{x}}_{\mathit{ref}}\right)+{\mathbf{x}}_{\mathit{ref}},$$ E41

where ${\mathbf{x}}_{\mathit{ref}}$ is the optimal transfer trajectory for $\mathbf{x}$ obtained by solving the system of equations for the states, parameters and co-states defined in Section 3, which together constitute a TPBVP. By solving the system of equations for the states, parameters and co-states defined in Section 3, one also seeks the trajectory coordinates of the reference trajectory and the total transfer time. It is important that the spacecraft’s attitude controller is able to track the reference trajectory. Consequently the control input is expressed as

$$\mathit{u}=-r^{-1}{\mathbf{B}}^T\mathbf{P}\left(\mathbf{x}-{\mathbf{x}}_{\mathit{ref}}\right)+r^{-1}{\mathbf{B}}^T{\mathit{\lambda }}_{\mathit{ref}},{\mathit{\lambda }}_{\mathit{ref}}={\mathbf{Px}}_{\mathit{ref}}.$$ E42

The second of Eq. (42) is used to obtain the matrix, P. Like in Eq. (34), the choice of the sign in the first of Eq. (42) will depend on the direction of the desired torque, forward or reverse torque.

5. Simplified formulation of the optimal angular rate trajectory synthesis problem

Some authors (for example, see Aghili [4]) have formulated the attitude rate acquisition problem without including the gravity gradient torques. Thus the quaternion kinematics could be ignored. Thus in this case one could set ${\mathit{\lambda }}_{\mathbf{q}}\left(t\right)\equiv 0$ , which results in considerable simplification of the trajectory synthesis problem. The downside of the approach is that the quaternion kinematics is ignored and consequently the quaternion attitude could be quite arbitrary. Quite often after a de-tumbling manoeuver, a precise orientation must be acquired. The required attitude could be acquired in an independent manoeuver and the methodology for this is developed in the next section. The associated tracking problem, which involves tracking the complete state vector set point, must then be separately addressed. Typically this is done by using a barrier Lyapunov function as illustrated by Vepa [17].

6. Optimal attitude orientation acquisition trajectory synthesis

To begin with the quaternion kinematics is given by Eqs. (9) and (11), and can be expressed in one of two alternate forms as.

$$\frac{d\mathbf{q}}{\mathit{dt}}=\frac{1}{2}{\mathbf{A}}_{\mathbf{\omega }}\left(\mathbf{\omega }\right)\mathbf{q}=\frac{1}{2}\mathbf{\Gamma }\left(\mathbf{q}\right)\mathbf{\omega }.$$ E43

In Eq. (43), the angular velocity vector is treated as a control variable and expressed as

$$\mathbf{\omega }\mathbf{=}{\left|\mathbf{\omega }\right|}_{\max }\mathit{u},$$ E44

where the direction vector $\mathit{u}$ is parametrized by an equation similar to Eq. (21). Thus,

$$\mathit{u}={\left[\begin{array}{lll}sin\alpha cos\beta & cos\alpha cos\beta & sin\beta \end{array}\right]}^T.$$ E45

When one is interested in the problem of finding the directional control

$$\mathit{u}=\mathit{u}\left(t\right),t_0\le t\le t_f,$$ E46

the angular velocity direction time history is sought, such that it minimizes the cost functional:

$$J=0.5{\left.{\left(\mathbf{q}\left(t\right)-{\mathbf{q}}_d\right)}^T{\mathbf{Q}}_f\left(\mathbf{q}\left(t\right)-{\mathbf{q}}_d\right)\right|}_{t=t_f}={\left.\mathrm{\Phi }\left\{\mathbf{q}\left(t\right)\right\}\right|}_{t=t_f},$$ E47

subject to, Eqs. (43), (44) and (45). The corresponding Hamiltonian function is

$$H={\mathit{\lambda }}_{\mathbf{q}}^T\left(\frac{1}{2}{\mathbf{A}}_{\mathbf{\omega }}\left(\mathbf{\omega }\right)\mathbf{q}\right)\mathit{=}{\mathit{\lambda }}_{\mathbf{q}}^T\left(\frac{1}{2}\mathbf{\Gamma }\left(\mathbf{q}\right)\mathbf{\omega }\right)\mathit{=}\frac{{\left|\mathbf{\omega }\right|}_{\max }}{2}\left({\mathit{\lambda }}_{\mathbf{q}}^T\mathbf{\Gamma }\left(\mathbf{q}\right)\right)\mathit{u}.$$ E48

The corresponding co-state differential equations are

$$\frac{d}{\mathit{dt}}{\mathit{\lambda }}_{\mathbf{q}}^T\left(t\right)=-\frac{\partial H}{\partial \mathbf{q}}=-\frac{1}{2}{\mathit{\lambda }}_{\mathbf{q}}^T{\mathbf{A}}_{\mathbf{\omega }}\left(\mathbf{\omega }\right).$$ E49

By using an argument similar to the one used in developing Eqs. (31)–(34), the optimal control is given by

$$\mathit{u}=\mathbf{-}{\mathbf{\Gamma }}^T\left(\mathbf{q}\right){\mathit{\lambda }}_{\mathbf{q}}/\left|{\mathbf{\Gamma }}^T\left(\mathbf{q}\right){\mathit{\lambda }}_{\mathbf{q}}\right|=\mathbf{\Gamma }\left(\mathbf{q}\right){\mathit{\lambda }}_{\mathbf{q}}/\left|\mathbf{\Gamma }\left(\mathbf{q}\right){\mathit{\lambda }}_{\mathbf{q}}\right|.$$ E50

For the co-state boundary conditions one has

$${\mathit{\lambda }}_{\mathbf{q}}\left(t_f\right)={\left.\frac{\mathrm{\partial \Phi }\left\{\mathbf{q}\left(t\right)\right\}}{\partial \mathbf{q}}\right|}_{t=t_f}={\mathbf{Q}}_f\left(\mathbf{q}\left(t_f\right)-{\mathbf{q}}_d\right).$$ E51

Once the control is found from Eqs. (50) and (44) is used to define the angular velocity vector and Eqs. (19) and (20) to define the optimal input control torque.

7. Shape based optimal trajectory synthesis

An alternative approach to the optimization based on the integration of co-states is to use a shape based approach as outlined by Caubet and Biggs [18, 19]. For purposes of comparison the shape based approach serves as a useful alternative. In a shape based approach, each of the quaternion components are expressed as a summation of polynomials in terms of a time variable, multiplied by coefficients which may be determined by applying the relevant boundary conditions at the initial and final values of the time variable over a finite time frame. Thus, for example, the quaternion components are expressed as

$$\begin{array}{ll}{q}_i& =q_{i0}+{\dot{q}}_{i0}{t}_f\left(\frac{t}{t_f}\right)+\left({\left.{\dot{q}}_{i0}-{\dot{q}}_i\right|}_{t=t_f}\right)t_f{\left(\frac{t}{t_f}\right)}^2\left(1-\frac{t}{t_f}\right)\\ & +\left(q_{\mathit{\text{if}}}-q_{i0}-{\dot{q}}_{i0}{t}_f\right)\frac{t^2}{t_f^2}\left(3-2\frac{t}{t_f}\right)+\sum _{j=5}^{N+4}{\overline{e}}_{\mathit{ij}-4}\left(\left(j-4\right){\left(\frac{t}{t_f}\right)}^2-\left(j-3\right){\left(\frac{t}{t_f}\right)}^3+{\left(\frac{t}{t_f}\right)}^{j-1}\right),\end{array}$$ E52

where the coefficients ${\overline{e}}_{\mathit{ij}-4}$ are yet to be determined. They are determined by minimizing the cost function

$$J={\left(\underset{0}{\overset{1}{\int }}\left(q_1^2+q_2^2+q_3^2+q_4^2-1\right)d\left(\frac{t}{t_f}\right)\right)}^2.$$ E53

Once all the coefficients of the quaternion components $q_i$ are determined, the angular velocity vector is defined by the inverse of the relation given by Eq. (43) which is

$$\mathbf{\omega }=\left[\begin{array}{llll}{q}_4& q_3& -q_2& -q_1\\ -q_3& q_4& q_1& -q_2\\ q_2& -q_1& q_4& -q_3\end{array}\right]\left[\begin{array}{l}{\dot{q}}_1\\ {\dot{q}}_2\\ {\dot{q}}_3\\ {\dot{q}}_4\end{array}\right].$$ E54

The angular velocity vector $\mathbf{\omega }$ is evaluated for a range of non-dimensional time values between 0 and 1. From the ratio of the maximum of this set, defining the angular velocity time history in terms of the non-dimensional time variable, and the maximum allowable angular velocity magnitude, the length of the time frame $t_f$ over which the control torques must be applied may be found. From the angular velocity vector $\mathbf{\omega }$ the torques that must be applied to the satellite including the gravity gradient torques may be found. From a range of choices for $N$ (say $1\le N\le 6$ ) in the Eq. (52) defining the quaternion, the one that gives the lowest value for $t_f$ is selected.

8. Typical simulation examples

The first example considered the attitude dynamics is defined by Eqs. (18) and (19). Thus the gravity gradient torques acting on the spacecraft are included in the dynamic model and they are responsible for coupling the attitude quaternion kinematics and the angular velocity dynamics. The objective is to spin the spacecraft so the final angular velocity vector is given by $\mathbf{\omega }={\left[\begin{array}{lll}1& 1& 1\end{array}\right]}^T\mathit{\text{rads}}/s$ . The initial angular velocity vector is $\mathbf{\omega }={\left[\begin{array}{lll}0& 0& 0\end{array}\right]}^T$ . The spacecraft is fitted with magnetic torque actuators and the maximum three axis torques are limited to ${\mathbf{T}}_c={\left[\begin{array}{lll}0.62& 1& 1\end{array}\right]}^T\mathit{Nm}$ . The diagonal non-zero elements of the weighting matrix $\mathbf{W}$ in Eq. (34) and the principal moment of inertia of the spacecraft are respectively given by

$${\mathbf{W}}_{\mathit{diag}}=\left[\begin{array}{lll}0.24& 0.9& 1\end{array}\right],\mathbf{I}=\left[\begin{array}{lll}2.27293& 3.27331& 0.3989\end{array}\right]{\mathit{kgm}}^2.$$ E55

In all examples the solution of the TPBVP is done using MATLAB’s function, bvp4c.m. Whenever there was a need to solve an initial value problem, the equations were integrated using MATLAB’s ode45.m.

In all cases, the time variable was made non-dimensional so it raged from [0–1]. The integration time step was chosen to be relatively small initially ( $\mathrm{\Delta }t=0.0001$ ), and automatically and iteratively reduced linearly as the final time was approached. The iterations were terminated when no further improvement in the accuracy of the predicted final time was feasible.

Figure 1 shows the time history of the reference quaternion components and Figure 2 shows the corresponding angular velocities (p, q, r). Figures 3 and 4 show the corresponding, actual, quaternion components and the corresponding angular velocity components, where an approximate optimal linear feedback law based on Eq. (42) is used to track the reference trajectory.

Figure 1.

Time history of the reference quaternion components for the first example.

Figure 2.

Time history of the corresponding reference velocity components for the first example.

Figure 3.

Time history of the actual quaternion components tracked by the feedback controller for the first example.

Figure 4.

Time history of the actual velocity components tracked by the feedback controller for the first example.

In the next example, the simplified attitude dynamics is used with the gravity gradient torques neglected. This decouples the angular rate dynamics from the attitude quaternion dynamics, which need to be considered for synthesizing the reference trajectories. The reference angular velocities are then integrated to obtain the spacecraft’s quaternion attitude time history. Figure 5 shows the reference angular velocity components.

Figure 5.

Time history of the reference velocity components for the second example.

Figure 6 shows the errors in the actual angular velocity components when compared with corresponding reference values and Figure 7 the corresponding quaternion components. Figure 8 shows the attitude in terms of the Euler axis and the Euler principal angle components.

Figure 6.

Time history of the errors in the actual velocity components corresponding to Figure 5, tracked by the feedback controller.

Figure 7.

Time history of the corresponding quaternion components tracked by the feedback controller for the second example.

Figure 8.

Time history of the attitude components in terms of the Euler axis and the Euler principal angle components for the second example.

In the final example it is desired to alter the attitude quaternion of the spacecraft, so as to point the spacecraft in a desired direction. In this case on the quaternion kinematics defined by Eq. (44) are used. The maximum angular velocity of the spacecraft is assumed to be limited to 0.001 rad/s. The desired pointing direction is specified as a yaw, roll pitch Euler angle sequence given by $\left[\begin{array}{lll}\psi & \theta & \varphi \end{array}\right]=\left[\begin{array}{lll}22°& 25°& 30°\end{array}\right]$ corresponding to the components of the quaternion ${\mathbf{q}}_d=\left[\begin{array}{llll}0.2068& 0.2518& 0.1682& 0.9303\end{array}\right]$ .

In Figure 9 are shown the reference optimal quaternion components and in Figure 10 are shown the corresponding angular velocity components. The optimum torque components required to affect the attitude change are shown in Figure 11. These include compensation for the gravity gradient torques. They show that they could be easily achieved by low thrust electric actuators such as electro-spray thrusters.

Figure 9.

Time history of the reference quaternion components for the final example.

Figure 10.

Time history of the corresponding reference velocity components for the final example.

Figure 11.

Time history of the corresponding reference torque components for the final example.

When small reaction wheels are used Eq. (19) may be modified to include the momentum of the wheels and the control inputs to the wheels could also be estimated. If reaction wheels are used much larger torques are possible and the time over which they are used could be shortened. In Figure 12 are shown the attitude time history components in terms of the Euler axis and Euler principal angle components.

Figure 12.

Time history of the attitude components in terms of the Euler axis and the Euler principal angle components for the final example.

This example is also solved using the shape based approach discussed briefly in Section 7. In Figure 13 are shown the required applied torque components obtained by the shape based approach with $N=1$ in Eq. (52). The time frame over which the control must be applied is $t_f=733.3s$ , which is the lowest for all $N$ considered and is about the same as the time required by the approach based on the integration of co-states. However the reference torque components shown in Figure 13, are much larger than those plotted in Figure 11, and for this reason, they are not referred to as ‘optimum torques.’ The corresponding gravity gradient torques, acting on a satellite orbiting the Earth at the geostationary orbit radius, are also shown in Figure 14. Although the gravity gradient torques are of the same orders of magnitude as the reference optimum torque components in Figure 11, they are much smaller than the corresponding torque components obtained by the shape based approach and shown in Figure 13. It must be recognized that the gravity gradient torques become much larger as the spacecraft orbits the Earth at a much closer orbit radius.

Figure 13.

Time history of the applied torque components for the final example, obtained by the shape based approach.

Figure 14.

Time history of the gravity gradient torque components for the final example, obtained by the shape based approach.

9. Discussion and conclusions

A close examination of the results in Figures 9 and 12 shows that while the quaternion component time histories are not linearly varying, the Euler axis and Euler principal angle components are almost linear. This allows for linear extrapolation of the trajectories if when desired. It is also observed that the acquisition of the Euler axis is relative fast in comparison with growth rate of the Euler angle which is relatively slower. The kinematics of the Euler axis seems to represent a fast sub-system while the growth of the Euler angle represents the slow sub-system. This observation, facilitates the construction of approximate sub-optimal trajectories where in the Euler axis is acquired instantly and the Euler angle increases or decreases linearly with time. Once a sub-optimal solution can be defined in terms of interpolating polynomial it can also be further optimized by using Eq. (13) and shape based methods as those proposed by Caubet and Biggs [18, 19], quite rapidly and if need be, by the pseudo spectral method, or other direct collocation methods. The advantage of further optimization using shape based methods is that the precise shape of the desired output could be achieved avoiding overshoot. However depending on the choice of the output shape function, the control could be restrictive and so the magnitudes of the torques required could be much larger in comparison with the co-states approach.

For the preceding example, where a set angular velocity components were desired, the velocity components and the corresponding quaternion components are shown in Figures 5, 6, 7 as they vary with time. Figure 8 also shows that the Euler axis and Euler principal angle components are varying as quadratic functions of time. A similar conclusion cannot be drawn as far as the angular velocity components and the components of the quaternion. This again is extremely useful in applying low order polynomials for developing formulae for extrapolating the optimal trajectories, by converting the quaternion components to the domain of the Euler axis and Euler principal angle components. It also facilitates the integration of various optimal segments into a single trajectory over an extended time frame.

From the first example, comparing Figures 1 and 2 with Figures 3 and 4, it is seen that the optimum tracking feedback control law obtained by linearly approximating the relationship between the states and co-states by Eq. (40), performs well. The errors between these two sets of trajectories, the reference trajectory and the actual tracked trajectory, are always within 5% of the corresponding reference value, over the time frame of the plots.

In this paper, either the simplest form of the attitude dynamics or the basic kinematic equations alone are used to construct the optimal trajectories. The required control torques are obtained from the inverse dynamic relations. The usefulness of transforming the attitude representation to the Euler axis and Euler principal angle components, as it facilitates the application of low order polynomials for the construction of approximate sub-optimal trajectories, is demonstrated. Furthermore it is shown how optimal feedback control laws may be constructed from the solution for the optimal trajectories, for tracking the reference trajectories.

Conflict of interest

There are no conflicts of interest to declare.