ModelingToolkitStandardLibrary SampledData Component

I’ve been trying to understand how to use the SampledData component from ModelingToolkitStandardLibrary, but I’ve been having trouble getting the example from documentation working. I’ve modified the code from the example to get it in a copy/paste runnable state up until the solve, since the code in the documentation didn’t seem to work if I just copied/pasted on my machine.

using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks
using DataInterpolations
using DifferentialEquations

function System(; name)
    @named src = SampledData(Float64)
            @variables t
    vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
    pars = @parameters m=10 k=1000 d=1
    D = Differential(t)

    eqs = [f ~ src.output.u
        ddx * 10 ~ k * x + d * dx + f
        D(x) ~ dx
        D(dx) ~ ddx]

    ODESystem(eqs, t, vars, pars; systems = [src], name)
end


dt = 4e-4
times = 0:dt:0.1

@named system = System()
sys = structural_simplify(system)
s = complete(system)
prob = ODEProblem(sys, [], (0, times[end]); tofloat = false)
defs = ModelingToolkit.defaults(sys)

function get_prob(data)
    defs[s.src.buffer] = Parameter(data, dt)
    # ensure p is a uniform type of Vector{Parameter{Float64}} (converting from Vector{Any})
    p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false))
    remake(prob; p)
end

data1 = 3*sin.(2 * pi * (0.5*times) * 100)
data2 = 5*sin.(2 * pi * times * 100)


prob1 = get_prob(data1)
prob2 = get_prob(data2)

#parameter values look okay
prob.p
prob1.p
prob2.p


solve(prob1)


#removed because a single solve should work first
#=
sol1 = Ref{ODESolution}()
sol2 = Ref{ODESolution}()
@sync begin
    @async sol1[] = solve(prob1, ImplicitEuler())
    @async sol2[] = solve(prob2, ImplicitEuler())
end
=#

solve(prob1) gives the following error:

ERROR: MethodError: no method matching getindex(::Parameter{Float64}, ::Int64)

The stacktrace is too long to include in this post since the max char length is about half what it needs to be for this error. I’ve put it here instead, hopefully that’s okay. MTSL SampledData Error · GitHub

I feel like it should be trying to get something from an array or vector of ::Parameter{Float64}, but I’m not sure how I’d achieve that since the parameter looks right when inspected in the problem.

It seems like remake always turns a parameter tuple into a parameter array, this breaks the generated code in this example. There must be an automated way to find that parameter index, but I don’t see it yet. Here is a workaround:

using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks
using DifferentialEquations

function System(; name)
    # Provide an empty array since it will be filled in later
    @named src = SampledData(Float64[], 1.0)
    @variables t
    vars = @variables f(t) = 0 x(t) = 0 dx(t) = 0 ddx(t) = 0
    pars = @parameters m = 10 k = 1000 d = 1
    D = Differential(t)

    eqs = [
        f ~ src.output.u
        ddx * 10 ~ k * x + d * dx + f
        D(x) ~ dx
        D(dx) ~ ddx
    ]

    ODESystem(eqs, t, vars, pars; systems = [src], name)
end

function main()

    dt = 4e-4
    times = 0:dt:0.1
    data1 = 3 * sin.(2 * pi * (0.5 * times) * 100)
    data2 = 5 * sin.(2 * pi * times * 100)

    @named system = System()
    sys = structural_simplify(system)
    s = complete(sys)
    prob = ODEProblem(s, [s.src₊sample_time=>dt], (0.0, times[end]); tofloat=false)

    # make a copy of the problem and change the parameter of interest.
    prob1 = deepcopy(prob)
    prob1.p[2][1] = data1
    prob2 = deepcopy(prob)
    prob2.p[2][1] = data2

    #parameter values look okay
    @show prob.p
    @show prob1.p
    # @show prob2.p

    solve(prob1)
end
1 Like

This looks like it’s an issue with remake after the heterogenous tuple change, @YingboMa

Thanks for the clarification! I assume for other systems I would just manually find the parameter of interest and redefine just that, avoiding remakes for now.