[MTK] Interpolating observed variables

Specific Domains Modelling & Simulations
1

I am toying around with ModelingToolkit, I have a very simple algebraic model representing the voltage/current relationship in a ideal resistor:

using ModelingToolkit, OrdinaryDiffEq

@parameters R
@variables t v(t) i(t)
D = Differential(t)

eqs = [
   i ~ v/R
   D(v) ~ 0
]

@named circuit = ODESystem(eqs, t)

u0 = [
   v => 5.0,
]

p = [
   R => 10,
]

tspan = (0.0, 0.01)
# convert the DAE problem to an ODE problem
circuit = structural_simplify(circuit)
prob = ODEProblem(circuit, u0, tspan, p)
sol = solve(prob, Tsit5())

# Interpolation can give the solution inside tspan and also further
#sol(0.009)

solve() solves the ODEProblem in the v variable, and I can either v or i as sol[v] or sol[i]. For v(t) I can also find the interpolated values using sol(t). Can the same be done easily for the observed variable i(t)?

2

sol(0.009;idxs=v)

1 Like
3

Thanks, that works indeed, even using idxs=i

4

Thanks. I was having a hard time finding this as well.

Is there a reason, why it’s idxs in some places and vars in others?

OK: just seen that idxs=[statevar1, statevar2] works in the plot recipe as well. So is idxs the preferred way?

1 Like
5

BTW, will the dump function create a useable form of the ODE solution that can be reread at a later time?

Since it is ASCII, it may be a way to keep the information even if the binary format changes.

6

We should fix that. vars is in the plot recipe, everything else is idxs.

The two-argument show method should?

7

What do you mean by that? Is there a way to read the output from dump, e.g. out of a text file, back into a julia object? Like a text-based serializer, so to speak?

8

Two argument show is supposed to be roundtrippable.

9

Looks like it (although I don’t know yet how to read it back in).

However, if I look at the output it seems that it does only include the solution vector to the actual ODE problem and the metadata for the solver run, but not the symbolic part that allows calculation of the other state variables.

Looking at the output of dump it has a lot of things like

var"##arg#408" = begin
                    #= /home/thor/.julia/packages/ModelingToolkit/9DDjU/src/structural_transformation/codegen.jl:154 =#
                    (ModelingToolkit.StructuralTransformations.numerical_nlsolve)(var"##fun#405", 0.0, (var"C3₊v(t)", R2₊R, var"C2₊v(t)"))
                end, var"R2₊i(t)" = #= /home/thor/.julia/packages/SymbolicUtils/Kp6RG/src/code.jl:169 =# @inbounds(var"##arg#408"[1]), var"##fun#409" = RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#410"), Symbol("##arg#411")), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0xb9251361, 0xd947528e, 0x1c66eb19, 0xd3949734, 0x4cd5ad53)}(quote
    #= /home/thor/.julia/packages/SymbolicUtils/Kp6RG/src/code.jl:282 =#
    #= /home/thor/.julia/packages/SymbolicUtils/Kp6RG/src/code.jl:283 =#
    let var"R3₊i(t)" = #= /home/thor/.julia/packages/SymbolicUtils/Kp6RG/src/code.jl:169 =# @inbounds(var"##arg#410"[1]), var"C3₊v(t)" = #= /home/thor/.julia/packages/SymbolicUtils/Kp6RG/src/code.jl:169 =# @inbounds(var"##arg#411"[1]), R3₊R = #= /home/thor/.julia/packages/SymbolicUtils/Kp6RG/src/code.jl:169 =# @inbounds(var"##arg#411"[2])
        (+)((*)(R3₊R, var"R3₊i(t)"), (*)(-1, var"C3₊v(t)"))
    end

in it, which - to the uninitiated, i.e. me - looks like the additional algebraic equations for the other states.

But this gets too far of topic on a closed thread. So maybe better to move the discussion to new topic.

10

Yes, new topic. You’re looking for prob.f.observed, which is also a RuntimeGeneratedFunction so it will similarly print out it’s code.