Quaternion

Preface

Let’s now face another central argument to solve the problem to find the attitude of an aircraft in a clean, consistent, reliable way on a digital computer. To follow this post is essential some basic knowledge of complex number (see “Complex numbers” post)

Quaternion, this amazing mathematical entity

Hamilton in 1843 introduced this new mathematical entity: here we don’t want to enter into specific mathematical important issues, like the consistency and the completeness of the theory. We assume its validity and we want only understand its meaning and its use.

The definition I like is my own interpretation of this entity. Let’s make this observation: let’s assume that a three dimensional vector is an extension of one dimensional representation of a real number. The quaternion has the same meaning, for complex number: they are a three dimensional extension of the imagery part of a complex number, that is:

$q=q_0+i{q}_1+j{q}_2+k{q}_3$   1)

where i, j, k are three unitary pure imaginary number (in analogy with i, j, k versors of a real three dimensional reference system).

q₀ is called the real or scalar part of the quaternion, while $i{q}_1+j{q}_2+k{q}_3$ is called the imaginary or vector part of the quaternion.

Furthermore i, j, k obey to the following rules:

$i^2=-1$ $j^2=-1$ $k^2=-1$

but they are non commutative respect to multiplication. In fact:

$\mathit{ij}=-\mathit{ji}$ $\mathit{kj}=-\mathit{jk}$ $\mathit{ki}=-\mathit{ik}$ It is easy to show that

$\mathit{ij}=k$ $\mathit{jk}=i$ $\mathit{ki}=j$ $\mathit{ji}=-k$ $\mathit{kj}=-i$ $\mathit{ik}=-j$

The conjugate of quaternion 1) is

$\tilde{q}=q_0-i{q}_1-j{q}_2-k{q}_3$

Using the shown quaternion properties

$q\ast \tilde{q}=q_0^2+q_1^2+q_2^2+q_3^2$

which is called the length or norm of the quaternion

Now we know almost all we need for using quaternions: one additional rules will be demonstrated later concerning the derivative of a quaternion.

Let’s now consider a vector V in the real world whose components are X, Y, Z. Let’s now “transport” our vector from the “real world”, in the “quaternion world”. In this new world our vector assume the form of a pure imaginary quaternion, that is

$V=iX+jY+kZ$

Let’s now examine the operation

$V\text{'}=\tilde{q}Vq$    2)

where q is a versor ( $q_0^2+q_1^2+q_2^2+q_3^2=1$ )

Equation 2) can be written as:

$V\text{'}=(q_0-i{q}_1-j{q}_2-k{q}_3)(\mathit{iX}+\mathit{jY}+\mathit{KZ})(q_0+i{q}_1+j{q}_2+k{q}_3)$   3)

Adopting the just learned quaternion rules, relation 3) becomes:

$V\text{'}=i[X(q_0^2+q_1^2-q_2^2-q_3^2)+Y(2q_0{q}_3+2q_1{q}_2)+Z(2q_1{q}_3-2q_0{q}_2)]\text{+}$

$\text{+}j[X(2q_1{q}_2-2q_0{q}_3)+Y(q_0^2-q_1^2+q_2^2-q_3^2)+Z(2q_0{q}_1+2q_2{q}_3)]\text{+}$

$\text{+}k[X(2q_0{q}_2+2q_1{q}_3)+Y(2q_2{q}_3-q_0{q}_1)+Z(q_0^2-q_1^2-q_2^2+q_3^2)]$   4)

Relation 4) can be written as:

$\left(\begin{array}{c}X\text{'}\\ Y\text{'}\\ Z\text{'}\end{array}\right)=\left(\begin{array}{ccc}(q_0^2+q_1^2-q_2^2-q_3^2)& (2q_0{q}_3+2q_1{q}_2)& (2q_1{q}_3-2q_0{q}_2)\\ (2q_1{q}_2-2q_0{q}_3)& (q_0^2-q_1^2+q_2^2-q_3^2)& (2q_0{q}_1+2q_2{q}_3)\\ (2q_0{q}_2+2q_1{q}_3)& (2q_2{q}_3-q_0{q}_1)& (q_0^2-q_1^2-q_2^2+q_3^2)\end{array}\right)\left(\begin{array}{c}X\\ Y\\ Z\end{array}\right)$   5)

which simply represents a coordinate transformation.

It is extraordinary important to note that relation 5) is identical to relation 3) of the post “Euler parameters”, when we pose:

$q_0=a_1=\cos \frac{\mu }{2}$    $q_1=a_2=\cos \alpha \sin \frac{\mu }{2}$   6)

$q_2=a_3=\cos \beta \sin \frac{\mu }{2}$   $q_3=a_4=\cos \gamma \sin \frac{\mu }{2}$

This means that operation 2) is equivalent to the rotation of the reference system around a vector represented by the imaginary (or vector) part of the quaternion, of an angle μ.

Furthermore, as anticipated in the post “Euler parameters”, it establishes an equivalence between the Euler parameters and the component of the quaternion representing the identical rotation.

Double rotations

The matter of two successive rotations may be easily handled. Consider the following transformation just seen above:

$V\text{'}=\widetilde{q_1}V{q}_1$   7)

and apply to $V\text{'}$ another rotation by means of quaternion q₂.

$V\text{'}\text{'}=\widetilde{q_2}V\text{'}{q}_2=\widetilde{q_2}\widetilde{q_1}V{q}_1{q}_2$   8)

Let’s call $q_3=q_1{q}_2$ and $q_4=\widetilde{q_2}\widetilde{q_1}$

Expression 8) takes the form

$V\text{'}\text{'}=q_{4}V{q}_3$     9)

It’s easy to demonstrate that both q₃ and q₄ are versors and that $q_4=\widetilde{q_3}$ and thus expression 9)

$V\text{'}\text{'}=\widetilde{q_3}V{q}_3$     10)

This means that two subsequent rotation, the first described by quaternion q₁, the second described by quaternion q₂, is equivalent to a rotation described by quaternion q₃ = q₁ q₂. This rule will be useful in the next paragraph

Quaternion differentiation

In this part we face the problem of quaternion differentiation which plays a fundamental role to find our solution. In the following I use a less formal mathematical approach, maybe not loved by mathematicians, but very effective for an engineer like me. Furthermore I’ve never seen in literature a similar approach. Let’s define the derivative of a quaternion as the limit of the incremental ratio when the increment going to zero, i.e.

$\frac{\mathit{dq}}{\mathit{dt}}=\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{q(t+\mathrm{\Delta }t)-q(t)}{\mathrm{\Delta }t}$    11)

Rotation represented by quaternion q(t+Δt) can be thought as a rotation until time t, and a subsequent infinitesimal rotation during time Δt, that we call $\mathrm{\Delta }\mu$ .Remembering what just learned in the previous paragraph this means:

$q(t+\mathrm{\Delta }t)=q(t)\ast q(\mathrm{\Delta }t)$     12)

Let’s consider now the quaternion $q(\mathrm{\Delta }t)$ related to the infinitesimal rotation $\mathrm{\Delta }\mu$ . According to the definition

$q(\mathrm{\Delta }t)=\cos \frac{\mathrm{\Delta }\mu }{2}+i\ast \sin \frac{\mathrm{\Delta }\mu }{2}\ast \cos \alpha +j\ast \sin \frac{\mathrm{\Delta }\mu }{2}\ast \cos \beta +k\ast \sin \frac{\mathrm{\Delta }\mu }{2}\ast \cos \gamma$  13)

We recall that $\alpha ,\beta$ and $\gamma$ represent the three angles the rotation axis makes with the three axis X, Y and Z.

Being $\mathrm{\Delta }\mu$ infinitesimal we can assume $\cos \frac{\mathrm{\Delta }\mu }{2}\simeq 1$ and $\sin \frac{\mathrm{\Delta }\mu }{2}\simeq \frac{\mathrm{\Delta }\mu }{2}$

Relation 13) becomes

$q(\mathrm{\Delta }t)=1+i\ast \frac{\mathrm{\Delta }\mu }{2}\ast \cos \alpha +j\ast \frac{\mathrm{\Delta }\mu }{2}\ast \cos \beta +k\ast \frac{\mathrm{\Delta }\mu }{2}\ast \cos \gamma =1+ϵ$   14)

where $ϵ=i\ast \frac{\mathrm{\Delta }\mu }{2}\ast \cos \alpha +j\ast \frac{\mathrm{\Delta }\mu }{2}\ast \cos \beta +k\ast \frac{\mathrm{\Delta }\mu }{2}\ast \cos \gamma$  15)

using 14), definition 12) becomes

$\frac{\mathit{dq}}{\mathit{dt}}=\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{q(t)q(\mathrm{\Delta }t)-q(t)}{\mathrm{\Delta }t}=\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{q(t)(1+ϵ)-q(t)}{\mathrm{\Delta }t}$  and finally

$\frac{\mathit{dq}}{\mathit{dt}}=\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{q(t)ϵ}{\mathrm{\Delta }t}=q(t)\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{ϵ}{\mathrm{\Delta }t}$  16)

but $\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{ϵ}{\mathrm{\Delta }t}$ represent the derivative of ε respect to time and relation 16) can be written as

$\frac{\mathit{dq}}{\mathit{dt}}=q(t)\frac{dϵ}{dt}$   17)

From 15) let’s calculate dϵ/dt

$\frac{dϵ}{\mathit{dt}}=\frac{1}{2}(ì\cos \alpha +J\cos \beta +k\cos \gamma )d\frac{\mu }{\mathit{dt}}=\frac{1}{2}(ì\cos \alpha +J\cos \beta +k\cos \gamma )\omega$

being the derivative of the angle of rotation around the axis of rotation equal to the angular velocity. Recalling that $\cos \alpha ,\cos \beta ,\cos \gamma$ are the cosines of the angles the axis of rotation forms with the X, Y and Z axis, we can affirm that

$\omega \cos \alpha ={\omega }_x=P$  component of the angular velocity along X axis

$\omega \cos \beta ={\omega }_y=Q$  component of the angular velocity along Y axis

$\omega \cos \gamma ={\omega }_z=R$  component of the angular velocity along Z axis

where we have used the classical aeronautical notation to identify the three angular velocities of roll, pitch and yaw around the X, Y, Z aircraft axis (P, Q, R).

Finally, adopting a synthetic notation

$\frac{\mathit{dq}}{\mathit{dt}}=\frac{1}{2}q(t)\overline{\omega }$    $\dot{q}=\frac{1}{2}q(t)\overline{\omega }$  18)

where

$\overline{\omega }=ì{\omega }_x+j{\omega }_y+k{\omega }_z=ìP+jQ+kR$

is the pure imaginary quaternion representing the overall angular velocity around the rotation axis

Relation 18) represents the key to calculate the quaternion components. Solving the differential equation 18) we know the $q_0,q_1,q_2,q_3$ components of the quaternion and from these components we will show how to calculate the attitude angles psi, theta and fi of the aircraft.

We now need to open a parenthesis over relation 18). In all text related to quaternions the derivative of a quaternion respect to time can assume the following expression:

$\dot{q}=\frac{1}{2}\overline{\omega }q(t)$

but in this case the considered rotation is related to a point around an axis of rotation identified by the quaternion, remaining in the same reference system, while expression 18) represent the rotation of the reference system around an axis of rotation identified by the quaternion, which is the case study we are interested on, as explained in the post “Euler parameters”.

Developing relation 18) we obtain:

$\dot{q_0}+i\dot{q_1}+j\dot{q_2}+k\dot{q_3}=\frac{1}{2}(q_0+i{q}_1+j{q}_2+k{q}_3)(ìP+JQ+kR)$

$\begin{array}{c}\dot{q_0}+i\dot{q_1}+j\dot{q_2}+k\dot{q_3}=\frac{1}{2}\text{(}iP{q}_0+JQ{q}_0+kR{q}_0-P{q}_1+kQ{q}_1-jR{q}_1\text{-}\\ \text{-}kP{q}_2-Q{q}_2+iR{q}_2+jP{q}_3-iQ{q}_3-R{q}_3\text{)}\end{array}$

$\begin{array}{c}\dot{q_0}+i\dot{q_1}+j\dot{q_2}+k\dot{q_3}=\frac{1}{2}(-P{q}_1-Q{q}_2-R{q}_3)+\frac{1}{2}i(P{q}_0+R{q}_2-Q{q}_3)\text{+}\\ \text{+}\frac{1}{2}j(Q{q}_0-R{q}_1+P{q}_3)+\frac{1}{2}k(R{q}_0+Q{q}_1-P{q}_2)\end{array}$   19)

if we adopt a matrix notation we can write 19) as

$\left(\begin{array}{c}\dot{q_0}\\ \dot{q_1}\\ \dot{q_2}\\ \dot{q_3}\end{array}\right)=\frac{1}{2}\left(\begin{array}{cccc}0& -P& -Q& -R\\ P& 0& R& -Q\\ Q& -R& 0& P\\ R& Q& -P& 0\end{array}\right)\left(\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right)$   20)

Equation 20) is the key mathematical relation that allows us, once solved, to determine the component of the quaternion when P, Q and R angular velocities are known. In the post LL_algo we focalize our effort in solving differential equation 20) with an optimize algorithm with the paramount characteristic to be immune from the “gimbal lock” problem.

Quaternion-Euler Angles relation

One more notion we will need in solving equation 20), is the definition of the initial value of the quaternion, giving the initial values of the Euler angles of the aircraft.

Let’s suppose that the initial attitude of the Aircraft is defined by the angles ψ, ϑ and φ.

Let’s start with ψ. We learned in the post “Euler Angles”, that ψ is the first rotation around the Z axis. Let’s try now to apply what we learned in this post about quaternions: how to express this rotation in terms of quaternions?

First af all we know the axis of rotation, the Z axis: this means that ɑ and β in 6) are $\frac{\pi }{2}$ while γ is equal 0. In addition we know that the angle of rotation is ψ: so the generic quaternion

$\left(\begin{array}{c}\cos (\frac{\mu }{2})\\ \cos (\alpha )\sin (\frac{\mu }{2})\\ \cos (\beta )\sin (\frac{\mu }{2})\\ \cos (\gamma )\sin (\frac{\mu }{2})\end{array}\right)$  assumes the form  $q_{\psi }=\left(\begin{array}{c}\cos (\frac{\psi }{2})\\ 0\\ 0\\ \sin (\frac{\psi }{2})\end{array}\right)$  21)

Repeating the same reasoning for , ϑ and φ we reach to the following :

$q_{\theta }=\left(\begin{array}{c}\cos (\frac{\theta }{2})\\ 0\\ \sin (\frac{\theta }{2})\\ 0\end{array}\right)$     $q_{\varphi }=\left(\begin{array}{c}\cos (\frac{\varphi }{2})\\ \sin (\frac{\varphi }{2})\\ 0\\ 0\end{array}\right)$   22)

The 3 quaternions above can be written as:

$q_{\varphi }=\cos (\frac{\varphi }{2})+i\sin (\frac{\varphi }{2})$    23)

$q_{\theta }=\cos (\frac{\theta }{2})+j\sin (\frac{\theta }{2})$    24)

$q_{\psi }=\cos (\frac{\psi }{2})+k\sin (\frac{\psi }{2})$   25)

In this post we have already learned how to face the issue of multiple rotation (see relation 9)). Applying the same concept, we can say that the resulting quaternion q_r which takes into account the 3 rotations expressed by 23), 24) and 25) is

$q_r=q_{\psi }{q}_{\theta }{q}_{\varphi }$    26)

substituting 23), 24), and 25) in 26)

$q_r=(\cos (\frac{\psi }{2})+k\sin (\frac{\psi }{2}))(\cos (\frac{\theta }{2})+j\sin (\frac{\theta }{2}))(\cos (\frac{\varphi }{2})+i\sin (\frac{\varphi }{2}))$

The result of this operation (taking into account all the rules related to multiplication among i, j and k, as discussed at the beginning of this post) is:

$q_r=\cos (\frac{\varphi }{2})\cos (\frac{\theta }{2})\cos (\frac{\psi }{2})+\sin (\frac{\varphi }{2})\sin (\frac{\theta }{2})\sin (\frac{\psi }{2})\text{+}$

$\text{+}i⟦\sin (\frac{\varphi }{2})\cos (\frac{\theta }{2})\cos (\frac{\psi }{2})-\cos (\frac{\varphi }{2})\sin (\frac{\theta }{2})\sin (\frac{\psi }{2})⟧\text{+}$

$\text{+}j⟦\cos (\frac{\varphi }{2})\sin (\frac{\theta }{2})\cos (\frac{\psi }{2})+\sin (\frac{\varphi }{2})\cos (\frac{\theta }{2})\sin (\frac{\psi }{2})⟧\text{+}$

$\text{+}k⟦\cos (\frac{\varphi }{2})\cos (\frac{\theta }{2})\sin (\frac{\psi }{2})-\sin (\frac{\varphi }{2})\sin (\frac{\theta }{2})\cos (\frac{\psi }{2})⟧$

which can be written in quaternion notation as

$q_r=\left(\begin{array}{c}\cos (\frac{\varphi }{2})\cos (\frac{\theta }{2})\cos (\frac{\psi }{2})+\sin (\frac{\varphi }{2})\sin (\frac{\theta }{2})\sin (\frac{\psi }{2})\\ \sin (\frac{\varphi }{2})\cos (\frac{\theta }{2})\cos (\frac{\psi }{2})-\cos (\frac{\varphi }{2})\sin (\frac{\theta }{2})\sin (\frac{\psi }{2})\\ \cos (\frac{\varphi }{2})\sin (\frac{\theta }{2})\cos (\frac{\psi }{2})+\sin (\frac{\varphi }{2})\cos (\frac{\theta }{2})\sin (\frac{\psi }{2})\\ \cos (\frac{\varphi }{2})\cos (\frac{\theta }{2})\sin (\frac{\psi }{2})-\sin (\frac{\varphi }{2})\sin (\frac{\theta }{2})\cos (\frac{\psi }{2})\end{array}\right)$   27)

Relation 27) will be used to initialize quaternion values once the 3 attitude angles are known