By Michal Jagodzinski - May 10th, 2023
Hello and welcome back to Star Coffee. We return to more spacecraft GN&C related topics, this time using
ModelingToolkitStandardLibrary.jl for acausal spacecraft attitude dynamics and control modelling. For a great overview of the differences between causal and acausal modelling, check out this talk from Dr. Chris Rackauckas: Causal vs Acausal Modeling By Example: Why Julia ModelingToolkit.jl Scales.
using CairoMakie, AlgebraOfGraphics, Latexify using ModelingToolkit, ModelingToolkitStandardLibrary using DifferentialEquations set_aog_theme!() t const Rot = ModelingToolkitStandardLibrary.Mechanical.Rotational const B = ModelingToolkitStandardLibrary.Blocks
ModelingToolkitStandardLibrary.jl is a standard library of reusable components used for acausal modelling with
ModelingToolkit.jl. It has plenty of components from various types of systems, including mechanical, electrical, and thermal systems. Using this library allows for easily creating complex models from the smaller parts.
Let's do a simple example to illustrate how to build systems with
ModelingToolkitStandardLibrary.jl. Let's simulate a small torque being applied to a body with some inertia. We define the components used in this system as:
0.1) torque = Rot.Torque() inertia = Rot.Inertia(J=100.0)torque_input = B.Constant(k=
Next we need to connect the components to form the system. The block diagram for the system is below:
In this diagram, each output from a component is connected to an input on another component, e.g.,
torque_input.output connects to
torque.tau. Components have distinct connectors, which themselves are other components, and these connector objects are used to connect different components together. The connections for this system can be defined as:
eqs = [ connect(torque_input.output, torque.tau), connect(torque.flange, inertia.flange_a) ]
Next we can create the
ODESystem object and simplify the system equations:
model = ODESystem( eqs, t; systems = [torque_input, torque, inertia] ) sys = structural_simplify(model)
As can be seen, using
structural_simplify reduces the system of equations down to just two. For reference, below is the full system of equations:
Finally we can define the
ODEProblem, solve it, and plot the results:
prob = ODEProblem(sys, , (0, 5*60.0), ) sol = solve(prob, Tsit5()) fig1 = Figure() ax1 = Axis(fig1[1,1], xlabel="Time (s)", ylabel="Angle (rad)") lines!(ax1, 0:1:5*60, sol(0:1:5*60)[inertia.phi]) fig1
Next, let's define a custom component to model the attitude dynamics of a spacecraft. This will be a linear model, utilizing the
Inertia component we used for the previous example. To create the more accurate nonlinear model, we would need to define the system of differential equations ourselves.
This component is very similar to the previous system, except it now has three distinct axes, each with their own inertias.
function LinearSpacecraftAttitude( ; name, Jx=100.0, Jy=100.0, Jz=100.0, u0=zeros(3), ω0=zeros(3), ω̇0=zeros(3) ) Ix = Rot.Inertia( J=Jx, phi_start=u0, w_start=ω0, a_start=ω̇0 ) Iy = Rot.Inertia( J=Jy, phi_start=u0, w_start=ω0, a_start=ω̇0 ) Iz = Rot.Inertia( J=Jz, phi_start=u0, w_start=ω0, a_start=ω̇0 ) x_flange_a = Rot.Flange() y_flange_a = Rot.Flange() z_flange_a = Rot.Flange() x_flange_b = Rot.Flange() y_flange_b = Rot.Flange() z_flange_b = Rot.Flange() ϕ_sensor = Rot.AngleSensor() θ_sensor = Rot.AngleSensor() ψ_sensor = Rot.AngleSensor() ps = Jx=Jx Jy=Jy Jz=Jz u0=u0 ω0=ω0 ω̇0=ω̇0 eqs = [ connect(x_flange_a, Ix.flange_a), connect(y_flange_a, Iy.flange_a), connect(z_flange_a, Iz.flange_a), connect(Ix.flange_b, x_flange_b), connect(Iy.flange_b, y_flange_b), connect(Iz.flange_b, z_flange_b), connect(x_flange_b, ϕ_sensor.flange), connect(y_flange_b, θ_sensor.flange), connect(z_flange_b, ψ_sensor.flange), ] compose( ODESystem(eqs, t, , ps; name = name), Ix, Iy, Iz, x_flange_a, y_flange_a, z_flange_a, x_flange_b, y_flange_b, z_flange_b, ϕ_sensor, θ_sensor, ψ_sensor ) end
With our custom component defined, let's recreate the previous example, but this time using our custom component we just defined.
spacecraft = LinearSpacecraftAttitude()
Next we can reuse the previously defined components to create this simulation. This simulation is essentially the same as the previous:
Creating the connections:
simple_eqs = [ connect(torque_input.output, torque.tau), connect(torque.flange, spacecraft.x_flange_a) ]
Defining the model and running the simulation:
0, 5*60.0), ) simple_sol = solve(simple_prob, Tsit5())simple_model = ODESystem(simple_eqs, t; systems = [ torque_input, torque, spacecraft ]) simple_sys = structural_simplify(simple_model) simple_prob = ODEProblem(simple_sys, , (
It's as simple as that. We can now visualize the angle:
times = 0:1:5*60 ϕ_vec = simple_sol(times)[spacecraft.Ix.phi] function handle_angle(angle) if abs(angle) > pi a = angle while a > pi a -= 2*pi end return a else return angle end end fig2 = Figure() ax2 = Axis(fig2[1,1], xlabel="Time (s)", ylabel="Angle (°)") lines!(ax2, times, rad2deg.(handle_angle.(ϕ_vec))) fig2
Looks pretty good to me. This is pretty much the same solution I got in a previous article of mine, so it's good to see my results should be correct.
Alright, now let's do a more complex example. Since we're using
ModelingToolkit.jl for acausal modelling, building upon our model to create more complex simulations is incredibly easy. From the
ModelingToolkit.jl standard library we get PID controllers for free, no additional setup is required. I'm sure as the ecosystem matures, we'll get things like Kalman filters and other useful components as well.
Let's now define our next spacecraft attitude simulation, this time using PID controllers to control all three axes. Here is the block diagram for this simulation setup:
This setup should be able to drive the spacecraft's attitude from an initial condition to a final setpoint. If you've taken a course on control systems this should look pretty familiar, it's just a simple feedback system. Let's define the components for this simulation:
150.0, Jy=100.0, Jz=100.0, u0=[0.5, 0.25, -0.5] ) setpoint = B.Constant(k=0.0) feedback_ϕ = B.Feedback() feedback_θ = B.Feedback() feedback_ψ = B.Feedback() ctrl_ϕ = B.PID(k=10.0, Td=32.0, Ti=100) torque_ϕ = Rot.Torque() ctrl_θ = B.PID(k=10.0, Td=32.0, Ti=100) torque_θ = Rot.Torque() ctrl_ψ = B.PID(k=10.0, Td=32.0, Ti=100) torque_ψ = Rot.Torque()ctrl_spacecraft = LinearSpacecraftAttitude( Jx=
Now let's define the connections:
ctrl_eqs = [ connect(setpoint.output, feedback_ϕ.input1), connect(setpoint.output, feedback_θ.input1), connect(setpoint.output, feedback_ψ.input1), connect(feedback_ϕ.output, ctrl_ϕ.err_input), connect(ctrl_ϕ.ctr_output, torque_ϕ.tau), connect(torque_ϕ.flange, ctrl_spacecraft.x_flange_a), connect(ctrl_spacecraft.ϕ_sensor.phi, feedback_ϕ.input2), connect(feedback_θ.output, ctrl_θ.err_input), connect(ctrl_θ.ctr_output, torque_θ.tau), connect(torque_θ.flange, ctrl_spacecraft.y_flange_a), connect(ctrl_spacecraft.θ_sensor.phi, feedback_θ.input2), connect(feedback_ψ.output, ctrl_ψ.err_input), connect(ctrl_ψ.ctr_output, torque_ψ.tau), connect(torque_ψ.flange, ctrl_spacecraft.z_flange_a), connect(ctrl_spacecraft.ψ_sensor.phi, feedback_ψ.input2), ]
Defining the connections is a little tedious, but it's not a hard task. GUI tools such as ModelingToolkitDesigner.jl are already in development, and I'm excited to see how this and other tools built on
ModelingToolkit.jl evolve. If we are able to get an open-source version of Simulink implemented in Julia I'll be incredibly happy.
Let's now simulate this system and plot the resulting angles over time of each axis:
0, 2.5), ) ctrl_sol = solve(ctrl_prob, Tsit5())ctrl_model = ODESystem(ctrl_eqs, t; systems = [ ctrl_spacecraft, setpoint, feedback_ϕ, ctrl_ϕ, torque_ϕ, feedback_θ, ctrl_θ, torque_θ, feedback_ψ, ctrl_ψ, torque_ψ, ]) ctrl_sys = structural_simplify(ctrl_model) ctrl_prob = ODEProblem(ctrl_sys, , (
times_ctrl = 0:0.01:2.5 ctrl_sol_interp = ctrl_sol(times_ctrl) fig3 = Figure() ax3 = Axis(fig3[1,1], xlabel="Time (s)", ylabel="Angle (°)") lines!(ax3, times_ctrl, rad2deg.(ctrl_sol_interp[ctrl_spacecraft.Ix.phi]), label="ϕ (Roll)") lines!(ax3, times_ctrl, rad2deg.(ctrl_sol_interp[ctrl_spacecraft.Iy.phi]), label="θ (Pitch)") lines!(ax3, times_ctrl, rad2deg.(ctrl_sol_interp[ctrl_spacecraft.Iz.phi]), label="ψ (Yaw)") hlines!(ax3, [0.0]; label="Setpoint", linestyle=:dash) axislegend(ax3) fig3
Looks good. As can be seen, the controllers successfully changed the attitude of the spacecraft to the setpoint. Obviously the PID controllers could be tuned to have better responses, but it's a good start. In addition to the spacecraft's attitude over time, we can also visualize the various states in our simulation. For example, let's see the torque inputs to the spacecraft from the controllers:
fig4 = Figure() ax4 = Axis(fig4[1,1], xlabel="Time (s)", ylabel="Torque (N m)") lines!( ax4, times_ctrl, ctrl_sol_interp[ctrl_spacecraft.x_flange_a.tau], label="ϕ Controller Output" ) lines!( ax4, times_ctrl, ctrl_sol_interp[ctrl_spacecraft.y_flange_a.tau], label="θ Controller Output" ) lines!( ax4, times_ctrl, ctrl_sol_interp[ctrl_spacecraft.z_flange_a.tau], label="ψ Controller Output" ) axislegend(ax4) fig4
Sick. We get these additional states for analysis for free again. I did not code in this state or anything, I simply accessed the states from the
ctrl_spacecraft component. Specific analysis points in the system can also be explicitly defined to conduct these kinds of analyses (see this documentation for more info).
Thanks for reading! Hope you learned something new and got inspired to give acausal modelling a try. I will continue learning about using
ModelingToolkit.jl for acausal modelling and I am planning on creating these kinds of models within SAT. More posts coming soon. Until next time.