Setu is ignored by solve

[1-Bart-1]

When using setu to update the problem state, this new state is simply ignored by solve. Using setp to update the parameters works fine.

MWE:

using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: D_nounits as D, t_nounits as t, setu, setp, getu, getp

Ts = 0.1
@mtkmodel Pendulum begin
    @parameters begin
        g = 9.8
        L = 0.4
        K = 1.2
        m = 0.3
        τ = 0.0 # input
    end
    @variables begin
        θ(t) = 0.0 # state
        ω(t) = 0.0 # state
        y(t) # output
    end
    @equations begin
        D(θ)    ~ ω
        D(ω)    ~ -g/L*sin(θ) - K/m*ω + τ/m/L^2
        y       ~ θ * 180 / π
    end
end

@mtkbuild mtk_model = Pendulum()
prob = ODEProblem(mtk_model, nothing, (0.0, Ts))
set_x = setu(mtk_model, unknowns(mtk_model))
get_x = getu(mtk_model, unknowns(mtk_model))
set_u = setp(mtk_model, [mtk_model.τ])
get_u = getp(mtk_model, [mtk_model.τ])
get_h = getu(mtk_model, [mtk_model.y])
p = (prob, set_x, set_u, get_h)

function f!(xnext, x, u, _, p)
    (prob, set_x, set_u, _) = p
    set_x(prob, x)
    set_u(prob, u)
    sol = solve(prob, Tsit5(); save_on=false, save_start=false)
    xnext .= sol.u[end] # sol.u is the state, called x in the function
    nothing
end

xnext = zeros(2)

println("Initial")
f!(xnext, zeros(2), 0.0, nothing, p)
@show xnext

println("Just changing state x - doesn't work, xnext stays at zero")
f!(xnext, ones(2), 0.0, nothing, p)
@show xnext

println("Just changing input u - does work, xnext is non-zero")
f!(xnext, zeros(2), 1.0, nothing, p)
@show xnext
nothing

From this discussion

[AayushSabharwal]

See https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/#Constant-constraints-in-initialization. Constant constraints in initialization have to be updated by setting Initial(x) not x.

[1-Bart-1]

Thanks, that works!

[1-Bart-1]

@AayushSabharwal is it correct that setp and setu are basically undocumented? It would be great if this Initial(x) trick was written in the documentation specifically for setp and setu.

[AayushSabharwal]

Initial(x) is an MTK concept. setp and setu are in SymbolicIndexingInterface and documented there.

[ChrisRackauckas]

I don't think this is a sufficient answer. Users using setu would expect u0 to change. If setu needs to update Initial(x) to do that, then setu should do that, not tell users that they need know that we need Initial(x) as a parameter and set that and then later that will set u0. That detail is not public API so users shouldn't have to know or understand it.

[AayushSabharwal]

Initial(x) is public API though? And setu is too low-level of a function to do this. It's basically setindex!. setu does change prob.u0, but MTK specifically doesn't respect prob.u0 ever since initialization was a thing.

[1-Bart-1]

But what if I want MTK to respect prob.u0? I don't want to run an initialization algorithm in my use case, I just want the simulation to start instantly from prob.u0, as I know this is a viable start state.

[1-Bart-1]

Initial(x) is an MTK concept. setp and setu are in SymbolicIndexingInterface and documented there.

I really like the Linearization documentation, which explains clearly that setu is from SymbolicIndexingInterface.jl, and how it can be used to update parameters. This is what I meant with there not being documentation on setu. It would be great if there would be documentation on how to use setu with an ODEProblem.

[AayushSabharwal]

But what if I want MTK to respect prob.u0? I don't want to run an initialization algorithm in my use case, I just want the simulation to start instantly from prob.u0, as I know this is a viable start state.

Sorry for the delay. In such a case, you could pass initializealg = CheckInit() to ODEProblem/remake/solve (either one) to ignore MTK's initialization. In this case, you would have to set x and not Initial(x). Alternatively, see #3570 and its fix.

[AayushSabharwal]

This is also now resolved in MTK@9.76. You should always set Initial(x), even if you just want to run CheckInit or BrownBasic. We will handle updating u0 appropriately from the initials.

[1-Bart-1]

This works with the original mwe I provided, using a prob and solve. But using the integrator interface, I still get the same problem:

using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: D_nounits as D, t_nounits as t, setu, setp, getu, getp

Ts = 0.1
@mtkmodel Pendulum begin
    @parameters begin
        g = 9.8
        L = 0.4
        K = 1.2
        m = 0.3
        τ = 0.0 # input
    end
    @variables begin
        θ(t) = 0.0 # state
        ω(t) = 0.0 # state
        y(t) # output
    end
    @equations begin
        D(θ)    ~ ω
        D(ω)    ~ -g/L*sin(θ) - K/m*ω + τ/m/L^2
        y       ~ θ * 180 / π
    end
end

@mtkbuild mtk_model = Pendulum()
prob = ODEProblem(mtk_model, nothing, (0.0, Ts); initializealg=CheckInit())
integ = init(prob, Tsit5())
set_x = setu(mtk_model, Initial.(unknowns(mtk_model)))
get_x = getu(mtk_model, unknowns(mtk_model))
set_u = setp(mtk_model, [mtk_model.τ])
get_u = getp(mtk_model, [mtk_model.τ])
get_h = getu(mtk_model, [mtk_model.y])
p = (integ, set_x, set_u, get_h)

function f!(xnext, x, u, _, p)
    (integ, set_x, set_u, _) = p
    set_x(integ, x)
    set_u(integ, u)
    reinit!(integ; reinit_dae=true)
    step!(integ)
    xnext .= integ.u # sol.u is the state, called x in the function
    nothing
end

xnext = zeros(2)

println("Initial")
f!(xnext, zeros(2), 0.0, nothing, p)
@show xnext

println("Just changing state x - doesn't work, xnext stays at zero")
f!(xnext, ones(2), 0.0, nothing, p)
@show xnext

println("Just changing input u - does work, xnext is non-zero")
f!(xnext, zeros(2), 1.0, nothing, p)
@show xnext
nothing

output:

Initial
xnext = [0.0, 0.0]
Just changing state x - doesn't work, xnext stays at zero
xnext = [0.0, 0.0]
Just changing input u - does work, xnext is non-zero
xnext = [1.0415277704015974e-7, 0.0020829166371642348]