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