doc: updating MTK example to v11 in manual.

[franckgaga]

Hello @baggepinnen,

I'm trying to update the MTK example to v11 but i'm stuck at the varmat_to_vars function here:

using ModelPredictiveControl, ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t, varmap_to_vars
@parameters g=9.8 L=0.4 K=1.2 m=0.3
@variables θ(t)=0 ω(t)=0 τ(t)=0 y(t)
eqs = [
    D(θ) ~ ω
    D(ω) ~ -g/L*sin(θ) - K/m*ω + τ/m/L^2
    y ~ θ * 180 / π
]
@named mtk_model = System(eqs, t)
function generate_f_h(model, inputs, outputs)
    (_, f_ip), x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function(
        model, inputs, split=false, simplify=true
    )
    if any(ModelingToolkit.is_alg_equation, equations(io_sys)) 
        error("Systems with algebraic equations are not supported")
    end
    nu, nx, ny = length(inputs), length(x_sym), length(outputs)
    function f!(ẋ, x, u, _ , p)
        try
            f_ip(ẋ, x, u, p, nothing)
        catch err
            if err isa MethodError
                error("NonLinModel does not support a time argument t in the f function, "*
                      "see the constructor docstring for a workaround.")
            else
                rethrow()
            end
        end
        return nothing
    end
    (_, h_ip) = ModelingToolkit.build_explicit_observed_function(
        io_sys, outputs; inputs, return_inplace = true
    )
    u_nothing = fill(nothing, nu)
    function h!(y, x, _ , p)
        try
            # MTK.jl supports a `u` argument in `h_ip` function but not this package. We set
            # `u` as a vector of nothing and `h_ip` function will presumably throw an
            # MethodError if this argument is used inside the function
            h_ip(y, x, u_nothing, p, nothing)
        catch err
            if err isa MethodError
                error("NonLinModel only support strictly proper systems (no manipulated "*
                      "input argument u in the output function h)")
            else
                rethrow()
            end
        end
        return nothing
    end
    println(bindings(io_sys))
    p = varmap_to_vars(bindings(io_sys), p_sym) # crashes here
    return f!, h!, p, x_sym, nu, nx, ny
end
inputs, outputs = [τ], [y]
f!, h!, p, x_sym, nu, nx, ny = generate_f_h(mtk_model, inputs, outputs)
x_sym

this code prints:

ReadOnlyDicts.ReadOnlyDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, ModelingToolkitBase.AtomicArrayDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}}}(τ(t) => missing)
ERROR: Initial condition underdefined. Some are missing from the variable map.
Please provide a default (`u0`), initialization equation, or guess
for the following variables:

L, K, m, g

Stacktrace:
 [1] varmap_to_vars(varmap::ReadOnlyDicts.ReadOnlyDict{…}, vars::Vector{…}; tofloat::Bool, use_union::Bool, container_type::Type, buffer_eltype::Type, toterm::Function, check::Bool, allow_symbolic::Bool, is_initializeprob::Bool, substitution_limit::Int64, missing_values::ModelingToolkitBase.MissingGuessValue.var"typeof(MissingGuessValue)")
   @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/uIKoY/src/systems/problem_utils.jl:378
 [2] varmap_to_vars
   @ ~/.julia/packages/ModelingToolkitBase/uIKoY/src/systems/problem_utils.jl:331 [inlined]
 [3] generate_f_h(model::System, inputs::Vector{Num}, outputs::Vector{Num})
   @ Main ~/.julia/dev/ModelPredictiveControl/docs/src/manual/mtk.md:89
 [4] top-level scope
   @ ~/.julia/dev/ModelPredictiveControl/docs/src/manual/mtk.md:93
Some type information was truncated. Use `show(err)` to see complete types.

Do you know why bindings(io_sys) does not contain the parameter values anymore ?

[codecov-commenter]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.57%. Comparing base (b9c5c77) to head (53c94fa).
⚠️ Report is 8 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #340   +/-   ##
=======================================
  Coverage   98.57%   98.57%           
=======================================
  Files          27       27           
  Lines        5190     5190           
=======================================
  Hits         5116     5116           
  Misses         74       74           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[baggepinnen]

Try

p = ModelingToolkit.get_p(iosys, op_dict)

[franckgaga]

Thanks for the answer but using p = ModelingToolkit.get_p(io_sys, bindings(io_sys)) returned the error:

ERROR: AssertionError: `promote_to_concrete` can't make type Missing uniform with Float64
Stacktrace:
 [1] promote_to_concrete(vs::Vector{Union{Missing, Real}}; tofloat::Bool, use_union::Bool)
   @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/uIKoY/src/utils.jl:1091
 [2] varmap_to_vars(varmap::ModelingToolkitBase.AtomicArrayDict{…}, vars::Vector{…}; tofloat::Bool, use_union::Bool, container_type::Type, buffer_eltype::Type, toterm::Function, check::Bool, allow_symbolic::Bool, is_initializeprob::Bool, substitution_limit::Int64, missing_values::ModelingToolkitBase.MissingGuessValue.var"typeof(MissingGuessValue)")
   @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/uIKoY/src/systems/problem_utils.jl:418
 [3] varmap_to_vars
   @ ~/.julia/packages/ModelingToolkitBase/uIKoY/src/systems/problem_utils.jl:331 [inlined]
 [4] get_p(sys::System, varmap::ReadOnlyDicts.ReadOnlyDict{…}; split::Bool, kwargs::@Kwargs{})
   @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/uIKoY/src/systems/problem_utils.jl:2059
 [5] get_p
   @ ~/.julia/packages/ModelingToolkitBase/uIKoY/src/systems/problem_utils.jl:2048 [inlined]
 [6] generate_f_h(model::System, inputs::Vector{Num}, outputs::Vector{Num})
   @ Main ~/.julia/dev/ModelPredictiveControl/docs/src/manual/mtk.md:91
 [7] top-level scope
   @ ~/.julia/dev/ModelPredictiveControl/docs/src/manual/mtk.md:96
Some type information was truncated. Use `show(err)` to see complete types.

If I do print(bindings(io_sys)), it prints a dict with only the τ key and the value is missing:

ReadOnlyDicts.ReadOnlyDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, ModelingToolkitBase.AtomicArrayDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}}}(τ(t) => missing)

[baggepinnen]

You should pass a varmap as the second argument, bindings should not be used

[franckgaga]

And how do I get the varmap of the 4 parameters ? I tried passing p_sym directly and it does not work, the error is ArgumentError: The operating_point passed to the problem constructor must be a symbolic map.

Nothing about this in the documentation also.

[baggepinnen]

a varmap is just

op = Dict(p1 => val1, p2 => val2, ...)

[franckgaga]

Alright thanks, but I still need to access the current value that I set at the line @parameters g=9.8 L=0.4 K=1.2 m=0.3, how do I access them ?

[baggepinnen]

This is how I typically do it
https://github.com/baggepinnen/LowLevelParticleFiltersMTK.jl/blob/main/src/LowLevelParticleFiltersMTK.jl#L122C10-L124C20

[franckgaga]

Not sure this is the good way of proceeding, constructing an ODEProblem just to extract the parameter vector.

I read the wall of text about this in the documentation here : https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/

Two observations:

1. I'm completely lost, way too long and completely unreadable by someone who is not familiar with MTK internals. Pretty sure a good portion of it was written by LLMs...
2. It seems that the only example about the initialization of parameters defined with @parameters p1=1 p2=4.3 is simply disabled and not executed during the doc deployment. My feeling is this feature is not ready. So I think the right move is to wait before updating my pendulum example to v11.

[baggepinnen]

The ODEProblem constructor is a thin wrapper around initialization so it isn't really much more expensive than just solving the initialization problem. The issue with your old approach is that parameters stored in the system may not be limited to simple values only, they may be the solution of a nonlinear problem (referred to as the "initialization problem")

[franckgaga]

I'm not sure I'm following you, how it could be a nonlinear root solving problem for parameters (as in non-differential and non-algebraic states) ?

[baggepinnen]

It's possible in MTK to specify initialization equations that must hold at initialization and solve for parameters that satisfy those equations. For example

@parameters constant_input=missing

equations = [
    D(x) ~ x + p
]

initialization_equations = [
    D(x) ~ 0
]

to solve for the value of p that satisfies the steady state equation D(x) ~ 0

[franckgaga]

It's cool and all that it is now supported in MTK but as I wrote in the disclaimer, the goal of this example is to provide a starting template to merge both packages, not something that cover all the corner cases. And this example focus on parameters that are well-defined and known, as this is the most common cases. I will try to come up to something that work well for this specific case.

[franckgaga]

Okay this what I do now:

function generate_f_h(model, inputs, outputs)
    (_, f_ip), x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function(
        model, inputs, split=false, simplify=true
    )
    if any(ModelingToolkit.is_alg_equation, equations(io_sys)) 
        error("Systems with algebraic equations are not supported")
    end
    nu, nx, ny = length(inputs), length(x_sym), length(outputs)
    function f!(ẋ, x, u, _ , p)
        try
            f_ip(ẋ, x, u, p, nothing)
        catch err
            if err isa MethodError
                error("NonLinModel does not support a time argument t in the f function, "*
                      "see the constructor docstring for a workaround.")
            else
                rethrow()
            end
        end
        return nothing
    end
    (_, h_ip) = ModelingToolkit.build_explicit_observed_function(
        io_sys, outputs; inputs, return_inplace = true
    )
    u_nothing = fill(nothing, nu)
    function h!(y, x, _ , p)
        try
            # MTK.jl supports a `u` argument in `h_ip` function but not this package. We set
            # `u` as a vector of nothing and `h_ip` function will presumably throw an
            # MethodError if this argument is used inside the function
            h_ip(y, x, u_nothing, p, nothing)
        catch err
            if err isa MethodError
                error("NonLinModel only support strictly proper systems (no manipulated "*
                      "input argument u in the output function h)")
            else
                rethrow()
            end
        end
        return nothing
    end
    # this is how I extract the correct `p` vector:
    ic = initial_conditions(io_sys)
    p_map = Dict(sym => ic[sym] for sym in p_sym)
    p = ModelingToolkit.varmap_to_vars(p_map, p_sym)
    return f!, h!, p, x_sym, p_sym, nu, nx, ny
end

It work well on MTK.jl v11. The initial_conditions method returns the values defined at the @parameters line, as described in MTK.jl NEWS.md:

"Defaults" specified via variable metadata are now translated into either initial_conditions or bindings depending on the value. If the value is a constant, it is part of initial_conditions. If it is an expression involving other variables/parameters, it is part of bindings. For example, the following are initial_conditions:

@variables x(t) = 1 y(t)[1:3] = zeros(3)
@parameters f(::Real) = sin

[...]

As expected, the p vector only include the 4 parameters, in the same order as defined by p_sym. It will presumably crashes if the parameter are not well defined, as explained above, but I'm okay with this.

edit: I will add a clear error message in such cases.

[franckgaga]

Oupsie there's a bug in the version that I just released. The order in p_sym will be random from run to another, thus $K$ is not always the third element. I will debug this tomorrow.

[baggepinnen]

In my experience, using a non-supported way of obtaining the parameter object leads to a constant maintenance burden.

Your three lines of code

    ic = initial_conditions(io_sys)
    p_map = Dict(sym => ic[sym] for sym in p_sym)
    p = ModelingToolkit.varmap_to_vars(p_map, p_sym)

could be just

p = ModelingToolkit.ODEProblem(iosys, [], (0.0, 1.0)).p

and it would work in more cases and be less likely to break on you.

[franckgaga]

How do you know this is the official way of doing this ? I did not found anything on this in the documentation.

Also, the returned p with this solution include 12 elements instead of 4 elements:

    p = ModelingToolkit.ODEProblem(io_sys, [], (0.0, 1.0)).p
    display(p)

displaying:

12-element Vector{Float64}:
 0.4
 0.3
 1.2
 9.8
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

It does work but this is misleading for the user, since the p_sym object informs the user on what is the order of the parameter in model.p (in case he want to modify them later). But p_sym includes 4 elements and model.p includes 12 elements.

[franckgaga]

How do you know this is the official way of doing this ? I did not found anything on this in the documentation.

Okay I see that's how initialization problems are solved in this section https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/#bindings_and_ics, so that's presumably the official way.

Do you know why p has 12 elements instead of 4 ?

[baggepinnen]

Do you know why p has 12 elements instead of 4 ?

MTK automatically adds parameters for initial conditions and inputs IIRC, possibly some additional things

[franckgaga]

Okay thanks. As I said above, if the user want to know how the returned p vector is sorted, to modify it later, p_sym has 4 elements so this is misleading for him. Except if you can confirm that I can safely do:

np = length(p_sym)
p = ModelingToolkit.ODEProblem(io_sys, [], (0.0, 1.0)).p[1:np]

I will stick to my unofficial and limited-in-scope solution.

[baggepinnen]

MTK typically does not guaerantee the order of variables, to modify values, symbolic indexing is supposed to be used. There are some tools listed here https://docs.sciml.ai/ModelingToolkit/stable/examples/remake/#Re-creating-the-problem

for example

remake(odeprob; p = [α => x[1], β => x[2], γ => x[3], δ => x[4]])

this does not rely on the order, it uses a symbolic map sym => val

[franckgaga]

Yeah I get this, but when I call (f_oop, f_ip), x_sym, p_sym, io_sys = generate_control_function( ... ), the order of the variables in p_sym (arbitrarily chosen by MTK) will specify the order of the p vector in the returned f_ip function:

f_ip(xout,x,u,p,t) -> nothing

Am I right ?

If yes, doing the 3 lines of code mentioned in #340 (comment) should work in most situation in which the parameters are constant.