*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 `ModelingToolkit.jl`

and `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!()
@parameters 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:

```
@named torque_input = B.Constant(k=0.1)
@named torque = Rot.Torque()
@named inertia = Rot.Inertia(J=100.0)
```

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:

```
@named model = ODESystem(
eqs, t; systems = [torque_input, torque, inertia]
)
sys = structural_simplify(model)
```

UndefVarError: `torque_input`

not defined 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:

`equations(model)`

UndefVarError: `model`

not defined
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.

```
@component function LinearSpacecraftAttitude(
; name, Jx=100.0, Jy=100.0, Jz=100.0, u0=zeros(3), ω0=zeros(3), ω̇0=zeros(3)
)
@named Ix = Rot.Inertia(
J=Jx, phi_start=u0[1], w_start=ω0[1], a_start=ω̇0[1]
)
@named Iy = Rot.Inertia(
J=Jy, phi_start=u0[2], w_start=ω0[2], a_start=ω̇0[2]
)
@named Iz = Rot.Inertia(
J=Jz, phi_start=u0[3], w_start=ω0[3], a_start=ω̇0[3]
)
@named x_flange_a = Rot.Flange()
@named y_flange_a = Rot.Flange()
@named z_flange_a = Rot.Flange()
@named x_flange_b = Rot.Flange()
@named y_flange_b = Rot.Flange()
@named z_flange_b = Rot.Flange()
@named ϕ_sensor = Rot.AngleSensor()
@named θ_sensor = Rot.AngleSensor()
@named ψ_sensor = Rot.AngleSensor()
ps = @parameters 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.

`@named 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:

```
@named simple_model = ODESystem(simple_eqs, t; systems = [
torque_input, torque, spacecraft
])
simple_sys = structural_simplify(simple_model)
simple_prob = ODEProblem(simple_sys, [], (0, 5*60.0), [])
simple_sol = solve(simple_prob, Tsit5())
```

It's as simple as that. We can now visualize the $\phi$ 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:

```
@named ctrl_spacecraft = LinearSpacecraftAttitude(
Jx=150.0, Jy=100.0, Jz=100.0, u0=[0.5, 0.25, -0.5]
)
@named setpoint = B.Constant(k=0.0)
@named feedback_ϕ = B.Feedback()
@named feedback_θ = B.Feedback()
@named feedback_ψ = B.Feedback()
@named ctrl_ϕ = B.PID(k=10.0, Td=32.0, Ti=100)
@named torque_ϕ = Rot.Torque()
@named ctrl_θ = B.PID(k=10.0, Td=32.0, Ti=100)
@named torque_θ = Rot.Torque()
@named ctrl_ψ = B.PID(k=10.0, Td=32.0, Ti=100)
@named torque_ψ = Rot.Torque()
```

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:

```
@named 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, [], (0, 2.5), [])
ctrl_sol = solve(ctrl_prob, Tsit5())
```

```
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.

CC BY-SA 4.0 Michal Jagodzinski. Last modified: May 09, 2024.

Website built with Franklin.jl and the Julia programming language.

Website built with Franklin.jl and the Julia programming language.