How to use the quaternion derivative

Question

I'm having troubles using the Quaternion derivative.

So I have:

q(t) = ... <-- my current attitude

The derivative (using w(t) as body rotation rate) is:

dq(t)/dt = 1/2 * q(t) * w(t)

So how do I use that? When doing stepwise integration (like in a time stepped simulation), I would have expected:

q(t+∆t) = q(t) + 1/2*q(t)*w(t)*∆t

(Assuming, of course, that the step size is small enough to warrant the linear approximation).

However, when I do that, the resultant vector (q(t+∆t)) no longer has a norm2 of 1. If I normalize it, it doesn't appear to do anything.

If I instead create a quaternion rotation with the angles in w(t)*∆t and apply that:

q(t+∆t) = q(t) * omega2quaternion(w(t)*∆t)

Then it works.

So, how do I use the quaternion derivative?

Answer 1

(assume "quaternion" implies unit-magnitude quaternion)

The thing to understand is that quaternions are not closed under elementwise-addition like vectors are. The composition operation for quaternions is the "quaternion multiplication" which you have been denoting as *.

For example, if you wanted to first rotate by q1 and then rotate by q2, the composite rotation would be

q_total = compose(q1, q2) = q2*q1

Note that order matters, as * is not commutative. You can read up more about why this is on Wikipedia and such, but let me explain something else that might help you understand what's going on with your integration scheme.

When we let q define the orientation of our body relative to the world, what we really mean is that q represents a change of basis from body-fixed coordinates to world-fixed coordinates. At time t we might say,

q = body(t)_to_world

noting that the world frame does not vary with time by definition.

The key is to realize that a differential dq on the quaternion manifold is really just the orientation of our body at time t+dt relative to our body at time t. That is,

dq = body(t+dt)_to_body(t)

It then makes sense that

body(t+dt)_to_world = compose(body(t+dt)_to_body(t), body(t)_to_world)

or

q(t+dt) = q(t)*dq

which is the integration scheme you have found to work.

To compute dq, the most accurate way is to find the quaternion associated with the small rotation w*dt. I call this a small rotation because the magnitude of w*dt is the small amount of angular change during dt and the direction of w is the axis about which that rotation takes place. We can compute this as such.

One way to derive this dq computation is to first derive the quaternion derivative, which can either be done rigorously with complex analysis or simply with a bit of hand-waving as is done here (the hand-waving is how they apply without derivation the "quaternion version of Euler’s formula" but it's easy to believe if you already understand basic complex numbers). The result is,

dq/dt = (1/2)*W*q

where

W = 0 + wx*i + wy*j + wz*k.

If we applied first-order integration at this point,

q(t+dt) = q(t) + (dq/dt)*dt

the result would be very inaccurate for any reasonable dt; quaternions are not closed under addition unless the added value is truly differential (this is the nature of being a manifold).

Instead, if we solve this differential equation from t=t0 to t=t0+dt over which we assume W is constant, we have,

q(t) = q(t0)*exp((1/2)*W*(t-t0))

(accurate while W is constant from t0 to t0+dt).

There is that exponential form of a quaternion again, which we will continue to assume makes sense based on Euler identity intuition. Thus at time t0+dt,

q(t0+dt) = q(t0)*exp((1/2)*W*dt)

which implies

dq = exp((1/2)*W*dt).

Performing the above operation is exactly the same as expressing w*dt as a quaternion through the axis-angle conversion.

Check out some example code here for doing a bunch of stuff with these concepts.