How can I define a vector in ModelingToolkit

Specific Domains Modelling & Simulations
, ,
1

The following code fails:

# Example one:
# Falling mass.

using ModelingToolkit, OrdinaryDiffEq

G_EARTH  = Float64[0.0, 0.0, -9.81]    # gravitational acceleration

# model
@variables t pos(t)[1:3]=[0.0, 0.0,  0.0]
@variables   vel(t)[1:3]=[0.0, 0.0, 50.0] 
@variables   acc(t)[1:3]=[0.0, 0.0, -9.81] 
D = Differential(t)

eqs = [D(pos) ~ vel,
       D(vel) ~ acc,
       acc ~ G_EARTH]

@named sys = ODESystem(eqs, t)

with the message:

ERROR: LoadError: axes of Differential(t)((pos(t))[1:3]) not known
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] axes(A::Symbolics.Arr{Num, 1})
   @ Symbolics ~/.julia/packages/Symbolics/5BibL/src/arrays.jl:537 [inlined]
 [3] combine_axes
   @ ./broadcast.jl:524 [inlined]
 [4] instantiate
   @ ./broadcast.jl:306 [inlined]
 [5] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(~), Tuple{Num, Num}})
   @ Base.Broadcast ./broadcast.jl:903 [inlined]
 [6] ~(lhs::SymbolicUtils.BasicSymbolic{Symbolics.Arr{Num, 1}}, rhs::Symbolics.Arr{Num, 1})
   @ Symbolics ~/.julia/packages/Symbolics/5BibL/src/equations.jl:104
 [7] top-level scope
   @ ~/repos/Tether/Tether_01.jl:17
 [8] include(fname::String)
   @ Base.MainInclude ./client.jl:489
 [9] top-level scope
   @ REPL[1]:1
in expression starting at /home/ufechner/repos/Tether/Tether_01.jl:17

julia> include("Tether_01.jl")
ERROR: LoadError: axes of Differential(t)((pos(t))[1:3]) not known
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] axes(A::Symbolics.Arr{Num, 1})
   @ Symbolics ~/.julia/packages/Symbolics/5BibL/src/arrays.jl:537 [inlined]
 [3] combine_axes
   @ ./broadcast.jl:524 [inlined]
 [4] instantiate
   @ ./broadcast.jl:306 [inlined]
 [5] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(~), Tuple{Num, Num}})
   @ Base.Broadcast ./broadcast.jl:903 [inlined]
 [6] ~(lhs::SymbolicUtils.BasicSymbolic{Symbolics.Arr{Num, 1}}, rhs::Symbolics.Arr{Num, 1})
   @ Symbolics ~/.julia/packages/Symbolics/5BibL/src/equations.jl:104
 [7] top-level scope
   @ ~/repos/Tether/Tether_01.jl:14
 [8] include(fname::String)
   @ Base.MainInclude ./client.jl:489
 [9] top-level scope
   @ REPL[1]:1
in expression starting at /home/ufechner/repos/Tether/Tether_01.jl:14

What am I doing wrong?

2

This does not give a syntax error:

# Example one:
# Falling mass.

using ModelingToolkit, OrdinaryDiffEq

G_EARTH  = Float64[0.0, 0.0, -9.81]    # gravitational acceleration

# model
@variables t pos(t)[1:3]=[0.0, 0.0,  0.0]
@variables   vel(t)[1:3]=[0.0, 0.0, 50.0] 
@variables   acc(t)[1:3]=[0.0, 0.0, -9.81] 
D = Differential(t)

eqs = vcat(D.(pos) ~ vel,
           D.(vel) ~ acc,
           acc    .~ G_EARTH)

@named sys = ODESystem(eqs, t)

Lets see if it works…

3

Works fine:

# Example one: Falling mass.
using ModelingToolkit, OrdinaryDiffEq, PyPlot

G_EARTH  = Float64[0.0, 0.0, -9.81]    # gravitational acceleration

# model
@variables t pos(t)[1:3]=[0.0, 0.0,  0.0]
@variables   vel(t)[1:3]=[0.0, 0.0, 50.0] 
@variables   acc(t)[1:3]=[0.0, 0.0, -9.81] 
D = Differential(t)

eqs = vcat(D.(pos) ~ vel,
           D.(vel) ~ acc,
           acc    .~ G_EARTH)

@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)

duration = 10.0
dt = 0.02
tspan = (0.0, duration)
ts    = 0:dt:duration
u0 = zeros(6)
u0[6] = 50.0

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

X = sol.t
POS_Z = sol(X, idxs=pos[3])
VEL_Z = sol(X, idxs=vel[3])

plot(X, POS_Z, color="green")
xlabel("time [s]")
grid(true)
twinx()
ylabel("vel_z [m/s]") 
plot(X, VEL_Z, color="red") 
nothing

Only remaining question:
Why does POS = sol(X, idxs=pos[1:3]) not work?

4

Symbolic arrays have a number of problems, you might need to do idxs = collect(pos) for now

1 Like
5

you can do


julia> sol[pos]
501-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0]
 [0.0, 0.0, 0.9980380000000015]
 [0.0, 0.0, 1.9921519999999922]
 [0.0, 0.0, 2.982341999999982]
 [0.0, 0.0, 3.96860799999997]
 [0.0, 0.0, 4.950949999999957]
 ⋮
 [0.0, 0.0, 14.260949999999628]
 [0.0, 0.0, 13.316607999999672]
 [0.0, 0.0, 12.368341999999714]
 [0.0, 0.0, 11.416151999999672]
 [0.0, 0.0, 10.460037999999674]
 [0.0, 0.0, 9.499999999999718]
1 Like
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