Attitude propagation stability and accuracy

[D H]

I am verifying and validating a rotational state propagator used in a dynamic simulation package. I have found some problems and solutions to them. This post outlines the propagation, the problems, and the solutions.

Some questions before I start:

- Has anyone else analyzed stability and accuracy for rotational state propagation techniques?
- Is there a "gold-standard" for propagating the rotational state?1. The propagation technique
The propagator currently integrates attitude and attitude rate as a quaternion and angular velocity using
$$\begin{align*} {\cal Q} \; &\text{is the inertial-to-body left transformation unit quaterion} \\ \vec\omega \; &\text{is the angular velocity wrt inertial expressed in the body frame} \\ \vec\tau \; &\text{is the external torque expressed in the body frame} \\ \mathbf I \; &\text{is the body inertia tensor expressed in the body frame} \\ \frac{d\cal Q}{dt} &= \left[\begin{matrix} 0 \\ -1/2 \vec\omega\end{matrix}\right]{\cal Q} \\ \frac{d\vec\omega}{dt} &= \mathbf{I}^{-1}\left(\tau - \vec\omega \times \mathbf I \vec\omega\right) \\ {\cal Q} (t+\Delta t) &= {\cal Q}(t) + \frac{d{\cal Q}(t)}{dt}\Delta t \\ \vec\omega(t+\Delta t) &= \vec\omega(t) + \frac{d \vec\omega(t)}{dt} \Delta t \end{align*}$$

We use several integrators, ranging including Euler, the standard RK-4, Adams-Bashford-Moulton, and others. All use the above as the basis for propagation.

The quaternion propagation does not result in a unit quaternion. To overcome this problem, we normalize the quaternion after each integration step.

2. Accuracy and stability analysis
I have found some apparently significant problems with this technique.

- The quaternion step is mathematically invalid for all but the trivial case, omega=0. Unit quaternions multiply; they do not add. The normalization hack fixes the problem to some extent. However, with some very simple test cases, I can make the attitude propagation fail (error grows unbounded) or even blow up (floating point underflow exception). For a simple Euler integrator, problems can appear when $\omega \Delta t$ is as small as 1e-5.

- The angular velocity propagation can fail catastrophically (floating point overflow). Numerical precision problems and poor quaternion propagation can couple with the inertial torque term $\vec\omega \times \mathbf I \vec\omega$ to result in a torque-free system that does not conserve angular momentum, making the angular velocity undergo secular growth.

3. New approach
I am doing two things to overcome these problems:

1. Fixing the quaternion step algorithm.
The quaternion can be propagated analytically for a constant rotation rate vector via
$${\cal Q} (t+\Delta t) = \exp\left(\left[\begin{matrix} 0 \\ -1/2 \vec\omega\end{matrix}\right]\right) {\cal Q}(t)$$
I am using this as the basis for propagation. The quaternion exponential of a pure vector quaternion expands to transcendental functions. To avoid calls to sin and cos, I am using the second-order Padé approximant for cos(x) and the third-order Padé approximant for sin(x).

2. Propagating angular momentum in the inertial frame rather than the angular velocity in the body frame.
This appears to work well on my test cases with analytic solutions (torque-free spherical body, torque-free symmetric top, and simple spherical body torsional oscillator).

 

[Aaron Bex]

To propagate angular momentum, I am using

\frac{d\vec L}{dt} = \vec\tau
\vec L(t+\Delta t) = \vec L(t) + \frac{d\vec L(t)}{dt} \Delta t

The angular velocity is determined from the angular momentum and the body inertia tensor.

\vec\omega(t+\Delta t) = \mathbf I^{-1}\vec L(t+\Delta t)

This is a valid algorithm for any external torque and does not suffer from numerical precision problems or instability issues.

Conclusion
I have outlined an improved rotational state propagator and tested it against some simple test cases. I believe this propagator will provide more accurate and stable propagation of the rotational state than the current technique in use. I am working on incorporating this new approach into the dynamic simulation package.

 

[ravenprp]

To answer your questions, yes, there have been previous analyses on stability and accuracy for rotational state propagation techniques. One well-known approach is the Rodriguez formula, which uses a rotation matrix instead of a quaternion. As for a "gold-standard" for propagating the rotational state, it really depends on the specific application and requirements. However, the Rodriguez formula is often used as a benchmark for comparison.

Your new approach seems promising and addresses the issues you have identified with the previous technique. By using the analytical solution for quaternion propagation and propagating angular momentum in the inertial frame, you are likely to achieve better stability and accuracy. It would be helpful to conduct further testing and comparisons with other techniques to validate your new approach.

Overall, it is important to thoroughly analyze and validate any propagator used in a dynamic simulation package, as it can greatly impact the accuracy and reliability of the simulation results. Your efforts in identifying and addressing these issues are commendable and will contribute to the improvement of the overall simulation package.