Sinusoidal Trajectory Generation Methods for Spacecraft Feedforward Control

Abstract

The following is a brief walkthrough of material related to the modeling of spacecraft dynamics with feedforward control as the self-awareness declaration for deterministic artificial intelligence. Specifically, the focus will be on the analysis of various sinusoidal trajectory methods. The methods utilized are the basic MATLAB sine generation function, a Taylor series implementation, and two alternate algorithms for higher speed, lower precision and lower speed, higher precision implementations. The chapter features a brief summary of previous work investigating the impact of step size on Euler and Body angles. This is followed by a high level overview of Euler angle theory, quaternions, direction cosine matrices, kinematics, and dynamics to form a mathematical basis for the core material. With the numerical basis for the modeling efforts outlined, the results of running a SIMULINK model of spacecraft dynamics with feedforward control will be briefly analyzed and explored. The analysis will cover the impacts of varying step size with various sinusoidal trajectory generation methodologies.

1. Introduction

The study of spacecraft rotational mechanics includes three core functional areas: kinematics, dynamics, and disturbances. This chapter will be primarily focused on the trajectory generation methodology used to drive the space vehicle dynamics. In the model shown in Figure 1, a commanded spacecraft body movement is taken in and translated into both an input torque (T) and an angular velocity (ω) for the spacecraft body. A quaternion and direction cosine matrix are then calculated and employed to find Euler Angles for the spacecraft’s attitude with respect to an inertial frame.

Figure 1.

Simulink model topology.

Previously unpublished academic work, based on the topology found in Figure 1, delved into the effects of varying time steps and maneuver time in feedforward control (everything to the left of dynamics) but considered only a single method of sinusoid generation. Eq. (1) is the idealized feedforward control for the simulation and it can additionally function as the self-awareness declaration for deterministic artificial intelligence [1, 2, 3, 4, 5, 6, 7, 8, 9]. In this chapter, the Torque generator driven by Eq. (1) will remain unchanged but we will investigate the impact of four different sinusoid generation methods in the Trajectory generator, namely:

• MATLAB’s sine function

• 4th order Taylor series approximation

• Low precision approximation algorithm

• High precision approximation algorithm

These methodologies require examination due to the inherent errors in an actual spacecraft’s measurement apparatus. No onboard system can ever measure a spacecraft’s angular position, velocity, or acceleration at infinite precision. Additionally, any disturbances accounted for in the coarse control from feedforward design implementations will still require additional feedback control for robust, fine tuning. As such, it may be more prudent to embrace a small amount of course control error in order to speed up overall computation time. In order to assess these methodologies, a SIMULINK model for each Trajectory generator method was created and their outputs were compared for error and computation time.

$$\left[\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right]=\left[\begin{array}{c}{J}_{\mathit{xx}}{\dot{\omega }}_x+J_{\mathit{xy}}{\dot{\omega }}_y+J_{\mathit{xz}}{\dot{\omega }}_z-J_{\mathit{xy}}{\omega }_x{\omega }_z-J_{\mathit{yy}}{\omega }_y{\omega }_z-J_{\mathit{yz}}{\omega }_z^2+J_{\mathit{xz}}{\omega }_x{\omega }_y+J_{\mathit{zz}}{\omega }_z{\omega }_y+J_{\mathit{yz}}{\omega }_y^2\\ J_{\mathit{yx}}{\dot{\omega }}_x+J_{\mathit{yy}}{\dot{\omega }}_y+J_{\mathit{yz}}{\dot{\omega }}_z-J_{\mathit{yz}}{\omega }_x{\omega }_y-J_{\mathit{zz}}{\omega }_x{\omega }_z-J_{\mathit{xz}}{\omega }_x^2+J_{\mathit{xx}}{\omega }_x{\omega }_z+J_{\mathit{xy}}{\omega }_z{\omega }_y+J_{\mathit{xz}}{\omega }_z^2\\ J_{\mathit{zx}}{\dot{\omega }}_x+J_{\mathit{zy}}{\dot{\omega }}_y+J_{\mathit{zz}}{\dot{\omega }}_z-J_{\mathit{xx}}{\omega }_x{\omega }_y-J_{\mathit{xz}}{\omega }_y{\omega }_z-J_{\mathit{xy}}{\omega }_y^2+J_{\mathit{yy}}{\omega }_x{\omega }_y+J_{\mathit{yz}}{\omega }_z{\omega }_x+J_{\mathit{xy}}{\omega }_x^2\end{array}\right]$$E1

2. Sinusoidal trajectory generation

In the previous work referred to in Section 1 of this chapter it was found that a sufficiently small time step size must be utilized in order to adequately model the commanded input as a sinusoid. This can be seen with a snapshot of the experimental results as shown in Figure 2.

Figure 2.

Position vs. time graphs as a function of time step.

These graphs were generated when the model shown in Figure 1 was tested with a commanded input of [0, 0, 30]^T which corresponds to zero roll, zero pitch, and a 30° yaw. This allowed analysis to focus on time step impacts to one axis of rotation of MATLAB sinusoidal generation methodology. The quiescent and post-maneuver timeframes were both set to 5 s and the maneuver was conducted over a 10 s interval. Figure 2 shows the results of a time step of 1 s and a time step of 0.01 s, both using the Runge-Kutta solver. With the larger step size of 1 s, the Euler angles lost tracking and went slightly above the commanded angle while the smaller step size of 0.01 s afforded better tracking between Euler angle and Body angle. The main takeaway was that with finer step size resolution one could essentially equate the Euler and Body angles on a given model even though the former is in the inertial reference frame and the latter is in the spacecraft body reference frame.

Additionally, 100 iterations were run over a range of step sizes between 1 and 0.01 s to see what impact a chosen step size would have on error between Body and Euler Angles at the end of the maneuver phase. The results of this study can be seen below in Figure 3.

Figure 3.

Error vs. step size.

From this graph you can clearly see the benefit of reducing step size. It should be noted, however, that once step size was reduced to 0.01119 s the model reached a point where CPU and MATLAB precision comes into play and shrinking step size did not necessarily result in reducing error. This is noticeable with the erratic nature of the curve on the left side of the graph when approaching a much small step size.

Since it was found that we could adequately model our system with a small enough step size this allowed further investigation into other areas within the model which could be changed and optimized. This line of thought led us to create a new model which could be utilized to compare and contrast various sinusoidal generation methodologies. This new model’s numerical basis and simulation results will be shown in the following two subsections.

2.1 Model creation

Before adding changes for various sinusoidal generation methodologies the mathematical basis and assumptions used in the previous model were re-verified. The first step used in model creation (and re-verification) was the implementation of a direction cosine matrix to numerically represent rotations about a set of axes to project a starting frame onto a desired reference frame in order to outline the system kinematics [10, 11, 12, 13, 14]. Figures 4 and 5 provide visual depictions of the process driven by Eq. (2) through Eq. (4). Each direction cosine matrix equation takes an axial rotation as depicted by the series of rotations in Figure 4 and represents it as an orthonormal matrix consisting of sines, cosines, zeroes, and ones per the trigonometry rules shown in Figure 5.

Figure 4.

Simple 3-2-1 rotation.

Figure 5.

Depiction of a Φ rotation about the X axis.

$$1\mathit{\text{RotationDCM}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathit{cos}\mathrm{\Phi }& \mathit{sin}\mathrm{\Phi }\\ 0& -\mathit{sin}\mathrm{\Phi }& \mathit{cos}\mathrm{\Phi }\end{array}\right]$$E2

$$2\mathit{\text{RotationDCM}}=\left[\begin{array}{ccc}\mathit{cos}\mathrm{\theta }& 0& -\mathit{sin}\mathrm{\theta }\\ 0& 1& 0\\ \mathit{sin}\mathrm{\theta }& 0& \mathit{cos}\mathrm{\theta }\end{array}\right]$$E3

$$3\mathit{\text{RotationDCM}}=\left[\begin{array}{ccc}\mathit{cos}\mathrm{\Psi }& \mathit{sin}\mathrm{\Psi }& 0\\ -\mathit{sin}\mathrm{\Psi }& \mathit{cos}\mathrm{\Psi }& 0\\ 0& 0& 1\end{array}\right]$$E4

For more complicated movements, direction cosine matrices are multiplied together resulting in a more intricate matrix. Eq. (5) is an example of this for a 3-2-1 Rotation. This direction cosine matrix sequence and others are provided a more in depth treatment in Refs. [14, 15, 16] and as well as other chapters in this book. Note that in Eq. (5), the “C” and “S” characters are utilized as shorthand for the sine and cosine trigonometric functions.

$$321\mathit{DCM}\mathit{\text{Rotation}}=\left[\begin{array}{ccc}\mathit{C\theta C}\mathrm{\Psi }& \mathit{C\Phi S}\mathrm{\Psi }& -\mathit{S\theta }\\ \mathit{\text{SΦSθC}}\mathrm{\Psi }-\mathit{C\Phi S}\mathrm{\Psi }& \mathit{\text{SΦSθS}}\mathrm{\Psi }+\mathit{C\Phi C}\mathrm{\Psi }& \mathit{\text{SΦCθ}}\\ \mathit{\text{CΦSθC}}\mathrm{\Psi }+\mathit{S\Phi S}\mathrm{\Psi }& \mathit{\text{CΦSθS}}\mathrm{\Psi }-\mathit{S\Phi }\mathrm{C}0& \mathit{\text{CΦCθ}}\end{array}\right]$$E5

The other major workhorse of the model’s kinematics is the orthonormal quaternion matrix [17, 18, 19], which accomplishes the same feat as the direction cosine matrix but with intrinsic divide by zero protection. The particular quaternion calculation format implemented in the model can be seen in Eq. (6).

$$\dot{q}=\left[\begin{array}{ccc}0& \mathrm{\omega }3& -\mathrm{\omega }2\\ -\mathrm{\omega }3& 0& \mathrm{\omega }1\\ \begin{array}{c}\mathrm{\omega }2\\ -\mathrm{\omega }1\end{array}& \begin{array}{c}-\mathrm{\omega }1\\ -\mathrm{\omega }2\end{array}& \begin{array}{c}0\\ -\mathrm{\omega }3\end{array}\end{array}\begin{array}{c}\mathrm{\omega }1\\ \mathrm{\omega }2\\ \begin{array}{c}\mathrm{\omega }3\\ 0\end{array}\end{array}\right]\left[\begin{array}{c}q1\\ q2\\ \begin{array}{c}q3\\ q4\end{array}\end{array}\right]$$E6

To find the quaternion from Eq. (6), it is a simple matter of initializing the quaternion vectors with an orthonormal set (i.e., [0,0,0,1]^T) and integrating the output.

The final matrix utilized, shown in Eq. (7), is an equation equivalent to the previously shown 3-2-1 rotational matrix Eq. (5) but written in terms of quaternion elements. This quaternion based 3-2-1 rotation is depicted below in Eq. (7).

$$321\mathit{\text{DCMRotation}}=\left[\begin{array}{ccc}1-2\left({\mathrm{q}}_2^2+{\mathrm{q}}_3^2\right)& 2\left({\mathrm{q}}_1{\mathrm{q}}_2+{\mathrm{q}}_3{\mathrm{q}}_4\right)& 2\left({\mathrm{q}}_1{\mathrm{q}}_3-{\mathrm{q}}_2{\mathrm{q}}_4\right)\\ 2\left({\mathrm{q}}_1{\mathrm{q}}_2-{\mathrm{q}}_3{\mathrm{q}}_4\right)& 1-2\left({\mathrm{q}}_1^2+{\mathrm{q}}_3^2\right)& 2\left({\mathrm{q}}_2{\mathrm{q}}_3+{\mathrm{q}}_1{\mathrm{q}}_4\right)\\ 2\left({\mathrm{q}}_1{\mathrm{q}}_3+{\mathrm{q}}_2{\mathrm{q}}_4\right)& 2\left({\mathrm{q}}_3{\mathrm{q}}_2-{\mathrm{q}}_1{\mathrm{q}}_4\right)& 1-2\left({\mathrm{q}}_1^2+{\mathrm{q}}_2^2\right)\end{array}\right]$$E7

In actual implementation, the SIMULINK model uses Eq. (7) but for calculation purposes takes advantage of Eq. (5). Combining terms leaves us with Eq. (8) through Eq. (10) for Euler Angles:

$$\mathrm{\theta }={\mathit{sin}}^{-1}\left(1-2\left({\mathrm{q}}_2^2+{\mathrm{q}}_3^2\right)\right)$$E8

$$\mathrm{\Phi }=\text{atan}2\left(2\left({\mathrm{q}}_2{\mathrm{q}}_3+{\mathrm{q}}_1{\mathrm{q}}_4\right)/1-2\left({\mathrm{q}}_1^2+{\mathrm{q}}_2^2\right)\right)$$E9

$$\mathrm{\Psi }=\text{atan}2\left(2\left({\mathrm{q}}_1{\mathrm{q}}_2+{\mathrm{q}}_3{\mathrm{q}}_4\right)/1-2\left({\mathrm{q}}_2^2+{\mathrm{q}}_3^2\right)\right)$$E10

Referring back to Figure 1, the Dynamics block is driven by Eq. (11) (the Euler moment equation) when a torque is provided by the Feedforward Control section’s Torque generator.

$$\sum T=\dot{H}=J\dot{\omega }+\omega \times \mathit{J\omega }$$E11

To derive the angular velocity of the body (ω_BODY) fed from dynamics to kinematics, we multiply moment of inertia (J) by the time derivative of angular velocity, i.e., angular acceleration ($\dot{\mathrm{\omega }}$). This J$\dot{\mathrm{\omega }}$ term is multiplied by the inverse of J and then integrated as per Eq. (12).

$$\int \left(J^{-1}\ast J\dot{\omega }\right)\mathit{dt}={\omega }_{\mathit{\text{BODY}}}$$E12

Afterwards, the ω_BODY is used to obtain the spacecraft attitude’s Euler angles.

The feedforward control blocks in Figure 1 contain the Trajectory generator and the Torque generator. The Trajectory generator takes in a commanded body angle and then uses a sine wave in order to approximate the commanded maneuver. To help elaborate on this point, Figure 6 shows a square wave and a sine wave, both shifted up in amplitude such that they operate from zero to two instead of between negative and positive one.

Figure 6.

Visual depiction of square and sine waves.

We can see that if a square wave was used to approximate the commanded maneuver, the spacecraft would essentially be expected to transition from the zero position up to its max amplitude at the two positions instantaneously. This is a physical impossibility. With the sine wave in the same figure we can see that the spacecraft machinery is afforded a ramp up and slow down period along with a relatively linearly slope in the middle. Basically, the half period of sinusoid within the square pulse of Figure 6 is used to provide a smooth, achievable input to the system.

Our original model’s Trajectory generator follows Eq. (13) through Eq. (15) to approximate the commanded maneuver for angular position, angular velocity (ω), and angular acceleration ($\dot{\mathrm{\omega }}$) via the MATLAB sine wave function.

$$\mathit{\text{AngularPosition}}=\frac{1}{2}(A+A\ast \mathit{sin}\left(\left(\frac{\pi }{\mathit{\Delta t}}\right)\left(t-t_{\mathit{\text{wait}}}\right)-\frac{\pi }{2}\right)$$E13

$$\mathit{\text{AngularVelocity}}=\frac{1}{2}A\ast \left(\frac{\pi }{\mathit{\Delta t}}\right)\mathit{cos}\left(\left(\frac{\pi }{\mathit{\Delta t}}\right)\left(t-t_{\mathit{\text{wait}}}\right)-\frac{\pi }{2}\right)$$E14

$$\mathit{\text{AngularAcceleration}}=-\frac{1}{2}A\ast {\left(\frac{\pi }{\mathit{\Delta t}}\right)}^2\mathit{sin}\left(\left(\frac{\pi }{\mathit{\Delta t}}\right)\left(t-t_{\mathit{\text{wait}}}\right)-\frac{\pi }{2}\right)$$E15

Eq. (13) models the input command, where “A” is the maneuver’s commanded angle. The base frequency of the sinusoid (${\omega }_{\mathit{\text{sine}}}$) is $\left(\frac{\mathrm{\pi }}{\mathrm{\Delta }t}\right)$ where $\left(\mathrm{\Delta }t\right)$ is the desired maneuver time. The $t_{\mathit{\text{wait}}}$ term allows for a quiescent period and $-\frac{\mathrm{\pi }}{2}$ term allows for a proper phase shift to implement the sinusoidal half period of Figure 4. The effect of this can be seen later on in Figure 7. Eqs. (14) and (15) are just successive derivatives of Eq. (13) used to generate angular velocity and acceleration which are fed into Eq. (1) in the Torque generator from Figure 1. This produces an output torque which drives the dynamics.

Figure 7.

Feedforward control for 30° yaw maneuver.

From Eq. (13) through Eq. (15) it can be clearly seen that the argument of the sine and cosine terms always follows the form of Eq. (16). We can use this to implement our Taylor series and the other two other algorithms on equal footing.

$$\mathit{Arg}=\left(\frac{\pi }{\mathit{\Delta t}}\right)\left(t-t_{\mathit{\text{wait}}}\right)-\frac{\pi }{2}$$E16

The Taylor series, as detailed in Ref. [20], is a numerical method that can be used to approximate other functions. In our model we substitute the sine and cosine in Eq. (13) through Eq. (15) with Eqs. (17) and (18).

$$\mathit{\text{Taylor}}\mathit{Sin}=\mathit{Arg}-\frac{{\mathit{Arg}}^3}{3!}-\frac{{\mathit{Arg}}^5}{5!}-\frac{{\mathit{Arg}}^7}{7!}$$E17

$$\mathit{\text{Taylor}}\mathit{Cos}=1-\frac{{\mathit{Arg}}^2}{2!}-\frac{{\mathit{Arg}}^4}{4!}-\frac{{\mathit{Arg}}^6}{6!}$$E18

The Taylor series is a power series and additional terms could be included ad infinitum for greater precision; however, initial testing found that four series terms provided a reasonable approximation while maintaining a viable runtime. Additionally, it was found that pre-calculating the factorial terms sped up computation time.

The last set of sinusoidal generation methodologies tested was the low precision (LP) and high precision (HP) algorithms. These were found to be used in many applications, especially ones which limited by lower processing power such as mobile gaming. Refs. [21, 22] provided the baseline for adaptation of the baseline Eqs. (19) and (20). Note that the same equation is used for sine and cosine except that the argument term is given a positive $\frac{\pi }{2}$ phase shift when applied to the cosine in Eq. (14). In Eq. (19), $\left(+/-\right)$indicates that a plus is used if Arg is less than zero and a minus is used otherwise.

$$\mathit{LP}=1.27323954\ast \mathit{Arg}\left(+/-\right)0.405284735\ast {\mathit{Arg}}^2$$E19

$$\mathit{HP}=.225\ast \left[\mathit{LP}\ast \mathit{abs}\left(\mathit{LP}\right)-\mathit{LP}\right]+\mathit{LP}$$E20

Eq. (20) provides additional smoothing to Eq. (16) at the cost of computation time to implement a high precision mode. This equation could be implemented such that “$\mathit{LP}\ast \mathit{abs}\left(\mathit{LP}\right)$” essentially become the magnitude of ${\mathit{LP}}^2$. While this analysis used “if then” statements, there may be more efficient ways to implement Eqs. (19) and (20) depending on the software package being utilized.

2.2 Model simulation and analysis

For the purpose of model verification, commands are initially held constant at zero resulting in zero torque to allow us to evaluate model operation during a quiescent period. With zero input torque, Eq. (11) shows us that the change in angular momentum should be zero as well; therefore the model should not change. After a 5 s quiescent period, a 30° yaw maneuver was input into each model and the maneuver was conducted over a 10 s period at which point the commanded body angle went to zero (indicating no more change required). Afterwards the output Euler angles should remain constant with the spacecraft’s new attitude. The results of this test were plotted in Figure 7. Additionally, note that the model was setup with an arbitrary diagonalized inertia matrix (J) per Eq. (21). The roll and pitch Euler angles remained at zero, as they should have, and as such are not shown.

$$J=\left[\begin{array}{ccc}150& 0& 0\\ 0& 90& 0\\ 0& 0& 35\end{array}\right]$$E21

From a qualitative standpoint you can see that each sine generation algorithm created the intended Euler angle movement shown earlier in Figure 4. Upon closer examination the only algorithm that appeared to lose tracking on the body angle was the low precision algorithm. Since we achieved the basic shape from each method, we now turn our attention to analyzing the error and computation time associated with each for a given time step interval.

In the previous work addressed in the beginning of Section 2 of this we found that for a Runge-Kutta solver with step size of 0.01119, the MATLAB sine function reached a point where CPU and MATLAB precision started to affect the gains offered by shrinking step size. With this in mind, step sizes of 0.1, 0.05, and 0.1 s were chosen for analysis, results of which can be seen in Figure 8 and Table 1.

Figure 8.

Computation time results.

Method	Step size		Step size		Step size
	0.01 s		0.5 s		0.1 s
	Time(s)	Error	Time(s)	Error	Time(s)	Error
1-MATLAB	0.2423	1.02E-9	0.2459	6.36E-11	0.2685	1.02E-13
2-Taylor	0.2434	3.53E-5	0.2471	3.53E-5	0.2720	3.53E-5
3-LP	0.2416	2.81E-2	0.2459	2.81E-2	0.2682	2.81E-2
4-HP	0.2417	3.25E-4	0.2456	3.25E-4	0.2666	3.25E-4

Table 1.

Timing and error results.

Time and error results averaged over 500 iterations per model.

One major inference to be drawn from Table 1 is that the normalized error columns show Taylor, LP, and HP algorithms are roughly step-size invariant within our test range; only the MATLAB function’s error reduces with step size in our model. Interestingly, the LP and HP functions can be more accurate than Taylor but we are more concerned with the final angle in coarse control (corresponding to the peak of sine curve) rather than the intermediate steps—so in this case, the Taylor series is considered more accurate.

The other major item of note is the difference in average computation time. While the MATLAB command is ultimately the most accurate at all step sizes, we can see that the LP and HP algorithms generally run faster. It should be noted that results were relatively inconsistent when run below 200 iterations but at 500 iterations they stabilized with HP and LP usually running faster during a given test. Since HP would run faster than LP sometimes and vice versa, even at greater iterations; it is believed that this may come down to specific CPU architecture used and how a dual core processor handles calculations.

3. Conclusions

When it comes to sinusoidal trajectory generation, if your spacecraft design has accurate fine control (e.g., robust feedback mechanisms) and can handle coarse control error on the order of 10⁻⁴, then it may be more prudent to utilize the HP algorithm. This would allow you to reduce processing power and time resulting in lower power requirements and faster onboard calculations. A fourth order Taylor series would not be beneficial due to longer computation times and less accuracy than the MATLAB function; extending the order of the series would increase accuracy but also extend its runtime, thereby, losing viability. Additionally, simulation results reveal that the fastest methodology may vary with CPU architecture indicating that there may not be a definitive answer for “best” trajectory generation method. This should be evaluated within the specific spacecraft design trade space. Very little information is available on the proprietary MATLAB sine calculations but future testing will investigate the LP and HP algorithms against MATLAB on a variety of CPU architectures and organic sinusoidal generation on other platforms such as PYTHON.