Integrating body angular velocity

I've been reading over some very comprehensive notes on attitude representation, which were compiled by James Diebel, a Stanford student: http://www.swarthmore.edu/NatSci/mzucker1/e27/diebel2006attitude.pdf

What is of particular interest to me is equation $266$, which states that the rotation vector representation of an attitude is the integral of the body angular velocities over the time frame of interest (assuming the body and inertial frames start out coincident)

I see no proof of this anywhere in the paper, can someone help me understand how this is possible?

Edit

To clarify, this is my issue. Say I want to represent that attitude of the difference between two coordinate systems (say inertial and body) using an angle/axis vector that rotates a vector in the inertial frame to one in the body frame: $$v_{bi}(t)$$ I have measurements of body angular rate (from, lets say a gyroscope) $$\omega_b(t)$$ I'm curious if the following is generally true: $$\dot{v}_{bi}(t)=\omega_b(t)$$

Equation $266$ suggests that it is, equation $265$ seems to suggest otherwise.

Solution 1:

TL;DR: Ignore what you read in that paper; use either Eq. (2) or (3) in this post if you want to "integrate" angular velocity.

There is a solution to your problem, in the sense that it is possible to "integrate" the angular-velocity vector $\omega$ to obtain the rotation taking the inertial "world" frame onto the rotating "body" frame. In fact, I've written a paper titled "The integration of angular velocity". It turns out that it's pretty simple to get the frame from the angular velocity, but it's not quite as simple as just integrating $\omega$ as a vector.

First, to clear up the confusion about that paper: Diebel uses some unfortunate notation. The equation (266) you point to says \begin{equation} \mathbf{v}_{\omega'} := \int_{t_0}^{t_1} \boldsymbol{\omega}' dt \tag{1} \end{equation} There are two important things to note here. First, $\boldsymbol{\omega}$ is primed, meaning (as you seem to understand) that this is measured with respect to the body. That is, the basis with respect to which it is defined is different at each instant. So this is not what we normally mean when we're talking about the angular-velocity vector. Second, the quantity on the left is most definitely not the rotation vector required to go from the inertial frame to the body frame. Instead, it is some kind of rotation-rate vector which, if you multiply it by some infinitesimal time, would generate that rotation that takes the body frame at this instant into the body frame at that infinitesimal time in the future. Basically, Diebel can define any vector he wants—and he did so in Eq. (1). But it's only any use if $t_1$ and $t_0$ are infinitesimally close. Of course, this is all he uses it for in the following paragraph, so it's okay for him. But in general, I suggest you ignore that equation entirely.

The reason that Eq. (1) is not useful in general is because, at different instants, $\boldsymbol{\omega}'$ is defined with respect to different frames. So it doesn't make a whole lot of sense to add vectors from different instants together, which is just what that integral tries to do, unless the rotation is about a fixed axis. However, it is possible to rotate that body-fixed vector into the inertial frame. (You also probably don't want to write it as an integral equation; it's better as a standard differential equation.) The best way to integrate angular velocity is to use quaternions. Of course, quaternions are basically always the right way to do anything with rotations anyway. I don't understand why Diebel insists on converting quaternions into matrices, so I'll just get away from him, and write things like any normal person would.

Let's suppose $R(t)$ is a unit quaternion function of time that takes your inertial frame onto the body frame. For example: \begin{equation} \mathbf{x}'(t) = R(t)\, \mathbf{x}\, \bar{R}(t), \end{equation} and similarly for $\mathbf{y}$ and $\mathbf{z}$. Using the usual definition of angular velocity as that thing which, when crossed with a vector in the body frame, gives the rate of change of that vector. In particular, we have: \begin{equation} \frac{d \mathbf{x}'} {dt}(t) = \boldsymbol{\omega}(t) \times \mathbf{x}', \end{equation} and similarly for $\mathbf{y}$ and $\mathbf{z}$. Now, we can combine these equations (and use the fact that they are generally true for any vectors, not just the basis vectors) to find a relation between the quaternion and the angular velocity: \begin{equation} \boldsymbol{\omega}(t) = 2\, \dot{R}(t)\, \bar{R}(t). \end{equation} I'm sure that this has been done elsewhere, but my favorite derivation is in Section 3 of my paper.

There are a few ways to go from here, but the most obvious is to just use this as the differential equation you were looking for. Just rearrange a little to find \begin{equation} \dot{R}(t) = \frac{1}{2}\, \boldsymbol{\omega}(t)\, R(t). \tag{2} \end{equation} The right-hand side is a valid quaternion product, so you can just integrate the four components that result to get your frame quaternion $R(t)$. You also need an initial condition to solve this, of course, but that's to be expected.

If, instead, you have $\boldsymbol{\omega}'$ (measured with respect to the body), you can simply rotate that back into the inertial frame as $\boldsymbol{\omega}(t) = R(t)\, \boldsymbol{\omega}'(t)\, \bar{R}(t)$, in which case you have \begin{equation} \dot{R}(t) = \frac{1}{2}\, R(t)\, \boldsymbol{\omega}'(t). \tag{3} \end{equation} This is the more relevant equation if, for example, you measure angular velocity using some gyroscope embedded in the body.

Finally, I'll point out that it's also possible to integrate to find the generator of the rotation, which I'll label $\boldsymbol{r}(t)$. This is precisely the angle/axis vector, where the axis is given by the direction of the generator and the angle is proportional to the length of the generator. It is related to the quaternion by $R(t) = \exp[\boldsymbol{r}(t)]$. In this case, the differential equation is a bit more complicated: \begin{equation} \dot{\boldsymbol{r}} = \frac{1}{2} \boldsymbol{\omega} \times \boldsymbol{r} + \boldsymbol{\omega} \frac{\lvert\boldsymbol{r}\rvert \cot \lvert\boldsymbol{r}\rvert} {2} + \boldsymbol{r} \frac{\boldsymbol{r} \cdot \boldsymbol{\omega}} {2 \lvert\boldsymbol{r}\rvert^2} \left(1 - \lvert\boldsymbol{r}\rvert \cot \lvert\boldsymbol{r}\rvert \right). \end{equation} Besides looking uglier, this equation also requires some careful handling, because $\boldsymbol{r}$ typically goes through some weird values that are precisely analogous to branch cuts of the complex logarithm, but which cause serious numerical problems. If you're careful about handling this behavior, you can make this system as good as (and occasionally even better than) the systems of Eq. (2) and (3), but it's probably not worth the effort. These issues are all discussed at length in my paper.

Solution 2:

First let me quickly explain what angular velocity is.

Consider two right-handed, orthonormal coordinate systems with coincident origins. We will consider one of these coordinate systems to be the "world" (some kind of reference) and the other to be the "body" (the thing who's orientation we are concerned with).

A point $p$ can be expressed in either world coordinates or body coordinates, related by a rotation matrix. $$p_w = R\ p_b$$

In this case, I have defined $R$ to convert body coordinates to world coordinates, but of course its inverse will perform the reverse transformation.

When we say that the body is "rotating" relative to the world, what we mean is that $R = R(t)$ is really a function of a scalar called time. So now lets take derivatives and employ the chain rule. Let $\dot{x} := \frac{dx}{dt}$. $$\dot{p}_w = R\ \dot{p}_b + \dot{R}\ p_b$$

What can be said about $\dot{R}$? This is where we have to make use of our choice of consistently-handed orthonormal coordinate systems, which is practically always the case. Under this assumption, $R$ is an orthonormal matrix, meaning $R^{-1}=R^T$ or, $$R^TR = I$$

Taking the derivative and employing the chain rule again, $$R^T\dot{R} + \dot{R}^TR = 0$$

Which implies, $$R^T\dot{R} = -\dot{R}^TR$$

The derivative and transpose obviously commute, so we realize that the right-hand-side is equal to the negative-transpose of the left-hand-side. This means that the matrix $R^T\dot{R}$ is skew-symmetric. Let's name it, $$\Omega := R^T\dot{R}$$

This is called the angular velocity matrix (in body coordinates), and being as it is a skew-symmetric 3x3 matrix, it only has 3 unique elements, so we can extract them and call it the angular velocity vector, $$\omega := \begin{bmatrix} \Omega_{2,1} \\ \Omega_{3,1} \\ \Omega_{3,2} \end{bmatrix}$$

for which the cross-product operation is equivalent to matrix multiplication, $$\omega \times v \equiv \Omega v,\ \ \forall v \in \mathbb{R}^3$$

Coming back to our original problem, we can substitute in an $RR^T=I$ to see, \begin{align*} \dot{p}_w &= R\ \dot{p}_b + RR^T\dot{R}\ p_b\\ &= R\ \dot{p}_b + R\ \Omega\ p_b\\ &= R(\dot{p}_b + \omega \times p_b)\\ \end{align*}

which is sometimes called the transport theorem, but is more just the very recognizable equation for relating velocities in different reference frames. You must "augment" $\dot{p}_b$ with its rotational velocity $\omega \times p_b$ before you can convert it to world coordinates. This is my "proof" to you that the $\omega$ I have derived here is the same $\omega$ you are familiar with.

Now then, we are ready to talk about the problems of "integrating" angular velocity. Considering, $$\Omega := R^T\dot{R}$$ we are reminded that we don't really care about the integral of $\Omega$. We only care about the integral of $\dot{R}$. Worded another way, the orientation, $R$, of the body evolves as, $$\dot{R} = R \Omega$$

Since $\Omega$ itself is also a function of time, this linear differential equation has a solution in terms of a (non-commutative) product integral, $$R(t) = R(t_0) \prod_{t_0}^t e^{\Omega(\tau) d\tau}$$

Notice that the matrix exponential of skew-symmetric matrices is always orthogonal, and the product of orthogonal matrices is always orthogonal (plenty of proofs out there). Thus the above equation is elegantly consistent.

Anyway, the key point here is that the angular velocity integral (in the exponent) is meaningless on its own, and as many people have pointed out is not on its own equal to change in orientation. Instead, it is related to change in orientation through that above equation.

Quaternions can come into play if you change $\Omega$ to $\omega$ in the exponent and consider the quaternion version of Euler's identity.

If your two coordinate systems do not have coincident origins, or rather are translating with respect to each other, the very first equation I wrote will have an "affine" offset, which you can deal with in many ways, but I'd leave that to you. It doesn't affect the evolution of $R$.

The angle-axis representation that you mention in your question can be introduced as follows. If we consider $\Omega$ to be constant over some (likely short) timestep, $\Delta t$, then, $$R(t_0+\Delta t) = R(t_0)e^{\Omega \Delta t}$$

The orthogonal marix $e^{\Omega \Delta t}$ can be considered a "rotation caused by $\Omega$" and it will have a central axis aligned with $\omega$ and an angle of rotation equal to $||\omega||\Delta t$. This is the axis-angle interpretation of $\omega$ your book refers to. It can be used to "step forward" $R$ by $\Delta t$ seconds, where $\Delta t$ is the period overwhich $\Omega$ is constant.

Final thing to mention: $e^{\Omega \Delta t}$ is, in practice, computed using Rodriguez's formula instead of taking an actual matrix exponential, but both do work.

Solution 3:

Integral of angular velocity vector has no physical meaning. It could be meaningful if only the rotation takes place through a single axis.

Let me explain it for the most general case. Say the frame {T} is rotated for about a degrees through x-axis (roll), then rotated for about b degrees through y-axis (pitch), then rotated for about c degrees through z-axis (yaw). Basically, a, b and c can be considered as 3 Euler angles that represent the successive rotations take place in the frame {T}. (All rotations took place successively)

There is no such thing as Rotation vector. You can stack [a b c], but this does not constitute a vector. Orientation of a frame can be configured with respect to another frame. Thus, a vector cannot represent the orientation quantity.

Since [a b c] vector DOES NOT exist, its time derivative [da db dc] CANNOT represent the related angular velocity. In order to find the angular velocity, let me use MATLAB Symbolic Toolbox.

Let's consider that Ra, Rb and Rc are the rotation matrices that are respectively associated with a b and c rotations.

$$ Ra = \begin{bmatrix} \cos(a) & 0 & sin(a) \\ 0 & 1 & 0 \\ -\sin(a) & 0 & cos(a) \end{bmatrix} $$ $$ Rb = \begin{bmatrix} 1& 0& 0 \\ 0 & \cos(b)& -\sin(b) \\ 0 & \sin(b)& \cos(b) \end{bmatrix}$$

$$Rc = \begin{bmatrix} \cos(c)& -\sin(c)& 0 \\ \sin(c)& \cos(c)& 0 \\ 0& 0& 1 \end{bmatrix}$$

The overall Rotation matrix, R is:

$$R = Ra Rb Rc = \begin{bmatrix} \cos(a)\cos(c) + \sin(a)\sin(b)\sin(c) & \cos(c)\sin(a)\sin(b) - \cos(a)\sin(c) & \cos(b) \sin(a) \\ \cos(b)\sin(c) & \cos(b)\cos(c) & -\sin(b) \\ \cos(a)\sin(b)\sin(c) - \cos(c)\sin(a) & \sin(a)\sin(c) + \cos(a)\cos(c)\sin(b) & \cos(a)\cos(b) \end{bmatrix}$$

In order to find the angular velocity, you can use the following equation:

$$Wsk = dR\cdot R'$$

Here, Wsk is the skew-symmetric form of angular velocity.
dR is the element-wise time derivative of R
R' is the transpose of R

The result is

$$W_{skew} = \begin{bmatrix} 0& db\sin(a) - dc\cos(a)\cos(b) & da - dc\sin(b) \\ dc\cos(a)\cos(b) - db\sin(a)& 0& - db\cos(a) - dc\cos(b)\sin(a) \\ dc\sin(b) - da & db\cos(a) + dc\cos(b)\sin(a) & 0 \end{bmatrix}$$

Using the tensor-vector identity that only applies to "angular velocity"

• $Wx = -Wskew(2,3)$ % Roll axis angular velocity, NOT EQUAL TO db
• $Wy = Wskew(1,3)$ % Pitch axis angular velocity, NOT EQUAL TO da
• $Wz = -Wskew(1,2)$ % Yaw axis angular velocity, NOT EQUAL TO dc

$$Wx =db\cos(a) + dc \cos(b) \sin(a)$$

$$Wy =da - dc \sin(b)$$

$$Wz =dc \cos(a) \cos(b) - db \sin(a)$$

If you want to go from the angular velocity to Euler angles, you need to solve the couple differential equations.

It is obvious that the time integral of Wx Wy Wz do not reveal a b c angles.

Over the years, I even saw professors who still confuse about this stuff.

Good luck, hope this helps.

Husn