Support for Inputs

[bradcarman]

There should be official support for the concept of "inputs" to a model. I have hacked together the following code to provide a start for the following features/requirements:

• provide a way to define model inputs so they are distinguishable from normal "unknowns". For example, inputs(sys::System) should list the inputs. Also an @inputs macro should be available to define model inputs. Input variables should be set directly and shouldn't change with time unless directed to, i.e. an input x has the equation D(x) ~ 0
• provide an API so that inputs are separate from a problem, such that one can call solve(prob, inputs). This way, if inputs change, the ODEProblem can remain the same.
• provide an API so that inputs can be set in an indeterminate form with set_input(integrator, var::Num, value::Number)

In the below example is a start to this feature request, the solve(prob, inputs) and set_input(integrator, var, value) are demonstrated. What's missing is the ability to influence initialization of the inputs from data. As can be seen, this provides the ability to run the same model in determinate and indeterminate forms and get the same result.

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using Plots
using DiffEqBase: reeval_internals_due_to_modification!, savevalues!
using Test

struct Input
    var::Num
    data::Vector{Float64}
    time::Vector{Float64}
end

function DiffEqBase.solve(prob::SciMLBase.AbstractDEProblem, inputs::Union{Input, Vector{Input}}, args...; kwargs...)

    tstops = Float64[]
    callbacks = DiscreteCallback[]
    if !isa(inputs, Vector)
        inputs = [inputs]
    end

    for input::Input in inputs
        tstops = union(tstops, input.time)

        condition(u,t,integrator) = any(t .== input.time)
        function affect!(integrator) 
            i = findfirst(integrator.t .== input.time)
            integrator[input.var] = input.data[i]
        end
        callback = DiscreteCallback(condition, affect!)

        push!(callbacks, callback)
    end


    return solve(prob, args...; tstops, callback=CallbackSet(callbacks...), kwargs...)
end

function set_input!(integrator, var, value)
    savevalues!(integrator)
    integrator[var] = value
    reeval_internals_due_to_modification!(integrator, false) 
    savevalues!(integrator, true)
end

# -----------------------------------------
# ----- example ---------------------------
# -----------------------------------------

vars = @variables begin
    # inputs
    x(t) = 1

    # states
    y(t) = 0
end

eqs = [
    # inputs    
    D(x) ~ 0
   
    # equations
    D(y) ~ x
]

@mtkcompile sys = System(eqs, t, vars, [])
prob = ODEProblem(sys, [], (0, 4))
input = Input(sys.x, [2,3,4], [1,2,3])

sol = solve(prob, input)


plot(sol; idxs=[sys.x, sys.y])


integrator = init(prob)

for i=1:4
    # set inputs
    set_input!(integrator, sys.x, i)

    # step
    step!(integrator, 1.0, true)

    # outputs
    integrator[sys.y] # |> f
end



p1= plot(integrator.sol; idxs=sys.x, label="indeterminate (x)")
plot!(sol; idxs=sys.x, label="determinate (x)")

p2 = plot(integrator.sol; idxs=sys.y, label="indeterminate (y)")
plot!(sol; idxs=sys.y, label="determinate (y)")

plot(p1, p2; layout=(2,1))

[AayushSabharwal]

provide a way to define model inputs so they are distinguishable from normal "unknowns". For example, inputs(sys::System) should list the inputs.

This exists:

@variables x(t) [input = true]

any such variables will be listed in inputs(sys). To simplify the system with these variables as inputs, they must be provided to structural_simplify(v9)/mtkcompile(v10) appropriately.

• structural_simplify(sys, (input_vars, output_vars))
• mtkcompile(sys; inputs = input_vars, outputs = output_vars)

Also an @inputs macro should be available to define model inputs.

@inputs is a nice idea.

Input variables should be set directly and shouldn't change with time unless directed to, i.e. an input x has the equation D(x) ~ 0

That is how they behave currently.

provide an API so that inputs are separate from a problem, such that one can call solve(prob, inputs). This way, if inputs change, the ODEProblem can remain the same.

The solve(prob, inputs) method in your code snippet can be implemented using MTK's SymbolicDiscreteCallback.

provide an API so that inputs can be set in an indeterminate form with set_input(integrator, var::Num, value::Number)

This I don't fully understand. Why does set_input need to exist? Why can't the symbolic affect function "pull" the input values from a data source rather than having to manually "push" them to the integrator via a set_input function.

[bradcarman]

Not every model will know what the input data is or how long it needs to run before the ODEProblem is compiled. A classic example is a digital twin that needs to mirror what is happening with it's physical counterpart. For a digital twin that is running continuously there needs to exist a way to update the inputs to the model from physical sensors. Furthermore, there needs to exist a way to input data, step the model, and then take action on the model outputs. In this way a model from MTK can be integrated into digital twin algorithms, or any other algorithms that need to intertwine the model inputs and outputs with other functions sending and receiving information.

Note: the ability to support both determinate and indeterminate forms of the model means (I think) that we don't want to implement SymbolicDiscreteCallback because in the case of indeterminate we don't use a callback. My code snippet is showing how we can implement the callback at the solve step so the callback is not intertwined with the model itself.

[bradcarman]

I attempted to use the inputs functionality you mentioned already existed. However I run into a problem...

vars = @variables begin
    # inputs
    x(t) = 1, [input=true]

    # states
    y(t) = 0
end

eqs = [  
    # equations
    D(y) ~ x
]

@mtkcompile sys = System(eqs, t, vars, []) inputs = [x]
prob = ODEProblem(sys, [], (0, 4))
input = Input(sys.x, [2,3,4], [1,2,3])

sol = solve(prob, input)

I can see, the variable x(t) is moved to be a parameter. So the following function is changed...

function DiffEqBase.solve(prob::SciMLBase.AbstractDEProblem, inputs::Union{Input, Vector{Input}}, args...; kwargs...)

    tstops = Float64[]
    callbacks = DiscreteCallback[]
    if !isa(inputs, Vector)
        inputs = [inputs]
    end

    for input::Input in inputs
        tstops = union(tstops, input.time)

        condition(u,t,integrator) = any(t .== input.time)
        function affect!(integrator) 
            i = findfirst(integrator.t .== input.time)
            integrator.ps[input.var] = input.data[i] #<--- here now using `ps`
        end
        callback = DiscreteCallback(condition, affect!)

        push!(callbacks, callback)
    end


    return solve(prob, args...; tstops, callback=CallbackSet(callbacks...), kwargs...)
end

But I get the following error:

julia> sol = solve(prob, input)
ERROR: Indexing with parameters is deprecated. Use `integrator.ps[x(t)] = 2.0` for parameter indexing.

[AayushSabharwal]

Not every model will know what the input data is or how long it needs to run before the ODEProblem is compiled. A classic example is a digital twin that needs to mirror what is happening with it's physical counterpart. For a digital twin that is running continuously there needs to exist a way to update the inputs to the model from physical sensors.

That still works - the callback affect needs to grab the data from the sensor instead of from a struct. Just [x ~ get_sensor_data()] in the callback.

My code snippet is showing how we can implement the callback at the solve step so the callback is not intertwined with the model itself.

Not tying the callback to the model is interesting, but the set_input! function in your example is effectively manually implementing a callback. It also doesn't account for reinitialization after the callback. What if updating the input causes an algebraic equation to not be satisfied?

But I get the following error:

What is the stacktrace for the error?

[bradcarman]

That still works - the callback affect needs to grab the data from the sensor instead of from a struct. Just [x ~ get_sensor_data()] in the callback.

The callback is applied at t=time+dt, not at t=time. I have studied this extensively. If the sensor knowledge exists for the current time, then a callback cannot be used. You can only use the callback when you already know what the input will be for all times that will be solved in advanced. By definition, with a digital twin for example, we don't know the future.

What if updating the input causes an algebraic equation to not be satisfied?

Maybe I haven't caught every edge case here, but by definition an input should be unconstrainted. If an input causes a model to crash because algebraic equations are unsatisfied, this is essentially a user error. A model crash from bad inputs is actually a useful thing, it informs the user that either they have put in a bad input, or the model was improperly constructed. For example, if a sensor gives a bad reading, maybe because the sensor itself broke, then a model error from the input would be informative, highlighting the fact that something is wrong, which would indeed be the case.

Stack trace: Note: step3 (inputs.jl:28) is line: integrator.ps[input.var] = input.data[i]

  [1] error(s::String)
    @ Base .\error.jl:35
  [2] setindex!(A::OrdinaryDiffEqCore.ODEIntegrator{…}, val::Float64, sym::Num)
    @ SciMLBase C:\Users\bradl\.julia\packages\SciMLBase\zk34N\src\integrator_interface.jl:523
  [3] (::var"#affect!#16"{…})(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ Main c:\Work\Packages\ModelingToolkit.jl\inputs.jl:28
  [4] apply_discrete_callback!
    @ C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\callbacks.jl:616 [inlined]
  [5] handle_callbacks!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore C:\Users\bradl\.julia\packages\OrdinaryDiffEqCore\zs1s7\src\integrators\integrator_utils.jl:396
  [6] _loopfooter!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore C:\Users\bradl\.julia\packages\OrdinaryDiffEqCore\zs1s7\src\integrators\integrator_utils.jl:284
  [7] loopfooter!
    @ C:\Users\bradl\.julia\packages\OrdinaryDiffEqCore\zs1s7\src\integrators\integrator_utils.jl:248 [inlined]
  [8] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore C:\Users\bradl\.julia\packages\OrdinaryDiffEqCore\zs1s7\src\solve.jl:621
  [9] #__solve#62
    @ C:\Users\bradl\.julia\packages\OrdinaryDiffEqCore\zs1s7\src\solve.jl:7 [inlined]
 [10] __solve
    @ C:\Users\bradl\.julia\packages\OrdinaryDiffEqCore\zs1s7\src\solve.jl:1 [inlined]
 [11] solve_call(_prob::ODEProblem{…}, args::CompositeAlgorithm{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:670
 [12] solve_call
    @ C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:627 [inlined]
 [13] #solve_up#45
    @ C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:1224 [inlined]
 [14] solve_up
    @ C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:1201 [inlined]
 [15] #solve#43
    @ C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:1100 [inlined]
 [16] solve
    @ C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:1086 [inlined]
 [17] #__solve#3
    @ C:\Users\bradl\.julia\packages\OrdinaryDiffEqDefault\kOV75\src\default_alg.jl:48 [inlined]
 [18] __solve
    @ C:\Users\bradl\.julia\packages\OrdinaryDiffEqDefault\kOV75\src\default_alg.jl:47 [inlined]
 [19] #__solve#64
    @ C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:1563 [inlined]
 [20] __solve
    @ C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:1554 [inlined]
 [21] solve_call(::ODEProblem{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:670
 [22] solve_up(::ODEProblem{…}, ::Nothing, ::Vector{…}, ::MTKParameters{…}; originator::SciMLBase.ChainRulesOriginator, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:1208
 [23] solve_up
    @ C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:1201 [inlined]
 [24] solve(::ODEProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\bradl\.julia\packages\DiffEqBase\dZqFN\src\solve.jl:1096
 [25] solve(::ODEProblem{…}, ::Input; kwargs::@Kwargs{})
    @ Main c:\Work\Packages\ModelingToolkit.jl\inputs.jl:36
 [26] solve(::ODEProblem{…}, ::Input)
    @ Main c:\Work\Packages\ModelingToolkit.jl\inputs.jl:14
 [27] top-level scope
    @ c:\Work\Packages\ModelingToolkit.jl\inputs.jl:67

[ChrisRackauckas]

The callback is applied at t=time+dt, not at t=time. I have studied this extensively. If the sensor knowledge exists for the current time, then a callback cannot be used. You can only use the callback when you already know what the input will be for all times that will be solved in advanced. By definition, with a digital twin for example, we don't know the future.

Yes but that's what initialize is for... just call it... You literally just do one call to up the pointer lol

[cstjean]

I've made it work here too with the integrator interface + pushing onto a DataInterpolations.ConstantInterpolation, but it's been challenging to piece together everything I needed.

A classic example is a digital twin that needs to mirror what is happening with it's physical counterpart. For a digital twin that is running continuously there needs to exist a way to update the inputs to the model from physical sensors.

If there's no higher-level interface planned for this, then a digital twin example would be a very helpful addition to the docs.

[ChrisRackauckas]

I think the simplest way to expose this would be pre-made components that are setup with certain channel constructors. Most of the features being asked for in this domain are things about blocking vs non-blocking, timing / clocking, getting data from different sources (files, HTTP, SSH, etc.), and these are all covered by channels, it's just not a normal thing for modelers to dig into. But if you use them it does work out:

using OrdinaryDiffEqTsit5, Setfield
function f1(m, observed, ctx, integrator)
    @show integrator.t
    input = parse(Float64,readline())
    if input == -Inf
        terminate!(integrator)
        @set! m.x = 0.0
    else
        @set! m.x = input
    end
end
@register_symbolic f1(t)

pars = @parameters begin
    # inputs
    x(t) = 1
end

vars = @variables begin
    # states
    y(t) = 0
end

eqs = [  
    # equations
    D(y) ~ x
]

data_read = ModelingToolkit.SymbolicDiscreteCallback(1.0, 
    ModelingToolkit.ImperativeAffect(f1, modified=(;x)))
sys = System(eqs, t; name=:sys, discrete_events=[data_read])
ssys = structural_simplify(sys)
prob = ODEProblem(ssys, [y=>0.0, x=>0.0], (0.0, Inf))
sol = solve(prob, Tsit5())

[bradcarman]

Thank you Chris for this example, unfortunately this does not satisfy the requirement. It seems to, but under close examination, it provides the following execution order:

1. integrator steps to time=t+dt
2. integrator waits for data input at time=t+dt
3. step is finalized and computed outputs are available for analysis at time=t+dt

What this means is I'm not able to do:

1. provide data input at time=t
2. step to time=t+dt
3. analyze computed outputs for analysis at time=t+dt

As a proof point, take the following...

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using Setfield

input_ref = Ref(1.0)
function fref(m, observed, ctx, integrator)
    @set! m.x = input_ref[]
    return m
end

pars = @parameters begin
    # inputs
    x(t) = 1
end

vars = @variables begin
    # states
    y(t) = 0
end

eqs = [  
    # equations
    D(y) ~ x
]

data_read_ref = ModelingToolkit.SymbolicDiscreteCallback(1.0, 
    ModelingToolkit.ImperativeAffect(fref, modified=(;x)))

@mtkcompile sys = System(eqs, t; name=:sys, discrete_events=[data_read_ref])
prob = ODEProblem(sys, [], (0.0, Inf))
integrator = init(prob)

# time = 0
input_ref[] = 1.0
step!(integrator, 1.0, true)
integrator.sol(0.0; idxs=x) #1.0
integrator.sol(0.5; idxs=x) #1.0
integrator.sol(1.0; idxs=x) #1.0

# time = 1.0
input_ref[] = 2.0
step!(integrator, 1.0, true)
integrator.sol(1.0; idxs=x) #1.0 <--- need this to be 2.0!
integrator.sol(1.5; idxs=x) #1.0 <--- need this to be 2.0!
integrator.sol(2.0; idxs=x) #2.0 <--- too late!!!

[baggepinnen]

This is the interface we use for inputs when doing state estimation and control. It uses generate_control_function to generate dynamics on the form $\dot x = f(x, u, p, t)$ and this can then be discretized with any integrator of choice, the example below uses SeeToDee

using ModelingToolkit, Test
using ModelingToolkit: t_nounits as t, D_nounits as D
using SeeToDee

pars = @parameters begin
end

vars = @variables begin
    u(t) = 1
    x(t) = 0
end

eqs = [  
    # equations
    D(x) ~ u
]

@named sys = System(eqs, t)

inputs = [u]
f, xs, ps, iosys = ModelingToolkit.generate_control_function(sys, inputs)
Ts = 1.0 # Sample interval
f_disc = SeeToDee.Rk4(f[1].f_oop, Ts) # This is now a function that integrates Ts=1.0s forwards in time

X, _ = ModelingToolkit.get_u0(iosys, [])
p = ModelingToolkit.get_ps(iosys)

# time = 0
input = 1.0
T = 0.0
X = f_disc(X, input, p, T)
@test X ≈ [1.0] # Integrate the input 1 for 1 second to get 1

# time = 1.0
input = 2.0
T += Ts
X = f_disc(X, input, p, T)
@test X ≈ [3.0] # Additionally integrate the input 2 for 1 second to get 3

[bradcarman]

I'm happy to report I have a workaround solution that is working for me for now. Thanks to @AayushSabharwal and @ChrisRackauckas for the help...

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using Plots
using SymbolicIndexingInterface
using Setfield
using Test

struct Input
    var::Num
    data::Vector{Float64}
    time::Vector{Float64}
end

function DiffEqBase.solve(prob::SciMLBase.AbstractDEProblem, inputs::Union{Input, Vector{Input}}, args...; input_funs, kwargs...)

    set_input!, finalize! = input_funs

    tstops = Float64[]
    callbacks = DiscreteCallback[]
    if !isa(inputs, Vector)
        inputs = [inputs]
    end

    for input::Input in inputs
        tstops = union(tstops, input.time)
        
        condition = (u,t,integrator) -> any(t .== input.time)
        affect! = function (integrator) 
            i = findfirst(integrator.t .== input.time)
            set_input!(integrator, input.var, input.data[i])
        end
        callback = DiscreteCallback(condition, affect!)

        push!(callbacks, callback)
    end

    # finalize!
    t_end = prob.tspan[2]
    condition = (u,t,integrator) -> (t == t_end)
    affect! = (integrator) -> finalize!(integrator)
    callback = DiscreteCallback(condition, affect!)

    push!(callbacks, callback)
    push!(tstops, t_end)

    return solve(prob, args...; tstops, callback=CallbackSet(callbacks...), kwargs...)
end

function setup_inputs(sys)

    inputs = ModelingToolkit.unbound_inputs(sys)
    setters = Dict{Num, Function}()

    if !isempty(inputs)
        sdcs = ModelingToolkit.SymbolicDiscreteCallback[]
        for x in inputs
            affect = ModelingToolkit.ImperativeAffect((m, o, c, i)->m, modified=(;x))
            sdc = ModelingToolkit.SymbolicDiscreteCallback(Inf, affect)

            push!(sdcs, sdc)
        end

        @set! sys.discrete_events = sdcs
        sys = complete(sys)  # @set! sys.index_cache = ModelingToolkit.IndexCache(sys)

        for (i,x) in enumerate(inputs)
            setter = SymbolicIndexingInterface.setsym(sys, x)
            sdc = sdcs[i]

            setval = function (integrator, set, val=NaN)
                if set
                    println("::setting $x to $val @ $(integrator.t)s")
                    setter(integrator, val)
                else
                    println("::saving $x @ $(integrator.t)s")
                end
                ModelingToolkit.save_callback_discretes!(integrator, sdc)
            end

            setters[x] = setval
        end
    end

    set_input! = function (integrator, var, value) 
        setters[var](integrator, true, value)
        u_modified!(integrator, true)
    end

    finalize! = function (integrator) 
        for ky in keys(setters)
            setters[ky](integrator, false)
        end
    end

    return sys, set_input!, finalize!
end

# -----------------------------------------
# ----- example ---------------------------
# -----------------------------------------

vars = @variables begin
    x(t)=1, [input=true]

    # states
    y(t) = 0
end

eqs = [  
    # equations
    D(y) ~ x

]

@mtkcompile sys = System(eqs, t, vars, []) inputs=[x]
sys, set_input!, finalize! = setup_inputs(sys);
prob = ODEProblem(sys, [], (0, 4))

# indeterminate form -----------------------

integrator = init(prob, Tsit5())

set_input!(integrator, sys.x, 1.0)
step!(integrator, 1.0, true)

set_input!(integrator, sys.x, 2.0)
step!(integrator, 1.0, true)

set_input!(integrator, sys.x, 3.0)
step!(integrator, 1.0, true)

set_input!(integrator, sys.x, 4.0)
step!(integrator, 1.0, true)

finalize!(integrator)

@test integrator.sol(0.0; idxs=sys.x) == 1.0
@test integrator.sol(1.0; idxs=sys.x) == 2.0
@test integrator.sol(2.0; idxs=sys.x) == 3.0
@test integrator.sol(3.0; idxs=sys.x) == 4.0
@test integrator.sol(4.0; idxs=sys.y) ≈ 10.0

# determinate form -----------------------
input = Input(sys.x, [1,2,3,4], [0,1,2,3])
sol = solve(prob, input, Tsit5(); input_funs = (set_input!, finalize!));


@test sol(0.0; idxs=sys.x) == 1.0
@test sol(1.0; idxs=sys.x) == 2.0
@test sol(2.0; idxs=sys.x) == 3.0
@test sol(3.0; idxs=sys.x) == 4.0
@test sol(4.0; idxs=sys.y) ≈ 10.0

This may not be the best interface, but it's a start. What is required here is to call the setup_inputs inbetween building the MTK system and building the ODEProblem