# 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.

• position $= \int$ velocity. Orientation (i.e., angular position) $= \int$ angular velocity. – Eric Towers Apr 29 '14 at 7:04
• @EricTowers Big thumbs down on that. Integrating velocity to find position works because the generators of translation commute; integrating angular velocity to find angular position does not work because the generators of rotation do not commute (unless the rotation is always about a single axis). – Mike Dec 9 '16 at 15:12
• Actually @Mike if you have a time series of local (IMU) angular velocity measurements with small sample interval dt, then you can in fact integrate these angular velocities over time to get the global orientation, you just have to integrate them in time order (due to the non-commutativity of rotation, as you point out). – Luke Hutchison Nov 21 '19 at 8:58
• @LukeHutchison No, either you're misusing language or you'll get the wrong result. You can do things like the evaluating the time-ordered exponential, which involves a series of integrals. Or you can solve differential equations — which is sometimes called integration — like the ones I showed in my answer. But you cannot simply compute "∫ angular velocity", as Eric said, because you will get the wrong result. – Mike Nov 21 '19 at 14:37
• @Mike then how do you get the correct result? Are you saying it's actually in theory impossible to correctly integrate rotations over time? I am doing literally this exact thing, applying a 1st order Taylor series approximation of the time integral of angular velocities using quaternion exponentiation (0th order) with a 1st order correction term (see my separate answer on Joan Solà's 2017 paper, "Quaternion kinematics for the error-state Kalman filter"), and I am able to track absolute angular orientation with extremely high fidelity over time. (ctd..) – Luke Hutchison Nov 22 '19 at 15:28

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 $$\mathbf{v}_{\omega'} := \int_{t_0}^{t_1} \boldsymbol{\omega}' dt \tag{1}$$ 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: $$\mathbf{x}'(t) = R(t)\, \mathbf{x}\, \bar{R}(t),$$ 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: $$\frac{d \mathbf{x}'} {dt}(t) = \boldsymbol{\omega}(t) \times \mathbf{x}',$$ 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: $$\boldsymbol{\omega}(t) = 2\, \dot{R}(t)\, \bar{R}(t).$$ 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 $$\dot{R}(t) = \frac{1}{2}\, \boldsymbol{\omega}(t)\, R(t). \tag{2}$$ 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 $$\dot{R}(t) = \frac{1}{2}\, R(t)\, \boldsymbol{\omega}'(t). \tag{3}$$ 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: $$\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).$$ 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.

• I believe R(t) actually is exp[1/2 r(t)]. – Forrest Voight Mar 9 '17 at 22:16
• @ForrestVoight It really just depends on how you define it. For consistency with Lie algebra notation, you would use $\exp[r(t)]$. For consistency with the axis-angle notation usually used in physics, you might prefer to use $\exp[r(t)/2]$. I defined it using the former choice, because in my paper I also explained how to derive the analog of that last equation in more general Lie groups. In particular, that last equation would have things like $\cot \lvert \boldsymbol{r} \rvert/2$ if I used your preference. But either one is acceptable. – Mike Mar 9 '17 at 22:46

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.

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.

Good luck, hope this helps.

Husn

• thank you very much whoever edited! – Husnu Nov 26 '14 at 6:07
• You need to be careful when you say things like "There is no such thing as Rotation vector". It's actually possible to define several such things. The one the OP is talking about isn't very useful, but it is possible to define more useful ones. – Mike Feb 4 '16 at 20:07

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

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.

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$.

• then I don't understand why in the same paper equation 265 shows a complicated relationship between rotation vector rates and body angular velocity. – staple Apr 29 '14 at 15:29
• @staple I am not very much familiar with the notations used in that pdf book. It would be better if you could rephrase your question with a bit explanation of the terms used there. I may read it but I have to go now, its night, cya. – user103816 Apr 29 '14 at 15:39
• I tried to clarify, thanks for your help. – staple Apr 29 '14 at 16:59
• This is wrong. Angular displacement cannot be represented as a vector quantity except as a differential, because it does not obey the commutative law of addition (except as a differential). – Radon Rosborough May 10 '16 at 14:48
• @raxod502 I don't remember it now. – user103816 May 11 '16 at 5:16