2.1.1. Kinetics

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 *T* = *J* .*ω* when rotation is added, where *T* = *J* .*ω* is expressed in the inertial reference frame's coordinates, while *T* = *J* .*ω* + *ω* × *Jω* 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 (1) as the ideal controller with (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 (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 (4) are neglected, removing coupled motion and leaving only the principle moments of inertia. Then, the inertia matrix *J* is removed from *J* .*ω*, and the remaining .*ω* is integrated into [*<sup>ω</sup>x*, *<sup>ω</sup>y*, *<sup>ω</sup>z*], which is fed into the Kinematics block of the model to finally determine the outputted Euler Angles.

$$T\_d = J\dot{\omega}\_d + \omega\_d \times J\omega\_d \tag{1}$$

$$\theta = \frac{1}{2} \left( A + A \sin \left( \omega\_f t + \varphi \right) \right) \tag{2}$$

$$T = \dot{H}\_{\dot{i}} = \dot{J}\dot{\omega}\_{\dot{i}} + \omega\_{\dot{i}} \times \text{J}\omega\_{\dot{i}} \tag{3}$$

#### 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 twelve 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 twelve possible sets of Euler angles. Because the body-axis and space-axis rotations are intimately related, only twelve Euler angle possibilities need be investigated; and the twelve sets from the body-axis sequence are typically used [26]. 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 [22] 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) [22].

**Figure 2.** Execution of a 3-2-1 rotation from C<sup>A</sup> to C<sup>B</sup> (left to right); blue-dotted arrows denote angle rotations. A direct rotation from C<sup>A</sup> to C<sup>B</sup> can be made about the Euler Axis, q4 in red. The set of three rotations may be depicted as four rectangular parallelepipeds, where each contains the unit vectors of the corresponding reference frame [29].

Kinematics is the process of describing the motion of objects without focus on the forces involved. The [*<sup>ω</sup>x*, *<sup>ω</sup>y*, *<sup>ω</sup>z*] values from the Dynamics are fed into the Quaternion Calculator where (5) and (6) yield *q*, the Quaternion vector. The Quaternions define the Euler axis in three dimensional space using [*q*1, *q*2, *q*3]. About this axis, a single angle of rotation [*q*4] 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 C<sup>A</sup> to CB. Figure 2 depicts a 3-2-1 sequence to rotate from C<sup>A</sup> to CB, where the Euler Axis is annotated by the thickest line.

$$\begin{bmatrix} j\_{12}\dot{\omega}\_{x} + J\_{xy}\dot{\omega}\_{y} + J\_{xz}\dot{\omega}\_{z} - J\_{yy}\omega\_{x}\omega\_{z} - J\_{yy}\omega\_{y}\omega\_{z} - J\_{yz}\omega\_{z}^{2} + J\_{xz}\omega\_{x}\omega\_{y} + J\_{zz}\omega\_{z}\omega\_{y} + J\_{yz}\omega\_{y}^{2} \\\ J\_{yz}\dot{\omega}\_{x} + J\_{yy}\dot{\omega}\_{y} + J\_{yz}\dot{\omega}\_{z} - J\_{yz}\omega\_{x}\omega\_{y} - J\_{xz}\omega\_{x}\omega\_{z} - J\_{xz}\omega\_{x}^{2} + J\_{xx}\omega\_{x}\omega\_{z} + J\_{xy}\omega\_{z}\omega\_{y} + J\_{xz}\omega\_{z}^{2} \\\ J\_{xz}\dot{\omega}\_{x} + J\_{xy}\dot{\omega}\_{y} + J\_{xz}\dot{\omega}\_{z} - J\_{xz}\omega\_{x}\omega\_{y} - J\_{xy}\omega\_{y}\omega\_{z} - J\_{xy}\omega\_{y}^{2} + J\_{yy}\omega\_{x}\omega\_{y} + J\_{yz}\omega\_{z}\omega\_{x} + J\_{xy}\omega\_{x}^{2} \end{bmatrix} = \begin{bmatrix} T\_{x} \\\ T\_{y} \\\ T\_{z} \end{bmatrix} \tag{4}$$

$$
\begin{bmatrix}
\dot{q}\_1\\ \dot{q}\_2\\ \dot{q}\_3\\ \dot{q}\_4
\end{bmatrix} = \frac{1}{2} \begin{bmatrix}
0 & \omega\_3 & -\omega\_2 & \omega\_1\\ -\omega\_3 & 0 & \omega\_1 & \omega\_2\\ \omega\_2 & -\omega\_1 & 0 & \omega\_3\\ -\omega\_1 & -\omega\_2 & -\omega\_3 & 0
\end{bmatrix} \begin{bmatrix}
q\_1\\ q\_2\\ q\_3\\ q\_4
\end{bmatrix} = \frac{1}{2} \begin{bmatrix}
q\_4 & -q\_3 & q\_2 & q\_1\\ q\_3 & q\_4 & -q\_1 & q\_2\\ -q\_2 & q\_1 & q\_4 & q\_3\\ -q\_1 & -q\_2 & -q\_3 & q\_4
\end{bmatrix} \begin{bmatrix}
w\_1\\ w\_2\\ w\_3\\ 0
\end{bmatrix} \tag{5}
$$

$$
\begin{bmatrix}
1 - 2\left(q\_2^2 + q\_3^2\right) & 2\left(q\_1q\_2 + q\_3q\_4\right) & 2\left(q\_1q\_3 - q\_2q\_4\right) \\
2\left(q\_2q\_1 - q\_3q\_4\right) & 1 - 2\left(q\_1^2 + q\_3^2\right) & 2\left(q\_2q\_3 + q\_1q\_4\right) \\
2\left(q\_3q\_1 + q\_2q\_4\right) & 2\left(q\_3q\_2 - q\_1q\_4\right) & 1 - 2\left(q\_1^2 + q\_2^2\right)
\end{bmatrix} = \begin{bmatrix}
\mathbb{C}\_2\mathbb{C}\_3 & \mathbb{C}\_2\mathbb{S}\_3 & -\mathbb{S}\_2 \\
\mathbb{S}\_1\mathbb{S}\_2\mathbb{C}\_3 - \mathbb{C}\_1\mathbb{S}\_3 & \mathbb{S}\_1\mathbb{S}\_2\mathbb{S}\_3 + \mathbb{C}\_1\mathbb{C}\_3 & \mathbb{S}\_1\mathbb{C}\_2 \\
\mathbb{C}\_1\mathbb{S}\_2\mathbb{C}\_3 - \mathbb{S}\_1\mathbb{S}\_3 & \mathbb{C}\_1\mathbb{S}\_2\mathbb{S}\_3 - \mathbb{S}\_1\mathbb{C}\_3 & \mathbb{C}\_1\mathbb{C}\_2
\end{bmatrix} \tag{6}
$$

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 *<sup>ω</sup>NO*, the orbital velocity relative to the Inertial Frame. Using (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.

$$
\omega^{OB} = \omega^{NB} - \omega^{NO} \tag{7}
$$

2.1.4. Disturbances

Multiple disturbances 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 (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 (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 (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 *<sup>ω</sup>OB*, the velocity of the body relative to the Inertial Frame.

$$T\_{\mathcal{S}} = 3\frac{\mu}{R^3} \mathbf{2} \times J\mathbf{2} \tag{8}$$

$$T\_a = \mathbb{C}\_p \times f\_a = \mathbb{C}\_p \times \left[ \left( \rho\_a V\_{\mathbb{R}}{}^2 A\_p \right) \hat{V}\_{\mathbb{R}} \right] \tag{9}$$

#### *2.2. Experimental Setup*

This experiment implemented and compared the 12 DCM to Euler Angle rotations using a variable step size. An angle of [*ϕ*, *θ*, *ψ*] = [30, 0, 0] was commanded, with the quaternion and torques initialized as *T* = [0, 0, 0] and *q* = [0, 0, 0, 1]. The spacecraft had an inertia matrix of *J* = [2, 0.1, 0.1; 0.1, 2, 0.1; 0.1, 0.1, 2]. The orbital altitude was set at 150 km with a drag coefficient of 2.5. Both orbital motion and torque disturbances were turned off.

Each simulation executed over a 5 s quiescent period, 5 s maneuver time, and 5 s post maneuver observation period, totaling 15 s. The sinusoidal trajectory was calculated to have *<sup>ω</sup>f* = *π*/2 and *ϕ* = *π*/2.

The model was built in Matlab and Simulink, where integrations were calculated using the Runge–Kutta solver (ode4) with variable time steps of 0.1, 0.001, and 0.0001 s. Euler Angles were resolved using the 12 unique DCM rotation sequences with the atan2 function.

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*

Each of the Euler angles was derived using the DCM and rotation matrices, creating a relationship like (7), but unique to each rotation. *ϕ*, *θ*, and *ψ* were isolated in this relationship as a method to calculate the Euler Angles. Once calculated, the Euler Angles were implemented in the simulation. However, when a [30, 0, 0] maneuver was commanded, discontinuities due to trigonometric quadrant error manifested. Post-processing removed the error, but yielded output rotation did not match the input command. In order to correct this, the derivations for each Euler Angle were revised to correlate six of the 12 rotations, yielding the results in Figure 3. Therefore, the rotations in Figure 3 are classified into two groups: the upper six non-symmetric rotations and lower six symmetric rotations. An example of symmetric rotations is 121, while a 132 is non-symmetric rotation. The commanded input and output maneuvers were not correlated for the 6 symmetric rotations.

**Figure 3.** Corrected Euler Angles vs time for all 12 DCM rotations.

#### *3.2. Euler Angle to Body Angle Accuracy*

The output Euler Angles are not the same as the commanded Body Angles, but measuring this delta is a method of determining accuracy. Figure 4 depicts the deviation over time and Table 1 provides the associated mean values and standard deviations for each of the rotations.

**Figure 4.** Euler and Body Angle deviation, using a 0.1 step size.


**Table 1.** Mean and standard deviation for all 12 rotations, using a 0.1 step size.

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 has been run until steady state was reached.

#### *3.3. Step Size Versus Accuracy*

The simulation step size was altered to determine the effects on accuracy for the 12 rotations. If variable step-size were used, there would be no way to assure a certain level of accuracy, thus fixed step size was utilized and iterated smaller-and-smaller until no discernable accuracy improvement is noted. Figure 5 shows the results of reducing the step size from 0.1 to 0.001 s. Figures 4 and 5 remain comparable, with the 10<sup>2</sup> order of magnitude decrease in step size yielded a comparable 10<sup>2</sup> order of magnitude increase in accuracy. When a step size of 0.0001 s was used, the accuracy increased by another order of magnitude, denoting the trend. Comparing against the accuracy of the rotations, they maintained their relative accuracy; the *132 and 312 remained the most accurate rotations when the step size decreased*, and therefore has limited to no effect.

**Figure 5.** Euler and Body Angle deviation, using a 0.001 step size.

#### *3.4. DCM to Euler Angle Timing*

The execution time of each maneuver was standardized at 15 s across all scenarios. Therefore, runtime deviations for each of the 12 rotations are attributable to the complexity of the calculations. Table 2 shows the results for three different step sizes and the runtimes for each rotation. Because the step size affected the simulation timing, comparisons were only valid between rotations of a similar step size; however, relative comparisons between step sizes were valid.

Analyzing the results yields several observations. Firstly, the 123 rotation was significantly slower than all the other rotations. Secondly, the symmetric rotations were on average slower than the non-symmetric rotations, despite the same mathematical process and number of steps to solve for the Euler Angles. Lastly, the fastest non-symmetric rotation was the 321 and the fastest symmetric was the 232, slightly faster than the 121 rotation. Taking all DCM rotations into account, the 232 rotation was the fastest.


**Table 2.** Simulation run times for all 12 Direction Cosine Matrices (DCM) rotations for a 30◦ roll maneuver, using 0.1, 0.001, and 0.0001 step sizes.
