Mean and standard deviation for all 12 rotations, using a 0.1 step size.
Motion mechanics (dynamics) comprises kinetics to describe the implications of applied forces and torques; and also kinematics (phoronomics). Developed in the 1700s, kinematics describes mathematical translations from one basis of measurement to another using common kinematic measurement variables like quaternions, Euler angles, and direction cosine matrices. Two ubiquitous rotation sequences are unquestionably adopted for developing modern direction cosine matrices from among the 12 potential options, stemming from applicability to aerospace systems, accuracy, and computation burden. This chapter provides a comprehensive reevaluation of all 12 options yielding a menu of options for accuracy and computational burdens, with the results illustrated compared to the ubiquitous two modernly adopted choices, broken into two rotational groups: symmetric rotations and nonsymmetric rotations. Validation will be provided by critical analysis of integration using step size to illustrate correlated minimal accuracy. No single rotational sequence is universally superior with respect to all figures of merit, enabling trade-space analysis between rotational sequences. One interesting revelation of one of the two ubiquitous sequences (the 3-1-3 symmetric sequence) is illustrated to have relatively less accuracy but lower computational burden than the other (the 3-2-1 nonsymmetric sequence). Meanwhile, a relatively unknown “2-3-1” rotational sequence is shown to have similar computational burden and accuracy.
- direction cosines
- Euler angles
- space dynamics
- digital computation
- control systems
- control engineering
In 1775, Leonhard Euler developed motion phoronomics  which immediately blossomed in the next two centuries [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. The space race between the now-defunct Soviet Union and the United State of the last century gave substantial impetus to development and adoption of motion kinematics together with survival imperatives driven by the nuclear cold war. The resultant lineage of literature contains seemingly countless technical and non-technical [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62]. With this heritage the two most common rotational sequences used to calculate direction cosine matrices are referred to as “aerospace” sequences for nonsymmetric sequences (where the resulting angles are referred to as Tait-Bryan angles), while the symmetric sequences are oft referred to as “orbital” sequences (where the resulting angles are called proper “Euler Angles”) .
In light of continued improvement in computational capabilities, the focus of this research is to evaluate all 12 rotation sequences comparing by mean and standard deviation of accuracy reflecting roll, pitch, and yaw angles; and also comparing by computational burden embodies by time required to perform the calculations. The chapter questions whether the 3-2-1 rotational sequence truly the best with respect to either of these figures of merit (statistical accuracy and computational burden). The results illustrate the two standard sequences are indeed good (with relative weaknesses). In particular the standard asymmetric sequence is more accurate, but slower than the standard symmetric sequence. On average, despite fewer mathematical steps, the symmetric rotations are on average slower to calculate. The 3-2-1 sequence is quickest to calculate amongst the asymmetric rotations, meanwhile the 2-3-2 and 1-2-1 are the quickest amongst the symmetric rotations.
Artificial intelligence and machine learning has evidenced the need for rapid calculations, so as motion mechanics incorporate adopt these new learning algorithms, the impact of this chapter become increasingly relevant in that options revealed in here illustrate simultaneous accuracy and favorable rapidity of calculation . This chapter also complements other algorithmic advances [37, 38, 39, 40, 41, 42, 43, 44, 45] like system identification [55, 56, 57, 58, 59] including nonlinear adaptive forms and also control [46, 47, 48, 49, 50, 51, 52, 53, 54] for space guidance, navigation, and control (GNC) missions [35, 36, 60, 61, 62, 63, 64, 65] in a time when the United States is developing and relying upon more advanced Machine Learning and AI products than ever before.
2. Materials and methods
One application of motion mechanics is the control of the attitude of spacecraft rotational maneuvers or even maintenance of a specified attitude. The key reminder is that Euler’s moment equations governing rotational movement apply in a non-moving reference frame referred to as “inertial,” which is a reference frame that has no meaning as a basis for measurement (i.e., it is not possible to identify a truly non-moving reference frame that can be used for measurement of angular position). This section of the chapter illustrates the method to numerically evaluate the options for kinematic expressions of rotations between chosen frames of reference (e.g., the body frame) and the inertial frame. MATLAB/SIMULINK depicted in Figure 1 is used to create a simulation of rotation of spacecraft and necessary components of the simulation to make it relatively high-fidelity include aerodynamic and gravity gradient disturbances; kinematic expressions including quaternions, direction cosine matrices, and Euler angles; and even incorporation of the motion of an object in a specified orbit. The simulation is elaborated in Section 2, while Section 3 will describe the experiments, with concluding results in Section 4.
2.1 Theory of dynamics
Mechanics and dynamics are synonyms. Interestingly, kinematics (which is also currently called phoronomics ) was referred to as “statics” in the era of Newton  indicating the lack of motion which is included in kinetics. This will be elaborated further in Section 2.1.2. Michael Chasle’s theorems permit us to simply “invoke” Euler’s moment equation to describe three-degrees of rotation and Newton’s law to describe three-degrees of translation; together comprising a full mathematical description of so-called 6DOF motion, or motion in six-degrees. Euler describes rotational motion expressed in a moving body frame as , where [J] is a matrix of mass moments of inertia explained by Kane . Measurements of rotational maneuvers are expressed in inertial coordinates by establishing an arbitrarily placed inertial reference frame [XI, YI, ZI], while kinematics relate the inertial coordinates to those expressed in the body reference frame [XB, YB, ZB]. References in the literature use the nomenclature “direction cosine matrix” , since the matrix is composed of projection components, where the dot-product projection operation is defined by the cosine of the angle between the two reference frames [17, 25, 26]. Individual vector components elaborate the orientation angle between reference frames .
Kinetics, or Dynamics, is the process of describing the motion of objects with focus on the forces involved. In the inertial frame, Newton’s F = ma is applied but becomes Euler’s when rotation is added, where is expressed in the inertial reference frame’s coordinates, while from above is still measured in the inertial frame, but expressed in body coordinates.
Combining the Euler and Newton equations, we can account for all six degrees of freedom. In application, when an input angle [φd, θd, ψd] is commanded, the feedforward control uses Eq. (1) as the ideal controller with Eq. (2) as the sinusoidal trajectory to calculate the required torque [Tx, Ty, Tz] necessary to achieve the desired input angle. The Dynamics calculator then uses Eq. (3) to convert the torques into ωB values, where ωB is defined as the angular velocity of the body. In order to calculate this, the non-diagonal terms in Eq. (4) are neglected, removing coupled motion and leaving only the principle moments of inertia. Then, the inertia matrix J is removed from , and the remaining is integrated into [ωx, ωy, ωz], which is fed into the Kinematics block of the model to finally determine the outputted Euler Angles.
2.1.2 Kinematics, phoronomics, or “The Laws of Going”
Formulation of spacecraft attitude dynamics and control problems involves considerations of kinematics, especially as it pertains to the orientation of a rigid body that is in rotational motion. The subject of kinematics is mathematical in nature, because it does not involve any forces associated with motion. The kinematic representation of the orientation of one reference frame relative to another reference can also be expressed by introducing the time-dependence of Euler Angles. The so-called body-axis rotations involve successive rotations three times about the axes of the rotated body-fixed reference frame resulting in 12 possible sets of Euler angles. The so-called space-axis rotations instead involve three successive rotations using axes fixed in the inertial frame of reference, again producing 12 possible sets of Euler angles. Because the body-axis and space-axis rotations are intimately related, only 12 Euler angle possibilities need be investigated; and the 12 sets from the body-axis sequence are typically used . Consider a rigid body fixed at a stationary point whose inertia ellipsoid at the origin is an ellipsoid of revolution whose center of gravity lies on the axis of symmetry. Rotation around the axis of symmetry does not change the Lagrangian function, so there must-exist a first integral which is a projection of an angular momentum vector onto the axis of symmetry. Three coordinates in the configuration space special orthogonal group (3) may be used to form a local coordinate system, and these coordinates are called the Euler angles.
Key tools of kinematics from which the Euler angles may be derived include direction cosines which describe orientation of the body set of axes relative to an external set of axes. Euler’s angles may be defined by the following set of rotations: “rotation about x axis by angle and θ, rotation about z′ axis by an angle ψ, then rotation about the original z-axis by angle φ”. Eulerian angles have several “conventions: Goldstein uses  the “x-convention”: z-rotation followed by x′ rotation, followed by z′ rotation (essentially a 3-1-3 sequence). Quantum mechanics, nuclear physics, and particle physics the “y-convention” is used: essentially a 3-2-3 rotation. Both of these have drawbacks, that the primed coordinate system is only slightly different than the unprimed system, such that, φ and ψ become indistinguishable, since their respective axes of rotation (z and z′) are nearly coincident. The so-called Tait-Bryan convention in Figure 2 therefore gets around this problem by making each of the three rotations about different axes: (essentially a 3-2-1 sequence) .
Kinematics is the process of describing the motion of objects without focus on the forces involved. The [ωx, ωy, ωz] values from the Dynamics are fed into the Quaternion Calculator where Eqs. (5) and (6) yield q, the Quaternion vector. The Quaternions define the Euler axis in three dimensional space using [q1, q2, q3]. About this axis, a single angle of rotation [q4] can resolve an object aligned in reference frame A into reference frame B. The Direction Cosine Matrix (DCM) then relates the input ω values to the Euler Angles using one of 12 permutations of possible rotation sequences, where multiple rotations can be made in sequence. Therefore, the rows of the DCM show the axes of Frame A represented in Frame B, the columns show the axes of Frame B represented in Frame A, and φ, θ, and ψ are the angles of rotation that must occur in each axis sequentially to rotate from orientation A to orientation B, turning CA to CB. Figure 2 depicts a 3-2-1 sequence to rotate from CA to CB, where the Euler Axis is annotated by the thickest line.
2.1.3 The orbital frame
In order to more completely represent a maneuvering spacecraft, orbital motion must be included with the Kinematics. This relationship is represented in Figure 1, where the output of the DCM is fed into the Orbital Frame Calculator, and the second column of the DCM is multiplied against the orbital velocity of the spacecraft. The second column of the DCM represents the Y axis of Frame B projected in the X, Y, and Z axes of Frame A. This yields ωNO, the orbital velocity relative to the Inertial Frame. Using Eq. (7), this velocity is removed from the velocity of the body relative to the Inertial Frame, leaving only the velocity of the body relative to the Orbital Frame for further calculations.
Multiple disturbance torques exist that effect the motion of a spacecraft in orbit, two of which are addressed in this paper. The first is the disturbance due to gravity acting upon an object in orbit, where the force due to gravity decreases as the distance between objects increases. The force is applied as a scaling factor to the mass distribution around the Z axis of a spacecraft. This force applied to a mass offset from the center of gravity is calculated through the cross product found in Eq. (8) and yields an output torque about the Z axis.
The second disturbance is an aerodynamic torque due to the force of the atmosphere acting upon a spacecraft, which also decreases as the altitude increases. In Eq. (9), the force due to air resistance is calculated by scaling the direction of orbital velocity by the atmospheric density, drag coefficient, and magnitude of orbital velocity. This force then acts upon the center of pressure, which is offset from the center of gravity, and yields a torque about the Z axis, due to the cross product in Eq. (9).
The disturbances are additive and act upon the dynamics in Figure 1. Because the ideal feedforward controller is the dynamics, an offsetting component equal to the negative anticipated disturbances can be used to negate the disturbance torque. This results in nullifying the disturbances when the two are summed to produce ωOB, the velocity of the body relative to the Inertial Frame.
2.2 Experimental setup
A model of the 12 DCM to Euler Angle rotations was implemented in Matlab and Simulink for this experiment. A [30, 0, 0] maneuver was commanded in the [φ, θ, ψ] channels, respectively. The expected runtime of each scenario was 15 s, comprised of a 5 s quiescent period, a 5 s maneuver time, and a 5 s post maneuver period for observations. The maneuver was initiated using a sinusoidal trajectory, calculated with and .
The simulated spacecraft had an inertia matrix of J = [2, 0.1, 0.1; 0.1, 2, 0.1; 0.1, 0.1, 2], the torque was initialized as T = [0, 0, 0], and the quaternion was initialized as q = [0, 0, 0, 1]. The spacecraft was simulated to fly at an altitude of 150 km, and received a drag coefficient of 2.5. For this experiment, both orbital motion and torque disturbances were turned off in order to simply the simulation.
The Matlab and Simulink models utilized the Runge-Kutta solver, with an ode4 back-end. Multiple step sizes were tested to determine accuracy variations for each of the rotations: 0.1, 0.001, and 0.0001 s. The trigonometric function used to mathematically solve for the Euler Angles was the atan2 function in Matlab.
Three figures of Merit were used to assess performance. The first two were the mean and standard deviation between the Euler Angles and Body Angles. The third was the calculation time for each rotation as a measure of complexity.
3. Experimental results and analysis
3.1 Euler angle calculations and post-processing
A relationship like Eq. (6) was created mathematically to relate the DCM and rotation matrices for each of the 12 rotation sequences. Then, φ, θ, and ψ were solved for, resulting in a mathematical process to determine the Euler Angles. This process was then coded in Matlab and Simulink, but the process was not perfect. Trigonometric quadrant errors caused the appearance of discontinuities from a commanded [30, 0, 0] maneuver. This artifact was resolved using post processing and further refinement of the DCM to rotation matrix derivations that correlated the 12 rotations found in Figure 3 into two groups of six rotations: symmetric and non-symmetric rotations. To further define the groups, symmetric rotations would be 1-2-1 or 2-3-2, while non-symmetric rotations would include 1-3-2 or 3-1-2 rotations.
3.2 Euler angle to body angle accuracy
Accuracy was measured in the experiment by measuring the difference between the Body Angles and output Euler Angles. The expectation was that a perfectly accurate system would have a difference of zero. Figure 4 depicts the deviation over time and Table 1 provides the associated mean values and standard deviations for each of the rotations.
The six non-symmetric rotations show consistent error in φ, and only begin to deviate beyond the fifth decimal place in both mean error and standard deviation. While φ is commanded to change to 30°, θ, and ψ are expected to remain at zero, but show non-zero values due to error incurred by step size.
The six symmetric rotations are substantially harder to draw conclusions from because of the uncorrelated rotations. The mean error and standard deviation values are drastically different from each other in Table 1 and visibly deviate in Figure 4. Therefore, further correlation is required to analyze accuracy. Table 1 values were calculated over the 15 s simulation time, noting that some sequences had not reached steady-state values making their error values even larger compared to others in Table 1 if the simulations had been run until steady state was reached.
3.3 Step size versus accuracy
This experiment implemented a variable step size to determine the accuracy delta resultant from the different step sizes. Figure 5 depicts analysis using a step sizes of 0.001 s, which can be compared against Figure 4, which used a 0.1 s step size. The primary difference between Figures 4 and 5 is the 2-order of magnitude increase of accuracy accompanying the two order of magnitude reduction in step size. A further reduction to a time step of 0.0001 s was made, with an additional order of magnitude increase in accuracy. Further reductions below this required more time than was feasible, but the trend holds that decreasing the step size increases accuracy. Furthermore, the relative accuracies between rotations held when the step sizes decreases, meaning the 1-3-2 and 3-1-2 rotations remained the most accurate.
3.4 DCM to Euler angle timing
All 12 rotation scenarios executed a maneuver within 5 s, with a standard pre and post maneuver observation period. However, actual runtimes sometimes exceeded this 15 s period; this is attributed to the complexity of the calculations and additional processes that were running at the time of the simulation. The results of each of the 12 rotations for each of three time steps are shown in Table 2. The simulation timing is effected by step size; therefore, the results can only be compared between different rotations (vertically in the Table 2) and not between step sizes (horizontally in Table 2); however, relative comparisons between step sizes are valid.
|Execution time (s)|
|DCM||0.1 step size||0.001 step size||0.0001 step size|
Three observations can be made from the results in Table 2. The first is that the slowest rotation is the 1-2-3 rotation, by a significant amount depending upon the step size. The second is that on average, non-symmetric rotations were faster than symmetric rotations. This result is unique because the same algorithm with the same number of mathematical steps yielded different execution times. Lastly, the fastest overall rotation was the 2-3-2 rotation, with 3-2-1 as the fastest non-symmetric rotation.
This chapter on modern kinematics or motion phoronomics elaborated all 12 possible instantiations of direction cosine matrices with comparisons of numerical accuracy representing how accurately the chosen Euler angle represents the roll, pitch, and yaw expressions of rotations about x, y, z axes respectively. Additionally, comparison is made by using the figure of merit of computation burden expressed in time-necessary to perform calculations using each respective kinematic instantiation. The results were listed in a large table of options available for trade-offs, where symmetric sequences proved more difficult to compare and correlate, meanwhile the non-symmetric rotational sequences proved easier to correlate to roll, pitch, and yaw due to the ease of allocating independent angles.
The “trade-space” of options is a key elaboration, since none of the options were unanimously best using more than one figure of merit. If accuracy measured by mean error is most relatively important, 1-2-3, 1-3-2, 2-1-3, 2-3-1, 3-1-2, or 3-2-1 rotational sequences best represent roll, while the ubiquitous 3-2-1 sequence cannot best to represent pitch, where the 1-2-3 sequence is superior; while 2-1-3, 2-3-1, and 3-2-1 rotational sequences can most accurately reflect yaw. Instead if accuracy measured by standard deviation of errors was most important, the results were not identical. The most computationally efficient rotational sequence was the 2-3-2 rotation, while the 3-1-3 and 1-2-1 performed next in the list of best options. Oddly, the ubiquitous 3-2-1 sequence was merely the fifth fastest option.
The demonstration of relative inferiority of the standard ubiquitous options is a key novel development in the chapter, and the novelties were validated using a relatively high fidelity simulation of spacecraft attitude dynamics, but the novel development are valid for other forms of rotational motion mechanics like naval vessels, airplanes, and even robotics.
Future works will validate these results on laboratory spacecraft hardware simulators at the Naval Postgraduate School, and if successful flight in space is available on the international space station making the technology available to enhance the aforementioned applications of the technology [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45].
We would like to thank our teacher, Dr. Timothy Sands for his guidance in developing this work, as well as our families for their support while we spent our time away from them developing this research, so thank you.
Conflict of interest
The authors declare no conflict of interest.