How can I define a vector in ModelingToolkit

Specific Domains Modelling & Simulations
, ,
6

Above code block doesn’t work anymore under ModelingToolkit v9.61.0 (Julia 1.11.3) but works with a few adjustments.

duration = 5.0
tspan = (0.0, duration)
dt = 0.02
ts = 0:dt:duration
	
G_EARTH  = Float64[0.0, -0.0, -9.81]    # gravitational acceleration
G_MOON   = G_EARTH * 0.166
	
# model
@variables (pos(t))[1:3] = [0.0, 0.0, 0.0]   # defaults to 'at rest at origin'
@variables (vel(t))[1:3] = [0.0, 0.0, 0.0] 
@variables (acc(t))[1:3]                    # NO defaults as balanced by force!

@parameters (force)[1:3] # do not provide default or structural_simplify will remove eqs where f=0 !
	
eqs = [
	collect( D.(pos) .~ vel )...   #   .~ instead of just ~ !
	collect( D.(vel) .~ acc )...
    collect( acc .~ force )...
]
	
@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)
	
u0 = [ pos=>[0.,0.,0.], vel=>[0.,0.,+10.0] ]  # MTK unknows are pos & vel: unknowns(simple_sys) 
p0 = [[0.,0.,-9.81]] # [G_EARTH] 

prob = ODEProblem(simple_sys, u0, tspan, p0)
sol = solve(prob, Rodas5(), dt=dt, saveat=ts)
plot(sol)

When constructing as a component block it becomes as follows (I added a force as well as default gravity). Note the requirement to state the initial condition keys as symbolic :pos and :vel if using the Dict version of u0:

# model as component
@component function TestMass3D(;name, m, g=G_EARTH)
	
	vars = @variables begin
		(pos(t))[1:3] = [0.0, 0.0, 0.0]   # defaults to 'at rest at origin'
		(vel(t))[1:3] = [0.0, 0.0, 0.0] 
		(acc(t))[1:3]                    # NO defaults as balanced by force!
	end
		
	pars = @parameters (force)[1:3] # do not provide default or structural_simplify will remove eqs where f=0 !
	
	eqs = [
		collect( D.(pos) .~ vel )...
		collect( D.(vel) .~ acc )...
	    collect( acc .~ force./m .+ g )...
	]

	return ODESystem(eqs, t, vars, pars;name)
end

@named sys = TestMass3D(m=1, g= G_MOON)
simple_sys = structural_simplify(sys)
	
u0 = [ :pos=>[0.,0.,0.], :vel=>[0.,0.,+10.0] ]  # MTK unknows are pos & vel: check with unknowns(simple_sys) 
p0 = [ [0.,0.,+2.] ]  # note vector inside array as only 1 parameter 

prob = ODEProblem(simple_sys, u0, tspan, p0)

sol = solve(prob, Rodas5(), dt=dt, saveat=ts)
plot(sol)

(this also links with https://discourse.julialang.org/t/broadcast-inside-of-mtkmodel-macro-causes-precompilation-failures-when-put-into-a-package/107277/18