- #1
- 15,524
- 769
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
<br /> \begin{align*}<br /> {\cal Q} \;<br /> &\text{is the inertial-to-body left transformation unit quaterion} \\<br /> \vec\omega \;<br /> &\text{is the angular velocity wrt inertial expressed in the body frame} \\<br /> \vec\tau \;<br /> &\text{is the external torque expressed in the body frame} \\<br /> \mathbf I \;<br /> &\text{is the body inertia tensor expressed in the body frame} \\<br /> \frac{d\cal Q}{dt} &=<br /> \left[\begin{matrix} 0 \\ -1/2 \vec\omega\end{matrix}\right]{\cal Q} \\<br /> \frac{d\vec\omega}{dt} &=<br /> \mathbf{I}^{-1}\left(\tau - \vec\omega \times \mathbf I \vec\omega\right) \\<br /> {\cal Q} (t+\Delta t) &= {\cal Q}(t) + \frac{d{\cal Q}(t)}{dt}\Delta t \\<br /> \vec\omega(t+\Delta t) &=<br /> \vec\omega(t) + \frac{d \vec\omega(t)}{dt} \Delta t<br /> \end{align*}<br />
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
<br /> {\cal Q} (t+\Delta t) = \exp\left(\left[\begin{matrix} 0 \\ -1/2 \vec\omega\end{matrix}\right]\right) {\cal Q}(t)<br />
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).
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
<br /> \begin{align*}<br /> {\cal Q} \;<br /> &\text{is the inertial-to-body left transformation unit quaterion} \\<br /> \vec\omega \;<br /> &\text{is the angular velocity wrt inertial expressed in the body frame} \\<br /> \vec\tau \;<br /> &\text{is the external torque expressed in the body frame} \\<br /> \mathbf I \;<br /> &\text{is the body inertia tensor expressed in the body frame} \\<br /> \frac{d\cal Q}{dt} &=<br /> \left[\begin{matrix} 0 \\ -1/2 \vec\omega\end{matrix}\right]{\cal Q} \\<br /> \frac{d\vec\omega}{dt} &=<br /> \mathbf{I}^{-1}\left(\tau - \vec\omega \times \mathbf I \vec\omega\right) \\<br /> {\cal Q} (t+\Delta t) &= {\cal Q}(t) + \frac{d{\cal Q}(t)}{dt}\Delta t \\<br /> \vec\omega(t+\Delta t) &=<br /> \vec\omega(t) + \frac{d \vec\omega(t)}{dt} \Delta t<br /> \end{align*}<br />
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
<br /> {\cal Q} (t+\Delta t) = \exp\left(\left[\begin{matrix} 0 \\ -1/2 \vec\omega\end{matrix}\right]\right) {\cal Q}(t)<br />
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).
Last edited: