# 2D arrays with ModelingToolkit

I have the following code, representing a mass attached to a spring-damper element:

``````using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra

G_EARTH     = Float64[0.0, 0.0, -9.81]          # gravitational acceleration [m/s²]
L0::Float64 = -10.0                             # initial spring length      [m]
V0::Float64 = 4                                 # initial velocity           [m/s]

# model, Z component upwards
@parameters mass=1.0 c_spring0=50.0 damping=0.5 l0=L0
@variables t pos(t)[1:3] = [0.0, 0.0,  L0]
@variables   vel(t)[1:3] = [0.0, 0.0,  V0]
@variables   acc(t)[1:3] = G_EARTH
@variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
@variables c_spring(t) = c_spring0
@variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
@variables force(t) = 0.0 norm1(t) = abs(l0) spring_vel(t) = 0.0
D = Differential(t)

eqs = vcat(D.(pos)      ~ vel,
D.(vel)      ~ acc,
norm1        ~ norm(pos),
unit_vector  ~ -pos/norm1,         # direction from point mass to origin
spring_vel   ~ -unit_vector ⋅ vel,
c_spring     ~ c_spring0 * (norm1 > abs(l0)),
spring_force ~ (c_spring * (norm1 - abs(l0)) + damping * spring_vel) * unit_vector,
acc          ~ G_EARTH + spring_force/mass)
``````

Now I want to model n+1 masses connected by n spring-damper elements to simulate a tether.

Any ideas how to achieve this? Can I define 2D arrays for position, velocity and acceleration? Can I differentiate a 2D array? Can I loop over an array when creating the symbolic equations?

Background: GitHub - ufechner7/Tethers.jl: Tether models

OK, found a first hint: Symbolic arrays · Symbolics.jl

I’m working on seemingly closely related topics, so pay attention to Broadcast inside of `@mtkmodel` macro causes precompilation failures when put into a package if you intend to use the new syntax at some point

1 Like

Could you not try to create an array of single spring mass components and then incorporate the connections between them? Something like this?

``````using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra

G_EARTH     = Float64[0.0, 0.0, -9.81]          # gravitational acceleration [m/s²]
L0::Float64 = -10.0                             # initial spring length      [m]
V0::Float64 = 4                                 # initial velocity           [m/s]

# model, Z component upwards
function springmass(;name)
@parameters mass=1.0 c_spring0=50.0 damping=0.5 l0=L0
@variables t pos(t)[1:3] = [0.0, 0.0,  L0]
@variables   vel(t)[1:3] = [0.0, 0.0,  V0]
@variables   acc(t)[1:3] = G_EARTH
@variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
@variables c_spring(t) = c_spring0
@variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
@variables force(t) = 0.0 norm1(t) = abs(l0) spring_vel(t) = 0.0
D = Differential(t)

eqs = vcat(D.(pos)      ~ vel,
D.(vel)      ~ acc,
norm1        ~ norm(pos),
unit_vector  ~ -pos/norm1,         # direction from point mass to origin
spring_vel   ~ -unit_vector ⋅ vel,
c_spring     ~ c_spring0 * (norm1 > abs(l0)),
spring_force ~ (c_spring * (norm1 - abs(l0)) + damping * spring_vel) * unit_vector,
acc          ~ G_EARTH + spring_force/mass)
return ODESystem(eqs, t; name=name)
end

n=5
sys_array = []
for i in range(1, n)
push!(sys_array, springmass(name=Symbol("springmass\$i")))
end
for sys in sys_array
# Connect the springs in a chain
end
``````
1 Like

I have now a nicely working example:

``````# Tutorial example simulating a 3D mass-spring system with a nonlinear spring (no spring forces
# for l < l_0) and n tether segments.
using ModelingToolkit, OrdinaryDiffEq, Plots, LinearAlgebra

G_EARTH     = Float64[0.0, 0.0, -9.81]          # gravitational acceleration     [m/s²]
L0::Float64 = 10.0                              # initial segment length            [m]
V0::Float64 = 2                                 # initial velocity of lowest mass [m/s]
segments::Int64 = 4                             # number of tether segments         [-]
α0 = π/8                                        # initial angle                   [rad]
POS0 = zeros(3, segments+1)
VEL0 = zeros(3, segments+1)
ACC0 = zeros(3, segments+1)
SEGMENTS0 = zeros(3, segments)
UNIT_VECTORS0 = zeros(3, segments)
for i in 1:segments+1
l0 = -(i-1)*L0
v0 = (i-1)*V0/segments
POS0[:, i] .= [sin(α0) * l0, 0, cos(α0) * l0]
VEL0[:, i] .= [sin(α0) * v0, 0, cos(α0) * v0]
end
for i in 2:segments+1
ACC0[:, i] .= G_EARTH
end
for i in 1:segments
UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
end

# model, Z component upwards
@parameters mass=1.0 c_spring0=50.0 damping=0.5 l_seg=L0
@variables t
@variables pos(t)[1:3, 1:segments+1]  = POS0
@variables vel(t)[1:3, 1:segments+1]  = VEL0
@variables acc(t)[1:3, 1:segments+1]  = ACC0
@variables segment(t)[1:3, 1:segments]  = SEGMENTS0
@variables unit_vector(t)[1:3, 1:segments]  = UNIT_VECTORS0
@variables norm1(t)[1:segments] = l_seg * ones(segments)
@variables rel_vel(t)[1:3, 1:segments]  = zeros(3, segments)
@variables spring_vel(t)[1:segments] = zeros(segments)
@variables c_spring(t)[1:segments] = c_spring0 * ones(segments)
@variables spring_force(t)[1:3, 1:segments] = zeros(3, segments)
@variables total_force(t)[1:3, 1:segments] = zeros(3, segments)
D = Differential(t)

eqs1 = vcat(D.(pos) ~ vel,
D.(vel) ~ acc)
eqs2 = []
for i in segments:-1:1
global eqs2
eqs2 = vcat(eqs2, segment[:, i] ~ pos[:, i+1] - pos[:, i])
eqs2 = vcat(eqs2, norm1[i] ~ norm(segment[:, i]))
eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])
eqs2 = vcat(eqs2, spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i])
eqs2 = vcat(eqs2, c_spring[i] ~ c_spring0 * (norm1[i] > l_seg))
eqs2 = vcat(eqs2, spring_force[:, i] ~ (c_spring[i] * (norm1[i] - l_seg) + damping * spring_vel[i]) * unit_vector[:, i])
if i == segments
eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i])
else
eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i]- spring_force[:, i+1])
end
eqs2 = vcat(eqs2, acc[:, i+1] .~ G_EARTH + total_force[:, i] / mass)
end
eqs2 = vcat(eqs2, acc[:, 1] .~ zeros(3))
eqs = vcat(eqs1..., eqs2)

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

duration = 10.0
dt = 0.02
tol = 1e-6
tspan = (0.0, duration)
ts    = 0:dt:duration

u0 = Float64[]
for i in 1:segments+1
global u0
u0 = vcat(u0, POS0[:, i], VEL0[:, i])
end

prob = ODEProblem(simple_sys, u0, tspan)
@time sol = solve(prob, Rodas5(), dt=dt, abstol=tol, reltol=tol, saveat=ts)

function plot2d(sol, reltime, segments)
index = Int64(round(reltime*50+1))
x = Float64[]
z = Float64[]
for particle in 1:segments+1
push!(x, (sol(sol.t, idxs=pos[1, particle]))[index])
push!(z, (sol(sol.t, idxs=pos[3, particle]))[index])
end
x_max = maximum(x)
z_max = maximum(z)
plot(x,z, xlabel="x [m]", ylabel="z [m]", legend=false)
annotate!(15, z_max-3.0, "t=\$(round(reltime,digits=1)) s")
plot!(x, z, seriestype = :scatter)
ylims!((-segments*10-10, 0.5))
xlims!((-segments*5, segments*5))
end

dt = 0.04
for time in 0:dt:10
display(plot2d(sol, time, segments))
sleep(0.25*dt)
end
nothing
``````

Try it out, it looks cool!

Only one question left:
I use this code to initialize the state vector of the solver:

``````u0 = Float64[]
for i in 1:segments+1
global u0
u0 = vcat(u0, POS0[:, i], VEL0[:, i])
end
``````

This seams to be a bit fragile. It is not clear to me if the order of the states in the state vector is well defined. Is there a better way to initialize the state vector of the solver?

UPDATE: This is better:

``````# initial state
u0 = Dict(pos=>POS0, vel=>VEL0)
``````

Don’t do this, the order can (and often do) change with any patch release of MTK.

The initial condition is supposed to be supplied as a dict

``````Dict(var => value, ...)
``````

since symbolic arrays are a bit buggy, you probably have to do

``````Dict([
collect(array_var1) .=> values1;
collect(array_var2) .=> values2;
])
``````
1 Like

For now I don’t want to use components, just build one large model as I did it in the past numerically without modelingtoolkit. See: https://github.com/ufechner7/KiteModels.jl/blob/main/src/KPS3.jl

When this kind of 1 to 1 translation works, than I might re-consider to use components… But in the end I probably need a particle system where the springs are not only between two masses, but in a more complex graph for the bridle of the kite…

So treating the kite as a component separate from the tether probably does not make much sense…

What can make sense is to treat the winch, the aerodynamic and the structural model of the kite separate…

1 Like

This actually works:

``````# initial state
u0 = Dict(pos=>POS0, vel=>VEL0)
``````

But I could not find that in the documentation of ODEProblem…

Is `vcat` the best way to add an equation in a loop?

``````eqs = vcat(eqs, total_force[:, i] ~ spring_force[:, i])
``````

I also tried `push!()`, but it did not work.

`push!` only works when adding a single element. To add multiple element, use `append!`.

Thank you!

One more question: I use events to detect when the length of a tether segment becomes less than the unstretched length, because then the stiffness becomes zero.

For one segment I can write:

``````function condition(u, t, integrator) # Event when condition(u,t,integrator) == 0
norm(u[1:3]) - abs(L0)
end
function affect!(integrator)
println(integrator.t)            # Not needed, just to show that the callback works
end
cb = ContinuousCallback(condition, affect!)
``````

But how can I create events if I have an array of segments, and each of them might reach a length below the unstretched length and must generate an event? Can I create an array of conditions? Or an array of callbacks?

Guess I need to use a VectorContinuousCallback in this case…

The questions that remains is: I just have state vector `u` … How do I know which indices relate to the `pos` and which to the `vel` vector?

Have you read the MTK docs on events? Event Handling and Callback Functions · ModelingToolkit.jl

1 Like

No, I was reading only about events in the documentation of DifferentialEquations, see: Event Handling and Callback Functions · DifferentialEquations.jl

Strange that there are no cross-links between the two documentations…

The first two paragraphs of the MTK event docs contains two different links to the DiffEq docs. The DiffEq docs do not point to MTK because MTK is a downstream package, of which there are 100s.

OK, I tried:

``````# Example two: Falling mass, attached to non-linear spring without compression stiffness,
# initially moving upwards with 4 m/s, using a callback to precisely calculate the discontinuities
using ModelingToolkit, OrdinaryDiffEq, PyPlot, LinearAlgebra

G_EARTH  = Float64[0.0, 0.0, -9.81]             # gravitational acceleration [m/s²]
L0::Float64 = -10.0                             # initial spring length      [m]
V0::Float64 = 4                                 # initial velocity           [m/s]

# model, Z component upwards
@parameters mass=1.0 c_spring0=50.0 damping=0.5 l0=L0
@variables t pos(t)[1:3] = [0.0, 0.0,  L0]
@variables   vel(t)[1:3] = [0.0, 0.0,  V0]
@variables   acc(t)[1:3] = G_EARTH
@variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
@variables c_spring(t) = c_spring0
@variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
@variables force(t) = 0.0 norm1(t) = abs(l0) spring_vel(t) = 0.0
D = Differential(t)

eqs = vcat(D.(pos)      ~ vel,
D.(vel)      ~ acc,
norm1        ~ norm(pos),
unit_vector  ~ -pos/norm1,         # direction from point mass to origin
spring_vel   ~ -unit_vector ⋅ vel,
c_spring     ~ c_spring0 * (norm1 > abs(l0)),
spring_force ~ (c_spring * (norm1 - abs(l0)) + damping * spring_vel) * unit_vector,
acc          ~ G_EARTH + spring_force/mass)

@named sys = ODESystem(eqs, t; continuous_events = [norm(pos) ~ abs(L0)])
simple_sys = structural_simplify(sys)

duration = 10.0
dt = 0.02
tol = 1e-6
tspan = (0.0, duration)
ts    = 0:dt:duration
# initial state
u0 = Dict(pos=>[0,0,L0], vel=>[0,0,V0])

# function condition(u, t, integrator) # Event when condition(u,t,integrator) == 0
#     norm(u[1:3]) - abs(L0)
# end
# function affect!(integrator)
#     println(integrator.t)
# end
# cb = ContinuousCallback(condition, affect!)

prob = ODEProblem(simple_sys, u0, tspan)
sol = solve(prob, Rodas5(), dt=dt, abstol=tol, reltol=tol, saveat=ts, callback = cb)

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

lns1 = plot(X, POS_Z, color="green", label="pos_z")
xlabel("time [s]")
ylabel("pos_z [m]")
lns2 = plot(X, L0.+0.005 .* sol[c_spring], color="grey", label="c_spring")
grid(true)
twinx()
ylabel("vel_z [m/s]")
lns3 = plot(X, VEL_Z, color="red", label="vel_z")
lns = vcat(lns1, lns2, lns3)
labs = [l.get_label() for l in lns]
legend(lns, labs)
PyPlot.show(block=false)
nothing
``````

But this does not work, I get the error:

``````julia> include("src/Tether_03b.jl")
ERROR: LoadError: UndefVarError: `pos` not defined
Stacktrace:
[1] macro expansion
@ ~/.julia/packages/SymbolicUtils/ssQsQ/src/code.jl:418 [inlined]
[2] macro expansion
@ ~/.julia/packages/Symbolics/5BibL/src/build_function.jl:537 [inlined]
[3] macro expansion
@ ~/.julia/packages/SymbolicUtils/ssQsQ/src/code.jl:375 [inlined]
[4] macro expansion
@ ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:163 [inlined]
[5] macro expansion
@ ./none:0 [inlined]
[6] generated_callfunc(f::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…}, __args::Vararg{…}) where {…}
@ ModelingToolkit ./none:0 [inlined]
[7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::Vector{…}, ::Float64)
@ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:150
[8] #346
@ ~/.julia/packages/ModelingToolkit/dCa81/src/systems/callbacks.jl:403 [inlined]
[9] determine_event_occurance
@ ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:290 [inlined]
[10] find_callback_time(integrator::OrdinaryDiffEq.ODEIntegrator{…}, callback::ContinuousCallback{…}, counter::Int64)
@ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:397
[11] macro expansion
@ ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:131 [inlined]
[12] find_first_continuous_callback
@ ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:126 [inlined]
[13] find_first_continuous_callback
@ ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:124 [inlined]
[14] handle_callbacks!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/integrators/integrator_utils.jl:324
[15] _loopfooter!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/integrators/integrator_utils.jl:247
[16] loopfooter!
@ ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/integrators/integrator_utils.jl:200 [inlined]
[17] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/solve.jl:539
[18] __solve(::ODEProblem{…}, ::Rodas5{…}; kwargs::@Kwargs{…})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/solve.jl:6
[19] __solve
@ ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/solve.jl:1 [inlined]
[20] solve_call(_prob::ODEProblem{…}, args::Rodas5{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:561
[21] solve_call
@ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:527 [inlined]
[22] #solve_up#42
@ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:1010 [inlined]
[23] solve_up
@ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:996 [inlined]
[24] #solve#40
@ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:933 [inlined]
[25] top-level scope
@ ~/repos/Tethers.jl/src/Tether_03b.jl:49
[26] include(fname::String)
@ Base.MainInclude ./client.jl:489
[27] top-level scope
@ REPL[1]:1
in expression starting at /home/ufechner/repos/Tethers.jl/src/Tether_03b.jl:49
Some type information was truncated. Use `show(err)` to see complete types.
``````

Any idea?

1 Like

I see that the newest version of `Symbolics.jl` is not used, but this is due to a restriction in `ModelingToolkit.jl`

@ufechner7 Did you find an approach to add multiple tether forces to your `acc` equation ? With a for loop or else ?

Yes, the multi-segment tether simulation works fine: Examples · Tethers.jl

Not yet working is the continues callback for this model, but I know how to implement it.

``````for i in segments:-1:1
global eqs2
eqs2 = vcat(eqs2, segment[:, i] ~ pos[:, i+1] - pos[:, i])
eqs2 = vcat(eqs2, norm1[i] ~ norm(segment[:, i]))
eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])
eqs2 = vcat(eqs2, spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i])
eqs2 = vcat(eqs2, c_spring[i] ~ c_spring0 * (norm1[i] > l_seg))
eqs2 = vcat(eqs2, spring_force[:, i] ~ (c_spring[i] * (norm1[i] - l_seg) + damping * spring_vel[i]) * unit_vector[:, i])
if i == segments
eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i])
else
eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i]- spring_force[:, i+1])
end
eqs2 = vcat(eqs2, acc[:, i+1] .~ G_EARTH + total_force[:, i] / mass)
end
``````

This is working indeed when you have a chain of tethers. What I was wondering about was how to connect multiple tethers to the same point, leading to additional terms in the same equation. A bit like a

``````for i in tethers
total_force .+~ spring_force[i]
end
``````

I do not know if that could work