Spacecraft Engineering with Julia - Attitude Dynamics

Using Julia to model the attitude dynamics of spacecraft

Introduction

As of writing this, I just finished up my final year of studying aerospace engineering. I enrolled in the program's spacecraft stream, where one of the courses I chose for my final semester was spacecraft attitude dynamics and control. This has become one of the most interesting (and hardest) classes I've taken, and I want to write a bit about the stuff I learned in that course. Specifically, how to model a spacecraft's attitude and some ways to control it. I’m also interested in getting into GNC (guidance, navigation, and control) so writing these kinds of articles is good to brush up my skills and apply some knowledge.

The full code associated with this article can be found at this link: Spacecraft Attitude Control with Julia Part 1 - Attitude Dynamics.ipynb

Spacecraft Attitude

To begin, what even is attitude? Attitude is essentially the same as orientation, so when I'm talking about a spacecraft's attitude, I'm talking about how a spacecraft is orientated in space relative to some reference frame. Spacecraft attitude dynamics describes how a spacecraft's attitude changes over time under the influence of forces/torques.

To describe a spacecraft's attitude dynamics, we have Euler's equation:

Where M is the torques acting on the spacecraft, J is the spacecraft's inertia matrix, and ω is the spacecraft's angular velocity.

Assuming principal axis rotations, this one equation can be expanded into three separate differential equations for each axis (x, y, and z):

These equations describe how the angular velocity change according to the moments applied to the spacecraft.

To model the spacecraft attitude, a state-space representation is used. The state variable for this case is:

As can be seen, we are using quaternions in this case to represent the attitude of the spacecraft. Quaternions have many benefits, however I will not get into too much detail in this post regarding that (I've attempted using Euler angles, however they did not work well unfortunately).

For the derivative of the state variable, we can rearrange the equations above to get equations for the angular acceleration, and some quaternion kinematic equations to get the derivatives of each quaternion value:

For more details regarding the theory and derivation of these equations, Bong Wie’s Space Vehicle Dynamics and Control is a good textbook that dives deeper into spacecraft attitude dynamics.

Implementation in Julia

To implement attitude dynamics in Julia, the DifferentialEquations.jl package is used, which provides a number of differential equation solvers and an intuitive method of defining systems of differential equations. The following function is used to solve this system of differential equations:

function diffeq_euler(initial_conditions, time_span, params; solver_args...)
    function _differential_system(u, p, t)
        q0, q1, q2, q3, ω1, ω2, ω3 = u
        J1, J2, J3, M1, M2, M3 = p

        return SA[
            -0.5 * (q1*ω1 + q2*ω2 + q3*ω3),
            0.5 * (q0*ω1 + q2*ω3 - q3*ω2),
            0.5 * (q0*ω2 - q1*ω3 + q3*ω1),
            0.5 * (q0*ω3 + q1*ω2 - q2*ω1),
            (M1(u, p, t) + (J2 - J3) * ω2 * ω3) / J1,
            (M2(u, p, t) + (J3 - J1) * ω1 * ω3) / J2,
            (M3(u, p, t) + (J1 - J2) * ω1 * ω2) / J3
        ]
    end

    problem = ODEProblem(_differential_system, initial_conditions, time_span, params)
    solution = solve(problem; solver_args...)

    return solution
end

A couple of things to note with this function. Firstly, it acts as an interface with the DiffEq functions, without needing to call them directly. It allows for inputting the initial conditions, time span, and parameters, and returns the solution of the inputted values (the solution being the returned object of the solve function from DifferentialEquations.jl). For this function, I am also defining the returned vector from the internal _differential_system function as a static array from the StaticArrays package, which actually provides a significant increase in performance.

Another thing to note are the M variables, which are not actually constant values, instead they are functions that take in the same inputs as the parent function itself. These M functions are the functions that represent the torque applied on the spacecraft, and declaring them as functions allows us to input any arbitrary type of behavior for our spacecraft model. This lets us use different control laws without changing the diffeq_euler function, rather we just need to code the corresponding torque functions. I will go over coding the control laws more in depth in the next part of this series of articles.

Just to simply test the diffeq_euler function, the following torque functions are defined:

M0(u, p, t) = 0
M_const(u, p, t) = 0.1

These torque functions are simple to understand, one outputs no torque, the other outputs 0.1 Nm of torque.

Next, the arrays for the initial conditions, time span, and parameters are defined:

init = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
tspan = (0.0, 5*60.0)
params = [250.0, 110.0, 110.0, M_const, M0, M0]

With these arrays defined, the result of this example can be calculated by simply calling the diffeq_euler function:

sol = diffeq_euler(init, tspan, params)

The sol variable contains all of the results from the defined problem. It is a unique type of object, and more information can be found on the DiffEq Solution Handling Docs. With the solution variable, we can now plot the results:

Well, it’s something right? The problem with quaternions of course is that their values are hard to understand without converting to a different attitude representation. Instead, we can convert the quaternion values to Euler angles, which give a more intuitive result:

As can be seen, this result is much easier to understand. Above the torques applied to the spacecraft are all zero except for in the x axis, resulting in rotation about the x axis. To make the conversions from quaternions to Euler angles, I used the ReferenceFrameRotations.jl package, provides a number of attitude representations and associated functions.

Conclusion

In conclusion, using a simple function allows us to model the attitude dynamics of a spacecraft. Using DifferentialEquations.jl gives a convenient method of solving the system of differential equations associated with attitude dynamics. In the next article, I aim to introduce some basic attitude control.