Attitude quaternion derivatives from Euler angular velocities

[vicjun]

I'm struggling to understand what the derivative of an attitude quaternion really is and how to use it. I need it to solve a problem relating to a rotating frame of reference relative an inertial frame.

The information I have is a vector of Euler angular velocities (i.e for roll $\phi$, pitch $\theta$ and yaw $\psi$), with components $\dot{\phi}$, $\dot{\theta}$ and $\dot{\psi}$, for each point along the path. The points correspond to a date/time $t_i$, rather than a spatial location (where $i$ is the index). The angular velocities represent the rotational velocity of the frame relative to an inertial frame at a point with date $t_i$.

The only other piece of information I have is the inital attitude of the frame at time 0, i.e $t_0$, in the form of a quaternion $Q_0$ with components $q_0, q_1, q_2$ and $q_3$

What I want to do is to calculate the quaternion derivative at every point in time. To do this I am using the method desribed here:
http://www.euclideanspace.com/physics/kinematics/angularvelocity/QuaternionDifferentiation2.pdf

Which works for the first point, because I can create a quaternion of angular velocities and I have the initial quaternion $Q_0$. I have a ready-made script that allows me to calculate this quaternion derivative using this method. So far, so good. But what then?

Using this method, I can't calculate the derivative in the next point without first knowing the attitude quaternion in that point, say $Q_1$ (I have the angular velocities in each point, as mentioned above). I tried approximating this quaternion by adding the derivative of the previous point times the timestep, as follows (for each quaternion component):

$q_{0,i+1} = q_{0,i} + (t_{i+1} - t_i) \cdot \dot{q}_{0,i}$

This effectively gives me an attitude quaternion in every point, allowing me to calculate the quaternion derivative in that point (as I have the angular velocities). However, when comparing with reference values for the quaternion components, this method is not at all accurate enough, and does not give reliable results.

This has led me to wonder what the quaternion derivative really is. As I understand it (and as indicated by the document in the link) the quaternion derivative is the rotational velocity equivalent to the Euler angle attitude quaternion, i.e an indicator of how fast the frame is rotating relative an inertial frame at any given time. I have used quaternions in the past and I think I have a pretty good understanding of them, but quaternion differentiation is new to me.

Ideally, I would like to be able to calculate the quaternion derivative in every point based entirely on the angular velocities and initial attitude, without ever having to use attitude quaternions (as they are not reliable when calculated using the equation above). Is there a way to do this?

Intuitively, I think there should be a way. Euler angular velocities should be able to directly yield a corresponding velocity quaternion. Or have I completely misunderstood the problem?I often have trouble putting my questions into words, so I'd be glad to clarify the problem if needed. I tried to find the answer on Google, but most topics seem to relate to 3D programming, without enough "mathematical background" to help me understand the principles.Thanks.

 

[D H]

This is too complex a topic to delve into in an internet forum such as this.

As a starter, the time derivatives of your roll-pitch-yaw Euler angles are not angular velocity. They are related to angular velocity by a somewhat messy transformation, and the relation depends on the order of the rotation sequence (e.g., a roll pitch yaw sequence versus a yaw pitch roll sequence). Angular velocity as measured by a rate gyro is a vector (better said, a pseudovector). Those time derivatives do not form a vector.

Simply integrating angular velocity does not yield anything meaningful. Properly integrating the quaternion derivative does. However, that "properly" is full of caveats. Typical numerical integration schemes will quickly fail because the integrated quaternion stops being unitary. You need to do something special (and maybe hackish) to keep the integrated quaternion a unit quaternion.

There's some stuff on this on the 'net, but it's piecemeal. What you need is a book. A good starting point is Quaternions and Rotation Sequences: A Primer With Applications to Orbits, Aerospace, and Virtual Reality by Jack B. Kuipers.

 

[vicjun]

Thank you for your answer, it was very helpful.

I understand that Euler angle time derivatives are related to angular velocities by a matrix consisting of cosine and sine functions of the Euler angles (too much LaTeX work to write it here, but you probably know which one I'm referring to). In a similar fashion a quaternion derivative is related to "its" quaternion by a tensor consisting of angular velocities. The latter relationship is what I'm using to calculate the quaternion derivative in every point.

But is there no way to bypass using the attitude quaternion as well as the Euler angles entirely? That is, is it possible at all to calculate the quaternion derivative directly from the Euler angle time derivatives? Maybe through the principal axis?

Even if the Euler angle time derivatives aren't the same as the angular velocities $\omega_1, \omega_2$ and $\omega_3$ isn't there a more direct relationship to quaternion differentiation?

Thanks for the litterature tip, it looks like a very good book. It's a little too expensive and would take a bit too long to get here, though. I only have some course literature available to me, such as Analytical Mechanics of Aerospace Systems by Schaub and Junkins, and Spaceflight Dynamics by Wiesel. They both go over quaternions, but I haven't found anything so far that would help me. Is there a more advanced book than these I would find useful?

As for numerical integration methods, I managed to avoid inaccuracies a little by normalizing the quaternion after each step. This only works for shorter paths, though. With many points, the result drifts too much from the correct value. I'd much rather use analytical methods to solve this problem.

I think this is largely due to me misunderstanding the physical meaning of Euler angle time derivatives, and how they are different from rotational velocities. They are more of a "rate of change" I guess, and are not used to describe attitude changes at any given moment. As you wrote, they cannot be assembled into a vector.

 

[D H]

The Euler angle time derivatives don't have much physical meaning. Your attitude rate sensors aren't measuring those derivatives. They're measuring angular velocity as expressed in the rate sensor case frame. I'd drop the Euler angle derivatives ASAP. Euler angles are handy to give people an understanding of orientation, but of very little use beyond that. The derivatives of those Euler angles aren't good for much of anything at all.

The basic problem here is that orientation is a beast of a rather different color compared to translation. Translation in three dimensional is described by the mathematics of vector space ℝ³, rotation by the mathematics of the Lie group SO(3). The problem with using techniques such as RK4 to integrate orientation is that it is a hammer when what you need is a screwdriver. It's the wrong tool. It can be forced to fit by hackery such as renormalizing the quaternion, but that is still a forced fit.

For something that isn't a forced fit, you need to look to geometric integration techniques such as Lie group integrators. There's lots of recent literature on geometric integration and Lie group integrators, and some of it you will find online. I don't think you'll find any of this in books yet; it's too new.

 

[vicjun]

I don't mean to gravedig, sorry, but I have another question about quaternion derivatives that does not deserve its own thread, and it is within the topic of the discussion.

I have now figured out how to calculate the quaternion derivatives from the quaternions and the angular velocities. However, when I interpolate these quaternion derivatives, and integrate the resulting polynomials, I don't get values close to my original quaternions. Strangely, only the first component of the reintegrated quaternion q⁰ matches its original value, while q₁, q₂ and q₃ deviate with time. I thought it might be related to the fact that the three latter components are a function of the principal rotation vector, while q₀ is a function of the principal angle only.

Either way, it is difficult to check where the error is. But the quaternion derivatives are obviously wrong. I suspect the angular velocities are incorrect. These are calculated from the Euler angles and the Euler angle derivatives.

First, I want to check my quaternion derivative. As the corresponding attitude quaternion is unitary, i.e a norm of one, should its derivative also be unitary? To me, it seems intuitive that a quaternion derivative should be a versor, as both frames have the same origin. But I'm not sure. If the norm of the quaternion derivative changes with time, shouldn't the norm of its quaternion also be changing? I would appreciate a clarification.

Thanks in advance.

 

[D H]

First, I want to check my quaternion derivative. As the corresponding attitude quaternion is unitary, i.e a norm of one, should its derivative also be unitary?

No. There's no reason the quaternion derivative is unitary. For example, the derivative of a constant unit quaternion is zero.

What you can do to check your quaternion derivative is to verify that the product $q_s \dot q_s + \vec q_v\cdot\dot{\vec q}_v$ is zero. (Here, I am decomposing a quaternion q into a real scalar part $q_s$ and an imaginary part expressed a real 3-vector $\vec q_v$.) Another way to look at it is that if you treat the quaternion and its derivative as real 4-vectors, the inner product is zero. In other words, the quaternion derivative is either zero or normal to the quaternion (normal in a Euclidean 4-space sense).

 

[vicjun]

I have been trying to find evidence of this relationship in literature. Why should the inner product of the quaternion and its derivative be zero? In other words, why is the derivative normal to the quaternion? I do understand how this leads to the relationship though, since it is just the computation of the inner product between these quaternions. But I would like some theoretical background.

Does it have to do with the attitude quaternions being unitary, so that they represent a point on a unit 3-sphere? Is the derivative then necessarily normal?

 

[D H]

The inner product of any constant length vector in ℝⁿ and its time derivative is necessarily zero. In other words, the time derivative either is zero or is normal to that constant length vector. Just differentiate $\vec v\cdot \vec v = c^2$. Since the vector is of constant length, this derivative must be zero. Differentiating $\vec v\cdot \vec v$ yields $\dot{\vec v}\cdot \vec v + \vec v \cdot \dot{\vec v} = 2 \dot{\vec v}\cdot \vec v$.

The magnitude of a quaternion can be computed by treating the quaternion as a 4-vector and using the standard Euclidean norm for 4-vectors.

 

[vicjun]

You're absolutely right, sorry. I understand the proof through the dot product, I had forgotten the quaternions were unitary.

It seems my derivatives are correct at least in direction, as I get inner products which are either 0, or in the order of 10^-19 and 10^-20, which can be considered zero, since I'm doing numerical calculations.

This might mean my interpolation methods are wrong. I'm going to test these.

 

[vicjun]

I've tested the interpolation methods in two different ways, and they're unlikely to be the source of error.

It's obvious the $\dot{q_1}, \dot{q_2}$ and $\dot{q_3}$ components are wrong. If I look at how the quaternion changes with time by printing out the component values at every point, they do not match the magnitude of the derivative.

What I do is this:

I have three frames, an inertial frame IN, a local orbital frame LO, and the body-fixed frame BODY. I have the angular velocity of the local orbital frame around the inertial frame, $\Omega^{LO/IN}$, and I have the Euler angles as well as their derivatives of the body-fixed frame relative to the local orbital frame. I can use these to calculate the angular velocity of the body-fixed frame around the local orbital frame with the following relationship:

$\Omega^{BODY/LO} = \begin{bmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \\ \end{bmatrix} = \begin{bmatrix} cos \theta cos \psi & sin \psi & 0 \\ - cos \theta sin \psi & cos \psi & 0 \\ sin \theta & 0 & 1\\ \end{bmatrix} \cdot \begin{bmatrix} \dot{\varphi} \\ \dot{\theta} \\ \dot{\psi} \\ \end{bmatrix}$

for the (1,2,3) i.e (roll,pitch,yaw) set of Euler angles I'm using. Now, to get the angular velocity of the Body-fixed frame around the Inertial frame I simply add the velocities:

$\Omega^{BODY/IN} = \Omega^{BODY/LO}+\Omega^{LO/IN}$


I then convert the Euler angles to quaternions in every point to get the quaternion that describes the orientation of the Body-fixed frame relative the Local Orbital frame. Quaternion multiplication by the quaternion relating the Local Orbital frame to the Inertial frame gives the quaternion describing the orientation of the Body-fixed frame relative the Inertial frame:

$Q^{BODY/IN} = Q^{BODY/LO} \cdot Q^{LO/IN}$


I compute the quaternion derivative in every point using the relationship

$\dot{Q} = \begin{bmatrix} \dot{q_0} \\ \dot{q_1} \\ \dot{q_2} \\ \dot{q_3} \end{bmatrix} =\frac{1}{2}\begin{bmatrix} 0 & -\omega_1 & -\omega_2 & -\omega_3 \\ \omega_1 & 0 & \omega_3 & -\omega_2\\ \omega_2 & -\omega_3 & 0 & \omega_1\\ \omega_3 & \omega_2 & -\omega_1 & 0 \end{bmatrix} \cdot \begin{bmatrix} q_0 \\ q_1 \\ q_2 \\ q_3 \end{bmatrix}$


Finally, I interpolate each component independently over time to get a polynomial for every component. I then integrate these polynomials, and recalculate a quaternion component value at every point in time. These do not match the original quaternions I set out with. Oddly enough, only $\dot{q_0}$ is correct, and matches its starting value almost perfectly, accounting for the numerical error of interpolation. I have tested the interpolation independently, and it is doubtful the error comes from there.

It seems then that the last three rows of the skew symmetric $\Omega$ matrix, used to calculate the derivatives above, are wrong. But I am pretty sure I used the correct formula. Is it a correct method to calculate a quaternion from a quaternion derivative?

 

[D H]

I compute the quaternion derivative in every point using the relationship

$\dot{Q} = \begin{bmatrix} \dot{q_0} \\ \dot{q_1} \\ \dot{q_2} \\ \dot{q_3} \end{bmatrix} =\frac{1}{2}\begin{bmatrix} 0 & -\omega_1 & -\omega_2 & -\omega_3 \\ \omega_1 & 0 & \omega_3 & -\omega_2\\ \omega_2 & -\omega_3 & 0 & \omega_1\\ \omega_3 & \omega_2 & -\omega_1 & 0 \end{bmatrix} \cdot \begin{bmatrix} q_0 \\ q_1 \\ q_2 \\ q_3 \end{bmatrix}$


Finally, I interpolate each component independently over time to get a polynomial for every component. I then integrate these polynomials, and recalculate a quaternion component value at every point in time. These do not match the original quaternions I set out with. Oddly enough, only $\dot{q_0}$ is correct, and matches its starting value almost perfectly, accounting for the numerical error of interpolation. I have tested the interpolation independently, and it is doubtful the error comes from there.

It seems then that the last three rows of the skew symmetric $\Omega$ matrix, used to calculate the derivatives above, are wrong. But I am pretty sure I used the correct formula. Is it a correct method to calculate a quaternion from a quaternion derivative?

That matrix is wrong. You need to either change the signs of the first row and first column or change the signs of the lower right 3x3 submatrix. And even that won't guarantee that you are doing things consistently. Rather than using a matrix, I suggest using a pure imaginary quaternion with the imaginary part being the angular velocity vector. Then use quaternion multiplication to multiply this pure imaginary quaternion and your quaternion to get the quaternion derivative. Finally, scale by 1/2 or -1/2.

There's a problem here: which comes first, the angular velocity or the quaternion, and is the product scaled by +1/2 or -1/2? The answer depends on how you are using quaternions. There are multiple sources of ambiguity with regard to using quaternions to represent rotations in ℝ³. Left versus right quaternions, transformation versus rotation quaternions, and does the quaternion the rotation/transformation from frame A to frame B or from frame B to frame A?

 

[vicjun]

I calculated the quaternion product

$\dot{Q}=\frac{1}{2} Q \cdot \Omega$

where I assumed that $\Omega$ is a pure imaginary quaternion, so that

$\Omega = \begin{bmatrix} 0 \\ \omega_1 \\ \omega_2 \\ \omega_3 \\ \end{bmatrix}$

This does indeed give the matrix you suggested, i.e

\begin{bmatrix}
0 & -\omega_1 & -\omega_2 & -\omega_3 \\
\omega_1 & 0 & -\omega_3 & \omega_2\\
\omega_2 & \omega_3 & 0 & -\omega_1\\
\omega_3 & -\omega_2 & \omega_1 & 0
\end{bmatrix}

for finding the quaternion derivative. However, it seems this doesn't give the correct answer either... But there aren't that many options. If I instead try $\dot{Q}=\frac{1}{2} \Omega \cdot Q$ I don't get the correct quaternion derivative. These are the only two products I can think of. Unless I can use the conjugate of any quaternion. But that would roughly be the equivalent of changing the sign of the equation (at least for $\Omega$), and $\dot{q_0}$ would now be wrong.

I tried to think about what equation made sense, instead of just trying out different combinations. If I'm not mistaken, the angular velocity vector/quaternion I'm using here is the angular velocity of the Body-fixed frame around the inertial frame. I multiply this by the quaternion representing the orientation of the Body-fixed frame relative the inertial frame. But what does the quaternion derivative represent? The change in time of the orientation of the Body-fixed frame relative to the Inertial frame? If so, then I could write

$\dot{Q}^{BODY/IN}=\frac{1}{2} Q^{BODY/IN} \cdot \Omega^{BODY/IN}$


$Q^{BODY/IN}$ is clearly a versor, since it is a unitary quaternion. $\Omega^{BODY/IN}$ is not, however, since its norm is not one. Neither is the derivative. Does this then mean that the quaternion derivative both rotates and transforms? Does this really make sense if the frame only rotates?

I'm also a little confused by the different derivations in literature. I got the previous matrix from Jack Kuiper's book (the one you recommended) chapter 11.5, "Derivative of the quaternion". The same matrix is presented in Schaub & Junkin's book "Analytical Mechanics of Aerospace Systems", chapter 3.4.

As understand it, if a quaternion represents the orientation of frame A relative to frame B it also transforms frame B into frame A. Right?

 

[D H]

However, it seems this doesn't give the correct answer either... But there aren't that many options. If I instead try $\dot{Q}=\frac{1}{2} \Omega \cdot Q$ I don't get the correct quaternion derivative. These are the only two products I can think of.

You can also negate. There are four possibilities, pre-multiply vs post-multiply, and multiply by +1/2 or -1/2. Which one is right depends on the frame in which you are representing angular velocity, the sense of the angular velocity vector, whether you are using left or right quaternions, whether the quaternion represents rotation or transformation, and the direction sense of the quaternion. I don't know which you are using, so I don't know which of the four alternatives you should be using.

I use left transformation quaternions and represent angular velocity as the angular velocity of the body with respect to inertial but represented in body coordinates. With this, the time derivative of the inertial frame to body frame left transformation quaternion is
$$\dot{\mathcal Q}_{I \to B} = - \frac 1 2 \begin{bmatrix} 0 \\ \vec{\omega} \end{bmatrix} {\mathcal Q}_{I \to B}$$
For right transformation quaternions, (e.g., Kuiper), simply conjugate:
$$\dot{\mathcal Q}_{I \to B} = \frac 1 2 {\mathcal Q}_{I \to B} \begin{bmatrix} 0 \\ \vec{\omega} \end{bmatrix}$$

As understand it, if a quaternion represents the orientation of frame A relative to frame B it also transforms frame B into frame A. Right?

Correct.

 

[vicjun]

OK, so I got an idea on how to test the handedness of the quaternions.

First, a correction. The quaternion describing the orientation of the Body-fixed frame relative the Inertial frame is given by the quaternion multiplication

$\mathcal Q_{I \to B} = \mathcal Q_{I \to O} \cdot \mathcal Q_{O \to B}$ (1)


and NOT the other way around, as I wrote earlier, it was a typo. I've used the above equation in my calculations. Also, I hope you don't mind me using your notation, it's a bit handier (O is the Local Orbital frame).

Now, what I did was calculate the quaternion derivatives separately. In other words, I calculated $\dot{\mathcal Q}_{I \to O}$ and $\dot{\mathcal Q}_{O \to B}$ separately, using $\Omega^{O/I}$ and $\Omega^{B/O}$, respectively (where these are pure imaginary quaternions as before). Working with two frames at a time rather than three is always easier. I tested the multiplication order of each step, and I noticed that it's different. The transition from Inertial to Orbital frames is given by


$\dot{\mathcal Q}_{I \to O} = \frac{1}{2} \Omega^{O/I} \cdot \mathcal Q_{I \to O}$

while the transition from Orbital to Body-fixed frames is given by

$\dot{\mathcal Q}_{O \to B} = \frac{1}{2} \mathcal Q_{O \to B} \cdot \Omega^{B/O}$

So apparently, one quaternion is left-handed, while the other is right-handed, as you wrote. Unfortunately, in this case I'm given many of these values from functions that I have little insight in. But I'll try to draw up the frames and see how they're related.

And the fact that both equations have +1/2 suggest that the angular velocities are expressed in the target frame, right? For instance: the angular velocity of the body with respect to the orbital frame is expressed in the body-fixed frame.

So how do I calculate the "total" derivative, i.e the time derivative of the body-fixed frame with respect to the inertial frame? Equation (1) above would suggest

$\dot{\mathcal Q}_{I \to B} = \dot{\mathcal Q}_{I \to O} \cdot \dot{\mathcal Q}_{O \to B}$

However, this turns out to not be the case. I tried to change the order of multiplication, but this does not give correct results either. How should I reason here? Shouldn't it be the same as (1)? I'll check my code again, I wrote this last calculation a bit hastily, and the results are weird.

 

[D H]

Also, I hope you don't mind me using your notation, it's a bit handier (O is the Local Orbital frame).

No problem. Have at it.

The transition from Inertial to Orbital frames is given by
$$\dot{\mathcal Q}_{I \to O} = \frac{1}{2} \Omega^{O/I} \cdot \mathcal Q_{I \to O}$$
while the transition from Orbital to Body-fixed frames is given by
$$\dot{\mathcal Q}_{O \to B} = \frac{1}{2} \mathcal Q_{O \to B} \cdot \Omega^{B/O}$$
So apparently, one quaternion is left-handed, while the other is right-handed, as you wrote. Unfortunately, in this case I'm given many of these values from functions that I have little insight in.

Whoa! It's not left-handed versus right-handed. It's just left versus right. A left transformation quaternion transform vectors from frame A to frame B via
$$\begin{bmatrix} 0 \\ \vec x_B \end{bmatrix} = {\mathcal Q}_{A \to B} \begin{bmatrix} 0 \\ \vec x_A \end{bmatrix} {\mathcal Q}_{A \to B}^{\ast}$$
The quaternion is unconjugated on the left, conjugated on the right. Hence the name left quaternion. With right transformation quaternions it's the other way around:
$$\begin{bmatrix} 0 \\ \vec x_B \end{bmatrix} = {\mathcal Q}_{A \to B}^{\ast} \begin{bmatrix} 0 \\ \vec x_A \end{bmatrix} {\mathcal Q}_{A \to B}$$

Which is right? Both are. There is no right or wrong here. Just right or left.

I suspect the problem is with your angular velocity vectors. Two issues that can trip you up: The angular velocity of what frame with respect to what other frame, and in which frame's coordinate system is the vector represented? The choice is pretty standard for position and linear velocity. Position is almost invariably the vector from the origin of the inertial frame to the origin of the rotating frame, represented in inertial frame coordinates. Or more generally, from the origin of the parent frame to the origin of the child frame, represented in child frame coordinates. Linear velocity is almost always the element-by-element time derivative of that position vector. In other words, it's the linear velocity of the origin of the rotating frame as observed by an inertial observer, with coordinates expressed in terms of the inertial frame.

Angular velocity is not so clear. Typically it is the angular velocity of the rotating frame with respect to the inertial frame, with coordinates expressed in terms of the rotating frame. It's a bit convoluted, but this convention goes back to Euler. You have to be very, very careful with angular velocities, particularly so when there are more than two frames involved.

 

[vicjun]

I understand the difference now. I agree, it is probably about finding out which frame the angular velocities are expressed in.

But if I manage to calculate the quaternion derivative from one frame to another, as I did, i.e $\dot{\mathcal Q}_{I \to O}$ as well as $\dot{\mathcal Q}_{O \to B}$ separately, shouldn't there be a way to find $\dot{\mathcal Q}_{I \to B}$ from these? Or should I instead try to conjugate one or both of my angular velocity vectors?

 

[D H]

I understand the difference now. I agree, it is probably about finding out which frame the angular velocities are expressed in.

But if I manage to calculate the quaternion derivative from one frame to another, as I did, i.e $\dot{\mathcal Q}_{I \to O}$ as well as $\dot{\mathcal Q}_{O \to B}$ separately, shouldn't there be a way to find $\dot{\mathcal Q}_{I \to B}$ from these? Or should I instead try to conjugate one or both of my angular velocity vectors?

Sure. Just use the chain rule. With right quaternions, ${\mathcal Q}_{I \to B} = {\mathcal Q}_{I \to O} {\mathcal Q}_{O \to B}$, so $\dot{\mathcal Q}_{I \to B} =\dot{\mathcal Q}_{I \to O} {\mathcal Q}_{O \to B} + {\mathcal Q}_{I \to O} \dot{\mathcal Q}_{O \to B}$

 

[AndrewKw]

Sorry for interference, but I need to understand almost the same thing. The difference is I need to calculate angular velocity from quaternion orientation difference and time.

I am using Bullet Physic library to program some function, where I have difference between orientation from gyroscope given in quaternion and orientation of my object, and time between each frame in milisecond. All I want is set the orientation from my gyroscope to orientation of my object in 3D space. But all I can do is set angular velocity to my object. I have orientation difference and time, that is why I need vector of angular velocity [Wx,Wy,Wz] from that.

And I found here:
http://www.euclideanspace.com/physics/kinematics/angularvelocity/QuaternionDifferentiation2.pdf
and here:
http://www.euclideanspace.com/physics/kinematics/angularvelocity/

that formula:

$(W_{x}(t),W_{y}(t),W_{z}(t)) = 2\frac{dq(t)}{dt}q^{-1}(t)$

But, probably I misunderstand something, because I did something like:

btQuaternion diffQuater = gyroQuater - boxQuater;
btQuaternion diffConjQuater;

diffConjQuater.setX(-(diffQuater.x()));
diffConjQuater.setY(-(diffQuater.y()));
diffConjQuater.setZ(-(diffQuater.z()));

////////////////
//W(t) = 2 * dq(t) /dt * conj(q(t))

btQuaternion velQuater;

velQuater = ((diffQuater * 2) / d_time) * diffConjQuater;

But it is not working as I expect, I mean, there is written, vector part of the result quaternion should be vector of angular velocity, and scalar part should be 0, but my result is not like that.

Can someone help me with that?

Edit:
I know where was a problem, I should do:

velQuater = ((diffQuater * 2) / d_time) * conjBoxQuater;

But I get another problems, when my difference between quaternions is not so small (I mean when angular velocity should be high), and when I get 120 degrees I have some weird results. Can someone help me with that?


I know where was a problem, I should do:
velQuater = ((diffQuater * 2) / d_time) * conjBoxQuater;

 

[D H]

You just can't take two random packages that use quaternions to represent orientation and expect them to work together without some glue. There are 32 combinations of representational choices, 30 of which represent ways things won't work together without glue.

1. Do the quaternions represent transformations or rotations?
2. Do quaternions transform (rotate) via QvQ* or Q*vQ ?
3. Is the scalar part of the quaternion first or last?
4. Is angular velocity from the inertial frame to the rotating frame, or the reverse?
5. Is angular velocity represented in the inertial frame or in the rotating frame?

That's 2⁵=32 opportunities for Murphy's law to rear it's ugly head. You need to know how each package works. The problem is confounded by the fact that a lot of authors don't even know these distinctions exist.

You've added yet another issue to the mix in using (Q(t+Δt)-Q(t))/Δt to estimate a quaternion derivative. This implicitly makes the small angle approximation due to the nature of the quaternion exponential. This will create problems eventually.

 

[vicjun]

I need to revive this thread, yet again, for another related question. Sorry. At least I won't have to re-explain the problem and all the parameters.

Thanks for the help with the quaternion derivatives, I understand them much better now.

I now wish to integrate the angular velocity $\boldsymbol{\omega}_n$ at every point n, where my time variable is $t_n$. This would yield the angular displacement $\boldsymbol{\theta}_n$. Now, since I have a vector with known components in every timstep, I can integrate component-wise:

$\theta_{i,n}=\omega_{i,n} \cdot t_n + C$

where C is the integration constant and i=1,2,3. But how do I find the value of this constant? Should it be zero? This would imply that at t=0, there is no angular displacement. Right? So should I then set the constant to whichever initial angular displacement I have? This seems like such a simple problem and a resulting stupid question, but I'm just confused by the physical implications.


Thanks in advance.

 

[D H]

I now wish to integrate the angular velocity $\boldsymbol{\omega}_n$ at every point n, where my time variable is $t_n$. This would yield the angular displacement $\boldsymbol{\theta}_n$.

This is not possible. While this works in SO(2) (rotation in 2D space, where angular velocity is a scalar), it does not work in any higher dimension. It has do with the noncommutativity of SO(3) and higher.

 

[vicjun]

Thanks for the reply, I should have been clearer about the context. I found a method to numerically integrate quaternions, in a book. It basically involves integrating the kinematic equation

$\dot{Q}=\frac{1}{2} Q \, \Omega$

where $\Omega$ is the pure quaternion representation of the angular velocity vector. In the method, it is assumed that it commutates with its integral, i.e.

$\Omega \cdot \int_{t_n}^{t_{n+1}}\Omega(\tau)d\tau=(\int_{t_n}^{t_{n+1}}\Omega(\tau)d\tau) \cdot \Omega$

which allows for analytical integration of the kinematic equation using

$Q(t)=Q \cdot exp(\frac{1}{2} \int_{t_n}^{t_{n+1}}\Omega(\tau)d\tau)$

and if we write

$\theta_{i,n}=\int_{t_n}^{t_{n+1}}\omega_i(\tau) d\tau$

then we get

$Q(t)=Q \cdot exp(\frac{1}{2} \theta_{i,n})$

and this can then be expanded with Euler notation to give an expression of the type

$Q_{n+1}=Q_n \cdot expression(\theta_{i,n}\, \text{as a quaternion and its norm})$

However, this procedure implies a commutation error. I quote from the book:

However the commutation relation only holds in the very specific case of rotation around a fixed axis. Once the rotation axis is mobile, applying the method leads to an error known as a commutation error

The method is then built upon by compensating for this error. It is assumed that both the angular velocity quaternion and the attitude quaternion allow for Taylor expansions at any point $t_n$. This effectively does just that, and so the method can be used for mobile axes.

So is it possible to integrate the angular velocity vector, if the commutation error is compensated for? I'm still trying to understand this method, but it is from a reliable source.

 

[D H]

So is it possible to integrate the angular velocity vector, if the commutation error is compensated for? I'm still trying to understand this method, but it is from a reliable source.

Yeah, but all you are doing is the exact same thing described in that text.

There are a number of ways to represent rotation in 3D space. Quaternions are but one of many. Transformation or rotation matrices are another. Euler angles, yet another. Rodrigues parameter and modified Rodrigues parameters, yet another. There are a number of people who are moving toward MRPs to represent rotation. Yet another is what is variously called an Euler rotation, an eigenrotation, or a single axis rotation. Per Euler, any rotation in three space can be represented as a single rotation by some angle about some axis.

Euler angles, Rodrigues parameters, MRPs, and eigenrotations are represented by three parameters. That's a good thing, particularly when it comes to Kalman filters. I hope that you aren't talking about Euler angles. Euler angles are a mess. They're nice for human representation, sometimes nice for commands to robots, but they are extremely nasty when it comes to propagating state. I doubt you are talking about Rodrigues parameters or MRPs (but you may want to look into them). So that leaves eigenrotations.

Eigenrotations are just quaternions in disguise, but without the solid mathematical basis that quaternions describe. Try to integrate angular velocity to yield the eigenrotation and you'll be doing the equivalent of integrating a quaternion, using the exact same semi-valid tricks used to integrate the quaternion. You'll be doing quaternion integration but without the benefits that come with representation as a quaternion.

 

[vicjun]

Yeah, but all you are doing is the exact same thing described in that text.

Do you mean what I did previously, i.e. integrating the quaternion derivatives?

I doubt you are talking about Rodrigues parameters or MRPs (but you may want to look into them). So that leaves eigenrotations.

In this case I am constrained to quaternions as a final result (interpolation polynomials), starting from Euler angles and Euler angle rates. I know about Rodrigues parameters, though I've never used them in practice. One of the advantages as I see it, is that there are no redundant parameters, while they are more suited for numerical calculations than Euler angles, and are better for describing large and small rotations. But for now, I have to use quaternions.

Eigenrotations are just quaternions in disguise, but without the solid mathematical basis that quaternions describe. Try to integrate angular velocity to yield the eigenrotation and you'll be doing the equivalent of integrating a quaternion, using the exact same semi-valid tricks used to integrate the quaternion. You'll be doing quaternion integration but without the benefits that come with representation as a quaternion.

The method I tried earlier with quaternion derivatives works quite well. The purpose of my work is to find different methods to reach the same results, a polynomial describing each quaternion component. I then compare these methods to reference results to see which one is best. I guess integrating the angular velocity is very similar to integrating quaternion derivatives (both methods are based on the kinematic equation), but maybe there's a difference in accuracy? I'm not sure pursuing this will give interesting results, but it might be worth a try, if I get it to work. The challenge is finding the correct eigenrotation (thanks for explaining what it is!) from the angular velocity.

 

[RG123]

If you use MATLAB to do this you can utilize built in functions to transform between quaternion and euler angles and from there transform to the body axis to determine angular rates. It's not easy but it is totally possible. However, you do really need to know what you're doing.