Chapter 2 Reference Frames and Roto-translations Before we start - - PDF document

▶

Mar 12, 2024 569 likes •1.32k views

Chapter 2 Reference Frames and Roto-translations Before we start to the study the kinematics and dynamics of rigid bodies and multi- body systems, it is appropriate to recall some geometrical concepts used to describe the basic quantities and

SLIDE 1

Chapter 2 Reference Frames and Roto-translations

Before we start to the study the kinematics and dynamics of rigid bodies and multi- body systems, it is appropriate to recall some geometrical concepts used to describe the basic quantities and the associate transformations transformations that characterize the motion of a rigid body in space. We will start with a formal definition of reference frames and then we will introduce the translation, rotation and roto-translation operators, that are essential for the study of motion of rigid bodies.

2.1 Tridimensional space

For simplicity, from now on, we will assume to be confined in a tridimensional world, except when we will study two-dimensional problem, as in planar motion or in robotic computer vision; consequently vectors will be described as elements of the 3D space R3, or E3 if the Euclidean norm is implicit. With no intention to raise philosophical questions, we can assume that the physical world around us, including the geometric entities we perceive, exist independently

f any reference frame. On the contrary, for modelling purposes, it is very often nec-

essary to express vectors with respect to one or more reference frames; we can say that fixing a coordinate system and the related reference frame “gives substance” to vectors: these can now be compared and measured relative to a common ruler. Moreover, suitable operations acting on vectors allow to determine, represents and measure geometric entities as angles, distances, orthogonality, projections, or physical quantities, as fields, powers, angular and linear velocities, etc. In principle, there is a difference between a reference frame and a coordinate system: the latter indicates the abstract structure used to define vectors, while the former 9

SLIDE 2

10 Basilio Bona - Dynamic Modelling indicates a specific way to define the parameters that characterize the vectors; examples of coordinate systems can be the cylindrical, the spherical or the rectangular (cartesian) coordinate system, while examples of (cartesian) reference frames are Earth centered inertial frame, the International Celestial reference frame (ICRF) and others [25]. In the following we will use indifferently the two terms to indicate a cartesian reference coordinate system, as specified in the following paragraphs.

2.2 Reference frames

A reference frame is defined by the symbol R(O, i, j , k) or R for short, where O is a particular point in space, called the origin, and i, j , k are three unit normal vectors, defining the metric properties of the space. Since it creates a loophole to define vectors in relation to a reference frame that is essentially built on vectors, it is necessary to understand how to construct R without making use of vectors and operations not yet defined. The only things we need in order to build a reference frame is the concept of square angles, i.e., orthogonality between lines, and three geometric directed segments or directed numbers − − → A′A, − − → B′B, and − − → C′C of equal length, that have the function of rulers of unit length in the 3D space. That of directed segments is a basic concept described in any introductory physics textbook, and has been briefly recalled in Appendix A.3. Having fixed the origin in O, we place the directed segments so that the three starting points A′, B′, and C′ coincide with O. Furthermore, the three segments must not be aligned with each other or lying all on the same plane. Usually these conditions are sufficient to characterize the reference frame, but, in

rder to preserve the orthogonal angles and orthogonal projections when using vector

products, we orient the three segments at angles of π/2 with respect to each other. These constraints allow only two possibilities, that are illustrated in Figure 1a) and b). The reference system in Figure 1a) is called a right-handed reference frame, while that of Figure 1b) is a left-handed reference frame. These names derive from the right-hand rule or the left-hand rule, illustrated in Figure ?? The most commonly used reference frame (from now on simply called RHF) is the first one, and we will tacitly adopt this convention in the rest of the book. In Figure 2.1, the three directed segments have been identified with the symbols i = − → OA j = − − → OB k = − → OC (2.1) The set {i, j , k} is called base vector set, and i, j , k are called base or basis vectors. We also say that {i, j , k} form a right or left orthonormal basis in R3. We are now able to use R(O, i, j , k) to characterize any other vector in R3.

SLIDE 3

Basilio Bona - Dynamic Modelling 11 Figure 2.1: a) right orthonormal reference frame; b) left orthonormal reference frame. Given a reference frame R, a vector v ∈ R3 is represented by three components, each

ne being the orthogonal projection of the directed segment on the three elements

i, j and k. This can be written as the linear combination of the base elements, and give origin to the well-known algebraic definition of vectors: v = v1i + v2j + v3k (2.2) where vi ∈ R, i = 1, 2, 3, are the three vector components with respect to R. In order to be able to define the values of vi, it is necessary to introduce a binary

peration between two vectors a and b, called the scalar vector product or dot

product, and defined as: a · b = a1b1 + a2b2 + a3b3 = b · a (2.3) We have defined this product and its properties in Appendix A.4.1; for the moment, relation (2.3) is sufficient for our aims. The three components vi are defined as: v1 = i · v v2 = j · v (2.4) v3 = k · v

SLIDE 4

12 Basilio Bona - Dynamic Modelling As i, j , k are vectors themselves, each one can be expressed according to (2.2): i = i1i + i2j + i3k j = j1i + j2j + j3k (2.5) k = k1i + k2j + k3k Since i is orthogonal to j and k, and has unit length, its first component must be equal to one, while the other two must be zero; the same argument holds for j and k as well. Hence, if we define e1 =   1   e2 =   1   e3 =   1   , (2.6) it follows that i = e1 j = e2 k = e3 (2.7) Therefore the representation of the reference frame R into itself can be given by the identity matrix [ i j k ] = [ e1 e2 e3 ] = I =   1 1 1   . (2.8) According to (2.2), each vector can be represented by a column of real components that we call a column vector: v = [ i j k ]   v1 v2 v3   = I   v1 v2 v3   =   v1 v2 v3   (2.9) In (2.9) we have implicitly assumed to know how to use the matrix notation and to understand the row by column product rule. Although notation (2.8) may appear pedantic and superfluous, we will see in Section ?? that it is possible to represent a reference frame R2 with respect to another reference frame R1 replacing the identity matrix I with a square orthogonal matrix R, having the specific geometric meaning properties, described in Section ??; Relation (2.9) is the simplest form of a more general relation v 1 = Rv 2 that provides the representation in R1 of a vector v 2 with components given in R2. To indicate the single components of a vector v one can adopt numerical k = 1, 2, 3

r literal k ∈ {x, y, z} indexes. While the literal indexes have an immediate meaning,

COMPLETARE COMPLETARE

SLIDE 5

Basilio Bona - Dynamic Modelling 13 they are more difficult to use when computer algorithms or mathematical formulas are considered; for example, the norm of a vector can be written as follows ∥v∥2 = ∑

k∈{x,y,z}

∥v∥2 =

∑

k=1

The second one is much more immediate, and except that some particular case, the numerical indexes will be adopted throughout these notes.li. Often we omit to specify the origin of a reference frame that is simply indicated as R(i, j , k). Other times we use an index to specify a particular reference frame; for example we use the symbol Rm(i m, j m, k m) to indicate a “local” frame and R0(i 0, j 0, k 0) to specify an “inertial” frame1. Given a specific reference frame Rm and a geometrical point P, this one is described in Rm by the (geometrical) vector v m

p that represents the oriented segment −

→ OP: v m

p = vm 1 i m m + vm 2 j m m + vm 3 k m m

(2.10) This notation puts in evidence the reference frame m that is used to represent P. The index m is omitted whenever there are no ambiguities. Given two reference frames R1 e R2, the same geometrical point P has two representations: v 1

p =

[ v1

]T in R1 (2.11) v 2

p =

[ v2

]T in R2 (2.12) In alternative, we can use one of the following symbols to indicate the representation

f v or −

− → UV in Rm: [− − → UV ]

[v]Rm [⃗ v ]Rm   v1 v2 v3  

(2.13)

2.3 Vector types

We represents vectors with a graphical icon: the most used icon is an arrow, as in Figure 2.2. Unfortunately this icon is sometimes misleading, as we will see, considering that there are two types of vectors with different properties polar vectors

1 We use the term inertial for pseudo-inertial reference frames, as those fixed in the environment.

SLIDE 6

14 Basilio Bona - Dynamic Modelling and axial vectors also called pseudovectors (for details see also [2]). Examples

f polar vectors are those used for geometrical points, linear velocities and acceler-

ations, forces, gradients, normals, unit vectors, etc., while examples of axial vectors are angular velocities, torques, moments, cross products, etc. Figure 2.2: The vector icon. To make things clear one should associate a different icon for each type of vectors: the arrow for polar vectors, and a different icon, e.g., a segment with a small curl, to the axial vectors, as sketched in Figure 2.3. Figure 2.3: Two different icons: an arrow for polar vectors, a segment with a curl for axial vectors. Unfortunately this is not the case, and we have the same symbol for both; when we want to interpret the arrow icon for polar vectors we should perform the implicit transformation depicted in Figure 2.4. A formal characterization of the two types of vectors is possible considering the property of invariance with respect to reflections. Although reflections will be math- ematically defined only in Sections ?? and 2.8;, we can for the moment rely on the intuitive meaning. COMPLETARE

SLIDE 7

Basilio Bona - Dynamic Modelling 15 Figure 2.4: To interpret the arrow icon for an axial vector one shall perform the transformation depicted here, i.e., align the right hand thumb with the arrow axis to obtain the curl from the other fingers.

2.3.1 Polar vectors

Consider Figure 2.5, where a vector v 1 is reflected with respect to two orthogonal planes; the plane π1 parallel to v 1, and the plane π2 orthogonal to v 1. The resulting vectors are v 2 and v 3: since v 2 = v 1, the properties of the two vectors are the same, while, given that v 3 = −v 1, the reflection reverse the direction of the vector, and consequently its physical significance. Examples of physical polar vectors are displacements, linear velocities and forces, among others.

2.3.2 Axial vectors

Axial vectors are vectors that represent directed quantities that are antisymmet- rical with respect to a reflection through a parallel plane π1, and are symmetrical with respect to a reflection through a perpendicular plane π2, as shown in Figure 2.6. In this case the reflections have completely different effects: considering the vector v 3

btained by a reflection with respect to π2, the curl icon does not change direction,

so we can say that v 3 = v 1. Should we have used the arrow icon our vector would have changed direction and therefore its physical meaning. The result of a reflection with respect to π1 would change the curl direction, and, using the right-hand rule, this means a reverse in effect, i.e., v 2 = −v 1. Examples of axial vectors are angular velocities and torques. The distinction between axial and polar vectors is often omitted in mechanical textbooks, but, among other application, it will become important when using quaternions to represent rotations, as presented in Sections 2.12.6 and 2.12.7.

SLIDE 8

16 Basilio Bona - Dynamic Modelling Figure 2.5: The symmetry and anti-symmetry of polar vectors. It is important to add that the cross product of two polar vectors produce an axial vector a (polar) × b (polar) = c (axial) The interested reader can find furthere details in [2] and [51].

2.4 Geometrical and Physical vectors

In these notes vectors will be used to represent both geometrical points, as the set

f vertices in a rigid object, the barycenter of a plate, the location of a fixture, the

position of a joint, etc., and signed quantities having a physical significance, as linear

r angular velocities and accelerations, forces, moments, torques, etc.

We will indicate the former as geometrical vectors, while the latter are known as physical vectors.

2.4.1 Geometrical vectors

A geometric point P ∈ R3, is described by the associated vector p that contains the representation of the directed segment − → OP in R(O, i, j , k).

SLIDE 9

Basilio Bona - Dynamic Modelling 17 Figure 2.6: Axial vector symmetry. Given two different reference frames, R1(O1, i 1, j 1, k 1) and R2(O2, i 2, j 2, k 2), the same geometric point P is associated to two different vectors p1 and p2, respectively the representation of P in R1 and R2. As an example, a velocity defined in the plane by the vector ⃗ v in Figure 2.7, had a representation v a = (va1, va2) in the cartesian reference frame Ra that is different from the representation v b = (vb1, vb2) in Rb. The relation between p1 and p2 will be discussed in Section 2.9.

2.4.2 Physical vectors

A physical vector − → QP is an oriented segment, also known as directed segment

r directed quantity, that represents a physical quantity, such as linear or angular

velocity, gravitational acceleration, force, torque, etc. We cite some lines from [24] (1) A physical vector is a quantity with some physical origin. However general this may be, it already expresses some other interpretation of a vector, since the mathematical vector space axioms make no requirement as to the origin or characteristic of the vectors.

SLIDE 10

18 Basilio Bona - Dynamic Modelling Figure 2.7: A vector as an abstract entity and its representation in two different reference frames. The vector ⃗ v is the same, but the two representations have different components in Ra and Rb. (2) A physical vector has a magnitudo, which is not part of the initial description of a mathematical vector. If, however, one introduces magnitudes by means of the additional structure of an inner product, these are real numbers, and not, as in physics, dimension-laden physical scalars. (3) Finally, a physical vector has a direction in (physical) space, because the physical vector spaces described above have a close relation to the position vector space. There is no correspondence here with the mathematical concept of a vector, since the axioms make no mention of a physical space. A physical vectors has an application point Q, that can be free or constrained, a direction, and a magnitude, as in Figure 2.2. It is usual to represent the vector − → QP by a difference between two geometrical vectors, as in Figure 2.9. ⃗ v = v QP ∈ R3 = v Q − v P =   p1 − q1 p2 − q2 p3 − q3   ≡   v1 v2 v3   It should be noted that these vectors can be of two types: what are called free vectors and what are called applied or point vectors. With free vectors, the application point Q, also called application point, has no physical or geometric significance; one can translate v parallel to itself without producing a variation in any effect it could

SLIDE 11

Basilio Bona - Dynamic Modelling 19 Figure 2.8: The oriented segment is represented by the difference of two geometrical points: ⃗ vQP = v P − v Q.

generate. On the contrary, for applied vectors, the application point U cannot be

changed without affecting in some way the physical significance of the effects it

represents. As an example, a force acting on a rigid body is an applied vector, as
ne cannot change its application point without, in general, affecting the torque

acting on the mass center of the body itself; the linear velocity of a rigid object,

n the contrary, can be translated in any other point of the object without causing

much trouble.

2.4.3 Vector units

When it is necessary to assign physical units to vector components, two approaches are possible: 1) assign the units to the axes; 2) assign the units to each component. The first approach is more elegant, but often we have vector components with different units; for instance, the state vector of a particle moving on the plane has two components, namely its position and its velocity. In this case it is more convenient to have adimensional axes and assign units to the components; however this approach produces some interpretation problems when we make a scalar product

r we compute the norm of the vector, that will result in the square root of space

squared plus velocity squared. Nevertheless it is a common rule to apply the second approach, that will be used also in the present notes.

2.5 Rigid Bodies and Their Displacements

With the term rigid body we define any tridimensional object for which the distance between any couple of its points remain constant in time, independently from any motion or any force or torque applied to it. Similarly we call rigid displacement

SLIDE 12

20 Basilio Bona - Dynamic Modelling

r rigid motion the motion of a rigid body in space.

A rigid body is an abstraction, since in nature (or at least in the macro-world where mechatronic modelling is used) perfectly rigid bodies do not exist. Every object is distorted, warped or deformed under the effects of static or dynamic forces or torques; similarly, the rigid motion is an abstraction. Consider a steel plate that under the effect of its mass and gravity acceleration (a static force) bends in various different shapes according to its position with respect to the vertical. Also not considering gravity, if we take a steel plate and accelerate it, it will flex, and return into its original shape only after the resulting vibrations are damped. Nevertheless those abstractions will make the characterization of the tridimensional space geometry simpler and rich of theoretical and practical developments. So, as put in evidence in the introduction, we will make this approximating assumption and study the models of rigid displacements. It has been demonstrated by Chasles (1830) that any rigid displacement in space can always be decomposed into a translation and a rotation; more precisely the Chasles theorem states that the most general displacement of a rigid body in R3 is the composition of a rigid translation along a line and a rotation around an axis parallel to the same line [3, 4]. Given the characteristics of invariance of the distances between any two points in a rigid body, it is possible to describe a rigid body B with respect to a reference frame RB (i B, j B, k B), since once the origin and the orientation of the body reference frame is known, all B points are also known or can be determined.

2.6 Translations

Given a vector v = [ v1 v2 v3 ]T, a rigid translation is the operator Trasl(v, t) that displaces v parallel to itself of a given vector t: Trasl(v, t) ≡ v t def =   v1 + t1 v2 + t2 v3 + t3   = v + t (2.14) Since the oriented segment − → AB is the difference between − − → OB and − → OA ⃗ v AB = v B − v A (2.15) the translation of an oriented segment is simply ⃗ v t

AB = v t B − v t A = v B + t − (v A + t) = v B − v A = ⃗

v AB (2.16) from which we deduce that the representation of an oriented segment is invariant with respect to translations. Since oriented segments usually represent physical

SLIDE 13

Basilio Bona - Dynamic Modelling 21 quantities, as forces, torques, velocities, etc., we say that their representation is invariant to translations. Now assume to have two different reference frames Rk and Rm, having the same

rigin and the same basis {i k, j k, k k} = {i m, j m, k m} = {i, j , k}, and assume to

translate Rm with respect to Rk so that its new origin is now O′. We call tk

m the

riented segment −

− → OO′, that represents the translation from reference frame k to reference frame m. A geometrical point P, represented in Rm by − − → O′P = v m

P will be represented in Rk

by − → OP = v k

P = v m P + tk m

(2.17) and so we can conclude that, while the representation of physical vectors is unaf- fected by rigid translations as in (2.16), the representation of a geometrical point in a translated reference frame adds the translation vector, as in (2.17). The inverse operator of the translation is defined as Trasl(v, t)−1 = Trasl(v, −t) = −Trasl(v, t) = v m

p = v i p − tk m = v k p + tm k

and −tk

m = tm k

Since the transaction operator is represented by a vector sum, it is commutative Trasl(v, a)Trasl(v, b) = Trasl(v, b)Trasl(v, a) = Trasl(v, a + b) i.e., given n translations ti, i = 1, . . . , n, the total translation is t =

∑

i=1

ti Strange as it may appear, the translation of a quantity t is not a linear operator, since it does not obeys the axioms listed in (A.1). Indeed, the axioms will require Trasl((λ1a + λ2b), t) = λ1a + λ2b + t instead it results λ1Trasl(a, t) + λ2Trasl(b, t) = λ1(a + t) + λ2(b + t) = λ1a + λ2b + (λ1 + λ2)t. The two relations coincide only if λ1 + λ2 = 1; this fact characterize the properties

f an affine space and not of a linear space, as described in [14].

SLIDE 14

22 Basilio Bona - Dynamic Modelling

2.7 Rotations

A well known theorem by Euler, published in 1776, states that in three-dimensional space, any displacement of a rigid body such that a point on the rigid body remains fixed, is equivalent to a single rotation of a given angle about some axis that contains the fixed point. In modern symbols we will indicate the rotation axis with a unit vector u and the angle with θ, such that the rotation is characterized by the vector v = uθ; unfortunately, as we will see, the composition of two or more rotations cannot be reduced to the sum of the related vectors; a more complex characterization is necessary. We can describe the rotation assuming a rigid body with its reference frame and considering a rigid displacement that leaves the reference frame origin O fixed, while the basis unit vectors are changed under the rotation. A rotation is therefore characterized by the mathematical relation between these two reference frames. Figure 2.9: An example of rigid rotation of a frame with respect to the origin O. In order to represent a rotation we must consider two reference frames Rm(O, i m, j m, k m) and Rn(O, i n, j n, k n), with common origin O; initially they are coincident, i.e., each unit vector of the first base {i m, j m, k m} coincides with the corresponding unit vec-

SLIDE 15

Basilio Bona - Dynamic Modelling 23 tor of the second base {i n, j n, k n}. Now we take one frame, for instance Rn and rotate it around the common fixed

rigin O of a arbitrary angle θ. The rotation axis is itself arbitrary, and we call u

the unit vector that represents it; the only condition is that the origin O must lie

n the rotation axis.

At the end of the rotation Rn has taken a different “orientation” with respect to Rm, and the two basis will be no more completely coincident, although it is always possible that some of the unit vectors remain coincident. We can now define the representation of each basis unit vector of Rn in Rm, as follows i m

n def

= im

n1i m + im n2j m + im n3k m =

   im

   (2.18) j m

n def

= jm

n1i m + jm n2j m + jm n3k m =

   jm

   (2.19) k m

n def

= km

n1i m + km n2j m + km n3k m =

   km

   (2.20) Considering the above relations and recalling that the components of a generic vector v =   v1 v2 v3   are given by v1 = v · i = v Ti = i Tv v2 = v · j = v Tj = j Tv v3 = v · k = v Tk = k Tv it is possible to introduce a matrix Rm

n , whose columns are the representations of

the basis unit vectors of Rn in Rk: Rm

n =

[ i m

j m

k m

] =    im

   (2.21) Therefore this matrix Rm

n represents the transformation that describes the rotated

frame Rn with respect to the fixed frame Rm, with common origins. According to (2.21), Rm

n has the following structure:

SLIDE 16

24 Basilio Bona - Dynamic Modelling

the first column is the representation of the first unit basis vectors of Rn in

Rm,

the second column is the representation of the second unit basis vectors of Rn

in Rm,

the third column is the representation of the third unit basis vectors of Rn in

Rm, It follows that the matrix Rm

n can be interpreted as the representation of Rn in Rm.

By a convention adopted in many textbooks, the lower index n of R denotes the represented reference frame (the rotated one), while the upper index m denotes the reference frame where we represent it (the fixed one). Notice that the terms “rotated” or “fixed” are only used to give a meaning to a relative displacement and does not always imply a “real” rotation; we can imagine as well to leave Rn “fixed” and rotate Rm. We can build Rn

m, instead of Rm n using similar arguments.

In conclusions, we have Rm

n =

    i T

mi n

j T

mi n

k T

mi n

i T

mj n

j T

mj n

k T

mj n

i T

mk n

j T

mk n

k T

mk n

    and Rn

m =

    i T

ni m

j T

ni m

k T

ni m

i T

nj m

j T

nj m

k T

nj m

i T

nk m

j T

nk m

k T

nk m

    and, by inspection, Rm

n = (Rn m)T

(2.22) Since the matrix Rm

n describes a rigid rotation around an axis that goes through

the common origin of two reference frames, this matrix is therefore called rotation matrix. Equation (2.22) shows that the inverse rotation from n to m is represented by the transpose of the matrix that represents the rotation from m to n. A matrix whose inverse is equal to its transpose is called orthonormal (or, as in some textbooks,

rthogonal).

All rotation matrices are orthonormal matrices, whose properties have been described in Appendix B.7.1 Now we consider two vectors, the first one indicated as v n

P is the representation of a

geometrical point P in Rn, while v n

AB = v n B −v n A is the representation of an oriented

segment − → AB in Rn.

SLIDE 17

Basilio Bona - Dynamic Modelling 25 We want to know how to represent both in the reference frame Rm. It turns out that the formula is: v m

P = Rm n v n P

and v m

AB = Rm n v n AB

i.e., when the reference frame origins are the same, both types of vectors transform in the same way v m = Rm

n v n

(2.23) and conversely v n = Rn

mv m = (Rm n )Tv m

(2.24) We conclude that a rotation is represented by a square matrix and the rotation

perator applied to a vector is

Rot(v, R) = Rv In conclusion, when a generic (orthonormal) rotation matrix Rm

n is considered, all

f the three following characterizations are true
1. Rm

n represents a geometrical rotation Rot(u, θ) of an angle θ around an axis

whose unite vector is given by u (θ > 0 is given by the right-hand-rule), that brings Rm to overlap with Rn (from m “fixed” a n “rotated” or “mobile”). The value of the angle θ and the components of u do not appear immediately from the matrix elements, but we can compute them, as shown in Section 2.8.1;

2. Rm

n characterize the description of the unit basis of Rn in Rm (frame n “ro-

tated” or “mobile” in the frame m “fixed”);

3. Rm

n is the representation of the linear operator that transforms a vector from

the reference frame Rn (“rotated” or “mobile”) into the reference frame Rm (“fixed” ). The rotation operator is linear, since it is represented by a matrix that obeys tot the following property Rot(Rm

n , λ1v n 1 + λ2v n 2) = λ1Rm n v n 1 + λ2Rm n v n 2 = λ1Rot(Rm n , v 1) + λ2Rot(Rm n , v 2)

If we need to make more than one rotation, each one represented by a matrix Rot(ui, θi) = Ri−1

with respect to the reference frame obtained by the previous rotation, the total rotation is given by the ordered product of the rotation operators, as Rot(u1, θ1)Rot(u2, θ2) · · · Rot(uN, θN) = Rot(u, θ) (2.25)

SLIDE 18

26 Basilio Bona - Dynamic Modelling It is very important to note that the total angle θ is not the sum of the single angles θ ̸= θ1 + θ2 + · · · + θN and the same applies for the unit vector u u ̸= u1 + u2 + · · · + uN The matrix representing the global rotation is the ordered product of the single matrices Rot(u, θ) ⇔ R0

N def

= R0

1R1 2 · · · RN−1 N

∏

k=1

Rk−1

(2.26) Since matrix product does not commute (apart from some particular cases, that we will see later), the factors order is important. The upper index of a rotation matrix shall be equal to the lower index of the preceding matrix (on its left). We will see in Section 2.10.2 that the order of the product is connected to a precise geometrical concept and we will provide a simple mnemonic rule to remember it.

2.7.1 Elementary Rotations

We call elementary rotations or basic rotations the rotations that take place around some particular axis; the first one is the matrix that represents a rotation of an angle θ around a generic axis, defined by the unit norm vector u = [ u1 u2 u3 ]T with∥u∥ = 1: Rot(u, θ) ≡ R(u, θ) ≡ Ru,θ

def

=   u2

1 (1 − cθ) + cθ

u1u2 (1 − cθ) − u3sθ u1u3 (1 − cθ) + u2sθ u1u2 (1 − cθ) + u3sθ u2

2 (1 − cθ) + cθ

u2u3 (1 − cθ) − u1sθ u1u3 (1 − cθ) − u2sθ u2u3 (1 − cθ) + u1sθ u2

3 (1 − cθ) + cθ

  (2.27) where we have adopted the following conventions sθ

def

= sin θ cθ

def

= cos θ The determinant is det R(u, θ) = 1 while the trace is tr R(u, θ) = (1 − cθ)(u2

1 + u2 2 + u2 3) + 3cθ

(2.28) Other basic rotation are those around the x, y and z axes; in these case u is given by the basis unit vectors i, j e k, as follows:

SLIDE 19

Basilio Bona - Dynamic Modelling 27

Rotation around the x axis of an angle α:

Rot(x, α) ≡ Rot(i, α) ≡ R(i, α) ≡ Ri,α

def

=   1 cα −sα sα cα   (2.29)

Rotation around the y axis of an angle β:

Rot(y, β) ≡ Rot(j , β) ≡ R(j , β) ≡ Rj ,β

def

=   cβ sβ 1 −sβ cβ   (2.30)

Rotation around the z axis of an angle γ:

Rot(z, γ) ≡ Rot(k, γ) ≡ R(k, γ) ≡ Rk,γ

def

=   cγ −sγ sγ cγ 1   (2.31) where sα

def

= sin α; cα

def

= cos α; sβ

def

= sin β; cβ = cos β; sγ

def

= sin γ; cγ = cos γ

Example 2.7.1

We want to compute the matrix representing the basic rotation around an axis given by the vector v =   1 1   The corresponding unit vector is u = v ∥v∥ =    

1 √ 2 1 √ 2

    Now, using (2.27) we have: R(u, π/2) =          − 1 √ 2 1 √ 2 1 √ 2 1 2 1 2 − 1 √ 2 1 2 1 2          (2.32)

SLIDE 20

28 Basilio Bona - Dynamic Modelling We can check that the determinant of R is equal to +1 det R(u, π/2) = (−1)1+2 ( − 1 √ 2 ) ( 1 2 √ 2 − ( − 1 2 √ 2 )) + (−1)2+2 ( 1 √ 2 ) ( 1 2 √ 2 − ( − 1 2 √ 2 )) = 1 4 + 1 4 + 1 4 + 1 4 = 1 (2.33) and that tr R = 1 2 + 1 2 = 1 Considering relation (2.28), given later, we have tr R = 1 + 2 cos θ from which cos θ = 0 → θ = ±π/2

From a simple inspection of the matrix in (2.29) we notice that it obeys to the

following relations: R(u, θ) = R(−u, −θ) R(u, −θ) = R(u, 2π − θ) = R(−u, θ) R(u, θ) = [R(u, −θ)]T R(−u, −θ) = [R(−u, θ)]T R(u, θ1)R(u, θ2) = R(u, θ1 + θ2) (2.34) We recall again that the sign of the rotation angle is given by the right-hand-rule, When we change the direction of the vector u, if we want to represent the same physical rotation we have to change also the sign of the angle θ; on the contrary we will obtain the inverse rotation, represented by the transpose matrix.

Example 2.7.2

We want to compose two different rotations, represented by the following matrices, in the stated order R(j , 45◦) =   0.7071 0.7071 1.0000 −0.7071 0.7071  

SLIDE 21

Basilio Bona - Dynamic Modelling 29 and R(k, 60◦) =   0.5000 −0.8660 0.8660 0.5000 1.0000   We have Ra = R(j , 45◦)R(k, 60◦) =   0.3536 −0.6124 0.7071 0.8660 0.5000 −0.3536 0.6124 0.7071   (2.35) While, exchanging the order we have Rb = R(k, 60◦)R(j , 45◦) =   0.3536 −0.8660 0.3536 0.6124 0.5000 0.6124 −0.7071 0.7071   (2.36) with Ra ̸= Rb

2.7.2

Planar rotations

If the motion takes place in a plane, it is customary to define the plane with the common reference frame, where i is horizontal, j is vertical and k is orthogonal to the plane and points toward the eyes of the reader. In this context, all rotations are around the unit vector k (positive counterclockwise) as in Figure 2.10, and the rotation matrices can be written as   cγ −sγ sγ cγ 1  

r simply as

[cγ −sγ sγ cγ ] In this case the matrix product commutes, since, according to the last of (2.34) [cθ −sθ sθ cθ ] [cγ −sγ sγ cγ ] = [cγ −sγ sγ cγ ] [cθ −sθ sθ cθ ] = [cθ+γ −sθ+γ sθ+γ cθ+γ ]

SLIDE 22

30 Basilio Bona - Dynamic Modelling Figure 2.10: A planar rotation of an angle θ.

2.8 The Rotation Matrix

In the following Section we will define the main mathematical properties of the rotation matrices. A detailed description can be found in several textbooks, as [1, 3, 5, 13, 17, 31, 32, 34, 35, 39, 41, 48, 52]. A generic rotation matrix is always square and orthonormal, and its properties have been described in B.7.1. It is interesting to note that the translations, the rotations and the reflexions are all members of a the set of the so-called isometries, i.e., those transformations of the three-dimensional space that keep constant the Euclidean distance d(P, Q) = ∥xP − xQ∥ between every couple of points belonging to a rigid body. Among the isometries there is also the so-called glide, i.e., a reflexion followed by a translation. However, among the isometries that are represented by orthogonal matrices, we are interested only in the rotations, since the reflections are nota applicable to rigid bodies: your left hand cannot be transformed into your right hand.

2.8.1 Rotation Matrix Properties

All rotations, also called proper rotations, are expressed by square orthonormal matrices R with determinant det R = +1; if det R = −1 the matrix represents a reflexion, also called roto-reflections. According to the Euler theorem any rotation or composition of more rotations is always represented by a single matrix R(u, θ). The subspace represented by u has dimension 1, and is invariant with respect to the rotation (since the point on it do

SLIDE 23

Basilio Bona - Dynamic Modelling 31 not rotate); hence one can write Ru = u (2.37) This relation correspond to the classic definition of the eigenvalues of a matrix Ru = λu with λ = 1. A proper rotation R has the following canonical modal decomposition R = M ΛM H = [ u v v ∗]   1 ejθ e−jθ   [ u v v ∗]H (2.38) where j

def

= √−1 and M H is the hermitian matrix (i.e., conjugate-transpose) of M . In (2.38) the matrix M is the modal matrix, i.e., the matrix of the eigenvalues of R; as said above, u identifies the rotation axis, while v and v ∗ define the plane normal to u. On this plane we have a planar rotation that can be represented by the complex rotation operator (also called a phasor) ejθ def = cos θ + j sin θ The right-handed reference frame R(u, Re(v), Im(v)) has i = u, while the other two basis unit vectors are j = Re(v) and k = Im(v). The determinant is det R = 1 · ejθ · e−jθ = +1 as we already knew, while the trace is tr R = 1 + (cos θ + j sin θ) + (cos θ − j sin θ) = 1 + 2 cos θ that is equal to (2.28) for ∥u∥ = 1 If we introduce S(u) as the skew-symmetric matrix S(u) =   −u3 u2 u3 −u1 −u2 u1   (2.39) R obeys to the following relations: R(u, θ)

def

= eS(u)θ def = I + sin θ ∥u∥ S(u) + 1 − cos θ ∥u∥2 S 2(u) (2.40)

SLIDE 24

32 Basilio Bona - Dynamic Modelling R(u, θ)T def = e−S(u)θ def = I − sin θ ∥u∥ S(u) + 1 − cos θ ∥u∥2 S 2(u) (2.41) where S 2(u) =   −u2

2 − u2 3

u1u2 u1u3 u1u2 −u2

1 − u2 3

u2u3 u1u3 u2u3 −u2

1 − u2 2

  (2.42) The two formulas (2.40) and (2.41) are very important, since they define in the most general way the rotation matrix as an exponential of an anti-symmetric matrix; we omit the proof since it requires some notion of differential geometry, and in particular, of the Lie algebras; the interested reader can look at the textbook [34] for further elements. Notice that a planar rotation is represented by ejθ while a tridimensional rotation is represented by eS(u)θ

2.8.2 Rotation Matrix Parametrization

We are now able to solve the following two problems Problem 1 Given the rotation matrix R, compute the rotation angle θ and the unit vector u. Solution The angle θ is computed according to (2.28) θ = ± arccos (tr (R) − 1 2 ) (2.43) The sign ambiguity in arccos cannot be avoided, being implicit in the inverse trigonometric formulas. The unit vector u could be obtained from the eigenvector relative to the unit eigenvalue of R: Ru = 1u (2.44) The rotation axis, without the positive direction information, can be computed also from any non zero column of the following matrix ( R + RT) − (tr (R) − 1) I . (2.45)

SLIDE 25

Basilio Bona - Dynamic Modelling 33 However this procedure does not allow to link the choice of the sign in ± arccos(·) with the positive verse of θ, because of the ambiguity that comes from the solution

f the eigenvector computation; indeed (2.44) allows to compute only the subspace

spanned by the eigenvector, but not its positive direction. We cannot therefore be a-priori sure that this eigenvector is coherent with the sign of θ, and we need to make an a-posteriori check. To solve this ambiguity it is better to follow a more complex method, that makes sure that the angle sign and the positive direction of u are coherent. The angle θ is computed following (2.43), with any choice of the sign; at this point

ne defines a symbolic skew-symmetric matrix S(u) according to (2.39), and nu-

merically compute S(u) subtracting (2.41) from (2.40): S(u) = ∥u∥ 2 sin θ ( R − RT) . (2.46) the components ui are obtained equating term by term the elements of the symbolic matrix (2.39) with those obtained from (2.46). As one can see, when R is symmetric, S(u) in (2.46) is zero. One can proceed considering two cases: Case 1. R = RT = I S(u) is not determined since sin θ = 0 ± 2kπ; by convention θ is set to zero and u is undetermined, since no rotation occurs. Case 2. R = RT ̸= I From the orthogonality properties of R follows that RTR = RR = R2 = I . This means that after two rotations of the angle θ around u the orientation is the starting one; this gives 2θ = 2π, therefore one sets θ = π. To compute u one builds S 2(u) as the sum of (2.40) and (2.41): S 2 = ∥u∥2 1 − cos θ (R + RT 2 − I ) (2.47) and then applies S 2 + ∥u∥2 I = uuT =   u2

u1u2 u1u3 u2u1 u2

u2u3 u3u1 u3u2 u2

  (2.48) in order to compute the components of u, equalling the numerical values of (2.47) with the corresponding symbolic terms in (2.48).

SLIDE 26

34 Basilio Bona - Dynamic Modelling Problem 2 Given the rotation axis represented by u =   u1 u2 u3  

and the angle θ, compute R. Solution One computes the skew-symmetric matrix S(u), from the known components of u, then computes the norm ∥u∥ and after thatthe matrix R according to (2.40). Now we present some examples to illustrate the various procedures.

Example 2.8.1

Let’s take Rb obtained from Example 2.7.1. Rb = R(k, 60◦)R(j , 45◦) =   0.3536 −0.8660 0.3536 0.6124 0.5000 0.6124 −0.7071 0.7071   (2.49) and compute the corresponding u and θ. We start computing θ from (2.43) θ = ± arccos (tr (R) − 1 2 ) = ± arccos 1.5607 − 1 2 = ±73.7◦ then we obtain u from (2.45), as A = ( R + RT) − (tr (R) − 1) I It results A =   0.1464 −0.2537 −0.3536 −0.2537 0.4393 0.6124 −0.3536 0.6124 0.8536   Once normalized, the three columns are respectively v 1 =   0.3190 −0.5525 −0.7701   v 2 = v 3 =   −0.3190 0.5525 0.7701   Now, in order to decide if the rotation axis is v 1 or v 2 one should do the a-posteriori check, computing both R(v 1, θ) = R(v 2, −θ) =   0.3536 0.6124 −0.7071 −0.8660 0.5000 0.3536 0.6124 0.7071  

SLIDE 27

Basilio Bona - Dynamic Modelling 35 and R(v 1, −θ) = R(v 2, θ) = R(v 1, θ)T = R(v 2, −θ)T We observe that Rb = R(v 1, −θ) = R(v 2, θ) hence we have at the end u = v 2 =   −0.3190 0.5525 0.7701   and θ = +73.7◦ Now we use the other approach using the skew-symmetric matrix from (2.46); we have already computed the angle θ = ±73.7◦ and we choose, for instance, the negative sign, so that θ = −73.7◦. If we now assume ∥u∥ = 1, we have the equality between the symbolic (2.39) and the numerical form (2.46) S(u) =   −u3 u2 u3 −u1 −u2 u1   =   0.7701 −0.5525 −0.7701 −0.3190 0.5525 0.3190   By inspection of the elements of S(u) we have u1 = 0.3190; u2 = −0.5525; u3 = −0.7701 but, considering that Rot(u, θ) = Rot(−u, −θ), we obtain the same results as the previous method, with no need to perform the sign check between θ and u. Example 2.8.2 We want to compute the rotation matrix R representing a rotation around the axis u = [

1 √ 3 1 √ 3 1 √ 3

f θ = 120◦

We compute the two matrices S e S 2 S = 1 √ 3   −1 1 1 −1 −1 1   ; S 2 = 1 3   −2 1 1 1 −2 1 1 1 −2   e successivamente calcoliamo R secondo l’equazione (2.40) ottenendo R =   1 1 1  

SLIDE 28

36 Basilio Bona - Dynamic Modelling Example 2.8.3 Given R =   1 1 1   (2.50) we want to compute u and θ. The eigenvector (with non-unit norm) relative to the unit eigenvalue of R is u = λ [ 1 1 1 ]T, as one can obtain from the equation   1 1 1     u1 u2 u3   =   u1 u2 u3   ⇒ u3 = u1 u1 = u2 u2 = u3 ⇒ u = λ   1 1 1   Using equation (2.45) we would have obtained the same results, since (R + RT) − (tr (R) − 1)I =   1 1 1 1 1 1 1 1 1   The angle θ is computed as ± arccos ( −1 2 ) = ±2π 3 = ±120◦. If we assume, for example, θ = − arccos ( − 1

) = −120◦, from (2.27) we would have

btained

R =   1 1 1   that is not equal to the initial matrix (2.50). So we use the method of the skew-symmetric matrix, after choosing θ1 = 120◦, sin θ1 = √ 3 2 . From (2.46) we compute S(u1) = 1 √ 3 ( R − RT) = 1 √ 3   −1 1 1 −1 −1 1  

btaining u1 = λ

[ 1 1 1 ]T, with λ = ∥u1∥ = 1 √ 3. Should we have assumed θ2 = −120◦, sin θ2 = − √ 3 2 , we would have obtained S(u2) = − 1 √ 3 ( R − RT) = − 1 √ 3   1 −1 −1 1 1 −1   and consequently u2 = λ [ −1 −1 −1 ]T = −u1.

SLIDE 29

Basilio Bona - Dynamic Modelling 37 Example 2.8.4 Given the matrix R =   −1 −1 −1   we want to compute u and θ; we observe that R is symmetric and R = RT ̸= I therefore we use relation (2.47) and (2.48). We know that in this case θ = π and cos θ = −1. Taking a unit norm vector u we compute ∥u∥2 1 − cos θ (R + RT 2 − I ) = 1 2   −1 −1 −2 −1 −1   ≡ S 2 now using (2.47) we obtain S 2 + ∥u∥2 I = uuT =   u2

u1u2 u1u3 u2u1 u2

u2u3 u3u1 u3u2 u2

  =   0.5 −0.5 −0.5 0.5   from which one obtains u = [√ 2 2 − √ 2 2 ]T

[ − √ 2 2 √ 2 2 ]T

Roto-translations

In the previous Sections we have detailed the representation of both translations and rotations; now we are ready to combine them in a single operator, the so called rototranslation or roto-translation operator. We take two reference frames, one identified as R0 the other identified as Rm, initially with common origins O and coinciding axes. For simplicity, R0 identifies the “fixed” frame, while Rm identifies the “mobile” frame; as already noted above, these two terms are used only for ease of explanation, in order to apply the roto- translation to Rm and characterize it with respect to R0. In addition we consider a geometrical point P, attached to Rm and represented by the vector v m

P ; we consider also a physical vector −

→ AB represented by v m

AB.

Now we want to apply a roto-translation to Rm, defined as follows

SLIDE 30

38 Basilio Bona - Dynamic Modelling

a rotation of an angle θ around the axis u going thought the common origin

O; the rotation is represented by the matrix Rot(u, θ) = R0

a translation of the origin of Rm with respect to the origin of R0; the trans-

lation is represented in R0 by the vector t0

If we want to compose rotations and translations (and rotations with rotations, as well) in a proper order, we should give additional information, namely

The order of the various displacements: what do you apply first, the translation
r the rotation?
The reference frame with respect to which the displacement is performed: you

do it with respect to the fixed frame or with respect to the mobile frame? Indeed, because of the previous displacements, the two reference frames do not coincide any more. According to the two decision we have four possible choices, considering that the first displacement is made indifferently with respect to the fixed or mobile frame, since they coincide. a) First perform a rotation of the initial reference frame around the common

rigin, then perform a translation with respect to the axes of the fixed frame

as in Figure 2.11 a), for a simple planar case. b) First perform a rotation of the initial reference frame around the common

rigin, then perform a translation with respect to the axes of the mobile

frame as in Figure 2.11 b). c) First perform a translation of the initial reference frame, then perform a rotation around the origin of the mobile frame as in Figure 2.11 c). d) First perform a translation of the initial reference frame, then perform a rotation with respect around the origin of the fixed frame as in Figure 2.11 d). As illustrated in Figure 2.11, the two cases a) and c) produce the same final roto- translation, while case b) and d) produce a different one. Now we consider the two vectors v m

P and v m AB; since they are “attached”, to the

mobile frame Rm, before the displacement their representation is v 0

P = v m P

and v 0

AB = v m AB

while after the displacement they will assume a different representation in R0, and precisely

SLIDE 31

Basilio Bona - Dynamic Modelling 39 Figure 2.11: Diversi modi di effettuare una rototraslazione planare.

SLIDE 32

40 Basilio Bona - Dynamic Modelling

In the two cases a) and c) we have

v 0

P = R0 mv m P

Rot +t

Transl

(2.51) and v 0

AB = R0 mv m AB

where R0

m is the rotation matrix between R0 e Rm (i.e., Rm represented in

R0) and t is the translation vector from R0 to Rm, represented in R0.

In the two cases b) and d) we have

v 0

P = R0 m (v m P + t′)

Trasl
Rot

= R0

mv m P + R0 mt′

(2.52) and v 0

AB = R0 mv m AB

but now t′ is the translation vector from R0 to Rm, represented in Rm, that becomes R0

mt′ when represented in R0.

First you notice that the “physical” vector v m

AB does not translate, but only change

its representation and becomes R0

mv m AB due to the rotation of the reference frame;

this is correct, since a change of reference frames does not change a force, a torque, a velocity, etc. except for their representations. Second, you see that the results are different; therefore it is necessary to solve any ambiguity that may arise. We will see in Section 2.10.2 a simple mnemonic rule to perform roto-translations in correct order.

Esempio_4-12

Inserire esempio di rappresentazione vettoriale da Matlab see Figure xxx.

2.10

Homogeneous Coordinates

We have seen in the previous Section that a rotation operator corresponds to a matrix product, while a translation operator corresponds to a vector sum. It is possible to use a unique operator for translations and rotations if we introduce the COMPLETARE

SLIDE 33

Basilio Bona - Dynamic Modelling 41 Figure 2.12: An example of vector transformations between two rotated frames. so-called homogeneous coordinates or, with a less common term, perspective coordinates. They come from the projective and epipolar geometry and find their use in many different contexts as computer graphics [41] and 3D vision systems [47]. The interested reader can find a detailed description of projective/epipolar geometry and its applications to computer vision in [10, 11, 19, 37, 44]. Given a point P in a 3D space, its homogeneous representation is given by the associated 4 × 1 homogeneous vector x, defined as

def

=     λp1 λp2 λp3 λ     = λ     p1 p2 p3 1     (2.53) where λ ∈ R is a scale factor; in the following this factor will be set to 1, giving

rigin to the homogeneous coordinates representation commonly adopted to study

the rigid roto-translations in the 3D space,

def

=     p1 p2 p3 1     (2.54) Now, considering again Eqn. (2.51), written in a simpler form v 0 = R0

mv m + t

we notice that it can be written as ˜ v 0 = T 0

m˜

v m = [ R0

t 0T 1 ] [ v m 1 ] (2.55) FARE FIGURA

SLIDE 34

42 Basilio Bona - Dynamic Modelling As one can see from the previous relation, a roto-translation, i.e., a rotation plus a translation is represented by a 4 × 4 homogeneous matrix (HM) T 0

m def

= [ R0

t 0T 1 ] (2.56) where R0

m and t are the rotation and the translation, respectively, and

0T def = [ ]

2.10.1 Homogeneous Transformations

We can define two basic homogeneous transformations, namely, the pure rotation T(R) ≡ T R

def

= [ R 0T 1 ] (2.57) and the pure translation T(t) ≡ T t

def

= [ I t 0T 1 ] (2.58) The generic homogeneous transformation in (2.56) can be obtained as the product

f a pure translation and a pure rotations as in

T = T tT R = [ I t 0T 1 ] [R 0T 1 ] = [R t 0T 1 ] (2.59) Reversing the factors order one obtains a different HM T ′ = T RT t = [R 0T 1 ] [ I t 0T 1 ] = [R Rt 0T 1 ] (2.60) We have already given in (2.55) the rule for transforming the representation of a geometrical point from Rm to R0 using the homogeneous coordinates; in brief:

1. We build the HM T 0

m from R0 m and t0 m in the right order: (2.59) or (2.60);

2. We write v m

P in homogeneous coordinates

v m

P as in (2.54);

3. We compute

v 0

P as

P = T 0 m

v m

4. We transform back from homogeneous coordinates to vector form

v 0

P in v 0 P.

SLIDE 35

Basilio Bona - Dynamic Modelling 43 It is important to notice that homogeneous vectors do no follow the classical sum

r difference rules, since
x +

y ̸=     x1 + y1 x2 + y2 x3 + y3 2     For this reason we introduce the homogeneous sum operator, denoted by ⊕, and the homogeneous difference operator, denoted by ⊖, respectively defined as:

x ⊕

y = [x + y 1 ]

x ⊖

y = [x − y 1 ] (2.61) Moreover these operators are not distributive, since: T( x ⊕ y) ̸= T x ⊕ T y T( x ⊖ y) ̸= T x ⊖ T y (2.62) Relations (2.62) are important only for their use on geometrical vectors, as v P; indeed when we have to do with oriented segments (physical vectors) v m

AB = v m B−v m A,

ne can use the homogeneous difference, as in
v 0

AB = T

v m

B ⊖ T

v m

A = Rv m B + t − Rv m A − t = R(v m B − v m A)

(2.63)

r the usual difference, as in
v 0

AB = T(

v m

B −

v m

A) = T

[v m

B − v m A

] = R(v m

B − v m A)

Inverse Homogeneous Matrix The generic HM T consists of two factors: a 3 × 3 submatrix R (the rotation

perator) and a 3 × 1 vector t (the translation operator).

For this reason the inverse matrix is not simply its transpose, but (see also Figure 2.13): [T]−1 = [ RT −RTt 0T 1 ] (2.64) It follows that, while the inverse rotation matrix is its transform, the inverse translation is t′ = −RTt; this has an immediate geometrical meaning: the minus sign changes the orientation of the transaction vector and RTt is the representation in Rm of the original t, that was represented in R0).

Example 2.10.1

SLIDE 36

44 Basilio Bona - Dynamic Modelling Figure 2.13: La trasformazione omogenea diretta e inversa tra due sistemi di riferi- mento. Given the physical vector v 0

AB = v 0 B −v 0 A, represented in R0, with v 0 A =

[ 1 ]T and v 0

B =

[ 1 ]T, we want to compute its representation v m

AB = v m A − v m B in

Rm, where T 0

m =

    −1 1 1 1 1 1 1     . We have to find the inverse of T 0

m as

T m

0 =

( T 0

)−1 =     1 −1 −1 1 1 1 1     ; Now using (2.63) we have:

AB =

    1 −1 −1 1 1 1 1         1 1     ⊖     1 −1 −1 1 1 1 1         1 1     =     −1 1 2 1     ⊖     −1 1 1     =     1 1 1     hence v 0

AB =

  1 1  .

SLIDE 37

Basilio Bona - Dynamic Modelling 45 The same results can be obtained computing directly the rotation of the oriented segment: v 0

AB = Rm

( v 0

B − v 0 A

) =   1 −1 1     −1 1   =   1 1  

2.10.2

Composition Rule for Roto-Translations

It is possible to perform any number of roto-translations obtaining a final roto- translation as the ordered product of the single displacements; it is therefore necessary to organize the product terms in the correct sequence. We assume to start with two reference frames, for simplicity one called “fixed”, the other “mobile”, that are related by a homogeneous transformation T. The we want to execute n generic displacements, each one represented by the related matrix T i, i = 1, . . . , n. To perform in the correct order the various matrix products we must apply the following rules. We call T(i) the matrix product obtained after the first i-th displacements

1. Set i = 0 and initialize T(0) = T.
2. If the i-th roto-translation T i is defined with respect to the fixed reference

frame, one should pre-multiply the previous matrix T(i − 1) by T i; the result will be T(i) = T iT(i − 1) Mnemonic Rule: pre-fix(ed)

3. If the i-th roto-translation T i is defined with respect to the mobile reference

frame, one should post-multiply the previous matrix T(i − 1) by T i; the result will be T(i) = T(i − 1)T i Mnemonic Rule: post-mob(ile) These rules are valid also for the product composition of rotation matrices. For the composition of translations alone, since the sum is commutative, the order is not important. Observing again Figure 2.11 we can interpret the results as follows

SLIDE 38

46 Basilio Bona - Dynamic Modelling

Figure 2.11 a)

It represents a rotation around the common origin of the two reference frames, followed by a translation with respect to the fixed frame T(0) = I , T(1) = T R, T(2) =

T tT R
Figure 2.11 b)

It represents a rotation around the common origin of the two reference frames, followed by a translation with respect to the mobile frame T(0) = I , T(1) = T R, T(2) =

T RT t
Figure 2.11 c)

It represents a translation with respect to the original frame followed by a rotation around the origin of the mobile frame T(0) = I , T(1) = T t, T(2) =

T tT R
Figure 2.11 d)

It represents a translation with respect to the original frame followed by a rotation around the origin of the fixed frame; T(0) = I , T(1) = T t, T(2) =

T RT t

It is evident that we can interpret the final result T(2) in two different ways, that differs only by the verbal description associated with their description. For instance, considering case a) and c), i.e., T(2) = T tT R we can give two different linguistic description, that nevertheless produce the same final result

a rotation with respect to the original frame followed by a translation with

respect to the fixed frame (rule pre-fix)

a translation with respect to the original frame followed by a rotation with

respect to the mobile frame (rule post-mob) In conclusion, it is worth noticing again that a HM as T 0

m represnts a “mobile” frame

Rm with respect to a “fixed” frame R0; the 3 × 3 rotation submatrix represents the mobile reference frame with respect to the fixed one, while the last column part provides the translation of the origin. All components are expressed in R0.

SLIDE 39

Basilio Bona - Dynamic Modelling 47 Example 2.10.2 We want to compute the HM T 0

m representing the frame Rm obtained rotating the

fixed frame around the basis unit vector i of an angle +90◦, followed by a translation t = [ 2 ] along the y axis of the resulting (mobile) frame, followed by a rotation around the axis z of the fixed reference frame of an angle −90◦. We have three displacements defined by the following HM: First displacement = rotation: T 1 = T (R(i, 90◦)) = [R(i, 90◦) 0T 1 ] =       1 −1 1   0T 1     Second displacement = translation: T 2 = T (t) =     I   2   0T 1     Third displacement = rotation: T 3 = T (R(k, −90◦)) = [R(k, −90◦) 0T 1 ] =       1 −1 1   0T 1     The composition rule is: first T 1, then post-product by T 2, then pre-product by T 3: T 0

= T 3T 1T 2 =     1 −1 1 1         1 −1 1 1         1 1 2 1 1     = =     1 −1 1 1         1 −1 1 2 1     = =     −1 −1 1 2 1    

SLIDE 40

48 Basilio Bona - Dynamic Modelling Example 2.10.3 A geometrical point P is represented in Rm by v m = [ 4 3 2 ]T; Rm is translated by t = [ −3 7 ]T and then rotated according to R =   −1 1 −1   with respect to the new frame. Find the components of P in R0. T 0

m is obtained post-multiplying the translation operator by the rotation operator,

as T 0

m =

    I   −3 7   0T 1    

translation

      −1 −1 −1   0T 1    

rotation

=       −1 −1 −1     −3 7   0T 1     At this point one computes the homogeneous representation of the transformed vector as

v m =

      −1 −1 −1     −3 7   0T 1         4 3 2 1     =     −6 −4 5 1    

btaining

v m =   −6 −4 5  

2.11

Rigid Body Representation

A rigid body B in a 3D space can be simply represented by a right hand reference frame RB associated to it; it is not strictly necessary that the frame lies inside the body or on its surface, but it must be rigidly joined to it. We call pose of a rigid body the set parameters that uniquely define its position and its orientation in R3. As we will detail in a successive Section, the pose of a rigid body, not subject to any constraint in R3 is given by three position parameters and three orientation parameters. In R2 the pose is given by two position parameters and one orientation parameter.

SLIDE 41

Basilio Bona - Dynamic Modelling 49 The pose is computed as follows: once the frame RB is given, its pose with respect to an external reference frame, often called world or (pseudo-)inertial reference frame R0 can be obtained from the homogeneous transformation T 0

B (see Figure

2.14). Figure 2.14: A rigid body pose is characterized by the relation between the body frame RB and the world frame R0. There are sixteen elements in T 0

B =

[ R0

t 0T 1 ] , but only six are independent, in particular the three elements of the translation vector and only three elements of the rotation matrix (the rotation matrix is orthonormal, so it is subject to six constraints). We can now state that the pose of a free rigid body in space is defined by six degrees

f freedom (dof), three related to the position of the origin reference frame that

describes the body, and three related to it the orientation with respect to a world reference frame. The motion of a rigid body is described by the time equations of its pose; the study of the pose and its time history defines what we call the body

kinematics. In this Section we will only detail the different ways in which the pose

is described, postponing to successive Sections the characterization of the time laws involved. Usually the three position dof’s are chosen as the cartesian coordinates of the origin, although sometimes it is useful to adopt spherical coordinates (but we will never use them in these notes). On the contrary, the orientation coordinates are more complex to choose: we will devote the entire Section 2.12 to describe the various possible parameterizations of orientation; for the moment we will call refer to these

SLIDE 42

50 Basilio Bona - Dynamic Modelling three numbers as the “vector” α = [ α1 α2 α3 ]T. The body pose is therefore formally defined as p(t)

def

=          p1(t) p2(t) p3(t) p4(t) p5(t) p6(t)          =          x1(t) x2(t) x3(t) α1(t) α2(t) α3(t)          = [x(t) α(t) ] (2.65) The first three components arranged in the vector x, that is a true vector, i.e., a geometrical vector, that can be subject to all the vector operations, while the second three components that describe the orientation, arranged in α do not formally obey to the vector rules, since the “sum” of two orientations α and β is obtained by the product of two rotation matrices and not by the sum [ α1 + β1 α2 + β2 α3 + β3 ]T α → R(α) β → R(β) } ⇒ R(α)R(β) = R(γ) ⇒ γ So it is not possible to consider p a vector, and the operations as sum and scalar product have no sense. This fact come from the formal definition of the quantities involved: indeed p ∈ R3 × SO(3), with x ∈ R3, α ∈ SO(3) where R3 is the real vector space of dimension 3, and SO(3) is the special orthonormal group of dimension 3, also known as the rotation group; for further details on this group SO(3), see [34]. Nevertheless the habit to call p a vector is common, and we will continue to refer to p or α with the term “vector”, since they are organized as a column of real numbers. With this in mind, we can state that two are the possible representations of a rigid body in space

1. The vector p(t) defined in (2.65).
2. The homogeneous matrix T(t), defined as:

T(t)

def

= [R(t) t(t) 0T 1 ] (2.66) These two representations are equivalent and it is possible to compute one from the

ther, as we will show in the next Section.

SLIDE 43

Basilio Bona - Dynamic Modelling 51 Figure 2.15: Dato p calcolare T.

2.11.1 From Pose to Homogeneous matrix

Given p, we have to compute T, i.e., R and t (see Figure 2.15). The vector t is simply obtained taking the first three elements of p, i.e., x t ⇐   p1 p2 p3   ≡ x (2.67) We note that this is true only if (as we do) we give the position in cartesian coordinates and not in cylindrical or spherical coordinates; if this is the case it results in general that t = f (x): as an example, from spherical coordinates x sph = [ ρ θ ϕ ]T, to cartesian coordinates x car = [ xc yc xc ]T, we have xc = t1 = ρ sin θ cos ϕ = f1(x sph) yc = t2 = ρ sin θ sin ϕ = f2(x sph) zc = t3 = ρ cos θ = f3(x sph) To compute R from α it is necessary to define the physical meaning of the three components αi of α. We assume that they are three angles characterizing the

rientation of RB in R0. The most used angles are the so called Euler angles or

the RPY angles, that will be introduced in Section 2.12.2 and 2.12.3, respectively. It exists a nonlinear function R = g E(αE) that gives the elements rij of a rotation matrix R from the Euler angles (2.69), and another nonlinear function R = g RPY (αRPY ) that gives the elements rij of a rotation matrix R from the RPY angles (2.76).

2.11.2 From Homogeneous matrix to Pose

Given T, we have to compute p (see Figure 2.16).

SLIDE 44

52 Basilio Bona - Dynamic Modelling The vector x is simply obtained taking the first three elements of the last column

x ≡   p1 p2 p3   ⇐   t1 t2 t3   ≡ t (2.68) Figure 2.16: Dato T calcolare p. Notice that in this way the obtained coordinates are relative to the origin of the reference frame described by T; if we need to compute the coordinates of a different geometrical point, we should characterize this new point as a vector in RB and transform it using eqn. (2.55). To compute α from R it is necessary to define the physical meaning of the three components αi of α. We assume to use the Euler angles or the RPY angles. It exists the inverse nonlinear function αE = g −1

E (R) that gives the Euler angles

from the elements rij of a rotation matrix R (2.73), and another inverse nonlinear function αRPY = g −1

RPY (R) that gives the RPY angles from the elements rij of a

rotation matrix R (2.77). The inverse solution, if exists, could be non unique, since a finite number of angles α may exist that produce the same rotation matrix R. Since inverse trigonometric functions are involved, we consider as equaivlent solutions those that differ for integer multiples of 2π. In the following Section we will illustrate the most used ways to characterize the

rientation of a rigid body in space.

2.12 Orientation Parameters

We have seen that a rotation is a geometrical transformation acting on R3, defined by a orthonormal matrix R. At the same time the orientation of a rigid body, with a local reference frame attached to it, is represented by the rotation of the local frame with respect to the world or global frame. We conclude that it is equivalent to speak of the orientation of a rigid body or of its reference frame with respect to

SLIDE 45

Basilio Bona - Dynamic Modelling 53 some other “fixed” frame Body B Orientation ⇔ Rotation Matrix RB We recall that a rotation is characterized by only three parameters, and all the representation described in this section will be a different way to “organize” such parameters; in eqn. (2.65) these three parameters were indicated as α = [ α1 α2 α3 ]T, but there are other parameterizations that use more than three parameters, as the quaternions, introduced in Section 2.12.6. Another example comes from the Euler theorem, that states that any composition of rotations is always a rotation of an angle θ around an axis defined by a unit vector u; therefore we can represent a rotation by the set (u, θ), where u has two free parameters, since there is an implicit constraint ∥u∥ = 1), and the third one coincide with θ. Other authors use a non unit vector v, whose norm provides the angle value ∥v∥ = θ. Now we describe the most common forms used to represent the orientation of a rigid body in R3.

2.12.1 Direction Cosines

Direction cosines are nothing else that the rotation matrix R itself. Indeed R contains in its columns (or rows) the representation of the unit basis vectors of the local frame with respect to the world frame. This representation needs nine parameters that must obey to six unit norm orthogonality constraints. In this case the parameters [ α1 α2 α3 ]T are hidden in R but can be extracted using the equations (2.12.2) e (2.12.3).

2.12.2 Euler Angles

Historically this is the first representation of the orientation of a body: it associates to α1, α2, α3 three angles, called Euler angles and usually denoted by ϕ, θ, ψ. To understand how to build the Euler angles it is necessary to define them through an implicit procedure. The orientation of a mobile reference frame Rm described by the Euler angles is obtained by three successive rotations around the principal axes, following a precise rule First rotation, angle ϕ Rz,φ ≡ R(k, ϕ) =   cφ −sφ sφ cφ 1   rotation of an angle ϕ around the local (mobile) axis z.

SLIDE 46

54 Basilio Bona - Dynamic Modelling Second rotation, angle θ Rx,θ ≡ R(i, θ) =   1 cθ −sθ sθ cθ   rotation of an angle θ around the local (mobile) axis x. Third rotation, angle ψ Rz,ψ ≡ R(k, ψ) =   cψ −sψ sψ cψ 1   rotation of an angle ψ around the local (mobile) axis z. Using the rule “pre–fixed” “post–mobile” we compute the complete rotation matrix based on the Euler angles as R (ϕ, θ, ψ) ≡ Rz,φRx,θRz,ψ ≡ R(k, ϕ)R(i, θ)R(k, ψ) =   cφcψ − sφcθsψ −cφsψ − sφcθcψ sφsθ sφcψ + cφcθsψ −sφsψ + cφcθcψ −cφsθ sθsψ sθcψ cθ   (2.69) This composition rule is not unique in technical literature, since many textbooks adopt a slightly different convention: the second rotation around the x axis is re- placed by a rotation around the mobile axis y, producing a different “Euler” matrix:

R (ϕ, θ, ψ) ≡ Rz,φRy,θRz,ψ ≡ R(k, ϕ)R(j , θ)R(k, ψ) =

  −sφsψ + cφcθcψ −sφcψ − cφcθsψ cφsθ cφsψ + sφcθcψ cφcψ − sφcθsψ sφsθ −sθcψ sθsψ cθ   (2.70) In this notes we will always adopt the first form (2.69), but we recall that in many aerospace engineering textbooks the form (2.70) is widely used. The above formulas are also called direct relations, since, given the three Euler angles they compute the rotation matrix R(ϕ, θ, ψ). The inverse relation, i.e. how to compute the Euler angles given an generic R matrix is solved considering the generic elements of a matrix, provided it is orthonormal R =   r11 r12 r13 r21 r22 r23 r31 r32 r33   , (2.71) where the elements rij are known. To obtain the Euler angles it is necessary to solve

SLIDE 47

Basilio Bona - Dynamic Modelling 55 the following nonlinear equation system: r11 = cφcψ − sφcθsψ r12 = −cφsψ − sφcθcψ r13 = sφsθ r21 = sφcψ + cφcθsψ r22 = −sφsψ + cφcθcψ r23 = −cφsθ r31 = sθsψ r32 = sθcψ r33 = cθ (2.72) that has the following generic solution θ = ± arccos (r33) ± 2kπ ψ = ± arccos (r32 sθ ) ± 2kπ ϕ = ± arccos (−r23 sθ ) ± 2kπ (2.73) Unfortunately this solution presents some drawbacks

1. the inverse trigonometric function arccos(·) is not unique: indeed cos(θ) =

cos(−θ);

2. the solution becomes not definite for r33 = 1, i.e., when sθ = 0; in this case

the angles ϕ and ψ are not uniquely known; only their sum is given;

3. when θ → 0◦ or θ → ±180◦, the second and the third equations provide inac-

curate solutions, since the numeric accuracy of the arccos(·) function depends

n the angle value.

To solve these drawbacks, instead of arccos(·) it is customary to use the function atan2(y, x), that is available in all the mathematical libraries of the most used computer languages. It is definite as: θ = atan2(y, x) = tan−1 (y x ) = =        0◦ ≤ θ ≤ 90◦ if x ≥ 0; y ≥ 0 90◦ ≤ θ ≤ 180◦ if x ≤ 0; y ≥ 0 −180◦ ≤ θ ≤ −90◦ if x ≤ 0; y ≤ 0 −90◦ ≤ θ ≤ 0◦ if x ≥ 0; y ≤ 0 (2.74) Moreover, by default, atan2(0, 0) = 0.

SLIDE 48

56 Basilio Bona - Dynamic Modelling Using this function, the solution to (2.72) is: ϕ = atan2 (r13, −r23) ± 2kπ ψ = atan2 (−cφr12 − sφr22, cφr11 + sφr21) ± 2kπ θ = atan2 (sφr13 − cφr23, r33) ± 2kπ (2.75) Euler angles singularity We observe that when r33 = 1 the matrix R(ϕ, θ, ψ) is equal to the elementary matrix R(k, γ), that is function of a single angle. In this case we say that the Euler representation is singular; from the three possible angles we can obtain only two angles; from (2.73) we have θ = 0, and the product R(k, ϕ)R(i, θ)R(k, ψ) reduces to R(k, ϕ)R(k, ψ) = R(k, (ϕ + ψ)) from which we have γ = (ϕ + ψ); we cannot compute separately the two angles phi and ψ, but only their sum. This situation is described saying that θ does not decouple any more the other two rotations and so a singular configuration is produced.

2.12.3 RPY Angles

Also in the case of Roll-Pitch-Yaw angles (RPY for short) θx, θy, θz it is necessary to define them through an implicit procedure. The orientation of a mobile reference frame Rm described by the RPY angles is obtained by three successive rotations around the principal axes, following a precise rule First rotation, angle θx Rx,θx ≡ R(i, θx) =   1 cθx −sθx sθx cθx   rotation of an angle θx around the world (fixed) axis x. Second rotation, angle θy Ry,θy ≡ R(j , θy) =   cθy sθy 1 −sθy cθy   rotation of an angle θy around the world (fixed) axis y. Third rotation, angle θz Rz,θz ≡ R(k, θz) =   cθz −sθz sθz cθz 1  

SLIDE 49

Basilio Bona - Dynamic Modelling 57 rotation of an angle θz around the world (fixed) axis z. Using the rule “pre–fixed” “post–mobile” we compute the complete rotation matrix based on the RPY angles as R (θx, θy, θz) ≡ Rz,θzRy,θyRx,θx ≡ R(k, θz)R(j , θy)R(i, θx)

def

=   cθzcθy sθxsθycθz − cθxsθz cθxsθycθz + sθxsθz cθysθz sθxsθysθz + cθxcθz cθxsθysθz − sθxcθz −sθy sθxcθy cθxcθy   (2.76) The RPY angles are computed applying the same approach adopted for the Euler angles, yielding: θx = atan2 (r32, r33) ± 2kπ θz = atan2 (−cθxr12 + sθxr13, cθxr22 − sθxr23) ± 2kπ θy = atan2 (−r31, sθxr32 + cθxr33) ± 2kπ (2.77) It is interesting to notice that the product R(k, θz)R(j , θy)R(i, θx) may also be read in a different order, starting from left to right, i.e., applying the “post-mobile” rule: first apply a rotation R(k, θz) around axis z, then apply a rotation R(j , θy) around mobile axis y, then apply a rotation R(i, θx) around mobile axis x. The difference is only in the description of the rotations, not in the final result (2.76), that is the same. Other definitions can be found in various textbooks and are sometimes used, for instance, the following alternative sequence is quite common 1) A rotation R(j , θy) of θy around the fixed axis y; 2) A rotation R(k, θz) of θz around the fixed axis z; 3) A rotation R(i, θx) of θx around the fixed axis x. The resulting rotation matrix is: R (θz, θy, θx) ≡ Rx,θxRy,θyRz,θz ≡ R(i, θx)R(j , θy)R(k, θz)

def

=   cθycθz −cθysθz sθy sθxsθycθz + cθxsθz −sθxsθysθz + cθxcθz −sθxcθy −cθxsθycθz + sθxsθz cθxsθysθz + sθxcθz cθxcθy   (2.78) In this notes we will apply only definition (2.76). RPY angles singularity Roll-Pitch-Yaw angles are subject to singularity too, as the Euler angles, and in general to all three-angles representation.

SLIDE 50

58 Basilio Bona - Dynamic Modelling In particular the following relations hold R(k, θz)R(j , 90◦) = R(j , 90◦)R(i, θz) (2.79) R(i, θx)R(j , 90◦) = R(j , 90◦)R(k, θx) (2.80) When θy = 90◦ the RPY matrix in (2.76) becomes singular, since from (2.79) R(k, θz)R(j , 90◦)R(i, θx) = R(j , 90◦)R(i, θz)R(i, θx) = R(j , 90◦)R(i, (θx + θz)) and also the RPY matrix in (2.78) becomes singular, since from (2.80) R(i, θx)R(j , 90◦)R(k, θz) = R(j , 90◦)R(k, θx)R(k, θz) = R(j , 90◦)R(k, (θx + θz)) Both representation are no more function of three angles, but only of a combination

f two of them; this fact is the cause of the so called gimbal-lockproblem that
ccurs in gyroscopes and is well known since the incident on the Apollo 10 Manned

Lunar Spacecraft [22].

2.12.4 Cardan Angles

We have seen in the previous Section that many alternatives are possible in order to define three rotations. All possible angles are called generically Cardan angles. Two possible groups of rotations arise: the first group includes the product obtained by the product of three elementary rotation (all different); we call θ1, θ2 and θ3 the three generic Cardan angles involved; Table 2.12.4 list all possible combinations. The second group includes the product of three rotation matrices, but introducing

nly two angles; in this case the Cardan angles are called α, β and γ. Table 2.12.4

list all possible combinations. Cardan angles singularities As already noted in Section 2.12.2 and 2.12.3 the sequence of three angles produces a singularity when a certain pattern of angles appear. For the first group it is the presence of a 90◦ angles that may produce a singular behaviour, since the following identities hold R(i, θx)R(j , 90◦) = R(j , 90◦)R(k, θx) R(i, θx)R(k, 90◦) = R(k, 90◦)R(j , θx) R(j , θy)R(i, 90◦) = R(i, 90◦)R(k, θy) R(j , θy)R(k, 90◦) = R(j , 90◦)R(i, θy) R(k, θz)R(i, 90◦) = R(i, 90◦)R(j , θz) R(k, θz)R(j , 90◦) = R(j , 90◦)R(i, θz) (2.81) while for the second group it is sufficient that the middle rotation matrix is equal to the identity to give origin of a singular representation.

SLIDE 51

Basilio Bona - Dynamic Modelling 59 R(i, θ1)R(j , θ2)R(k, θ3)

Eqn. (2.78)

R(i, θ1)R(k, θ3)R(j , θ2)   cθ1cθ3 −sθ3 sθ2cθ3 cθ1cθ2sθ3 + sθ1sθ2 cθ1cθ3 cθ1sθ2sθ3 − sθ1cθ2 sθ1cθ2sθ3 − cθ1sθ2 sθ1cθ3 sθ1sθ2sθ3 + cθ1cθ2   R(j , θ2)R(i, θ1)R(k, θ3)   sθ1sθ2sθ3 + cθ2cθ3 sθ1sθ2cθ3 − cθ2sθ3 cθ1sθ2 cθ1sθ3 cθ1sθ3 −sθ1 sθ1cθ2sθ3 − sθ2cθ3 sθ1cθ2cθ3 + sθ2sθ3 cθ1cθ2   R(j , θ2)R(k, θ3)R(i, θ1)   cθ2cθ3 −cθ1cθ2sθ3 + sθ1sθ2 sθ1cθ2sθ3 + cθ1sθ2 sθ1 cθ1cθ3 −sθ1cθ3 −sθ2cθ3 cθ1sθ2cθ3 + sθ1cθ2 −sθ1sθ2sθ3 + cθ1cθ2   R(k, θ3)R(i, θ1)R(j , θ2)   −sθ1sθ2sθ3 + cθ2cθ3 −cθ1sθ3 sθ1cθ2sθ3 + sθ2cθ3 sθ1sθ2cθ3 + cθ2sθ3 cθ1cθ3 −sθ1cθ2cθ3 + sθ2sθ3 −cθ1sθ2 sθ1 cθ1cθ2   R(k, θ3)R(j , θ2)R(i, θ1)

Eqn. (2.76)

Table 2.1: Cardan angles obtained by three elementary rotations around three different angles.

2.12.5 Euler Parameters

Given the rotation described in (2.27) R(u, θ) of an angle θ around an axis represented by the unit vector u = [ u1 u2 u3 ]T, we introduce four parameters, called Euler parameters vi (not to be confused with the Euler angles), defined as: v1 = u1 sin θ 2, v2 = u2 sin θ 2, v3 = u3 sin θ 2, v4 = cos θ 2 (2.82) Only three out of four of these parameters are independent, since the following constraint holds

∑

i=1

i = 1

(2.83) The Euler parameters may be interpreted as the components of a unit quaternion

u. The quaternions will be introduced in Section 2.12.6.

SLIDE 52

60 Basilio Bona - Dynamic Modelling R(i, α)R(j , β)R(i, γ)   cβ sβsγ sβcγ sαsβ −sαcβsγ + cαcγ −sαcβcγ − cαsγ −cαsβ cαcβsγ + sαcγ cαcβcγ − sαsγ   R(i, α)R(k, β)R(i, γ)   cβ −sβcγ sβsγ cαsβ cαcβcγ − sαsγ −cαsβsγ − sαcγ sαsβ sαcβcγ + cαsγ −sαcβsγ + cαcγ   R(j , α)R(i, β)R(j , γ)   −sαcβsγ + cαcγ sαsβ sαcβcγ + cαsγ sβsγ cβ −sβcγ −cαcβsγ − sαcγ cαsβ cαcβcγ − sαsγ   R(j , α)R(k, β)R(j , γ)   cαcβcγ − sαsγ −cαsβ cαcβsγ + sαcγ sβcγ cβ sβsγ −sαcβcγ − cαsγ sαsβ −sαcβsγ + cαcγ   R(k, α)R(i, β)R(k, γ)

Eqn. (2.69)

R(k, α)R(j , β)R(k, γ)

Eqn. (2.70)

Table 2.2: Cardan angles obtained by two elementary rotations around three different angles. Given u, θ and the Euler parameters vi, the rotation matrix R(u, θ) can be computed as: R(u, θ) =    v2

1 − v2 2 − v2 3 + v2 4

2(v1v2 − v3v4) 2(v1v3 + v2v4) 2(v1v2 + v3v4) −v2

1 + v2 2 − v2 3 + v2 4

2(v2v3 − v1v4) 2(v1v3 − v2v4) 2(v2v3 + v1v4) −v2

1 − v2 2 + v2 3 + v2 4

   (2.84) Conversely, given R(u, θ), the Euler parameters are computed as: v4 = ±1 2 √ (1 + r11 + r22 + r33) v1 = 1 4v4 (r32 − r23) v2 = 1 4v4 (r13 − r31) v3 = 1 4v4 (r21 − r12) (2.85)

SLIDE 53

Basilio Bona - Dynamic Modelling 61 The sign ambiguity in v4 can be eliminated constraining the angle as −π 2 ≤ θ 2 ≤ π 2 ,

r −π ≤ θ ≤ π; in this way v4 can be only positive.

The rotation angle can be computed from cos θ = v2

4 − (v2 1 + v2 2 + v2 3)

and the unit vector u as u = 1 sin(θ/2)   v1 v2 v3   The Euler parameters can also be computed directly from the Euler angles (ϕ, θ, ψ) as: v1 = sin (ϕ − ψ 2 ) sin (θ 2 ) v2 = cos (ϕ − ψ 2 ) sin (θ 2 ) v3 = sin (ϕ + ψ 2 ) cos (θ 2 ) v4 = cos (ϕ + ψ 2 ) cos (θ 2 ) (2.86) Often the Euler parameters are denoted by the symbol v =     v1 v2 v3 v4     that is similar to a vector, but it is not, since the application of vector operators has no sense. The Euler parameter are a most convenient form for parameterizing the rotations: they are more compact than the rotation matrix R and more efficient from an algorithmic point of view, since to compute R = R(v) from them, as in (2.84), does not require the use of trigonometric functions. Moreover, given two rotations R(v a) and R(v b), with their relative Euler parameters v a and v b, one can compute directly the rotation product R(v c) = R(v a)R(v b) using the following matrix product: v c = F(v a)v b =     va4 −va3 va2 va1 va3 va4 −va1 va2 −va2 va1 va4 va3 −va1 −va2 −va3 va4     v b (2.87)

SLIDE 54

62 Basilio Bona - Dynamic Modelling Only 16 multiplications are required compared to the 27 that are involved in the product between two 3 × 3 rotation matrices.

2.12.6 Quaternions

Quaternions were introduced by Hamilton who discovered them in 1843, in order to generalize in three dimensions the complex numbers and their characteristic of being plane rotators. For an detailed history of the problems and discussion raised by the quaternions in the scientific community in the middle of XIX century, read the interesting monograph by Crowe [7]. The quaternion algebra, with its definitions and operators used for representing rotations, is described in details in Appendix D. The symbol for the generic quaternion is h.

2.12.7 Quaternions and Rotations

In order to use a quaternion u = (u0, u1, u2u3) to represent a rotation, we assign to each element ui, an Euler parameter vi u0 = v4 u1 = v1 u2 = v2 u3 = v3 (2.88) The quaternion so defined is therefore u = (u0, u1, u2, u3) = (cos θ 2, u1 sin θ 2, u2 sin θ 2, u3 sin θ 2) and has unit norm. Since a biunivocal correspondence exists between the Euler parameters and the unit quaternions, any unit quaternion represents a rotation in the 3D space, as any unit complex number represents a rotation in the 2D plane. We write R(u) to indicate that for every rotation there is a corresponding a unit quaternion u and vice versa. In order to convert a unit quaternion u = (u0, uv) into the corresponding matrix R(u), we use the following relation, that is equal to (2.84), valid for Euler parameters: R(u) = (u2

0 − uT v uv)I + 2uvuT v − 2u0S(uv) =

   u2

0 + u2 1 − u2 2 − u2 3

2(u1u2 − u3u0) 2(u1u3 + u2u0) 2(u1u2 + u3u0) u2

0 − u2 1 + u2 2 − u2 3

2(u2u3 − u1u0) 2(u1u3 − u2u0) 2(u2u3 + u1u0) u2

0 − u2 1 − u2 2 + u2 3

   (2.89)

SLIDE 55

Basilio Bona - Dynamic Modelling 63 Conversely, to compute the quaternion u starting from the elements rij of the rotation matrix R(h)we use a relatione that is equal to (2.85), valid for Euler parameters: u0 = ±1 2 √ (1 + r11 + r22 + r33) u1 = 1 4u0 (r32 − r23) u2 = 1 4u0 (r13 − r31) u3 = 1 4u0 (r21 − r12) (2.90) When u0 = 0, i.e., θ/2 = π/2 → θ = π = 180◦, we use a different formula u0 = 1 2 √ (1 + r11 + r22 + r33) u1 = 1 2 sgn(r32 − r23) √ (1 + r11 − r22 − r33) u2 = 1 2 sgn(r13 − r31) √ (1 − r11 + r22 − r33) u3 = 1 2 sgn(r21 − r12) √ (1 − r11 − r22 + r33) (2.91) where sgn(x) is the sign function of x sgn(x) = +1 for x > 0 sgn(x) = 0 for x = 0 sgn(x) = −1 for x < 0 According to (2.29) – (2.31), elementary rotations R(i, α), R(j , β) e R(k, γ) correspond to the following elementary quaternions: R(i, α) → u1 = (cos α 2 , sin α 2 , 0, 0) R(j , β) → u2 = (cos β 2 , 0, sin β 2 , 0) R(k, γ) → u3 = (cos γ 2, 0, 0, sin γ 2) (2.92) therefore the vectorial base of the quaternions corresponds to the three elementary rotations by 180◦ angles around the principal axes: i = (0, 1, 0, 0) → R(i, π) j = (0, 0, 1, 0) → R(j , π) k = (0, 0, 0, 1) → R(k, π) (2.93)

SLIDE 56

64 Basilio Bona - Dynamic Modelling Observe that, while the product of two equal elementary quaternions gives iı = jj = kk = ijk = (−1, 0, 0, 0), the analog product of two equal elementary rotation matrices gives the identity matrix, i.e., a rotation of 2kπ: R(i, π)R(i, π) = R(j , π)R(j , π) = R(k, π)R(k, π) = R(i, π)R(j , π)R(k, π) = I (2.94) corresponding to the quaternion (1, 0, 0, 0); this strange property is related to a new entity, called spinor, that will not be further discussed here; the interested reader can find more details in [1], [18] and [38, Ch. 11]. The quaternion operations that are related to the computation of rotations are the following:

1. Given n rotations R1, R2, · · · , Rn and the corresponding unit quaternions u1,

u2, · · · , un, the product of rotations R(u) = R(u1)Ru2) · · · R(un) corresponds to the quaternion product u = u1u2 · · · un in the same order;

2. given the rotation R(u) and the corresponding unit quaternion u, the transpose

matrix RT corresponds to the conjugate unit quaternion u∗ R(u) ⇔ RT(u∗)

3. given a generic vector x, that corresponds to a generic quaternion consisting
n the vectorial part only x = (0, x T) = (0, x1, x2, x3) and given the rotation

R(u) corresponding to the unit quaternion u, the rotated vector x ′ = R(u)x coincides with the quaternion product, that has always a zero real part x ′ = R(u)x ⇔ x′ = uxu∗

4. the product of two or more rotation matrices product R1(u1)R2(u2) is com-

puted adopting the following identity: R1(u1)R2(u2)x ⇔ u1(u2xu∗

2)u∗ 1 = (u1u2)x(u∗ 2u∗ 1) = (u1u2)x(u1u2)∗.

As a last comment, it should be taken into consideration the fact that in space applications and aerospace textbooks often the quaternions are organized in a list theta is different from the one adopted here: namely the real part is the fourth term

f the quaternion, not the first as in our notations.

SLIDE 57

Basilio Bona - Dynamic Modelling 65 Example 2.12.1 Given the elementary rotation R(j , 90◦), find the relative Euler parameters and the corresponding unit quaternion u Since R(j , 90◦) =   1 1 −1   applying relations (2.85) and (2.90), we have u0 = v4 = √ 2 2 u1 = v1 = 0 u2 = v2 = √ 2 2 u3 = v3 = 0 that are the same as those obtained applying (2.82). Example 2.12.2 Find the quaternion and the Euler parameters that represent the rotation R = R1 (i, 90◦) R2 (j , 90◦) , then compute the rotation axis u and the angle θ. The two quaternions are: R1 (i, 90◦) → u1 = (√ 2 2 , √ 2 2 , 0, 0 ) R2 (j , 90◦) → u2 = (√ 2 2 , 0 √ 2 2 , 0 ) from which we obtain the product quaternion R1 (i, 90◦) R2 (j , 90◦) → u1u2 =  1 2 − 1 2 i · j

1 2i + 1 2j + 1 2 (i × j )

  = (1 2, 1 2, 1 2, 1 2 ) Therefore the Euler parameters are v4 = h0 = 1 2; v1 = h1 = 1 2; v2 = h2 = 1 2; v3 = h3 = 1 2

SLIDE 58

66 Basilio Bona - Dynamic Modelling The angle is θ 2 = arccos (1 2 ) ⇒ θ = 120◦ and the axis is u = [ 1 1 1 ]T; since u norm is not unit we compute its norm ∥u∥ = √ 3 and after that, applying (2.82) we

btain the same result

1 2 = 1 √ 3 sin θ 2 ⇒ sin θ 2 = √ 3 2 ⇒ θ 2 = 60◦ ⇒ θ = 120◦

2.12.8

Cayley-Klein Parameters

Another parameterizations of rotations was introduced by Felix Klein (1849-1925), with the aim to make easier the integration of differential equations in complex gyroscopic problems [16]. This parameterizations, as well as that with quaternions or Euler parameters, has the advantage of not requiring the computation of trigonometric functions. The Cayley-Klein Parameters (CK parameters) are represented by complex 2×2 matrices Q = [α β γ δ ] (2.95) where { α β γ δ } are complex variables. Q must be unitary, hence the CK parameters obeys to the following constraints: α = δ∗ β = −γ∗ and the matrix Q can be written as Q = [ α β −β∗ α∗ ] with an additional constraint αα∗ + ββ∗ = 1 In Q there are three free parameters and they can be used to characterize the rotations: indeed from the CK parameters it is possible to compute the rotation matrix R =       1 2(α2 − β2 − γ2 + δ2) j 2(−α2 − β2 + γ2 + δ2) γδ − αβ j 2(α2 − β2 + γ2 − δ2) 1 2(α2 + β2 + γ2 + δ2) −j(αβ + γ + δ) βδ − αγ j(αγ + βδ) αδ + βγ       (2.96)

SLIDE 59

Basilio Bona - Dynamic Modelling 67 where j is the imaginary unit, j = √−1; the matrix R has real elements, also if it is function of several complex variables. Given the Euler angles (ϕ, θ, ψ), the CK parameters are computed as: α = ej(φ+ψ) cos (θ 2 ) β = jej(φ−ψ) cos (θ 2 ) γ = je−j(φ−ψ) cos (θ 2 ) δ = e−j(φ+ψ) cos (θ 2 ) (2.97) The relation between the CK parameters (α, β), the quaternion (h0, h1, h2, h3) or the Euler parameters (v1, v2, v3, v4) is the following: α = h0 + jh3 = v4 + jv3 β = h2 + jh1 = v2 + jv1 Matrix Q can also be written as: Q = h01 + j (h1σ1 + h2σ2 + h3σ3) (2.98) where 1 is the 2 × 2 identity matrix and σi are the so called Pauli Spin matrices 1 = [1 1 ] σ1 = [0 1 1 ] σ2 = [0 −j j ] σ3 = [1 −1 ] . (2.99) defined by the Pauli algebra σjσk = δjk1 + jϵjkmσm where δjk is the Kronecker delta and ϵjkm are the Levi-Civita symbols, defined in (2.127). For further details on the CK parameters one can refer to [16, pages 145-158].

2.12.9 Rotation Vectors – Rodrigues Vectors

Another way to represent the orientation of a rigid body consists in using only three parameters without incurring in the singularity problems encountered with Euler or RPY angles. Instead of using four parameters and a unit norm constraint, like in Euler parameters or in quaternions, one can introduce the so-called rotation vectors r, whose

SLIDE 60

68 Basilio Bona - Dynamic Modelling components ri describe the rotation axis, while its norm ∥r∥ provides the rotation angle or a trigonometric function of the rotation angle. In general the rotation vector r is defined as: r = f (θ) u (2.100) where u is the unit norm vector characterizing the rotation axis and ∥r∥ = f(θ) is an odd function2 of θ. The function f(θ) is usually chosen among the following list of odd trigonometric functions: a) f(θ) = θ b) f(θ) = sin θ c) f(θ) = sin θ 2 d) f(θ) = tan θ 2 Notice that the choice c) is related to quaternions, being the three components of r equal to the last three quaternion elements, or the first three Euler parameters. In this case the rotation vector r is called Euler (rotation) vector, not to be confused with the Euler parameters or the Euler angles. Choice d) brings to the so-called Rodrigues (or Gibbs) rotation vectors, quite used in theoretical kinematics. The relation between the Rodrigues vector r and the Euler parameters v is the following: r1 = v1 v4 , r2 = v2 v4 , r2 = v3 v4 (2.101) Notice that the Rodrigues vectors are undefined for odd angles θ = ±(2k + 1)π. Given two Rodrigues vectors r a and r b, their “product”, denoted by the symbol ⊙, is computed as: r a ⊙ r b = r a + r b − r b × r a 1 − r T

ar b

(2.102) The relations between rotation matrices and Rodrigues vectors, are established considering that:

1. the product of two rotation matrices RaRb is equivalent to the product r a⊙ r b
f the corresponding Rodrigues vectors;
2. in order to compute R given r, after setting the identities v1 = r1, v2 =

r2, v3 = r3, v4 = 1, one applies relation (2.84) and the elements found are divided by (1 + r2

1 + r2 2 + r2 3);

2An odd function is defined as −f(θ) = f(−θ).

SLIDE 61

Basilio Bona - Dynamic Modelling 69

3. in order to compute r given R one builds the skew-symmetric matrix S(r),

given by S(r) = R − RT 1 + tr (R) =   −r3 r2 r3 −r1 −r2 r1   = (2.103) and has immediately the elements ri.

Example 2.12.3

Considering again Example 2.12.1, we want to compute the Rodrigues vectors that represent the two rotations R1(i, 90◦) and R2(j , 90◦) and then compute the Ro- drigues vector of the product R = R1(i, 90◦)R2(j , 90◦). The Rodrigues vector relative to R1(i, 90◦) is r 1 = [ 1 ]T, while the Rodrigues vector relative to R2(j , 90◦) is r 2 = [ 1 ]T; The overall rotation corresponds to a Rodrigues vector computed as: R = R1(i, 90◦)R2(j , 90◦) ⇒ r = r 1 ⊙ r 2 =   1 1 1   Since ∥r∥ = √ 3, we have θ 2 = arctan √ 3 = 60◦ ⇒ θ = 120◦.

2.12.10

Orientation Error and Angle-Axis Representation

Often, in satellite or UAV orientation control problems it is necessary to compute the angular error between two different orientations: when Euler or RPY angles are used, the error between the desired αd and the measured ones αm is simply defined by the difference between the respective angles ∆α = αd − αm

∆α′ = αm − αd The first definition ∆α is commonly used by control designers that define as error the difference between the desired and the measured angles, while the second definition ∆α is used by measurement people that call error the opposite difference. Instead, when the rotation matrices are used, it is necessary to define the error in another way since it is a big error to compute the difference between two rotation matrices! Let us assume that Rd is the desired rotation matrix and Rm the measured one; the error may be defined in four different ways:

SLIDE 62

70 Basilio Bona - Dynamic Modelling 1) Re1: defined as RmRe1 = Rd → Re1 = RT

mRd.

2) Re2: defined as RdRe2 = Rm → Re2 = RT

d Rm.

3) Re3: defined as Re3Rm = Rd → Re3 = RdRT

4) Re3: defined as Re4Rd = Rm → Re4 = RmRT

d .

Notice that Re1 = RT

e2 and Re3 = RT e4.

Relations 1) and 2) define the orientation error as a post-product, and therefore are to be understood as a perturbation of the rotation with respect to the “mobile” reference frame, while relations 3) and 4) define the orientation error as a pre- product, and therefore are to be understood as a perturbation of the rotation with respect to the “fixed” or world reference frame. In any case the orientation error is the product of two rotation matrices, and in

rder to compute the axis unit vector ue and the angle ∆θ it would be necessary to

perform the products and then use the identities (2.43) and (2.46). Fortunately there exist another formula that one can use to compute the angle nd the axis of the rotation resulting from the product of rotation matrices, without actually doing the product. This relation is just used to compute the axis-angle of the orientation error. Assuming to have defined the orientation error as Re = RT

1 R2 – cases 1) and 2) –

then the following identity holds ue sin ∆θ = 1 2 (¯ r 2 × ¯ r 1 + ¯ g 2 × ¯ g 1 + ¯ b2 × ¯ b1 ) (2.104) where ¯ r i, ¯ g i and ¯ bi are the first, second and third rows of Ri, i = 1, 2. One can notice that the left-hand term in (2.104) is a type b) rotation vector, that goes to zero for angles equal to 0 or ±π, for which we know that identity (2.46) in not applicable. If we define the orientation error as Re = R1RT

2 – cases 3) and 4) – then the

following identity, also known as the unit vector lemma, holds u sin θ = 1 2 (r 2 × r 1 + g 2 × g 1 + b2 × b1) (2.105) where r i, g i and bi are the first, second and third columns of Ri, i = 1, 2.

2.12.11 Conclusions

In this Section we have seen that three parameters αi are sufficient to define the

rientation of a rigid body; the problem with three (angle) parameters is that they

are not singularity-free.

SLIDE 63

Basilio Bona - Dynamic Modelling 71 If we use the quaternions, the Euler parameters, the Cayley-Klein parameters or the Rodrigues vectors, it in unclear what are the parameters that take the place

f the αi in (2.65).

Using these alternative representations, we cannot say that the αi are “angles”, but it is necessary to speak in more general terms of “angular parameters”, whose knowledge nevertheless allows to compute the chosen (Euler or RPY or Cardan) angles. For instance, using the quaternions we can define: α1 = u1 sin θ 2, α2 = u2 sin θ 2, α3 = u3 sin θ 2 (2.106) and implicitly assume the unit norm. However it is a very common practice to associate to the αi the Euler or the RPY angles, although in the last few years quaternion representation has gained much attention, mainly in the aerospace and computer graphics applications.

2.13 Point Kinematics

Having introduced in the previous Sections the different representations of a geometrical point or of a vector, both polar or axial, in a given reference frame, and having also characterized the various types of rigid displacements in the Euclidean space, we are now ready to describe the motion of geometrical points, assumed massless. This description takes the name of kinematics and is distinct from dynamics since the former studies the motion of points or rigid bodies establishing the relations among positions, velocities and accelerations from a pure geometrical point of view, while the latter studies the influence of the forces and torques on the bodies on these quantities. In few words, kinematics study the motion without considering its causes and its effects, while dynamics studies how the external actions on the given bodies are related to their motion. Suppose to have a geometrical point P moving in the 3D Euclidean space with a time law specified by a function P(t); we can assume a vector representation of P in a given reference frame x P(t) =   x1P(t) x2P(t) x2P(t)  

r, to keep notation simple

x(t) =   x1(t) x2(t) x2(t)  

SLIDE 64

72 Basilio Bona - Dynamic Modelling Recalling the vector derivative in A.5, we write d dtx(t) =           d dtx1(t) d dtx2(t) d dtx2(t)           ≡ ˙ x(t) =   ˙ x1(t) ˙ x2(t) ˙ x3(t)   ≡ v(t) Since the velocity is the limit dP(t) dt = lim

∆T→0

− − → ∆P ∆T the signed segment − − → ∆P can be different in different reference frames. If the point P is fixed in a local reference frame Rm, we have [− − → ∆P ]

Rm = 0; but if

Rm moves with respect to a world frame R0, then [− − → ∆P ]

R0 ̸= 0. In the first case

we speak on “local” velocity, in the second case we speak of “global” or absolute velocity or total velocity. Usually when we do not specify otherwise, we will consider always the total velocity of a point. When the vector x(t) is subject to a time-varying rotation R(t)x(t), its derivate is computed according to the normal derivative product rule, i.e., d dt [R(t)x(t)] = [ d dtR(t) ] x(t) + R(t) [ d dtx(t) ] therefore it is necessary to derive the rotation matrix R(t). To understand the consequences of this operation, we have to study its properties.

2.13.1 Rotation Matrix Derivative and Angular Velocity

To start we assume that the generic rotation matrix R is function of a generic variable x, that is R = R(x). We recall that the orthogonality of the matrix implies R(x)RT(x) = I We derive this relation with respect to a generic variable x obtaining: dR(x) dx RT(x) + R(x)dRT(x) dx = O (2.107) It is evident that the first term is the transpose of the second term, therefore, recalling the definition of skew-symmetric matrix in Appendix B.6, we can define the following skew-symmetric matrix S (u(x))

def

= dR(x) dx RT(x), (2.108)

SLIDE 65

Basilio Bona - Dynamic Modelling 73 we have introduced a generic vector u(x) because we know that a skew-symmetric matrix in R3 embeds in its structure the components of a vector. We will later discuss its precise meaning. Taking both terms of (2.107) and post-multiplying them by R(x) we obtain dR (x) dx = S (u (x)) R (x) (2.109) and it results that the derivative of an orthonormal matrix is the matrix itself pre- multiplied by a suitable anti-symmetric matrix. S(u) in (2.108) is function of a generic vector u(x), itself a function of a generic scala variable x. In particular, when the rotation matrix R depends on an angle θ(t) around an axis given by the unit vector u, we can write: dR (u, θ(t)) dt ≡ ˙ R (u, θ(t))

def

= S (ω(t)) R (u, θ(t)) (2.110) where ω(t) is the instantaneous total angular velocity vector of the reference frame represented by the matrix R(t). The angular velocity is represented in the world reference frame, but we omit to write it as ω0 for notational simplicity. Now we take into consideration some simple cases of vector transformation R(t)r(t): we start with r(t) = r constant in time in a “mobile” reference frame Rm a we transform it n the world reference frame R0 r 0(t) = R0

m(t)r m

the time derivative is computed as: ˙ r 0(t) = ˙ R

m(t)r m

(2.111) the vectors in this expression are represented in two different reference frames, so it is necessary to represent both in R0; considering that r m = (R0

m)Tr 0

we have at the end, omitting for notational convenience the dependence of R on t, ˙ r 0(t) = ˙ R

1(R0 1)Tr 0(t)

(2.112) In mechanics textbooks one often encounters the following notation Ω(t) ≡ ˙ R

1(R0 1)T

that coincide with our previous skew-symmetric matrix S(ω); indeed Ω = ˙ R

1(R0 1)T = S(ω)R0 1(R0 1)T = S(ω)

SLIDE 66

74 Basilio Bona - Dynamic Modelling This matrix is often calle the angular velocity matrix; in some textbooks the skew-symmetric matrix Ω(t) = S(ω(t) is denoted with the symbol ω(t). Considering again the identity (2.111), we observe that it expresses the well known formula that gives the linear velocity of a point P fixed in a reference frame that rotates with an angular velocity ω: ˙ r 0(t) = Ω(t)r 0(t) = S(ω(t))r 0(t) = ω(t) × r 0(t) (2.113) Other interesting properties of the skew-symmetric matrices are related to the rigid body representation in 3D space. As shown in (2.40), it is possible to compute the rotation matrix R from the rotation axis represented by the non-unit vector v and the rotation angle θ = ∥v∥, considering the skew-symmetric matrix S(v). This equality is a theoretical consequence of the following property, in general valid for Lie groups and algebras (refer to [34] for details): R(v) ≡ R(u, θ) = eS(v) ≡ eS(θ, u) =

∞

∑

k=0

1 k!S k(θu) (2.114) that relates the rotation matrix R to the exponential of the skew-symmetric matrix S(v) = S(u, θ), function of the unit axis u and angle θ. Taking the Taylor series of the matrix exponential and considering that the rotations are elements of a cyclic group, and this reflects on the strusture of S i.e., S 3 = −S, S 4 = −S 2 we obtain: R = I + sin θ θ S(v) + 1 − cos θ θ2 S 2(v) (2.115) from which one gets the identities (2.40) and (2.41). Without entering into details, that can be found in the cited textbook, we write here the fundamental identity that relates the angular velocity vector ω(t), the unit vector u(t) and the angle θ(t): ω(t) = ˙ θ(t)u(t) + sin θ(t) ˙ u(t) + (1 − cos θ(t)) S (u(t)) ˙ u(t) (2.116) From (2.116) we see that the vector ω(t) is not the formal derivative of another vector, except from the simple case that u is constant, i.e., u(t) = c. In such a case we have: ω(t) = ˙ θ(t)c = ˙ r(t) (2.117) that represents the total time derivative of the vector r(t) = θ(t)c. Another useful formula to get ω(t) from a rotation matrix R = [ r g b ] is the following (see [17]): ω(t) = 1 2(r × ˙ r + g × ˙ g + b × ˙ b) (2.118) Notice the similarity of this equality with (2.105).

SLIDE 67

Basilio Bona - Dynamic Modelling 75

2.13.2 Infinitesimal Rotations

We have seen that, given two rotations respectively represented by two angular parameters α1 = [ α11 α12 α13 ]T and α2 = [ α21 α22 α23 ]T, the rotation resulting from their composition cannot be obtained (except particular cases) from the simple sum of the two angles: α ̸= α1 + α2 Now let us consider an infinitesimal rotation defined by the infinitesimal angle vector dα = [ dα1 dα2 dα3 ]T the associated rotation matrix can be written as R(dα) = R(u(t), dθ(t)) = Ru(t),dθ(t) (2.119) From (2.115), and taking into account that for dθ → 0 we have 1 − cos (dθ) (dθ)2 → 0 sin (dθ) (dθ) → 0 the following approximate relation holds R (dα) ≃ I + S (dα) . (2.120) The inverse of R(dα) is computed as R(dα)−1 = R(dα)T ≃ I − S(dα) = I + S(−dα) ≃ R(−dα) This relation can be verified writing R(dα)TR(dα) = (I + S)(I − S) = I + S − S + S 2 ≃ I where S 2(dα) ≃ O since it contains second order infinitesimals, we conclude that the inverse of R(dα) is R(−dα). Using (2.120) it is possible to show that dα behaves like a proper vector. For instance, the commutative property holds; indeed the following identity is true: R (dα1) R (dα2) = (I + S (dα1)) (I + S (dα2)) = I + S (dα1) + S (dα2) + S (dα1) S (dα2) Since the last term goes to zero, and recalling (B.3), we have R (dα1) R (dα2) = I + S (dα1) + S (dα2) = I + S (dα1 + dα2) = R (dα1 + dα2) = R (dα2 + dα1) = R (dα2) R (dα1) (2.121)

SLIDE 68

76 Basilio Bona - Dynamic Modelling We know that from (2.110)) it results ˙ R(t) = S (ω(t)) R (u, θ(t)) (2.122) Since R does not depend on the choice of the representation of the parameters α used to characterize the rotation, the formula (2.122) can be re-written as: ˙ R (α(t)) = S (ω(t)) R (α(t)) (2.123) Considering the incremental rate we have ∆R = R (α + dα) − R (α) ≃ ˙ R(α)dt = S (ω(t)) R (α) dt = S (ω(t)dt) R (α) = S (dα) R (α) from which it follows R (α + dα) ≃ R (α) + S (dα) R (α) = (I + S (dα)) R (α) = R (dα) R (α) (2.124) We recall that the last term in the above identity does not commute. We conclude this Section with an observation: the role played by the skew-symmetric matrix S(ω) is fundamental in the definition of the angular velocity of a reference frame, and the two relations (2.114) and (2.122) are an example of this role. Now, if we introduce the following three matrices M i: M 1 =   −1 1   M 2 =   1 −1   M 3 =   −1 1   (2.125) thta are commonly called infinitesimal rotation generators), we are able to build S(ω) as a weighted sum S(ω) = ∑

i=1,3

ωiM i and observe that the M i form a basis for S(ω). Moreover the M i have the property that the differences between their mutual products obey to M iM j − M jM i ≡ [M i, M j] = ϵijkM k (2.126) where the ϵijk are called permutation symbols or Levi-Civita symbols defined as ϵijk =    i = j, j = k, k = i; two symbols out of three are equal +1 (i, j, k) ∈ {(1, 2, 3), (2, 3, 1), (3, 1, 2)} −1 (i, j, k) ∈ {(1, 3, 2), (2, 1, 3), (3, 2, 1)} (2.127) The matrix operator in (2.126) is referred as a commutator or Lie bracket and it defines the Lie algebra on the orthonormal matrix group. The properties and the importance of such an algebra will not be detailed here; the interested reader can find more details in [34].

SLIDE 69

Basilio Bona - Dynamic Modelling 77

2.13.3 Infinitesimal Rotations and Quaternions

If the rotation angles are small it is possible to directly compute the corresponding quaternions and vice-versa. Let us assume to have a small rotation ∆θ represented by the (small) angle parameters ∆α, and the corresponding quaternion h = [ h0 h1 h2 h3 ] ; from (2.120) we have R(∆α) = I + S(∆α) =   1 −∆α3 ∆α2 ∆α3 1 −∆α1 −∆α2 ∆α1 1   =   r11 r12 r13 r21 r22 r23 r31 r32 r33   (2.128) comparing (2.128) with (2.89), we obtain: ∆α1 = 1 2(r32 − r23) = 2h1h0 ≃ 2h1 ∆α2 = 1 2(r13 − r31) = 2h2h0 ≃ 2h2 ∆α3 = 1 2(r21 − r12) = 2h3h0 ≃ 2h3 Taking into consideration definition (2.119), the approximation is reasonable; indeed, applying (2.114), we obtain h0 = cos ∆θ 2 ≃ 1 since ∆θ is a small angle. Hence, if we know the value of an infinitesimal rotation dα = [ dα1 dα2 dα3 ]T, we can easily compute the corresponding quaternions as h ≈ [ 1 dα1 2 dα2 2 dα3 2 ]

2.14 Total Velocity and Acceleration of Points

When we speak of linear or angular velocites of a geometrical point or a rigid body, we mean the relative velocity of the geometrical point with respect to a well defined reference frame, or the velocity of the reference frame describing the rigid body with respect to another well defined reference frame. To be precise, when we say that a rigid body has at the time t an angular velocity ω(t), we mean that at t the relation (2.116), that we report below for the sake of clarity, is true: ω(t) = ˙ θ(t)u(t) + sin θ(t) ˙ u(t) + (1 − cos θ(t)) S (u(t)) ˙ u(t) This means that we have to give non only ˙ θ(t), but also the instantaneous rotation axis u(t) at time t. If the rotation axis remains constant, then (2.116) reduces to (2.117).

SLIDE 70

78 Basilio Bona - Dynamic Modelling It is important to note that the vector ω(t) cannot be translated parallel to itself without modifying the physical significance of the described phenomenon, since the application point of ω(t) gives the position of the axis around which the rotation takes place; this is not valid for the linear velocity given by the vector v(t): this one is a free vector, while ω(t) is an applied vector. We recall that an applied vector is defined by the couple (P, w), where P is the geometrical point of application and w is defined such that − → PQ = w To be more precise, when we speak of angular (or linear) velocities we should use the following (cumbersome) notation: ωk

ij(P)

(2.129) where we denote the relative angular velocity of the reference frame (or of the associated rigid body) Rj with respect to the reference frame Ri, represented in the reference frame Rk, and where P is the application point of ω; when P is known

r superfluous, we only write ωk
ij. Usually we implicitly assume to represent the

angular velocity as ωi

ij(P).

Whit this notation it is always true that ωk

jj = 0

∀j, k (2.130) since no relative motion will take place, while, in general, it results ωj

ij ̸= ωi ij

∀i, j (2.131) provided that Ri and Rj are non-coincident. Moreover, given the vectorial nature

f the velocity, the additive property holds

ωk

im + ωk mj = ωk ij

from which we have the following property ωk

ij + ωk ji = ωk ii = 0

⇒ ωk

ij = −ωk ji

(2.132) that says that, if a reference frame rotates with an angular velocity with respect to another reference frame, this last rotates with an angular velocity that is the same in module, but opposite in direction to the previous one. Now assume that x 1 is the representation of the geometrical point A in R1 and that this frame roto-translates with respect to R0; assume also that R0

1(t) and t0 1(t) are

respectively the rotation matrix and the translation vector between the two frames, and that ω0

01(P) is the relative angular velocity of R1 with respect to R0, represented

in R0. We have still to define the application point P; it varies according to the fact that the rotation takes place with respect to the fixed frame R0 axes, or with respect to the translated local frame R1 axes. In the following we will consider the two cases.

SLIDE 71

Basilio Bona - Dynamic Modelling 79 Case A We consider here a roto-translation described by T 0

1 = T(d 0 1)T(R0 1)

(2.133) that represents a translation d 0

1, followed by a rotation R0 1, with instantaneous

rotational velocity ω0

01, performed with respect to the mobile frame (post-mobile

rule). The motion defined in (2.133) allows the following alternative description: “a rotation followed by a translation performed with respect to the fixed frame (pre-fixed rule)”. In this case we can state that the application point P of the applied vector(P, ω) is the origin of R1. In this case the vector r 1(t) represented in R1 is described in R0 by: r 0(t) = R0

1(t)r 1(t) + d 0 1(t);

(2.134) Computing the time derivative of this identity, recalling (2.113) and omitting for brevity to write the time dependence in the various terms, we have: ˙ r 0 = ˙ R

1r 1 + R0 1 ˙

r 1 + ˙ d

= S (ω0

01) R0 1r 1 + R0 1 ˙

r 1 + ˙ d

= ω0

01 ×

( R0

1r 1)

+ R0

1 ˙

r 1 + ˙ d

(2.135) If now we define ρ0(t) = R0

1r 1(t)

(2.136) we have ˙ r 0(t) = ω0

01(t) × ρ0(t)

+ R0

1 ˙

r 1(t)

+ ˙ d

1(t) C

(2.137) This is the classical relation that gives the total linear velocity in R0 of a point P that moves with a linear velocity ˙ r 1(t) in the reference system R1, itself moving at a linear velocity ˙ d

1(t) and angular velocity ω0 01(t) applied at the origin of R1.

The three terms are: A) the tangential velocity due to the rotational velocity of the frame; B) the proper velocity of the point, represented in R0; C) the translational velocity of the frame. The acceleration ¨ r 0(t) of the point P expressed in R0 can be computed considering the identity: d (a × b) dt = (da dt × b ) + ( a × db dt )

SLIDE 72

80 Basilio Bona - Dynamic Modelling Taking the time derivative of ˙ r 0(t) in (2.137), using (2.110) and (2.136), and omitting again to write the time dependence, we obtain: ¨ r 0(t) = d dt (ω0

01(t) × ρ0(t)) + d

dt ( R0

1 ˙

r 1(t) ) + d dt ˙ d

1(t)

= ˙ ω0

01 × ρ0 + ω0 01 × ˙

r 1 + ˙ R

1 ˙

r 1 + R0

1¨

r 1 + ¨ d

= ˙ ω0

1 × ρ0 + ω0 01 × d

dt ( R0

1r 1)

+ ˙ R

1 ˙

r 1 + R0

1¨

r 1 + ¨ d

= ˙ ω0

01 × ρ0 + ω0 01 ×

( ˙ R

1r 1 + R0 1 ˙

r 1) + ˙ R

1 ˙

r 1 + R0

1¨

r 1 + ¨ d

= ˙ ω0

01 × ρ0 + ω0 01 × (ω0 01 × ρ0) + ω0 01 ×

( R0

1 ˙

r 1) + ˙ R

1 ˙

r 1 + R0

1¨

r 1 + ¨ d

= ˙ ω0

01 × ρ0 + ω0 01 × (ω0 01 × ρ0) + ω0 01 ×

( R0

1 ˙

r 1) + ω0

01 ×

( R0

1 ˙

r 1) +R0

1¨

r 1 + ¨ d

= ˙ ω0

01 × ρ0

+ ω0

01 ×

( ω0

01 × ρ0)

+ 2ω0

01 ×

( R0

1 ˙

r 1)

+ R0

1¨

r 1 + ¨ d

(2.138) where A) is the tangential acceleration; B) is the centripetal acceleration; C) is the Coriolis acceleration; D) is the sum of the linear acceleration of the point represented in the fixed frame and the acceleration of the frame (already represented in R0). Case B We consider here a roto-translation described by

1(t) = T

( R0

) T ( d 0

) (2.139) that represents a rotation R0

1 with instantaneous rotational velocity ω0 01, followed

by a translation d 0

1 performed with respect to the mobile frame (post-mobile rule).

The motion defined in (2.139) allows the following alternative description: “a translation followed by a rotation performed with respect to the fixed frame (pre-fixed rule)”. In this case we can state that the application point P of the applied vector(P, ω) is the origin of R0. The point r 1 is described in R0 by: r 0(t) = R0

1(t)r 1(t) + R0 1(t)d 0 1(t);

(2.140) Taking the time derivative and omitting again to write the time dependence, we

btain:

˙ r 0 = ˙ R

1r 1 + R0 1 ˙

r 1 + ˙ R

1d 0 1 + R0 1 ˙

= S(ω0

01)R0 1

( r 1 + d 0

) + R0

1 ˙

r 1 + R0

1 ˙

= ω0

01 × (R0 1r 1 + R0 1d 0 1) + R0 1 ˙

r 1 + R0

1 ˙

(2.141) If now we define ρ0 = R0

1r 1 + R0 1d 0 1

(2.142)

SLIDE 73

Basilio Bona - Dynamic Modelling 81 we have ˙ r 0 = ω0

01 × ρ0 + R0 1 ˙

r 1 + R0

1 ˙

(2.143) that is identical to (2.137), with R0

1 ˙

1 instead of ˙

1. Taking the time derivative of

(2.143) we obtain the total acceleration: ¨ r 0(t) = ˙ ω0

1 × ρ0 A

+ ω0

01 ×

( ω0

01 × ρ0)

+ 2ω0

01 ×

( R0

1 ˙

r 1 + R0

1 ˙

)

+ R0

1¨

r 1 + R0

1¨

(2.144) The formula has the A), B), C) and D) terms identical with those in (2.138), the only difference being the translational velocity and acceleration of the frame R1 expressed in the fixed reference R0. In conclusions, the total acceleration of a point that moves in a local reference frame, that itself moves with respect to a fixed frame, is the sum of the following terms (time dependency omitted for brevity): Tangential acceleration Due to rotational acceleration of the local frame with respect to the fixed frame: A = ˙ ω0

01 × ρ0

where ρ0 is the distance of the point P from the origin of the fixed frame. Centripetal acceleration Due to the variation of the tangential acceleration, that produce an acceleration directed toward the instantaneous rotation center: B = ω0

01 ×

( ω0

01 × ρ0)

Coriolis acceleration Due to the linear velocity in a rotating frame: C = 2ω0

01 × v 0

where v 0 = R0

1 ˙

r 1

v 0 = R0

1( ˙

r 1 + ˙ d

is the total linear velocity; Linear acceleration Due to the linear acceleration of the point P and of the reference frame: a0 = R0

1¨

r 1 + ¨ d

a0 = R0

1(¨

r 1 + ¨ d

Notice that the symbol ρ0 has a different meaning in the two cases, as put in evidence by (2.136) and (2.142).