Integrating body angular velocity

Question

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.

Answer 1

The proof is by "Fundamental Theorem Of Calculus".
By definition: $\dfrac{d\vec \theta}{dt}=\vec\omega$. By Fundamental Theorem Of Calculus $\vec \theta=\int\vec \omega dt$.

Answer 2

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

Answer 3

You might find the answer you're looking for in Joan Solà's 2017 paper, Quaternion kinematics for the error-state Kalman filter, section 4.6: Time integration of rotation rates.

Angular velocities are instantaneous measures, sampled at some interval $\mathrm{d}t$. You can integrate these by treating the angular velocity vector as a pure quaternion, and using quaternion exponentiation to integrate the velocity over interval $\mathrm{d}t$. See Eq. 218 for a single integration step of zeroth order, and Eq. 226 for a single integration step of first order. Quaternion exponentiation is covered in Eq. 46 and the unnumbered equation between Eq. 213 and Eq. 214. Proofs are inline in the paper.

Answer 4

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.

Answer 5

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.

Answer 6

Integrating angular rates using Euler angles -

Let FI be an inertial frame, and FB the body frame of an airplane. Then the FI frame can be transformed to the FB frame by a sequence of three rotations in the usual order yaw (about z axis), pitch (about yawed y axis), and roll (about yawed and pitched x axis).

Let the yaw, pitch, and roll Euler angles be y, p, and r. Let Ry, Rp, and Rr be the rotations matrices for these rotations about the z, y, and x axes.

Let F1 = Ry * FI be the yawed FI frame, and let F2 = Rp * FI be the pitched F1 frame.

Then FB = Rr * Rp * Ry * FI

Now we need to transform the body rotation rates (or small angles) rb', pb', yb' to the Euler angle rates r', p', y'.

Note that y is measured in the FI frame, p in the F1 frame, and r in the F2frame. We can transform y, p, and r rates to the FB frame as follows:

[rb', rb', ry'] = [r', 0, 0] + Rr * [0, p', 0] + Rr * Rp * [0, 0, y'] = M * [r', p', y']

So we have [r', p', y'] = Minverse * [rb', rp', ry']

and the Euler angle rates can be integrated.

Because I don't know the language for formulas I can't expand the matrices, so this answer is only meant to be indicative of the explicit solution. To see the matrices written out I have two good references http://www.chrobotics.com/library/understanding-euler-angles and https://www.princeton.edu/~stengel/MAE331Lecture9.pdf

Integrating angular rates using quaternions -

Quaternions make it easy to integrate angular rates. The salient facts are that a unit quaternion represents a rotation about a unit vector, that is,

Q = [cos(a/2) sin(a/2)*u1 sin(a/2)*u2 sin(a/2)*u3]

represents a rotation of angle a about the unit vector

u = [ 0 u1 u2 u3]

in this sense, let v be a vector, then

Q * v * Q**

gives the coordinates of v rotated about the u axis a radians.

Also, the rotation corresponding to Qpr = Qp * Qr is the rotation corresponding to the rotation corresponding to Qr followed by the rotation corresponding to Qp. Just what the Dr. ordered for integrating angular rates.

Suppose we have an airplane with a spinning gyroscope which maintains its axes aligned with the inertial axes (the axes of a laser gyro rotate with the airplane). The the gyro rate vector

s = [0 s1 s2 s3]

can be converted to a rotation quaternion as follows:

r = norm(x)

a = r*dt

Qs = [cos(a/2) sin(a/2)*s1/r sin(a/2)*s2/r sin(a/2)*s3/r]

Let Q be the initial orientation of the airplane, with the airplane frame aligned with the inertial frame. Then after each time interval dt Q is updated as follows:

Q = Qs** * Q

Now Q transforms vectors in the airplane frame F to the inertial frame I.

The coordinates of a vector v with coordinates vF in the F frame are given in the I frame by

vI = Q * vF * Q**

Note: * denotes multiplication, ** denotes conjugate

It's late, and it's taken me 3 days to get this to work, so ... if this be error and upon me proved .... Great refs: https://www.lucidarme.me/quaternions-and-gyroscope/ and http://www.cs.iastate.edu/~cs577/handouts/quaternion.pdf