Derivation of All Attitude Error Governing Equations for Attitude Filtering and Control

This article presents the full analytical derivations of the attitude error kinematics equations. This is done for several attitude error representations, obtaining compact closed-forms expressions. Attitude error is defined as the rotation between true and estimated orientations. Two distinct approaches to attitude error kinematics are developed. In the first, the estimated angular velocity is defined in the true attitude axes frame, while in the second, it is defined in the estimated attitude axes frame. The first approach is of interest in simulations where the true attitude is known, while the second approach is for real estimation/control applications. Two nonlinear kinematic models are derived that are valid for arbitrarily large rotations and rotation rates. The results presented are expected to be broadly useful to nonlinear attitude estimation/control filtering formulations. A discussion of the benefits of the derived error kinematic models is included.


Introduction
Many attitude representations are available for modeling problems in science and engineering [1][2][3][4][5][6][7][8][9]. Nonlinearity of the representation of a given physical motion and location of geometric singularities are dependant on (1) the true motion, (2) the attitude representation selected, and (3) the particular choice of estimated axes. Selecting the appropriate representation is highly linked with the kind of problem being considered. The most popular attitude parameterizations are as follows [1,[10][11][12]: Definitions, characteristics, and transformations between these representations can be found in many references [1,[4][5][6]8,10,13]. For applications requiring large and rapid rotational dynamics, there exists a need for developing attitude error kinematic models that exactly describe arbitrary large rotational motions. In particular, for control problems, Markley [8] has considered different attitude error representations for estimating the state of a maneuvering spacecraft. He has clarified the relationship between the four-component quaternion and the Multiplicative Extended Kalman Filter. Cao et al. [14] developed an unscented predictive filter for satellite formation. Later, Cao et al. [15] proposed a huber-based Kalman filter for the attitude estimation problem of small satellites. Batch and Paielli [16] investigated the rigid-body attitude error using inversion techniques to obtain a robust linearized model of attitude control error dynamics. Sun et al. have derived relative rotational and transnational dynamics error for spacecraft rendezvous control [17]. Crassidis et al. [18] studied a variable-structure control strategy for maneuvering vehicles. In their work, they used a feedback linearizing technique and added an additional term to the spacecraft maneuvers to deal with model uncertainties, which they demonstrated always provides an optimal response. Ahmed et al. [19] extended previous work to consider adaptive asymptotic tracking during maneuvers while estimating inertia properties. They used a Lyapunov argument to generate an unconditionally robust control law with respect to their assumed parametric uncertainty. Bani Younes et al. [20,21] considered generalized optimal control formulations that handle nonlinear system dynamics and enable the development of control gain sensitivities to handle plant model uncertainties during maneuvers. Sharma and Tiwari [22] introduced MRP for parameterization of the orientation error. They defined the attitude error as an additive quantity. Their work is extended by retaining a rigorous nonlinear MRP-based error equation. Schaub et al. [23] developed a new penalty function for optimal control formulation of spacecraft attitude control problems. This function returns the same scalar penalty for a given physical attitude regardless of the attitude coordinate choice. A recent work by Tanygi investigated the projective geometry of vectorial attitude parameterizations with applications for control [24] and estimation [25]. Junkins [26] discussed the link between designing a good controller and the choice of coordinates to represent the attitude kinematics. He linearized the attitude error equations by defining the departure motion as an additive error from a nominal trajectory. Unfortunately, the error in orientation cannot be rigorously represented as additive (linear) because of the nonlinear behavior of underlying kinematical descriptions [2].
Error equations are challenging to define for quaternion variables because of the implied coupling effects between the quaternion components and the implicit norm constraint. The position error is linear and, therefore, can be described by the distance between the two vectors representing the true and the estimated position states. Unfortunately, the error in orientation in not linear and, therefore, must be described by a matrix product (corrective attitude matrix [2]). An exact quaternion error equation is defined as the quaternion product between the true quaternion and the inverse quaternion of the estimated (or reference) state. The significant advantage of this approach is that the norm constraint is preserved and, most importantly, that the quaternion error product remains valid for arbitrarily large relative rotations. This fact is exploited to develop equivalent arbitrarily large relative rotational representations for the vehicle attitude motion. The resulting expressions are compact, accurate, and computationally efficient.
The primary goal of this paper is to develop uniformly valid kinematics equations to describe the time-evolution of the attitude error. References [1,[3][4][5]10] contain good review for the attitude kinematics equations. This paper presents compact nonlinear attitude error kinematics equations that can be used in attitude control and/or estimation dynamics problems. Exact analytical (closed-form and no approximations) attitude error kinematic equations are derived for most popular and known attitude representations. This paper extends previous work originally initiated by the authors on developing attitude error kinematics [27][28][29][30], where the estimated angular velocity is defined in the true attitude axes frame. It builds on our initial findings and extends the formulations to address more detailed developments and to present two distinct error kinematics approaches treated with more compact depth and insights in the derivations. It also presents simulation examples that demonstrate the applications in attitude filtering and tracking control. The multiplicative Kalman filter is discussed to enable the use of these attitude error kinematics for two different coordinate choices. The resulting expressions have been optimized to obtain the most compact and computationally efficient forms.
Also, the applications of these formulations are discussed by solving the open-loop optimal spacecraft tracking control problem.
The major contribution of this paper is the complete development of attitude error kinematics equations where the estimated angular velocity is defined in either the true attitude axes frame or in the estimated attitude axes frame. These equations can be used for arbitrarily large relative rotations and rotation rates. These attitude errors represent rotations between estimated and true attitudes. Numerical integrations of these kinematic equations are performed to validate the machine error accuracy for each attitude representation. Singularities and constraints are discussed for minimum and non-minimum attitude parameter representations, respectively. Applications are expected in rotational dynamics problems for both nonlinear attitude estimation filtering and attitude tracking.

Attitude Error Kinematics
With reference to Figure 1, let C be the true attitude matrix (DCM) of which the axes are [x, y, z] (dashed, black) and letĈ be the estimated attitude matrix of which the axes are [x,ŷ,ẑ] (dashed, red). The relationships between the DCM with the corresponding Euler's angles associated with the rotation axes and sequences are given in Table A1. The attitude error is described by the corrective attitude matrix: In fact, the product between δC and the estimated attitudeĈ provides the true attitude: In particular, matrix δC is the transformation matrix for vectors from the estimated [x,ŷ,ẑ] frame to the true [x, y, z] frame.
Let ω be the true angular velocity vector (black solid in Figure 1) of which the elements are defined in the [x, y, z] axes of the true attitude. Now, considerω to be the estimate of the true angular velocity vector (red solid in Figure 1). Two distinct cases appears: (1) The estimatedω is defined in the true body axes frame (x, y, z). In this case, ω andω are defined in the same reference frame. Therefore, the angular velocity error vector is as follows: This case finds application in simulations when the true attitude is known.
(2) The estimatedω is defined in the estimated body axes frame (x,ŷ,ẑ). In this case, ω andω are defined in two different reference frames. Therefore, the angular velocity error vector is as follows: where δC is provided by Equation (1). This case is important in real applications as the true attitude is always unknown and the estimated angular velocity can be defined in the estimated attitude frame only.
The true attitude dynamics is defined by ω and C. They satisfy the DCM kinematic equations: The estimated angular velocity satisfies the dynamic equation: where I ∈ 3×3 is the moments of inertia tensor and u ∈ 3 is the torque acting on the system. We assume that the system parameters are deterministic.
To develop the attitude error kinematics equations, two distinct analysis must be performed. These two derivations are associated to the two angular velocity errors definitions provided in Equations (2) and (3), respectively. The following two sections derive the attitude error kinematics equations valid for arbitrarily large angles and rate errors.

Simulation Case: Estimated Angular Velocity Given in the True Attitude Frame
In this section, the estimated and the true angular velocities are expressed in same coordinate frames, δω = ω −ω. This approach is suitable in simulations where the true body attitude is known. When using the same reference frame for true and estimated angular velocities, the angular velocity error dynamics equation is written as follows [21,27,28]: In the following subsections, the attitude error kinematic equations are derived for the most important attitude representations and for arbitrarily large rotational angles and angular rates. These equations are mathematically simple and compact. They can be used, for instance, to validate novel control theories and/or attitude estimation filtering techniques.

Quaternion Error Kinematics
The quaternion error is a four-dimensional vector defined as follows: , e is the principal axis, and φ is the principal angle. The kinematic solution for the true quaternion trajectory defines the desired relative rotational motion for the spacecraft. Equation (1) written in term of quaternion is as follows: whereq −1 is the inverse of the estimated quaternion rotation and ⊗ represents the quaternion product. Note that the error, δq, is also a quaternion, that is, a unit vector representing the rotation from the estimated axes to the true axes. This expression is valid for arbitrarily large rational motions and provides a foundation for developing all other kinematic variable generalizations presented. We follow reference [1] in writing the quaternion product as follows: where q and q represent two arbitrary quaternions. The quaternion error rate is as follows: Quaternion kinematics evolves according to the kinematic equation: The derivative of the identityq ⊗q −1 = {0, 0, 0, 1} T leads toq ⊗q −1 +q ⊗q −1 = 0. Hence, the inverse quaternion evolves:q where Γ(ω) is the estimated angular velocity matrix. Substituting Equation (7) and Equation (8) into Equation (6), yields the following: This allows us to write the following: The quaternion error rate equation becomes the following: Now, by substituting the angular velocity error given in Equation (2) into Equation (9), the bilinear differential equation for the tracking error kinematics becomes the following: whereΓ(ω) is a matrix defined as follows: Equation (10) can be rewritten in the following compact form: This equation can be split into the scalar and vector part of the quaternion as follows [28]: These above equations are exact for arbitrarily large and rapid relative rotational motions.

Rodrigues Parameter Error Kinematics
Rodrigues parameters are a minimum attitude parametrization. The RP vector is defined in terms of quaternion parameters as follows [10]: where the inverse transformation is given as follows: where ρ 2 = ρ T ρ. Note that the attitude error given in Equation (12) is represented in the quaternion form. An exact nonlinear model for the quaternions with no approximation is used to preserve the quaternion unit norm. This implies that the quaternion error still represents a finite orientation that can be mapped to any other attitude representations. Here, the quaternion error to RP using Equation (13) is mapped. Thus, the RP error vector is simply expressed as follows: The inverse mapping for quaternion variables follows: where δρ 2 = δρ T δρ. The RP error differential equation is obtained by taking the derivative of Equation (14) and substituting Equations (12) and (15), which yields the following [28]: Equation (16) provides the desired RP attitude error nonlinear kinematic differential equation.

Modified Rodrigues Parameter Error Kinematics
Modified Rodrigues parameters are an elegant addition to the family of attitude representations. MRP vector is defined in terms of the quaternion parameters by the following [10] The inverse transformation is given by the following: Similarly, since the attitude error in Equation (12) is expressed by quaternion, the mapping into MRP can be performed using Equation (17). Thus, MRP error vector is expressed as follows: The inverse mapping for quaternion variables follows: The MRP error differential equation is obtained by taking the derivative of Equation (18) and substituting Equation (12) and Equation (19), which yields the following compact, nonlinear, third-order form [28]: Equation (20) provides the exact MRP attitude error kinematic differential equation.

Euler Angles Error Kinematics
The most famous attitude representation is described by three angles, known as Euler angles, (θ 1 , θ 2 , θ 3 ) associated with subsequent rotations about three coordinate axes. The variations of these angles represent the attitude error of δC. There are several conventions for Euler angles, depending on the axes about which the rotations are carried out. In the following, we assume the rotation is 3-1-3 (yaw-roll-yaw) sequences, then we generalize the expression for arbitrary rotation sequence. We start from the mapping equations from quaternion to Euler angles [6]. For the 3-1-3 set, the transformation is as follows: Differentiating Equation (21), we obtain the following: Thus, Euler angle rates can be written in the following least-squares solution: Substituting Equation (12) and making use of Equation (21), we obtain the following: In general, for the generic i-j-k Euler angles sequence, the formula is as follows [28]: Equation (22) represents the attitude error kinematic equation using any Euler angle sequence.

Principal Axis and Angle Error Kinematics
Any rigid-body rotation can be obtained by a single rotation about a principal axis, e, by a principal angle, φ. To derive the kinematics of the principal axis/angle for attitude error, we start from the definition of the quaternion: Taking the derivative of Equation (23), substituting Equation (12), and solving for δφ and δė, we obtain the following [28]: where δê = δe Equations (24) and (25) are the desired nonlinear kinematic differential equation of the attitude error using principal axis and principal angle.

Direction Cosine Matrix Error Kinematics
The DCM error can be written as follows: Performing the derivative of Equation (26) and making use of the attitude kinematics, C = − [ ω×] C, and the DCM inverse identity,˙Ĉ T = −Ĉ TĊĈ T , we obtain the following [28]:

Cayley-Klein Error Parameters Kinematics
Cayley-Klein parameters is an attitude representation provided by a 2 × 2 complex matrix. This matrix is as follows: where angle can be computed from the following: Rewriting Equation (28) in column form, we obtain the following: where Ψ 0 is a non-singular constant matrix. Differentiating Equation (29), substituting Equation (12), and using q v = Ψ −1 0 col(δK), we obtain the following [28]: where the col(δK) has to satisfy the constraint equation

Estimation/Control Case: Angular Velocity Estimated in the Estimated Attitude Frame
In this section, the estimated and the true angular velocities are expressed in different coordinate frames: δω = ω − δCω.
This equation explicitly accounts for the obvious truth that the estimated and true angular velocity vectors are expressed in different coordinate frames. This definition explicitly computes the angular velocity error in the current body frame. This kinematic variable definition leads to the following angular velocity error dynamics equation: It is clear that the angular velocity error dynamics equation is coupled with the attitude solution. This is expected as the estimated angular velocity has to be mapped to the true frame using the attitude corrective matrix, δC.

Direction Cosine Matrix Error Kinematics
To investigate the new attitude error parametrization for the estimation/control case, the direction cosine matrix is differentiated for the attitude error given in Equation (1): Recall the attitude kinematics of the current motionĊ = −[ ω×]C and the estimated motioṅĈ = −[ω×]Ĉ. Equation (32) can be written: The angular velocity error is δω = ω − δCω. Therefore, we obtain the following: Using the transformation of skew-symmetric tensor identity [4], we obtain the following: Also, since δC is an orthogonal matrix, Equation (34) reduces to the following: It is interesting that the attitude error kinematics equation is similar to the attitude kinematics equation. This result agrees with the demonstration presented in Tanygin's work [24]. Hence, the attitude error kinematics using other attitude parameterizations can be simply obtained.

Quaternion Error Kinematics
The quaternion error is a four-dimensional vector defined as follows: The transformation from quaternion error, δq, to DCM error, δC, is given by the following: The inverse transformations is as follows: This transformation shows the redundancy in selecting the sign of the scalar parameter, δq 4 . The positive sign is associated to δφ ≤ π, while the negative is associated to π < δφ ≤ 2π.
The kinematics differential equation of the quaternion error is derived by differentiating Equation (38): Substituting the expressions of the diagonal elements in Equation (40) yields the following: Similar derivations are applied on δq 1 , δq 2 , and δq 3 in Equation (39). The kinematic error differential equation for quaternion is then the following: where I 3×3 is the 3 × 3 identity matrix.

Rodrigues Parameter Error Kinematics
Rodrigues parameters error vector is simply expressed by the following: Also, the inverse mapping for quaternion variables follow: where δρ 2 = δρ T δρ. The differential equation for the RP error follows by taking the derivative of Equation (43): Substituting Equation (42) into Equation (45) yields the following: This equation can be further simplified by substituting Equation (44) to yield the compact quadratic form: Equation (46) is the RP attitude error kinematic differential equation. The above form approaches singularity when |δφ| → ±π.

Modified Rodrigues Parameter Error Kinematics
Modified Rodrigues parameters error vector is expressed as follows: The inverse mapping to quaternion is as follows: The differential equation for the MRP error follows from the derivative of Equation (47): Substituting Equation (42) into Equation (49) yields the following: This equation is simplified by substituting Equation (48) to obtain a compact quadratic form: Equation (50) provides the MRP attitude error kinematic differential equation. The above form approaches singular behavior as |δφ| → ±2π.

Euler Angles Error Kinematics
There are 12 distinct Euler angles sequences, depending on the axes about which the rotations are carried out. Review the development in Section 2.1.4; similar derivations are followed here that yield to the general i-j-k Euler angle error kinematics equations: Equation (51) is the desired kinematic differential equation using Euler angles to represent the attitude error.
Alternatively, the Euler angles error kinematics can be derived by writing the angular velocity error vector in the body frame in terms of the Euler angles error rates. For example, the angular velocity error vector for 3-1-3 rotation sequence is given by the following: where R 313 (δθ 1 , δθ 2 , Equation (52) can be written in the compact form: where The kinematic differential equation of 3-1-3 Euler angles error is the inverse of Equation (54): The complete set of 12 Euler angles error kinematics in terms of M −1 ijk and the angular velocity error vector is provided in Appendix A.

Principal Axis and Angle Error Kinematics
To derive the kinematics of the principal axis/angle for attitude error, we start from the definition of the quaternion: δq 4 = cos(δφ/2) and δq v = δe sin(δφ/2).

Cayley-Klein Error Parameters Kinematics
Cayley-Klein parameters is an attitude representation provided by a 2 × 2 complex matrix. This matrix is as follows: where σ 1 , σ 2 , and σ 3 are the three Pauli spin matrices. The principal angle can be computed from the following: It has been already shown in (Section 2.1.7) that the Cayley-Klein parameters can be written in the column form: Substituting Equation (42) into Equation (58) and making use of col(δK) = Ψ 0 δq, the kinematics equation of the Cayley-Klein error column vector yields the following: where the col(δK) has to satisfy the constraint equation col(δK) T col(δK) = 2. Appendix B summarizes the attitude error kinematics using the two definitions of the angular velocity error.

Numerical Validation
Numerical integrations of all the kinematic equations presented are performed for arbitrary initial conditions to validate the machine error accuracy for each attitude representation. For this particular example, the initial conditions are given in Table 1, where constant angular velocity is considered. The validation test is performed as follows: 1.
Integrate the true attitude kinematics for each attitude representation using the attitude kinematics equations; 2.
Transform the time history solution of the true attitude for each attitude representation to the direction cosine matrix C or quaternion q; 3.
Integrate the estimated attitude kinematics for each attitude representation using the attitude kinematics equations; 4.
Transform the time history solution of the estimated attitude for each attitude representation to the direction cosine matrixĈ or quaternionq; 5.
Calculate the attitude error of the two solutions in steps 2 and 4 using δC = CĈ T or δq = q ⊗q −1 ; then, calculate the principal angle of the attitude error, δφ * ; 6.
For approaches 1 and 2, integrate the attitude error kinematics for each attitude representation using the attitude error kinematics equations in Appendix B (Table A1); then, calculate the principal angle of the attitude error δφ for each approach; 7.

Kalman Filter
In this section, a sequential extended Kalman filter (EKF) formulation is developed and presented for the two approaches. The development of compact forms of nonlinear attitude error kinematics enables the applications of arbitrarily large relative rotations and rotation rates in different attitude coordinates. Since quaternions present no singularities, it is the most popular coordinate for attitude estimation. However, it must obey the norm constraint.

Estimated Angular Velocity Defined in the True Attitude Frame
The estimated angular velocity vectorω is defined in the true body axes frame (x, y, z). In this case, ω andω are defined in the same reference frame. Therefore, the angular velocity error vector is as follows: This case finds application during simulations when the true attitude is available. The quaternion error kinematics is given by Equation (12), which can be linearized for first-order approximation to the following [8]: Unfortunately, this linearization requires that the estimated quaternion is close to the true quaternion to preserve the normalization constraint. A common sensor model to measure the angular rate is the rate-integrating gyro, which is defined by the following: where η v ∼ N (0, σ 2 v I 3×3 ) and η u ∼ N (0, σ 2 u I 3×3 ) are zero-mean Gaussian noise with variances given by σ 2 v and σ 2 u , respectively; β is the gyro bias vector; andω is the measured angular velocity, which is given by the following:ω =ω +β, whereω andβ are the estimated angular velocity and the estimated gyro bias vectors, respectively. Hence, the angular velocity error is given by the following: where δβ = β −β. The multiplicative extended Kalman filter error model is now given by the following (note:˙β = 0):

Estimated Angular Velocity Defined in the Estimated Attitude Frame
In this section, the estimated and the true angular velocities are expressed in different coordinate frames, δω = ω − δCω. It has been discussed earlier that the attitude error kinematics equations for this approach take the same form as the attitude kinematics equations; see Figure A1. The quaternion error kinematics is linearized for first-order approximation to the following: where the first-order approximation of the attitude error matrix is given by the following: Considering the same rate-integrating gyro sensor model defined in Equation (61) and the observed angular velocity definition in Equation (62), the angular velocity error is written as follows: where δβ = β − δCβ. The multiplicative extended Kalman filter error model is now given by the following (note:˙β = 0):

Extended Kalman Filter Error Model
For a single sensor, we define the true and estimated body vectors as follows: where r = {x, y, z} T is the 3 × 1 vector of the vehicle position in the earth-centered, earth-fixed (ECEF) coordinates, the true attitude matrix is C(q) = C (δq)C(q − ), C(δq) is defined in Equation (37), and Thus, the estimation error of the body vector is defined as follows: The sensitivity matrix of n measurement sets is given by the following: The EKF error model is given by the following: where and F(x(t), t), G(t), and Q(t) are given by Table 2: [β×] Discrete-time attitude observation n × 1 model at time t k is given by the following: where v k ∼ N (0, R k ) is zero-mean gaussian measurements noise with covariance error matrix R k . Thus, the error state update is given by the following: where ∆x is the estimate output, and K k is Kalman gain, as given in Equation (75).
The EKF is implemented following the below sequential steps [30,31]: • Initialization: at t 0 and for given initial states x 0 = [q 0 , β 0 ] T and initial value of the covariance matrix P 0 , the initial values are given by the following: where the superscript (−) denotes priori values and E{ } is the expectation operator. Assume that x(0) ∼ N (0, P(0)).
• Gain: compute the Kalman gain matrix: where • Update: update the state estimatex + k and covariance P + k at each measurement: where the superscript (+) denotes posteriori values. • Propagation: propagate both the state estimatex k and covariance P k using the posteriori estimatê x + k and posteriori covariance P + k . The estimated angular velocity,ω =ω −β, is used to propagate the quaternion kinematics:˙q where the matrices F, G, and Q are given by Table 2.

Numerical Simulation
This section presents simulation results that utilize data provided by a star tracker to estimate the attitude of slow spinning spacecraft under the following scenario: • Optical axis aligned withẑ axis of the body reference frame. The directions captured by the sensor field of view during the whole trajectory is shown in Figure 3. Figure 4 shows the attitude estimation error, which is defined as the principal angle between estimated and true DCM and is small and within the predicted covariance (±3σ) error.

Optimal Tracking Control
A general scalar nonnegative attitude penalty function is utilized to formulate an optimal feedback control for the spacecraft tracking problem. This new variable yields identical performance index values, regardless of the attitude variables selected. The general final finite-time optimal control formulating is given by minimizing the following performance index [29]: which is subject toẋ = [δω T δζ T ] T = f (δω, δζ, u), where δω is the angular velocity error, δζ is an arbitrary attitude representation of the attitude error, and the penalty functions are (using δω(t f ) = δω t f and δζ(t f ) = δζ t f ): The weights Q 1 and Q 3 are scalars, and the weights Q 2 , Q 4 , and R are 3 × 3 matrices. The f (δω, δζ, u) is the spacecraft error dynamics. The scalar function g(δζ) is a general nonnegative attitude penalty function. The function is chosen to produce the same cost for a given physical orientation [23,29]: when using exact nonlinear attitude error kinematics, the orientation will work for large angles; δφ = ±180 • . Therefore, the function is bounded 0 ≤ g(δC) ≤ 1 for all possible rotations. Thus, the attitude cost reaches its highest value at the maximum rotation angle. Defining the attitude cost function in terms of the DCM makes it universally valid for arbitrary choice of attitude coordinates. It can be simply parameterized by any other attitude coordinate. The universal quadratic penalty function for arbitrary attitude error representations is given in Table 3. This penalty function returns the same cost for a given physical attitude, thereby eliminating the dependency of the optimal control solution on the choice of the attitude coordinate. CRPs g(δC(δρ)) = δρ T δρ 1 + δρ T δρ Principal angle/axis g(δC(δê, δφ)) = sin 2 (δφ/2) Cayley-Klein g(δC(δK)) = 1 − 1 4 tr(δK) 2 * For the 3-1-3 rotation sequence.

Reference Angular Velocity Defined in the Body Attitude Frame
The reference ω r and current ω angular velocities are expressed in the same coordinate frames; i.e., δω = ω − ω r . It can be shown that the general expression of the attitude error kinematics of this set can be given as follows: with initial state δζ(t 0 ) = δζ 0 . The Hamiltonian H for this system of equations is defined, for the given optimal control problem in Equations (79) and (80), as follows: where λ δω and λ δζ are the co-state variables for the angular velocity error and the attitude error, respectively. Invoking the standard necessary condition for optimality, the co-state differential equations are given by the following: The two co-state differential equations are integrated backward in time with given terminal The optimal control is given by the first-order necessary conditions for an extremum, H u = 0, leading to u = −(IR) −1 λ δω .

Reference Angular Velocity Defined in the Reference Attitude Frame
The reference and current angular velocities are expressed in different coordinate frames, i.e., δω = ω − δCω r . Therefore, this definition explicitly computes the angular velocity error in the current body frame. As an important result of this set, the expressions of the attitude error kinematics follow the attitude kinematics equation for any given attitude representation choice. Therefore, the general expression of the attitude error kinematics of this set can be given as follows: with initial state δζ(t 0 ) = δζ 0 . The Hamiltonian H for this system of equations is as follows: Invoking the standard necessary condition for optimality, the co-state differential equations are given by the following: The optimal control is given by the first-order necessary conditions for an extremum, H u = 0, leading to u = −(IR) −1 λ δω . Note the last term in the (λ δζ ) expression involves calculating the partial derivative of the angular velocity error kinematic, ∂ ∂(δζ) [δ˙ω], which obviously leads to difficult math. This step is also required when performing coordinate mapping between the DCM into other attitude parameters.

Numerical Simulation
This section presents simulation results of a fixed final-time and final-state open-loop optimal control solution for the spacecraft tracking problem. Modified Rodrigues parameters are used for the attitude motion. The initial and final state variable conditions for this example are given in Table 4. The spacecraft moment of inertia tensor is given in Table 5. For simplicity, the weighting matrices are Q 1 = 0, Q 3 = 1, Q 2 = 0 3×3 , and Q 4 = R = I 3×3 . However, one can sweep those penalties to obtain different solutions sets. This example stands to solve the open-loop optimal control problem for the spacecraft tracking problem using arbitrary attitude representation. We consider the universal performance index given in Table 3. The state and co-state differential equations are solved in a Boundary Value Problem (BVP) framework using a shooting method (MATLAB fsolve) [29,33]. The optimal open-loop solution is shown in Figure 5. The time history of the optimal trajectory and control is shown in Figure 5a for MRPs representation. The trajectory is controlled to drive the spacecraft for a given initial state error, δω(0) and δζ(0), to rest at zero attitude error after 25 seconds. It is noted that the optimal open-loop solution obtained for various attitude representation produces the same angular displacement δφ and total cost, as given in Figure 5b.

Conclusions
Full analytical derivations of attitude error kinematics equations have been presented. Compact forms of attitude error kinematics equations are derived for various attitude parameterizations. The attitude error is defined as the rotation error between the true and estimated orientations. Several attitude error representations are developed for describing the orientation error kinematics. Two approaches to attitude error kinematics are introduced. The first one considers the estimated angular velocity defined in the true body axes, while in the second one, the estimated angular velocity is defined in the estimated body axes. These two angular velocity definitions are usually adopted in simulations and in real estimation/control applications, respectively. These two nonlinear kinematic models are valid for arbitrarily large relative rotations and rotation rates. These results are expected to be broadly useful for generalizing extended Kalman filtering formulations and optimal control tracking problems.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Angular Velocity Error Vector and Euler Angle Error Rates Mapping
This appendix presents forward and inverse mapping between the angular velocity error vector δω and the Euler angle error rates δθ = δθ 1 δθ 2 δθ 3 T in the following form: Note that the shorthand notations c i = cos(δθ i ) and s i = sin(δθ i ) are used in Table A1. Table A1. Angular velocity error vector and Euler angle error rates mapping.
i-j-k M −1 Figure A1. Summary of the attitude error kinematics equations. The terms in red indicate the differences between the two forms