Hey @baggepinnen since I am a bit impatient I start experimenting with some custom multibody components I created and the compiler options you suggested and I had some intresting observations.
I created a custom Frame() connector with torque and force vectors as flow variables and euler angles and inertially expressed position vector as non flow variables ( I understand the limitations of using euler angles as part of the connector and for a proper solution that scales creating an orientation object that encodes both rotation matrix and angular velocity is propably beneficial but for the purposes of my toy example I think it should work) . I also created a simple Body() compoent that just writes the kinematics + Euler/Newton equations for dynamic equilibrium.
Depending on how I formulate the angular velocity kinematics in the model I get different response from the symbolic engine and compiler. I will try to list them into four scenarios, two of which contain only one body which should reduce to 12 unknowns, and another two which contain two bodies rigidly connected but should again reduce to 12 unknowns with summed inertia and mass.
1 . For the one_body system in which I explicitly wrote the kinematics that transform angular velocity to euler angle rates (use_w = true) the compile engine correctly takes the “simplest” approach and “roots” angular velocity picking them as the state which leads to symbolically lean equations.This creates 12 differential and 0 algebraic equations. Happy Days!
2 . For the one_body_alt system in which I explicitly wrote the inverse kinematics that transform euler angle rates to angular velocity (use_w = false) the compile engine picks the euler angle rates as states and then in order to calculate angular velocity components x and y creates two extra algebraic equations. I hypothesized this happens because of the singularity created at pitch angle equal to ± pi/2 but I could be wrong cause in the next example it does things that could lead to divisions by zero under certain circumastances. This creates 12 differential and 2 algebraic equations.
3 . For the two_bodies system in which use_w = true the compile engine strangely now picks the euler angle rates as states and then in order to calculate angular velocity components x and y creates two extra algebraic equations for each body similar to case 2. This creates 12 differential and 4 algebraic equations.
4 . For the two_bodies_alt system in which use_w = false the compile engine again picks the euler angle rates as states and manages to symbolically resolve the angular velocities without introducing any algebraic equations. However the expressions created are enormous, which for sure not lead to efficient code. This creates 12 differential and 0 algebraic equations.
It seems that depending one how you write the equations the same physics can lead to different sets of equations which is not always the simplest. It would be intresting if there was a way to force the compiler to pick certain variables as root/states just like you can do in Modelica. I tried using the state_priority metadata but it seems to be completely ignored by mtkcompile.
Below is the code showcasing all four cases.
using ModelingToolkit
using LinearAlgebra
@parameters t
D = Differential(t)
@connector function Frame(; name , r_0 = nothing)
@variables begin
r_0(t)[1:3] = r_0
θ(t)[1:3]
f(t)[1:3],[connect = Flow]
tau(t)[1:3],[connect = Flow]
end
(r_0, f, tau, θ) = collect.(( r_0, f, tau, θ))
return System(Equation[], t, [r_0; f; tau; θ],[]; name = name)
end
@component function Body(; name, use_w = false)
@named frame_a = Frame(r_0 = [0,0,0])
@parameters begin
m = 1, [description = "mass of the body (kg)"]
Ixx = 0.1, [description = "moment of inertia around the body x-axis (kg⋅m²)" ]
Iyy = 0.1, [description = "moment of inertia around the body y-axis (kg⋅m²)" ]
Izz = 0.1, [description = "moment of inertia around the body z-axis (kg⋅m²)" ]
g = 9.81, [description = "gravity vector, resolved in inertial frame"]
end
g_0 = [0, 0, g]
I = [Ixx 0 0; 0 Iyy 0; 0 0 Izz]
@variables begin
v_0(t)[1:3]
θ_t(t)[1:3]
w_a(t)[1:3]
end
(θ, r_0, f_a, tau_a, v_0, θ_t, w_a) = collect.((frame_a.θ, frame_a.r_0, frame_a.f, frame_a.tau, v_0, θ_t, w_a))
R_a = [
cos(θ[2])*cos(θ[3]) -cos(θ[1])*sin(θ[3]) + sin(θ[1])*sin(θ[2])*cos(θ[3]) sin(θ[1])*sin(θ[3]) + cos(θ[1])*sin(θ[2])*cos(θ[3]);
cos(θ[2])*sin(θ[3]) cos(θ[1])*cos(θ[3]) + sin(θ[1])*sin(θ[2])*sin(θ[3]) -sin(θ[1])*cos(θ[3]) + cos(θ[1])*sin(θ[2])*sin(θ[3]);
-sin(θ[2]) sin(θ[1])*cos(θ[2]) cos(θ[1])*cos(θ[2])
] # assuming ZYX Euler angles from local to global
θ_t_exp = [
w_a[1] + sin(θ[1]) * tan(θ[2]) * w_a[2] + cos(θ[1]) * tan(θ[2]) * w_a[3],
cos(θ[1]) * w_a[2] - sin(θ[1]) * w_a[3],
sin(θ[1]) / cos(θ[2]) * w_a[2] + cos(θ[1]) / cos(θ[2]) * w_a[3]
] # assuming ZYX Euler angles
w_a_exp = [
θ_t[1] - θ_t[3]*sin(θ[2]),
θ_t[2]*cos(θ[1]) + θ_t[3]*cos(θ[2])*sin(θ[1]),
-θ_t[2]*sin(θ[1]) + θ_t[3]*cos(θ[2])*cos(θ[1])
] # assuming ZYX Euler angles
ang_eqs = if use_w
θ_t .~ θ_t_exp
else
w_a .~ w_a_exp
end
eqs = [
ang_eqs
D.(θ) .~ θ_t
D.(r_0) .~ v_0
m*(D.(v_0) + g_0) .~ R_a*f_a # expressed in world frame
I*D.(w_a) + cross(w_a, I*w_a) .~ tau_a # expressed in frame_a (body frame)
]
return System(eqs, t; name = name, systems=[frame_a])
end
@component function TwoBodies(; name, use_w = false)
@named B1 = Body(use_w = use_w)
@named B2 = Body(use_w = use_w)
eqs = Equation[
connect(B1.frame_a, B2.frame_a)
]
return System(eqs, t; name = name, systems=[B1, B2])
end
@named one_body = Body(use_w = true) # uses angular velocity w_a as state so equations are symbolically simpler and correctly reduced to 12 unknowns
@named one_body_alt = Body(use_w = false) # uses Euler angle rates θ_t as state which leads to the creation of two extra algebraic equations 14 in total to calculate w_a from θ_t
@named two_bodies = TwoBodies(use_w = true) # produces similar results to one_body_alt picking θ_t as state variables for both bodies and creates algebraic equations to calculate w_a
@named two_bodies_alt = TwoBodies(use_w = false) # reduces to 12 unknowns per body by picking θ_t as state variables directly and since expressions are explicit no extra algebraic equations are created BUT the equations are more symbolically complex
one_body = mtkcompile(one_body; reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 3))
two_bodies = mtkcompile(two_bodies; reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 3))
one_body_alt = mtkcompile(one_body_alt; reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 3))
two_bodies_alt = mtkcompile(two_bodies_alt; reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 3))