This works (Julia 1.11.3, MTK 9.61.0 and 9.62.0)
This little test system comprises a mass-flange-inputforce system with gravity, a forced sinusoidal input function via the flange, and a constant force direct on the mass.
(edit: added the packages and constants used)
using ModelingToolkit
using DifferentialEquations
using LaTeXStrings
using Plots
import LinearAlgebra as LA
import ModelingToolkitStandardLibrary.Mechanical.Translational as T
import ModelingToolkitStandardLibrary.Blocks as B
using ModelingToolkit: t_nounits as t, D_nounits as D
# a few constants for easy test switching
const G_EARTH = [0.0, 0.0, -9.81 ] # m/s²
const G_MOON = [0.0, 0.0, -1.625 ]
const G_MARS = [0.0, 0.0, -3.72076]
const Earth="Earth"
const Moon="Moon"
const Mars="Mars"
const gravities=Dict(Earth => G_EARTH, Moon=>G_MOON, Mars=>G_MARS)
@connector function Flange3D(;name)
vars = @variables begin
f(t)[1:3], [connect=Flow, description = "Cut force resolved in connector"]
v(t)[1:3], [ description = "Cut velocities in connector"]
end
pars = []
v,f = collect.((v,f))
vars = [v;f]
eqs = [1~1] # This avoids "BoundsError: attempt to access 0-element Vector{Vector{Any}} at index [0]"
eq = eqs==[] ? Equation[] : eqs
ODESystem(eq, t, vars, pars; name)
end
@component function Mass3D(; name, m, g=G_EARTH)
pars = @parameters begin # this allows the params to be set in mtkbuild call OR as system parameters
m=m
g[1:3]=g
f0[1:3]=[0.,0.,0.]
end
vars = @variables begin
s(t)[1:3] # =[0.,0.,0.], [description = "Positions in [m] along X, Y and Z axes of the mass"]
v(t)[1:3] # =[0.,0.,0.], [description = "Velocities in [m/s] along X, Y and Z axes of the mass"]
a(t)[1:3], [description = "Acceleration in [m/s²] along the X,Y,Z axes of the mass"]
f(t)[1:3], [description = "Forces in [N] along X,Y,Z axes of the mass"]
end
systems = @named begin
flange = Flange3D() # connects v & f
end
eqs = [
collect(D.(s) .~ v)...
collect(D.(v) .~ a)...
collect(v .~ flange.v)...
collect(f .~ flange.f)...
collect(a .~ (f .+f0)./m .+ g)...
]
initialization_eqs=[]
return ODESystem(eqs, t, vars, pars; systems, name, initialization_eqs)
end
@component function ForceFunction3D(; name, a0=[0.,0.,0.])
pars = @parameters a0[1:3]=a0
vars = @variables f(t)
systems = @named begin
output = B.RealOutputArray(; nout=3)
end
@constants π=pi # this makes for 'pretty' equations in display
f = a0 .* sin(π*t) # this gets passed through!
eqs = [
collect(output.u .~ -f)...
]
return ODESystem(eqs, t, vars, pars; systems, name)
end
@component function Force3D(; name)
pars = []
vars = []
systems = @named begin
input = B.RealInputArray(; nin=3, guess=[0.,0.,0.])
flange = Flange3D()
end
eqs = [
collect(flange.f .~ input.u)...
]
return ODESystem(eqs, t, vars, pars; name, systems)
end
@component function ForcedMass3D(;name, m, g=G_EARTH, a0=[0.,0.,0.] )
pars = []
vars = []
systems = @named begin
force = Force3D()
force_func = ForceFunction3D(a0=a0)
mass = Mass3D(m=m,g=g)
end
eqs = [
ModelingToolkit.connect(force_func.output, force.input )
ModelingToolkit.connect(force.flange, mass.flange)
]
return ODESystem(eqs, t, vars, pars; systems, name)
end
@mtkbuild sys1 = ForcedMass3D(m=1.0); # = @named + structural_simplity()
parameters(sys1) # mass.f0, mass.g, mass.m, force_func.a0
Simulate it with the variable parameters available without recreating the ODESystem:
dt = 0.02
duration = 5.0
tspan, ts = time_span(duration, dt) # duration, dt
# MTK unknows are pos s & vel v: unknowns(simple_sys)
s00 = [0.,0.,0.]
v00 = [0.,0.,0.]
s0 = [10., 20., 10.]
v0 = [5., 2., 3.]
a0 = [0.,10.,+20] # amplitude of input force function
f0 = [1.,0.,10.] # constant forces on mass over and above gravity
z0 = [0.,0.,0.]
p = Earth
g = gravities[p]
m = 1.
u0 = []
p0 = Dict( :mass₊m=>m, :mass₊g=>g, :mass₊f0=>f0, :force_func₊a0=>a0) # ... this works :-)
forcedmass_init_eqs = [
collect(sys1.mass.s .~ s0 )...
collect(sys1.mass.v .~ v0 )...
]
prob = ODEProblem(sys1, u0, tspan, p0; initialization_eqs=forcedmass_init_eqs)
sol = solve(prob, Rodas5(), dt=dt, saveat=ts)
plot(sol, title="$p? trajectory")