# Optimization and Gimbal Lock Control via Shifted Decomposition of Rotations

^{*}

**Corresponding Author:**Brezov D, Department of Mathematics, University of Architecture, Civil Engineering and Geodesy, Hristo Smirnenski, Bulgaria, Tel: +359 2 963 5245, Email: [email protected]

*
Received Date: Jul 10, 2018 /
Accepted Date: Aug 07, 2018 /
Published Date: Aug 09, 2018 *

**Keywords:**
Rotations; Spacecraft attitude; Robot kinematics; Optimal control

#### Introduction

Euler angles have been a major tool in the description of attitude orientation of spacecraft and aircraft since these technologies exist. However, the method causes some troubles, such as the so-called gimbal lock problem related to the inevitable topological singularity of the parameterization map between the three-dimensional torus and the rotation group in which has the topology of projective space . This problem does not appear if we use quaternions or axis-angle decomposition [1] associated with the Hopf bration , but in practice this may be even less efficient in terms of attitude adjustment. It has been [2] the gimbal lock singularity may be utilized within the Euler angles' framework using a simple optimization procedure on the excessive free parameter at the singular point. The present paper suggests an analogous procedure in the regular case as well, by means of introducing a fourth “shift" parameter. Such an additional degree of freedom allows for optimizing and thus, saving a significant portion of the energy spent on spacecraft maneuvers as well as manipulating the gimbal lock condition. A similar idea for optimization [3,4] for sequences of arbitrary length in the two-axes setting. Here we allow three axes while restricting to four factors, which turns out to be sufficient in many cases [5]. Hence, the cost function may be expressed explicitly in terms of α, due to the closed form solutions obtained below. The derivation of these results is based on a recent study on the generalized Euler decomposition problem [6] considering arbitrarily oriented axes and the corresponding necessary and sufficient condition [7]. In this extended framework we may weaken the Davenport condition^{1} γ_{12}+γ_{23} ≥ π for the relative angles between the axes [8] to γ_{12}+γ_{23}+γ_{31} ≥ π, thus allowing for different geometries of attitude control mechanisms in vehicles and robots, beyond the classical orthogonal settings. The text is organized as follows: in the following section the quaternion based vector-parameter is introduced and utilized to solve the generalized three-axes decomposition problem. Section three explains in detail the shift parameter construction, while sections four and five focus on its applications in gimbal lock control and optimization of rotational sequences, respectively. Several examples are given to illustrate the way variations of the shift parameter influence the energy efficiency of attitude maneuvers. Here we consider only minimization of the rotation path, while the impact of specific asymmetries is being modelled by introducing weights in the cost function. A more detailed discussion from dynamical perspective [9,10].

#### Vector Decompositions of Rotations

The vector-parameter may be introduced as

is a representative of the spin group so c actually inhibits the projective space . Note that it may also be expressed via the Euler’s trigonometric substitution as** c=τn** where is usually referred to as the *scalar parameter* (ϕ being the angle of rotation) and stands for the unit vector along the invariant axis. This parameterization provides several major advantages. Firstly, it allows for expressing the rotation matrix entries as rational functions of **c**, namely

(1)

where stands for the identity, **cc ^{t}** denotes the tensor (dyadic) product and c

^{×}is an extension of the Hodge duality that de nes the cross product in as a×b=a×b. As the right-hand side in eqn. (1) is the well-known Cayley map, we have

Another advantage is the compact form of the group composition law

(2)

where stands for the scalar (dot) product in . Note also that formula in eqn. (2) is just the quaternion multiplication rule projected onto . Composing rotations in this manner appears far more efficient in terms of computational complexity compared to the usual matrix multiplication [5].

With the aid of formula in eqns. (1) and (2) one may derive [6] the explicit solutions to the generalized Euler decomposition problem where stands for the vector-parameter of the i-th rotation in the decomposition, i.e., More precisely, we have

(3)

where ε_{ijk} denotes the Levi-Civita symbol. Furthermore, we make use of the notation

as well as

and the necessary and sufficient condition for the existence of real solutions

(4)

Whenever any of the conditions γ_{ji}=g_{ji} (i≠j) holds, one may decompose into a pair of rotations about the i-th and the j-th axis as with

(5)

where the notation and is used. With i=1,j=2 in particular, these expressions apply also to the gimbal lock setting

(6)

in which the first equality in eqn. (5) yields the scalar parameter of ϕ1 ± ϕ3.[2,6].

#### The Shift Parameter Construction

Let us consider a specific type of decompositions involving four rotations about three initially given axes. Namely, we rotate twice about one of the axes introducing an additional parameter in one of the following ways

where subscripts label the invariant rotation axes used in the factorization and Greek letters denote the corresponding angles. With the aid of group conjugation it is straight-forward to [2] that

(7)

where stands for a rotation about the axis . This technique allows also for the second decomposition above to be written in a quite similar form as

(8)

With . Thus, both cases may effectively be reduced to decompositions into three factors with one of the axes “shifted” by an α-rotation.

The next step is to use the additional fourth parameter to guarantee the necessary and sufficient condition eqn. (4) for the corresponding decomposition. In the case eqn. (7) instead of eqn. (4) this procedure yields.

(9)

where we may express the α-dependent matrix entry as

with the standard notation a[_{i}b_{j}]=a_{i}b_{j}−a_{j}b_{i}. On the other hand, the condition in eqn. (9) may be written in the equivalent form [7].

Moreover, using a standard trigonometric substitution one may express

with

This allows for writing the above inequalities in the simpler form

where we make use of the notation

Note that since we have the restriction |cos(α−α0)|≤1, formula eqn. (15) determines a real interval for the shift angle α for arbitrary compound rotation if and only if

On the other hand, one has

so the above conditions may be written also as

(10)

which finally yields for the unknown parameter^{2}

(11)

The interval determined in this way is typically larger compared to

(12)

that guarantees decomposition into three consecutive rotations about the [7]. In particular, for γ_{12}=γ_{13}=γ the lower bound in eqn. (11) is trivial and we end up with . Moreover, one obviously has (the extremal values of correspond to )

(13)

and since the case also yields a trivial lower bound in eqn. (11). Hence, the sufficient condition for the compound rotation's angle is

(14)

Note that as long as (11) holds, one may choose an arbitrary

(15)

The endpoints of the interval in eqn. (15), for example, yield Δ'=0 and hence, rational expressions for the τ_{k}'s. More generally, using the quantities (otherwise g_{ij}'=g_{ij} and r_{ij}'=r_{ij} )

as well as

one obtains the explicit solutions with the aid of formula in eqn. (3) in the form

(16)

which gives for the remaining three rotation angles in the decomposition

If, on the other hand, is parallel to either or one cannot use as a shift parameter: in the former case one ends up with only two factors that imposes a far more restrictive condition r_{31}=g_{31} and in the latter, a decomposition of the type takes place. Note that formula in eqn. (14) is too restrictive as it always ensures such decompositions into three factors. The shift parameter, however, effectively weakens in eqn. (12) and may thus be utilized for engineering purposes as already discussed above.

Next, we consider in eqn. (8) with the discriminant condition

(17)

where the matrix entry r'_{31} is explicitly given in the form

and may be further simplified as

with

This yields the inequalities

which have real solutions for α as long as

Moreover, taking into account that

we end up with the system

(18)

and with the aid of formula in eqn. (13), we obtain a condition for ϕ in the form

(19)

In particular, in the symmetric case γ12=γ23=γone has while for r_{21}=±1, the interval in eqn. (15) for α is either empty, or there is no restriction, the condition being depending on the sign. Now, if in eqn. (18) holds and , one may choose arbitrary α in eqn. (15), where

and calculate (for all other values of i and j we have respectively )

in order to obtain the solutions based in eqn. (3) in the form

(20)

which finally provides the remaining three rotation angles

Note also that the decomposition in eqn. (8) is, in some sense, dual to in eqn. (7). Namely, using matrix transposition one easily shows that the former is equivalent to

(21)

Where stands for a rotation about the unit vector so the shift angle in this case is ϑ . It is not hard to see that it is real as long as

with

that leads to in eqn. (18). Then, if the latter holds and , we choose arbitrary

(22)

where and calculate

Finally, Δ" is obtained in eqn. (4) with replacing respectively g_{ij} and γ_{31}. With the above data the solutions are expressed via formula in eqn. (3) as

(23)

and provide the remaining three rotation angles in the form

Let us also consider a third type of decompositions, namely

(24)

in which it is no longer possible to conjugate any of the factors. However, one may work with the compound rotation instead, writing the above as

Similarly to the previous cases, we express the α-dependent parameter

with the notation in the form

where

Then, the necessary and sufficient condition

(25)

may be written for the shift angle α as

On the other hand, one has

so real solutions for α exist as long as

that yields for the interval

(26)

Note that the lower bound is trivial since defined as nonnegative, while the expression on the left is non-positive due to the triangle inequality. Then, the corresponding condition for the compound rotation’s angle reads

(27)

which allows for a much larger class of configurations compared to the classical Davenport setting. Consider for example the Pyraminx, in which the three axes are oriented along the edges of a proper tetrahedron. The existence of the decomposition for arbitrary demands hile the restriction imposed by formula in eqn. (27) namely is trivial. In the generic case one proceeds with the solution as follows: first, choose a value of α in the non-empty interval in eqn. (15), where^{2}

and use it to obtain the shifted rotation that is going to replace in formula in eqn. (3). Then, calculate

With as well as

and thus express

(28)

where Δ' is given by formula in eqn. (25). Hence, one has for the remaining angles

Taking into account in eqns. (11), (18) and (27) it is not hard to justify the following

**Theorem**

Given two or three distinct points on connected by a threesegment spherical geodesic path (with possible repetitions) of length Lγ ≥ π, any SO(3) transformation is decomposable into four consecutive rotations, each fixing one of those points.

In particular, the vertices of any spherical triangle on with perimeter at least π determine the invariant axes in for such a factorization. This triangle, however, may be degenerate, i.e., the points are allowed to lie on the same big circle. Moreover, two of them may coincide, in which case for instance γ_{12}=γ_{23}=γ and one has L_{γ}=3_{γ} that is an old result due to Lowenthal on the order of finite generation of the orthogonal group in . Note that the above theorem allows in principle for generalization to longer factorization sequences, but in this study we are interested in the case of four factors.

#### Gimbal Lock Control

Since condition in eqn. (6) determines a zero measure set, the above construction may be used to avoid it as long as the interval for the shift parameter contains more than one point. On the other hand, to inflict gimbal lock on purpose using the shift parameter, one needs another nontrivial condition. Let us consider once more the decomposition in eqn. (8), written as a three-factor product. The gimbal lock condition yields for the shifted third axis

which imposes the geometric restriction

The latter has the trivial solution α=0 for (the gimbal lock is already present without the shift) plus another one in the form . The interval for α in those cases is non-empty as long as the inequality holds for the positive sign and respectively for the negative one. Note that in this setting we have so the gimbal lock condition yields cos(α-α_{0})=±1 for the shift parameter, which has the solutions α=α0 and respectively, α=α_{0}-π. Thus, for the α-interval we have in the former case and in the latter, which has a solution only if γ_{12}=γ_{23} for the positive sign and γ_{12}+γ_{23}=π for the negative one. Note that in both cases this allows for factorization into a pair of rotations, since one has γ_{21}=g_{21}. In particular, in the Davenport setting g12=g23=0 one needs also γ_{21}=0, which yields^{4}. and , so we may set and thus obtain for the shift parameters respectively and . Therefore, avoiding the above values guarantees avoiding gimbal lock even if γ_{21}=0. Moreover, in this case both signs in eqn. (6) are possible since we have

One may approach the decomposition in eqn. (7) in a similar way. Note that the α-parameter only alters the middle axis, so it cannot affect the gimbal lock. Actually, γ_{31}=±1 yields A=0 and in eqn. (10) is possible only if g_{12}=0 or g_{13}=1. On the other hand, using matrix inversion, one may express in eqn. (7) in the form

(29)

where stands for a rotation about the unit vector , so φ plays the role of a shift angle. The discriminant condition for this setting is

(30)

Where

which may be written also as

With

Taking into account that

real solutions for φ are easily seen to be available as long as the inequalities

hold, which leads to in eqn. (11). Thus, in the generic case γ_{12}≠ γ_{13} namely formula in eqn. (11) determines the resolution set and the condition for the compound rotation’s angle in eqn. (14).

Then, for we set the value of ϕ in the interval

(31)

where

and proceed as before calculating

Using these quantities one easily obtains the scalar parameters

(32)

where Δ' is given by formula in eqn. (30) and for the remaining angles one has

The gimbal lock condition for this case may be expressed in the form

which naturally yields α=0 for γ_{32}=± 1 and the second solution is obviously γ_{31}= ± g_{21}. Identical arguments here show that in order to inflict gimbal lock one needs γ_{31}=g_{31} as well as either of the relations γ_{12}=γ_{13} and γ_{12}+γ_{13}=π for the two signs in formula in eqn. (6), respectively. In both cases, however, a decomposition of the type is possible. For instance, in the setting g_{12}=g_{13}=0 one easily obtains the values of the shift parameter ϕ_{+}=γ_{32} and ϕ_{-}=π+γ_{32} and . Since the admissible interval is

we see that both possibilities for gimbal lock are present in this case.

Finally, the gimbal lock condition in eqn. (24) may be written as

which yields, apart from the trivial solution , another one in the form . Since in this case is equivalent to cos(α−α_{0})= ± 1, one needs either γ_{12}=γ_{23} or γ_{12}+γ_{23}=π, respectively. In the Davenport setting for example, the corresponding values of α are once more given by while

Therefore, just as before, both signs in formula in eqn. (6) may occur.

Allowing permutation of axes, we systemize the above results in the following

**Theorem**

Avoiding gimbal lock via the shift parameter technique is always possible as long as there are two distinct angles γ_{ij}, whose sum is

Inflicting it, however, may be done only in one of the two cases:

1. γ_{ij}=g_{ij}= ± g_{jk}=0 that yields a decomposition into two factors

2. γ_{ii}= ± g_{ij} and g_{ik}= ± g_{jk}

where the values of i; j and k are assumed to be different.

#### Optimization

The additional fourth parameter introduced above (usually denoted by α) may be used for optimization of rotational sequences as long as it is defined over a set of non-zero measure. One may naively suggest that the total minimum of the cost function

(33)

is reached at α=0, but our numerical tests reveal quite a different picture. Consider for example the decomposition in eqn. (7) and let the unit vectors be oriented along the coordinate axes (the so-called Bryan setting). Then, we obviously have g_{ij}=δ_{ij}, where denotes the Kronecker symbol, and γij are simply the matrix entries of in the standard basis. Moreover, the α-shift of the vector is obviously a rotation in the YOZ plane and may thus be written as , which facilitates the calculation of the other primed quantities we need, namely

and finally Substituting the above expressions in formula in eqn. (16) yields for the scalar parameters

(34)

and the cost function may easily be obtained in the form

Moreover, the gimbal lock singularity in eqn. (6) persists for all values of α, but the condition in eqn. (9) demands α=0 or α=π, so we have a transfer of parameters α→ψ ∈ [−π,π] and minimize εγ(ψ). Our example is a rotation by an angle ϕ=−120° about the vector n ∼(3,4,5) t. The minimum of the cost function ε_{γ}(α_{0}) ≈ 179.81° obtained at α0 ≈ 13.49° yields about 36.49% efficiency increase compared to the standard Bryan setting ε_{γ}(0) ≈ 245.31° **(Figure 1a)** and the optimal sequence of rotation angles is easily derived as

Another relatively simple setting is provided by the classical Euler ZXZ decomposition.

We consider this time (8) with oriented along the z-coordinate axis aligned with OX. This yields which allows us to obtain the scalar parameters in the form

(35)

where we have denoted Δ'=sec^{2}α-(γ_{33}-γ_{23} tanα)^{2} and by abuse of notation, γ_{ij} are the matrix entries of in the standard basis. Moreover, the shift parameter α ∈ [−π,π] is unrestricted since one has . To illustrate this case, we consider a half-turn about n∼(5,4,3)^{t} and depict graphically the dependence **(Figure 1b)**.

εγ(α)=|α|+|2arctanτ_{1}|+|2arctan τ_{2} −α|+|2arctan τ_{3}|.

Note that here α=0 is a local maximum of the cost function εγ(α) that is somewhat surprising. The global minimum at α0≈−105.37° yields a total difference of ε_{γ}(0)− ε_{γ}(α_{0})≈50.65°, i.e., about 16.35% decrease in energy consumption, or approximately 19.54% increase in efficiency compared to the standard Euler decomposition. The two points of discontinuity α1≈−143.13° and α2≈36.87° are related to the condition γ'_{32}=g'_{32}, i.e., φ=0. The optimal sequence of rotation angles in this case has the form

Finally, considering the XY ZX setting in formula in eqn. (24) one easily finds

(36)

with Δ''=sec^{2}α − (γ_{31}-γ_{21} tanα)^{2}. Note that the condition γ_{11}=0 holds, so the endpoints of the interval for α correspond to gimbal lock singularity. Consider again the former rotation discussed above in this section **(Figure 1c)**. The minimum at α_{0}≈−74.97° yields about 35.82% increase in efficiency. The optimal sequence of rotation angles in the decomposition this time is

Note that the cost functions ε_{γ}(α) in all examples above imply either spherical or cylindrical symmetry of the spacecraft or aircraft. In the case of significant anisotropy, however, due to different moments of inertia or external forces, e.g., gravitational, magnetic, air flow etc., it is possible to adjust εγ(α) by introducing specific weight for each rotation axis. For instance, in the decomposition in eqn. (7) one may replace in eqn. (33) with

**Figure 2 **depicts the weighted cost functions in the examples considered above, the weights being κ_{1}=1, κ_{2}=1/2 and κ_{3}=1/3 in all three cases. Note that in the first setting the violation of the rotational symmetry does not alter much the graph of εγ(α) apart from re-scaling and a slight shift of the minimum. In the second example, on the other hand, the minimum has moved to the singular point, thus reducing the length of the optimal sequence from four to two, so we encounter a sort of phase transition. Finally, in the last case we observe an energy transfer between the two local minima.

#### Concluding Remarks

The construction proposed in this paper has several advantages that seem quite useful in various practical situations. To begin with, it allows us to weaken the classical conditions due to Davenport and Lowenthal and thus, consider a wide variety of decomposition geometries, e.g., tetrahedral or asymmetric, which satisfy a condition of the type Σγ_{ij} ≥ π. Using this technique, on the other hand, it is easy to avoid gimbal lock configurations and more importantly, increase (in some cases significantly) the energy efficiency of attitude adjustments used in spacecraft and aircraft maneuvers. Unlike other optimization methods (e.g. the famous Wahba problem, it provides explicit exact solutions: numerical procedures are involved only at the final stage of obtaining the minimum of the cost function ε_{γ}(α), which is an almost trivial operation. We have chosen orthogonal axes to illustrate the method just for the sake of simplicity and because this choice allows for a wider range of α. Note also that our construction naturally extends to the spin covering group SU(2) and therefore, may be applied (among other things) to two level systems [3]. Quite a similar procedure can be derived also for the group SO(2,1) and its spin cover based on analogous results [6] as long as one takes into account both the gimbal lock and isotropic singularities. Let us point out that in the hyperbolic case the former is far easier to control - one simply needs to choose the first and the last invariant axes in the sequence to be of different geometric types (like in the Iwasawa decomposition), while the latter is not affected by the shift parameter since in the singular setting the tangent to the light cone is an invariant subspace independent on the specific values of the τi’s. Moreover, the definition of the cost function in SO (2,1) needs to be modified as formula in eqn. (33) tends to mix angles with rapidity’s. Another possible generalization may be seen in the context of screw dynamics modelled with dual quaternions and respectively, dual vector-parameters. The cost function in this case would be a (weighted) path on the Galilean group of rather than one on as in the above considerations. The advantages of this approach in optimization and avoiding singularities may be utilized not only in motion control but also areas like quantum information [3] or computer vision and virtual reality [1], but this goes beyond the scope of the present article and will be left for further study.

#### Acknowledgements

I am grateful to Professor Ivaılo Mladenov at the Bulgarian Academy of Sciences for inspiring my interest in vectorial parameterization and decomposition of rotations in .

^{1}Here γ_{ij} denotes the minimal angle between the rotation axes, so orthogonality is implied.

^{2}We may always assume for two of the relative angles , while

^{3}We may always assume as the case does not impose any restriction on α.

^{4}In this case the unit vector is confined to the plane spanned by and

#### References

- Beckermann C (2002) Modeling of Solidification, in: Purdue Heat Transfer Celebration. April 3-5 University of Iowa, pp: 19-22.
- Rao PP, Chakrayerthi GACSK, B. Balakrishna (2011) Application of Casting Simulation for Sand Casting of a Crusher Plate. International Journal of Thermal Technologies.
- Guharaja S, Noorul HA, Karuppannan KM (2006) Optimization of Green Sand Casting Process Parameters by Using Taguchis Method. The International Journal of Advanced Manufacturing Technology 30: 1040-1048.
- Jun HH (2006) Application of the software ProCAST in the casting of solidification simulation. Journal of Materials Science & Technology 14: 292-295.
- Chokkalingam B, Mohamed NSS (2009) Analysis of Casting Defect Through Defect Diagnostic Study Approach. Journal of Engineering Annals of the Faculty of Engineering Hunedoara 2: 209-212.
- Trovant M, Argyropouos SA (1996) Mathematical modeling and experimental measurement of shrinkage in the casting of metals. Canadian Metallurgical Quarterly 35: 77-84.
- Bride M (1999) Numerical simulation of incompressible flow driven by density variation during phase change. International. Journal for Numerical Methods in Fluids 3.1: 787-800.
- Kim C, Ro S (1993) Shrinkage Formation during the Solidification Process in an Open Rectangular Cavity. ASME Journal of Heat Transfer 115: 1078-1081.
- Andreause U, Delliisola F (1996) On Thermo kinematic Analysis of Pipe Shaping in Cast Ingot: A Numerical Simulation via FDM. International Journal of Engineering Science 34: 1349-1367.
- Paszndiden-Fard M, Chandra S, Mostaghimi J (2002) A Three Dimensional Model of Droplet Impact and Solidification. Journal of Heat Mass Transfer 45: 2229-2242.

Citation: Brezov D (2018) Optimization and Gimbal Lock Control via Shifted Decomposition of Rotations. J Appl Computat Math 7: 410. DOI: 10.4172/2168-9679.1000410

Copyright: © 2018 Brezov D. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.