# Context

As discussed in Broadcast inside of `@mtkmodel` macro causes precompilation failures when put into a package - #11 by BambOoxX, I’m trying to define `mtkmodel`s to model a moving rigid body submitted to some actions.
Although I could do everything by hand the purpose is to test the `@mtkmodel` syntax and see where the trouble comes.

In Mechanical Components · ModelingToolkitStandardLibrary.jl there are nice simple components that show how to do most things, but I’ve reached a point where I understand the problems, but not the solutions.
Given that I’m using the latest syntax, it seems previous answers on related topics offer limited help.

# Actual problem

What I am trying to do first is to setup a simulation of a ball submitted to its weight, ballistics 101, in a 3D space.
For that I’ve re-implemented 3D versions of serveral components from `ModelingToolkitStandardLibrary`. In need in particular, 3D versions of `MechanicalPort`, `Force`, `Mass` and `Free`. These are implemented below

``````# Example to demonstrate the trajectory of a moving mass in 3D using ModelingToolkit and modified ModelingToolkitStandardLibrary components
using ModelingToolkit, LinearAlgebra, DifferentialEquations
using ModelingToolkitStandardLibrary.Blocks: RealInput
using GeometryBasics, GLMakie

@variables t
D = Differential(t)

@connector MechanicalPort3 begin
v(t)[1:3] = zeros(3), [description = "Translational velocity of the node in [m/s]"]
f(t)[1:3] = zeros(3), [description = "Force entering the node in [N]", connect = Flow]
end

@mtkmodel Force3 begin
@parameters begin
f_0[1:3] = zeros(3)
end
@components begin
node = MechanicalPort3()
f = RealInput(; nin=3, u_start=f_0)
end
@equations begin
collect(node.f .~ -f.u)...
end
end

@mtkmodel Mass3 begin
# Definition of the model icon
#@icon "..."

# Definition of the model parameters
@parameters begin
m, [description = "Mass of the moving rigid body in [kg]"]
v_0[1:3] = zeros(3), [description = "Initial velocities in [m/s] along X, Y and Z axes of the rigid body"]
s_0[1:3] = zeros(3), [description = "Initial positions in [m] along X, Y and Z axes of the rigid body"]
f_0[1:3] = zeros(3), [description = "Initial forces in [N] along X, Y and Z axes entering the rigid body"]
end

@components begin
node = MechanicalPort3(v=v)
end

@variables begin
(s(t))[1:3] = s_0, [description = "Positions in [m] along X, Y and Z axes of the rigid body"]
(v(t))[1:3] = v_0, [description = "Velocities in [m/s] along X, Y and Z axes of the rigid body"]
(f(t))[1:3] = f_0, [description = "Forces in [N] along X, Y and Z axes entering the rigid body"]
end

@equations begin
collect(node.v .~ v)...
collect(node.f .~ f)...
collect(D.(s) .~ v)...
collect(D.(v) .~ f ./ m)...
end
end

@mtkmodel Free3 begin
@components begin
node = MechanicalPort3()
end
@variables begin
f(t)[1:3] = zeros(3)
end
@equations begin
collect(node.f .~ f)...
end
end
``````

Now I can actually instantiate my components with

``````@named mass = Mass3(m=1, v_0=[0.0, 0.0, 0.0])
@named g = Force3(f_0=[0.0, 0.0, -9.81])
@named free = Free3()
``````

My problem comes when I try to make the connections to assemble my model.
I figured that the `g` and `mass` nodes should be connected together to obtain the proper results with

``````eqs = Equation[
ModelingToolkit.connect(g.node, mass.node)
]

@named model = ODESystem(eqs, t, [mass.s; mass.v; mass.f],[];systems = [mass,g])
sys = structural_simplify(model)
``````

However that leaves me with an `ExtraVariablesException`, so there is something I’m missing there, though I can’t figure what.
Nevertheless, if I add manually a relation between `mass.f` and `g.f_0` with

``````eqs = Equation[
collect(mass.f .~ g.f_0)...
ModelingToolkit.connect(g.node, mass.node)
]
``````

I get a seemingly correct behavior, while I would have though initially that `connect` would take care of this relation.
Could someone explain if I missed something in the components or in the way connections works ?

Ok so I think I managed to find the missing part. It does not seem to be enough to state a value for the external forcing directly, so I created an additional `Constant3` block to handle such forcings. This, with some corrections based on ModelingToolkitStandardLibrary.jl/src/Mechanical/TranslationalPosition/TranslationalPosition.jl at v2.3.4 · SciML/ModelingToolkitStandardLibrary.jl · GitHub, gives

``````# Example to demonstrate the trajectory of a moving mass in 3D using ModelingToolkit and modified ModelingToolkitStandardLibrary components
using ModelingToolkit, LinearAlgebra, DifferentialEquations
using ModelingToolkitStandardLibrary.Blocks: RealInput
using GeometryBasics, GLMakie

@variables t
D = Differential(t)

@connector Flange3 begin
s(t)[1:3] = zeros(3), [description = "Translational position of the node in [m/s]"]
f(t)[1:3] = zeros(3), [description = "Force entering the node in [N]", connect = Flow]
end

@mtkmodel Force3 begin
@parameters begin
f_0[1:3] = zeros(3)
end
@components begin
node = Flange3()
f = RealInput(; nin=3, u_start=f_0)
end
@equations begin
collect(node.f .~ -f.u)...
end
end

@mtkmodel Fixed3 begin
@parameters begin
s = zeros(3)
end
@components begin
node = Flange3(; s=s)
end
@equations begin
collect(node.s .~ s)...
end
end

@mtkmodel Free3 begin
@components begin
node = MechanicalPort3()
end
@variables begin
f(t)[1:3] = zeros(3)
end
@equations begin
collect(node.f .~ f)...
end
end

@mtkmodel Constant3 begin
@components begin
output = RealOutput(; nout=3)
end
@parameters begin
k[1:3] = zeros(3), [description = "Constant output value of block"]
end
@equations begin
collect(output.u .~ k)...
end
end

@mtkmodel Mass3 begin
# Definition of the model icon
#@icon "..."

# Definition of the model parameters
@parameters begin
m, [description = "Mass of the moving rigid body in [kg]"]
end

@variables begin
(s(t))[1:3] = zeros(3), [description = "Positions in [m] along X, Y and Z axes of the rigid body"]
(v(t))[1:3] = zeros(3), [description = "Velocities in [m/s] along X, Y and Z axes of the rigid body"]
(f(t))[1:3] = zeros(3), [description = "Forces in [N] along X, Y and Z axes entering the rigid body"]
end

@components begin
node = Flange3(; s=s)
end

@equations begin
collect(node.s .~ s)...
collect(node.f .~ f)...
collect(D.(s) .~ v)...
collect(D.(v) .~ f ./ m)...
end
end

@mtkmodel SpringDamper3 begin
@parameters begin
k, [description = "Spring stiffness in [N/m]"]
c, [description = "Damper viscous damping in [N/m]"]
l0, [description = "Spring length at rest in [m]"]
end

@components begin
origin_node = Flange3()
target_node = Flange3()
end

@variables begin
l(t), [description = "Spring length in [m]"]
v(t), [description = "Damper velocity in [m]"]
f(t)[1:3]
dir(t)[1:3], [description = "Spring direction"]
end

@equations begin
D(l) ~ v
l ~ norm(collect(target_node.s .- origin_node.s))
collect(dir .~ collect(target_node.s .- origin_node.s) ./ l)...
collect(f .~ k .* collect(dir) .* (l - l0) + c .* collect(dir) .* v)...
collect(origin_node.f .~ .+f)...
collect(target_node.f .~ .-f)...
end
end

@named mass = Mass3(m=1, s=[0.0, 0.0, 0.0])
@named springdamper = SpringDamper3(k=10, l0=0.5, c=2)
@named fix = Fixed3(s=[1.0; -1.0; 2.0])
@named g = Force3()
@named sig = Constant3(k=[0.0; 0.0; -9.81])

eqs = Equation[
ModelingToolkit.connect(sig.output, g.f)
ModelingToolkit.connect(g.node, mass.node)
ModelingToolkit.connect(springdamper.target_node, mass.node)
ModelingToolkit.connect(springdamper.origin_node, fix.node)
]

@named model = ODESystem(eqs, t, [], []; systems=[mass, fix, springdamper, g, sig])
sys = structural_simplify(model)
``````

At that point, the `structural_simplify` works and gives as expected 6 degrees of freedom, 3 positions and 3 velocities.

However, if I try to solve that, it fails on some sort of broadcasting.

``````initial_conditions = reduce(vcat, [collect(mass.s .=> [0.0, 1.0, -0.1]), collect(mass.v .=> [0.5, 1.0, 0.0])])
prob = ODEProblem(sys, initial_conditions, (0.0, 100.0))
sol = solve(prob, Rodas5(), abstol=1e-8, reltol=1e-8, saveat=0:0.1:200)
``````

ERROR: MethodError: no method matching +(::Vector{Float64}, ::Float64)
Closest candidates are:
+(::Any, ::Any, ::Any, ::Any…)
@ Base operators.jl:578
+(::T, ::T) where T<:Union{Float16, Float32, Float64}
@ Base float.jl:408
+(::AbstractAlgebra.MatrixElem{T}, ::Union{AbstractFloat, Integer, Rational}) where T<:AbstractAlgebra.NCRingElement

@ChrisRackauckas I’m using ModelingToolkit v8.73.2, do you think it is the same bug as in Broadcast inside of `@mtkmodel` macro causes precompilation failures when put into a package - #11 by BambOoxX or something else ?

@ChrisRackauckas Given that I was not sure if this was due to the lack of array support for such models, I also rewrote my case by composing three 1D block per 3D block but then I get a `ExtraEquationSystemException` so I must over constraining my 3D blocks. What strategy would you consider the best suited for such multidimensional analysis ?

For now, using them with `@mtkmacro` is just suspect. It’s early and doesn’t support all of the array stuff at this point. That coverage is increasing daily so any rule of thumb would be wrong too quickly other than, symbolic arrays need more work for now.

1 Like

Thanks @ChrisRackauckas, so for the moment I tried producing 3D blocks based on the 1D components in `ModelingToolkitStandardLibrary`. However my system blocks on `structural_simplif` seemingly due to the `Spring3` definition. After checking multiple times and testing different implementations, I can’t figure why `structural_simplify` finds extra equations. Could you explain ?
The `structural_simplify` suggests the following extra equations

0 ~ springdamper₊origin₊x₊f(t)
0 ~ springdamper₊origin₊y₊f(t)
0 ~ springdamper₊origin₊z₊f(t)
0 ~ springdamper₊target₊x₊f(t)
0 ~ springdamper₊target₊y₊f(t)
0 ~ springdamper₊target₊z₊f(t)

but none of them should hold. Is there a way to expand the `connect` calls to see exactly what relations are actually built ?

The MWE is the following

``````# Example to demonstrate the trajectory of a moving mass in 3D using ModelingToolkit and modified ModelingToolkitStandardLibrary components
using ModelingToolkit, LinearAlgebra, DifferentialEquations
using ModelingToolkitStandardLibrary
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition
import ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition: Fixed
using GeometryBasics, GLMakie

@variables t
D = Differential(t)

@connector Flange3 begin
@components begin
x = Flange()
y = Flange()
z = Flange()
end
end

@mtkmodel Force3 begin
@parameters begin
f_x_0 = 0.0
f_y_0 = 0.0
f_z_0 = 0.0
use_support = false
end
@components begin
x = Force(s=f_x_0)
y = Force(s=f_y_0)
z = Force(s=f_z_0)
end
end

@mtkmodel Fixed3 begin
@parameters begin
s_0_x = 0.0
s_0_y = 0.0
s_0_z = 0.0
end
@components begin
x = Fixed(; s_0=s_0_x)
y = Fixed(; s_0=s_0_y)
z = Fixed(; s_0=s_0_z)
end
end

@mtkmodel Free begin
@components begin
flange = Flange()
end
@variables begin
f(t) = 0.0
end
@equations begin
flange.f ~ f
end
end

@mtkmodel Free3 begin
@components begin
x = Free()
y = Free()
z = Free()
end
end

@mtkmodel Constant3 begin
@parameters begin
k_x = 0.0
k_y = 0.0
k_z = 0.0
end
@components begin
x = Constant(k=k_x)
y = Constant(k=k_y)
z = Constant(k=k_z)
end
end

@mtkmodel Mass3 begin
@parameters begin
m
s_0_x = 0.0
s_0_y = 0.0
s_0_z = 0.0
v_0_x = 0.0
v_0_y = 0.0
v_0_z = 0.0
f_0_x = 0.0
f_0_y = 0.0
f_0_z = 0.0
end
@components begin
x = Mass(; m=m, s=s_0_x, v=v_0_x, f=f_0_x)
y = Mass(; m=m, s=s_0_y, v=v_0_y, f=f_0_y)
z = Mass(; m=m, s=s_0_z, v=v_0_z, f=f_0_z)
end
end

@mtkmodel SpringDamper3 begin
@parameters begin
k, [description = "Spring stiffness in [N/m]"]
c, [description = "Damper viscous damping in [N/m]"]
l0, [description = "Spring length at rest in [m]"]
end

@components begin
origin = Flange3()
target = Flange3()
end

@variables begin
l(t)
f(t) = k * (l - l0)
end

@equations begin
l ~ sqrt(abs2(target.x.s - origin.x.s) + abs2(target.y.s - origin.y.s) + abs2(target.z.s - origin.z.s))
f ~ k * (l - l0)
origin.x.f ~ f * (target.x.s - origin.x.s) / l
origin.y.f ~ f * (target.y.s - origin.y.s) / l
origin.z.f ~ f * (target.z.s - origin.z.s) / l
target.x.f ~ -f * (target.x.s - origin.x.s) / l
target.y.f ~ -f * (target.y.s - origin.y.s) / l
target.z.f ~ -f * (target.z.s - origin.z.s) / l
end
end

@named mass = Mass3(m=1, s_0_x=0.0, s_0_y=0.0, s_0_z=0.0)
@named springdamper = SpringDamper3(k=10, l0=0.5, c=2)
@named fix = Fixed3(s_0_x=1.0, s_0_y=-1.0, s_0_z=2)
@named free = Free3()
@named g = Force3()
@named sig = Constant3(k_z=-9.81)

eqs = Equation[
ModelingToolkit.connect(sig.x.output, g.x.f)
ModelingToolkit.connect(sig.y.output, g.y.f)
ModelingToolkit.connect(sig.z.output, g.z.f)
ModelingToolkit.connect(springdamper.target.x, mass.x.flange, g.x.flange)
ModelingToolkit.connect(springdamper.target.y, mass.y.flange, g.y.flange)
ModelingToolkit.connect(springdamper.target.z, mass.z.flange, g.z.flange)
ModelingToolkit.connect(fix.x.flange, springdamper.origin.x)
ModelingToolkit.connect(fix.y.flange, springdamper.origin.y)
ModelingToolkit.connect(fix.z.flange, springdamper.origin.z)
]

@named model = ODESystem(eqs, t; systems=[mass, fix, springdamper, g, sig])
sys = structural_simplify(model)
``````

`expand_connections` IIRC

1 Like

Yes, that seems to be it, thanks. This is exactly what I needed. So now I use this to compare expanded and non-expanded connections. I’ve isolated the `springdamper` to make things clearer.

I have two questions here

• I suppose the additional equations added by `expand_connections` are supposed to handle the equilibrium of a `Flange` i.e. sum of all forces = 0. Can you confirm ?

• Why are there double equations ? I would have expected 3 equilibrium condition per flange, but here I get 6 of them. Doing the same comparison with s simple `Spring` component yields as expected only 2 equations with expansion…

`equations(expand_connections(springdamper))` gives

``````20-element Vector{Equation}:
springdamper₊l(t) ~ sqrt(abs2(springdamper₊target₊y₊s(t) - springdamper₊origin₊y₊s(t)) + abs2(-springdamper₊origin₊z₊s(t) + springdamper₊target₊z₊s(t)) + abs2(-springdamper₊origin₊x₊s(t) + springdamper₊target₊x₊s(t)))
springdamper₊f(t) ~ springdamper₊k*(-springdamper₊l0 + springdamper₊l(t))
springdamper₊origin₊x₊f(t) ~ ((-springdamper₊origin₊x₊s(t) + springdamper₊target₊x₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊origin₊y₊f(t) ~ ((springdamper₊target₊y₊s(t) - springdamper₊origin₊y₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊origin₊z₊f(t) ~ ((-springdamper₊origin₊z₊s(t) + springdamper₊target₊z₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊target₊x₊f(t) ~ ((springdamper₊origin₊x₊s(t) - springdamper₊target₊x₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊target₊y₊f(t) ~ ((-springdamper₊target₊y₊s(t) + springdamper₊origin₊y₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊target₊z₊f(t) ~ ((springdamper₊origin₊z₊s(t) - springdamper₊target₊z₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
0 ~ springdamper₊origin₊x₊f(t)
0 ~ springdamper₊origin₊y₊f(t)
0 ~ springdamper₊origin₊z₊f(t)
0 ~ springdamper₊target₊x₊f(t)
0 ~ springdamper₊target₊y₊f(t)
0 ~ springdamper₊target₊z₊f(t)
0 ~ springdamper₊origin₊x₊f(t)
0 ~ springdamper₊origin₊y₊f(t)
0 ~ springdamper₊origin₊z₊f(t)
0 ~ springdamper₊target₊x₊f(t)
0 ~ springdamper₊target₊y₊f(t)
0 ~ springdamper₊target₊z₊f(t)
``````

while `equations(springdamper)` gives

``````8-element Vector{Equation}:
springdamper₊l(t) ~ sqrt(abs2(springdamper₊target₊y₊s(t) - springdamper₊origin₊y₊s(t)) + abs2(-springdamper₊origin₊z₊s(t) + springdamper₊target₊z₊s(t)) + abs2(-springdamper₊origin₊x₊s(t) + springdamper₊target₊x₊s(t)))
springdamper₊f(t) ~ springdamper₊k*(-springdamper₊l0 + springdamper₊l(t))
springdamper₊origin₊x₊f(t) ~ ((-springdamper₊origin₊x₊s(t) + springdamper₊target₊x₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊origin₊y₊f(t) ~ ((springdamper₊target₊y₊s(t) - springdamper₊origin₊y₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊origin₊z₊f(t) ~ ((-springdamper₊origin₊z₊s(t) + springdamper₊target₊z₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊target₊x₊f(t) ~ ((springdamper₊origin₊x₊s(t) - springdamper₊target₊x₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊target₊y₊f(t) ~ ((-springdamper₊target₊y₊s(t) + springdamper₊origin₊y₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
springdamper₊target₊z₊f(t) ~ ((springdamper₊origin₊z₊s(t) - springdamper₊target₊z₊s(t))*springdamper₊f(t)) / springdamper₊l(t)
``````

Hi, I’ve found what I’ve been doing incorrectly, except from a sign error somewhere, I defined `Flange3` using `@connector` which doubles the equilibrium equations, changing that to `@mtkmodel` works just fine.
Thanks for the help @ChrisRackauckas, do you think this type of redundant equations could be avoided with `structural_simplify` ?

They should get eliminated. If not, you can open an issue that highlights what you think should be eliminated but isn’t and we can investigate why.

Ok, thanks. I’ll set up an issue.

1 Like