Closed-form quaternion representations for rigid body rotation: application to error assessment in orientation algorithms of strapdown inertial navigation systems

Closed-form analytical representations of the rigid body orientation quaternion, angular velocity vector and the external moment vector satisfying kinematic equations and equations of motion are derived. In order to analyze errors of orientation algorithms for strapdown inertial navigation systems, reference models for specific rigid body rotation cases are formulated. Based on solutions, analytical expressions for ideal signals of angular velocity sensors in the form of quasi-coordinates are derived. For several sets of parameters, numerical implementations of the reference models are performed and trajectories in the configuration space of orientation parameters are presented. Numerical analysis of the drift error for the third-order orientation algorithm is performed. The results show that the value of the accumulated drift error using the derived two-frequency models exceeds the value of the accumulated drift error in the conventional case of a regular precession.

where ω ω ω(t) is the angular velocity vector and θ θ θ n is the apparent rotation vector. A number of algorithms for determination of orientation quaternions have been developed and discussed [6,7,13,22,[25][26][27]35]. The approach, presented in [7], is based on the decomposition of the particular solution of the quaternionic kinematic equation in the Picard series by approximating of the apparent rotation vector θ θ θ n within an interval [t n−1 , t n ] by a fourth-order polynomial of finite differences. In [7], a series of algorithms for the determination of the rotation quaternion is developed, which have different orders up to the fourth inclusive. Another approach is based on the Taylor time series decomposition of sine and cosine functions of half true rotation angle in the representation of the orientation quaternion components and the use the rotation vector as an intermediate parameter [22]. This approach was first suggested in [6]. It was further developed in [7,13,22,25], among others.
The error estimation of orientation algorithms is usually based on the use of special benchmark tests, for which the current orientation and components of the angular velocity vector of a rigid body are presented in a closed analytical form. To analyze the accuracy of orientation algorithms, an analytical reference model of the conic motion is applied in [22]. The analytical SPIN-CONE model presented in [27] is widely used as a benchmark test.
Since for benchmark tests the local or accumulated error is presented in analytical form, it can be used to optimize the orientation algorithm for a specific test motion. Optimization of orientation algorithms for conic motion by special tuning of the algorithm coefficients without changing its structure was proposed in [22]. Optimization is based on obtaining an analytical expression for the algorithm error in the form of a stepped series with further determination of unknown coefficients, based on the condition of cutting off the highest order terms of the series. Another method, which implements essentially the same approach, was presented in [13]. The method of optimizing algorithms for regular precession and conic motion which is based on minimizing asymptotic estimations of computational drift error was presented in [25]. The papers [8][9][10]12,20,21,23,[31][32][33] discuss improved methods for optimization of orientation algorithms and present the results of research on algorithm optimization under conditions of generalized conical motion, regular precession, random angular motion, and practical issues of effective use of algorithms, including for a specific structure of SINS. Analytical reference models that differ from cases of regular precession and conic motion are presented in [14,15]. The improvements of orientation algorithms in conditions of vibrational conic motion are discussed in [17,18,37]. The issues of amplitude and frequency expansion of the application of orientation algorithms are presented in [28][29][30].
Since conic motion and regular precession are rather specific cases of angular motion of the base of SINS, it is of practical interest to analyze the orientation algorithms for a wide class of test motions. In this paper, new analytical two-frequency solutions to the equations of rotation motion of a rigid body are formulated, and reference models of rotation-based solutions are presented. The paper is organized as follows. In order to fix the notation, Sect. 2 presents classical governing equations of rotation motion of rigid bodies as well as different representations of orientation. The attention is placed on the quaternion representation. The relationships between the rotation tensor and quaternion components are discussed. Section 3 presents analytical solutions for several cases of rigid body motion. Special cases of the derived solutions are discussed in detail. Based on solutions, analytical expressions for ideal signals of angular velocity sensors in the form of quasi-coordinates are derived in Sect. 4. Numerical implementations of the reference models for several sets of parameters are discussed in Sect. 5. Furthermore, trajectories in the configuration space of orientation parameters are presented. Numerical analysis of the drift error for the third-order orientation algorithm is performed in Sect. 6. Finally, conclusions are drawn with respect to the value of the accumulated drift error based on the derived two-frequency models if compared to the conventional case of regular precession, discussed in the literature.

Governing equations and problem statement
In order to fix our notation, let us recall the governing equations describing rotation of a rigid body. To compute the orientation of the rigid body, the following Darboux problem should be solved [1][2][3]38] P P P = × P P P = P P P × ω ω ω, P P P(0) = P P P 0 (2) where P P P is the rotation tensor, P P P · · · P P P T = I I I , detP P P = 1, is the spatial angular velocity vector, also known as left angular velocity [38] , ω ω ω = P P P T · · · is the body angular velocity vector, also known as right angular velocity, P P P 0 is the initial condition and I I I is the second rank unit tensor.
In this paper, the direct tensor calculus in the sense of Gibbs [36] and Lagally [16] is applied. A secondorder tensor is a finite sum of dyads of vectors A A A = a a a ⊗ b b b + c c c ⊗ d d d + · · · + e e e ⊗ f f f . The basic operations for dyads can be summarized as follows: Operations (3)-(8) are generalized for finite sums of dyads and tetrads. The direct tensor calculus is widely used in continuum mechanics and rheology, see for example [5,19,24]. The angular momentum L L L of the rigid body with respect to the center of mass is defined as follows: where J J J ref is the tensor of inertia in the reference configuration and J J J is the corresponding tensor in the actual configuration defined as follows ⊗ e e e 1 + J 2 e e e 2 ⊗ e e e 2 + J 3 e e e 3 ⊗ e e e 3 , J J J = J 1 n n n 1 ⊗ n n n 1 + J 2 n n n 2 ⊗ n n n 2 + J 3 n n n 3 ⊗ n n n 3 n n n i = P P P · · · e e e i where J i , i = 1, 2, 3 are principal moments of inertia and the orthogonal unit vectors e e e i and n n n i are principal directions connected with the rigid body in the reference and actual configurations, respectively. The kinetic energy K is computed as follows: 2K = ω ω ω · · · J J J ref · · · ω ω ω = · · · J J J · · · (10) The rotational motion of a rigid body around its center of mass is described by the following vector equation: where M M M is the resultant moment vector acting on the body and 0 is the initial condition for the left angular velocity. With respect to the right angular velocity vector, Eq. (11) is formulated as follows: where M M M ref = P P P T · · · M M M. The scalar products of Eq. (12) with the vectors e e e i provide the following dynamic Euler equations where M i = M M M · · · n n n i are components of the resultant moment vector, ω i = ω ω ω · · · e e e i are components of the right angular velocity vector and (1,2,3) is the symbol of circular permutation of indices. In order to compute the orientation of a rigid body, equations of motion (11) or (12) together with Eq. (2) by taking into account the initial conditions must be solved. The key step for the efficient solution is the suitable representation of the rotation tensor. According to the Euler theorem, it can be represented as follows: where the unit vector m m m denotes the axis of rotation and ϕ is the angle of rotation about this axis. In many cases, the rotation tensor can be represented as a composition of three rotations (14) about three different axes [38]. Examples include the classical Euler representation and yaw-pitch-roll composition. Although compositions of three rotations provide clear interpretations of complex angular motions and closed-form solutions can be derived in many cases by specifying special axes of rotation [38], this representation suffers from drawbacks. First, nonlinear implicit relations between the angles of rotation exist [38]. Second, for any composition of three rotations, gimbal locks (the loss of one or two degrees of freedom) are unavoidable, leading to difficulties in numerical solution procedures [11]. To avoid these drawbacks, quaternion representation of orientation is widely in development of orientation algorithms for SINS. To introduce the orientation quartenion, let us consider the following abbreviations Equation (14) takes the form The angular velocity vectors corresponding to the rotation tensor (15) are computed as follows: With λ 0 and the components of the vector λ λ λ, λ i = λ λ λ · · · e e e i the orientation quartenion is introduced as follows In addition, the following kinematic quaternion equation can be deriveḋ is the angular velocity quaternion, and • is the quaternion multiplication symbol. For operations of the quartenion algebra we refer to [4,34]. Once the components of the orientation quartenion are given, the rotation tensor can be computed from Eq. (15). Vice versa, for the given rotation tensor P P P, λ 0 and λ λ λ can be computed from the following relations: From Eqs. (9), (10) and (18) 2 , the following integral can be derived [38] 2K = L L L · · · P P P · · · J J J −1 ref · · · P P P T · · · L L L = const (19) In general, the rotation tensor can be parameterized with three angles of rotation about three fixed axes.
In the case of free rotation, the relationship between the angles exists according to Eq. (19). Closed form two-parametric solutions for free rotations are presented in [38].
In what follows, we derive two-frequency solutions of the system of Eqs. (13) and (17) and formulate reference models of rigid body rotation in the form of analytical expressions for quasi-coordinates where the apparent rotation vector θ is defined by Eq. (1). Furthermore, we present analytical expressions for the components of the orientation quaternion (t), which corresponds to the angular motions. For the numerical implementation of the reference model, it is necessary to estimate the accumulated orientation error obtained by the algorithm for calculating the third-order orientation quaternions. The technique for constructing analytical reference models is based on a special representation of the orientation quaternion by trigonometric functions of the angles of rotations, which simultaneously change in time. This representation provides automatic fulfillment of the normalization condition for the quaternion (t) = 1, corresponding to Eq. (15) 2 .

First solution
Let us specify the components of the orientation quaternion as follows: where k α and k β are frequencies, k α = k β . η and ξ are parameters satisfying the condition η 2 + ξ 2 = 1. For the given orientation quaternion, the angular velocity quaternion is computed from the inverse kinematic equation where˜ (t) is the conjugate orientation quaternion. The components of the angular velocity vector are The proposed analytical solutions (21) and (23) to Eqs. (13) and (17) satisfy the initial In this case, the first and second components of the angular velocity vector change according to a harmonic law, and the third component is constant, as is the case with the classical regular precession of a transversely isotropic rigid body with (I 1 = I 2 = I 3 ). Indeed, in the case of (I 1 = I 2 = I 3 ), the magnitude of the angular momentum vector L is constant, i.e., The kinetic energy K is constant too However, these conditions generally do not correspond to the free rotation of a rigid body. Indeed, the proposed solutions (21) and (23) describes a more general case of rotational motion. To show this, we will solve the inverse problem of dynamics, that is, we will find the external moments that provide the motion of a rigid body with components of angular velocity (23). From Eq. (13), we obtain Therefore, in the general case of an asymmetric rigid body, its rotation under the action of moments (24) is not free. Since the kinetic energy is conserved, the power of external moments must be zero, i.e., One may verify that the moment vector with the components (24) is orthogonal to the right angular velocity vector with the components (23). For a transversally isotropic body in the case when I 1 = I 2 , we find the following components of the external moment vector that supports regular precession It is obvious that only under the conditions , that is, the motion of a rigid body with components of the angular velocity vector (23) is the classical regular precession. In this case, the parameter ξ must take the following value Consider the case of the solution (21) when η = 0, ξ = ±1. Then, we have a plane-free rotational motion with a constant component of the angular velocity along the third principal axis, i.e., The orientation quaternion takes the following form: and (0) = [0, 0, 1, 0] Consider the case of solution (21) when ξ = 0 and η = 1. The angular velocity vector has the following components: under the action of the moment vector The orientation quaternion (21) then takes the following form: For the magnitude of the angular momentum vector, we obtain L 2 = 4k 2 α (I 2 1 cos 2 (2k β t) + I 2 2 sin 2 (2k β t) + 4I 2 3 k 2 β ) = const Therefore, the regular precession is not observed. For a transversally isotropic rigid body (I 1 = I 2 ), we have that is, in this case of a forced regular precession.

Second solution
Let us formulate the analytical solution to the kinematic equation (17) as follows: where k α and k β , k α = k β are frequencies. The parameters η and ξ must satisfy the condition η 2 + ξ 2 = 1.
The corresponding solution to the inverse kinematic equation (22) takes the form Let us note that the first and second components of the angular velocity vector contain constant parts, while the third component changes according to the harmonic law and does not depend on the parameters ξ and η. This solution satisfies the initial conditions (0) = [η, 0, 0, ξ], ω i (0) = [2ξ k α + 2ηk β , 2ηk α − 2ξ k β , 0]. Since none of the components of the angular velocity vector (29) in the general case of the parameters ξ and η is constant, the rotation is significantly different from the cases of regular precession and conical motion.
Let us find the moment vector under the action of which the motion of the rigid body has the components of the angular velocity vector (29). According to equations of motion, we obtain Consider a particular case of solution (28) for η = 0 and ξ = 1. The orientation quaternion takes the following form: For the angular velocity vector, we obtain This rotary motion is provided by the action of the moment vector with the components The initial conditions are . For a transversely isotropic rigid body with I 1 = I 3 , we obtain that L 2 = 4I 2 1 k 2 α + 4I 2 2 k 2 β = const and 2K = 4I 1 k 2 α + 4I 2 k 2 β = const. Therefore, this motion is a forced regular precession. If in the expression for the model quaternion (28) the parameters η = 1 and ξ = 0 are specified, the components of the orientation quaternion are The corresponding angular velocity vector has the following components: The rotational motion with the angular velocity vector (32) is provided by the action of an external moment with the following components The initial conditions are (0) = [1, 0, 0, 0], ω i (0) = [2k β , 2k α , 0]. For transversally isotropic rigid body with I 2 = I 3 , this motion is a forced regular precession, since the magnitude of the angular momentum vector has the form The kinetic energy is computed as follows: The moment vector is not zero vector. Again, one may verify that the moment vector is orthogonal to the angular velocity vector.

Third solution
Let us define the orientation quaternion as a solution to the kinematic equation (17) as follows: where k α and k β , k α = k β are frequencies. ν and μ are parameters satisfying the condition ν 2 +μ 2 = 1. We find that according to Eq. (22), this quaternion corresponds to the following angular velocity vector components The initial conditions are In this case, the vector components do not depend on the parameters μ and ν. Therefore, there is no relationship between (0) and ω i (0), as it was in the cases of the first and second two-frequency solutions. Note that the angular velocity vector (35), which corresponds to quaternion (34), is exactly the same as obtained earlier by (32) for the case of quaternion (28) with the parameters η = 1 and ξ = 0. Rotational motion with the orientation (34) and the angular velocity vector components (35) is provided by the action of an external moment with components (33). Since for μ = 1 and ν = 0, the components of the quaternion (35) coincide with the components of the quaternion (31), the third two-frequency solution can be considered as an independent case of integrability in quadratures to the equations of motion, which is valid for arbitrary values of the parameters μ and ν in Eq. (34), different from the case μ = 1 and ν = 0. The considered solution is a special case of forced regular precession according to Eqs. (28) and (29) for η = 1 and ξ = 0.

Fourth solution
The solution to Eq. (17) is represented in the following form The components of the angular velocity vector which correspond to the quaternion with components (36), according to Eq. (22) take the following form: The initial conditions of motion are This solution describes a motion with the following time-variable magnitude of the angular velocity vector If we set k β = 0 in (36), then (37) provides ω i = [0, 0, 2k α ], which corresponds to the plane rotation with constant velocity around the third axis and orientation quaternion (t) = [cos(k α t), 0, 0, sin(k α t)].
Let us consider in more detail the case of a transversally isotropic rigid body I 2 = I 3 under the conditions when k α = 0. From (37), we obtain ω i (t) = [2k β , 2k β cos(2k β t), −2k β sin(2k β t)], and the orientation Since the external moment vector depends on time M i (t) = [0, −4k 2 2 I 1 sin(2k 2 t), −4k 2 2 I 1 cos(2k 2 t)] = const, the motion can be interpreted as forced regular precession around the first axis. In the general case of a rigid body, we will obviously have a more complex motion.
where t = t n − t n−1 .

Second reference model
This model includes analytical solution for the orientation quaternion (28) with given values of η, ξ , k α and k β . The analytical expression for quasi-coordinates takes the following form:

Third reference model
This model is based on the quaternion representation (34), when the numerical values of the parameters, and frequencies, are given. The reference values of the quasi-coordinates are defined by following analytical expressions:

Fourth reference model
This reference rotation model is based on Eq. (36). The quasi-coordinates are computed as follows:

Numerical implementation of two-frequency reference rotation models
The proposed two-frequency models are evaluated numerically for specific values of the frequencies and  Fig. 8. For this implementation, the magnitude of the angular velocity vector of a rigid body is not constant, but its average value over the time interval was 0.583 rad/s. Figure 9 shows the trajectories in the configuration space of the orientation parameters obtained for the model of regular precession on the time interval with the initial conditions (0) = [1, 0, 0, 0] and ω i (0) = [0, 15, −0, 356, 0, 437] rad/s. For this implementation, the constant magnitude of the angular velocity vector with the value 0.583 rad/s was assumed. The trajectories in the configuration space of orientation parameters for the case of regular precession have features typical for a regular motion including symmetry, repeatability, and specific form. Comparison of the results based on the reference models of two-frequency solutions with the results of the regular precession model allows us to conclude that the trajectories obtained for the proposed two-frequency models in the general case of the frequency values and parameters differ significantly from the trajectories obtained for the regular precession, and also are more complex.

Accuracy analysis of the third-order orientation algorithm based on reference models
Let us apply the developed implementations of the two-frequency reference models to obtain an estimate of the accuracy of the third-order algorithm, for which the rotation quaternion n = [ λ n:0 , λ λ λ n ] within the where θ 2 n = θ θ θ n · · ·θ θ θ n . The final orientation of the rigid body the time t n is calculated using the following formula: where n = [θ θ θ n ], n−1 = (θ θ θ n−1 ). To evaluate the accuracy of the algorithm, let us define the fatal orientation error-the accumulated small angle of rotation (drift) δ of the calculated tripple of axes relative to its true position, which is specified by the reference model. To this end, we apply the method for determining the drift as proposed in [7]. , 437] rad/s, for which the condition |ω ω ω(t)| = 0.583 rad/s is also satisfied. Analysis of the given dependences in Fig. 10 allows us to draw the following conclusions -The drift error under the same magnitude of the angular velocity vector for all test motions increases with time, but its rate of increase for the derived reference models exceeds the rate of increase in the case of a regular precession -It is found that the worst test motion from the point of view of the drift error for the selected algorithm is the reference model based on the first two-frequency solution -The test motion in the form of a regular precession of a rigid body is not the worst case for analyzing the errors of the SINS orientation algorithms.

Conclusions
Two-frequency solutions of the system of dynamic and kinematic equations of rotation of a rigid body are derived and test motions on their basis in the form of analytical reference models are developed. The numerical implementation of the reference models and the analysis of the obtained trajectories in the configuration space of the orientation parameters show that the corresponding motions differ significantly from the case of regular precession and are more complex. The implementations of the reference models are used as test motions to estimate the drift error of the selected third-order orientation algorithm for SINS. It was found that the drift error, which occurs in the case of the proposed models, exceeds the drift error, which is obtained for the regular precession with the same magnitude of the angular velocity vector. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by Projekt DEAL.