# Why Runge-Kutta for Quaternion integration in Kalman filter?

I'm reading up on Kalman filtering at the moment. In particular, I'm interested in using the "extended" and "unscented" variants for IMU sensor fusion and calibration.

In A comparison of unscented and extended Kalman filtering for estimating quaternion motion, quaternions are used to represent 3d rotation.

I understand unit quaternions can be used to represent a 3d rotation. They suit for representing absolute attitude (a rotation from a universal reference), relative rotation, or angular velocity (a rotation representing rate per second, or some other fixed time period).

However, this papers discuss using Runge-Kutta integration, specifically RK4. It uses RK4 with the quaternions but doesn't seem to furnish details of what this involves or why it is necessary. Here's the part of the paper that mentions it…

Given the state vector at step k − 1, we first perform the prediction step by finding the a priori state estimate xˆ−k by integrating equation 1 [f = dq/dt = qω/2] through time by ∆t (i.e., 1.0 divided by the current sampling rate) using a 4th Order Runge-Kutta scheme.

I've encountered Runge Kutta before for integrating positions in kinematics. I don't really understand how or why it would be needed here.

My naive approach would be to simply multiply the existing attitude q by the angular velocity ω to get the expected new q – I don't see why numerical integration is necessary here? Perhaps it's to "scale" the unit time ω to the change that occurs in in ∆t, but surely that can be done very simply by directly manipulating of ω (raising it to fractional power ∆t)?

Anyone know?

• The paper is poorly written from a KF / EKF viewpoint. It doesn't specify what the signal model is, though it alludes to it when referring to equation 1. Also, equation 1 itself isn't clear to me: $q$ is usually a 4-vector, and $\omega$ is usually a 3-vector... so the term $q\omega$ is dimensionally unsound. – Peter K. Jan 18 '16 at 13:34
• Hi @PeterK. It mentions ω "is a pure vector quaternion representing angular velocity" – I presume this means that it has no real part, or that the real part is 0? – Benjohn Jan 18 '16 at 14:26
• Yes, you may be correct... though I think it's the $w$ component that's zero... "real" when it comes to quaternions is... hard to fathom. :-) – Peter K. Jan 18 '16 at 14:35
• :-) I've not yet worked out how an angular velocity quaternion is derived from the tri-axis rate gyros – I'm wondering if this could need some kind of integration step? Grasping-at-straws, though! – Benjohn Jan 18 '16 at 15:42

The main difference between normal vector algebra and geometric algebra is that we can multiply vectors. So of $e_x$, $e_y$, and $e_z$ are our (orthonormal) basis vectors, we also have $e_x e_y$ which is not a vector but a bivector which can be thought of as an oriented plane element. (By "pure vector quaternion" they mean a bivector; you are correct that that means $0$ $w$ part and it's reasonable to say it has a $0$ "real" part.) A key property is that the basis vectors anti-commute with each other, i.e. $e_x e_y = -e_y e_x$. This leads to $(e_x e_y)^2 = -1$. In general, any unit bivector squares to $-1$, which means we can apply Euler's formula: $$R = e^{\theta B} = \cos\theta + B\sin\theta$$ where $B$ is a unit bivector and we call $R$ a rotor. A unit quaternion is just a 3D rotor. (Complex numbers are just the even sub-algebra of the 2D geometric algebra. Here we are looking at quaternions, the even sub-algebra of 3D GA.)
If we want a rotation by $\theta$ in the plane of a unit bivector $B$ we use the rotor $R = e^{\frac{\theta}{2}B}$. (The half comes in because we rotate a vector $v$ via $Rv\widetilde{R}$; see the link above.) Of course, we want to allow both the plane and the angle to vary so we define the (non-unit) bivector $\Theta(t) = \theta(t)B(t)$ so if $$R(t) = e^{\frac{1}{2}\Theta(t)}$$ then $$\dot{R}(t) = \frac{1}{2}\Omega(t)R(t)$$ where $\Omega = \dot{\Theta}$. This is their $f = q\omega/2$; their $q$ is our $R$ and their $\omega$ is our $\Omega$. Now if we have $R(t_0)$ and we want $R(t_0 + \Delta{}t)$ we need to integrate $\dot{R}$ from $t_0$ to $t_0 + \Delta{}t$.
The scheme you describe is roughly analogous to doing forward Euler integration, which essentially assumes their $q$ and $\omega$ our $R$ and $\Omega$ is constant over $\Delta{}t$. RK4 is just a better integration method. I doubt there is any special reason they chose to use RK4 as opposed to other integration methods. It's just the default choice typically.
• My guess is that they keep $\omega$ constant while doing the integration. Since it doesn't sound like they get any angular velocity measurements the angular velocity only changes when it is "corrected" by the Kalman filter, i.e. the angular velocity is unobserved state that the Kalman filter is trying to reconstruct. – Derek Elkins left SE Jan 19 '16 at 14:20