MTK example no longer work

[franckgaga]

@baggepinnen The MTK example in the manual stop working and I'm not able to pinpoint why. It was working well when I released v1.0.1. It now crashes inside the h_ function (the function returned by ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs) with the stacktrace:

ERROR: BoundsError: attempt to access Float64 at index [2]
Stacktrace:
  [1] getindex
    @ ./number.jl:98 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/SymbolicUtils/ij6YM/src/code.jl:387 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [4] macro expansion
    @ ./none:0 [inlined]
  [5] generated_callfunc
    @ ./none:0 [inlined]
  [6] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::Vector{…}, ::Nothing)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
  [7] (::var"#h!#13"{…})(y::Vector{…}, x::Vector{…}, ::Vector{…}, p::Vector{…})
    @ Main ~/.julia/dev/ModelPredictiveControl/docs/src/manual/mtk.md:86
  [8] h!
    @ ~/.julia/dev/ModelPredictiveControl/src/model/nonlinmodel.jl:211 [inlined]
  [9] evaloutput(model::NonLinModel{…}, d::Vector{…})
    @ ModelPredictiveControl ~/.julia/dev/ModelPredictiveControl/src/sim_model.jl:283
 [10] sim_closedloop!(est_mpc::NonLinMPC{…}, estim::UnscentedKalmanFilter{…}, N::Int64, u_ry::Vector{…}, d::Vector{…}, ru::Vector{…}; plant::NonLinModel{…}, u_step::Vector{…}, u_noise::Vector{…}, y_step::Vector{…}, y_noise::Vector{…}, d_step::Vector{…}, d_noise::Vector{…}, x_noise::Vector{…}, x_0::Vector{…}, xhat_0::Nothing, lastu::Vector{…}, x̂_0::Vector{…})
    @ ModelPredictiveControl ~/.julia/dev/ModelPredictiveControl/src/plot_sim.jl:282
 [11] #sim!#157
    @ ~/.julia/dev/ModelPredictiveControl/src/plot_sim.jl:244 [inlined]
 [12] sim!
    @ ~/.julia/dev/ModelPredictiveControl/src/plot_sim.jl:236 [inlined]
 [13] top-level scope
    @ ~/.julia/dev/ModelPredictiveControl/docs/src/manual/mtk.md:139
Some type information was truncated. Use `show(err)` to see complete types.

Do you have any idea why ?

edit: Also, what's your trick to debug these kind of error? I looked inside the lines at [1], [2], and so on in the stacktrace but the code is "hidden" because of the macros and stuff.

[franckgaga]

@baggepinnen I just tested and the example work well with MTK v9.41. It stopped working with v9.42. Do you see anything v9.42 release that could have break this code ? Should I open an issue on MKT github ? For now I'll just restrict compat to 9.41 in the documentation Poject.toml.

[baggepinnen]

@AayushSabharwal did anything change in parameter handling or otherwise for observed functions recently?

[AayushSabharwal]

Yes, but it shouldn't cause that error. If you give me something I can run, I'll pinpoint what broke

[baggepinnen]

This is the example
https://juliacontrol.github.io/ModelPredictiveControl.jl/stable/manual/mtk/

[franckgaga]

Here's a copy-paste of the code in the manual:

using ModelPredictiveControl, ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t, varmap_to_vars
@mtkmodel Pendulum begin
    @parameters begin
        g = 9.8
        L = 0.4
        K = 1.2
        m = 0.3
    end
    @variables begin
        θ(t) # state
        ω(t) # state
        τ(t) # input
        y(t) # output
    end
    @equations begin
        D(θ)    ~ ω
        D(ω)    ~ -g/L*sin(θ) - K/m*ω + τ/m/L^2
        y       ~ θ * 180 / π
    end
end
@named mtk_model = Pendulum()
mtk_model = complete(mtk_model)

function generate_f_h(model, inputs, outputs)
    (_, f_ip), dvs, psym, io_sys = ModelingToolkit.generate_control_function(
        model, inputs, split=false; outputs
    )
    if any(ModelingToolkit.is_alg_equation, equations(io_sys))
        error("Systems with algebraic equations are not supported")
    end
    nu, nx, ny = length(inputs), length(dvs), length(outputs)
    vx = string.(dvs)
    p = varmap_to_vars(defaults(io_sys), psym)
    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_ = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs)
    u_nothing = fill(nothing, nu)
    function h!(y, x, _ , p)
        y .= try
            # MTK.jl supports a `u` argument in `h_` function but not this package. We set
            # `u` as a vector of nothing and `h_` function will presumably throw an
            # MethodError it this argument is used inside the function
            h_(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
    return f!, h!, p, nu, nx, ny, vx
end
inputs, outputs = [mtk_model.τ], [mtk_model.y]
f!, h!, p, nu, nx, ny, vx = generate_f_h(mtk_model, inputs, outputs)
Ts = 0.1
vu, vy = ["\$τ\$ (Nm)"], ["\$θ\$ (°)"]

model = setname!(NonLinModel(f!, h!, Ts, nu, nx, ny; p); u=vu, x=vx, y=vy)

mtk_model.K = defaults(mtk_model)[mtk_model.K] * 1.25
f_plant, h_plant, p = generate_f_h(mtk_model, inputs, outputs)
plant = setname!(NonLinModel(f_plant, h_plant, Ts, nu, nx, ny; p); u=vu, x=vx, y=vy)

α=0.01; σQ=[0.1, 1.0]; σR=[5.0]; nint_u=[1]; σQint_u=[0.1]
estim = UnscentedKalmanFilter(model; α, σQ, σR, nint_u, σQint_u)
Hp, Hc, Mwt, Nwt = 20, 2, [0.5], [2.5]
nmpc = NonLinMPC(estim; Hp, Hc, Mwt, Nwt, Cwt=Inf)
umin, umax = [-1.5], [+1.5]
nmpc = setconstraint!(nmpc; umin, umax)

using Plots
N = 35
res_ry = sim!(nmpc, N, [180.0], plant=plant, x_0=[0, 0], x̂_0=[0, 0, 0])
plot(res_ry)

[AayushSabharwal]

The observed function expects an MTKParameters, but the code uses a parameter vector

[baggepinnen]

That looks like a bug in the observed function building then, the system io_sys was simplified using split = false?

[AayushSabharwal]

Right, I didn't notice the split = false. Weird.

[AayushSabharwal]

This is because prior to calling generate_control_function, there's this:

mtk_model = complete(mtk_model)

Which doesn't have split = false. complete does get called with split = false inside generate_control_function, but that just prevents it from doing @set! sys.index_cache = IndexCache(sys). It doesn't remove the IndexCache already created.

[AayushSabharwal]

I can make complete(sys; split = false) always remove the IndexCache, but it is worth noting that these two will not be identical:

sys1 = ODESystem(...)
sys2 = deepcopy(sys1)
sys1 = complete(sys1; split = false)
sys2 = complete(complete(sys2); split = false)

This is because complete with split = true reorders parameters.

[baggepinnen]

I can make complete(sys; split = false) always remove the IndexCache

👍

but it is worth noting that these two will not be identical:

Is that likely to be a concern?

[AayushSabharwal]

Well a common workflow is parameters(sys) .=> values, which will break. Apart from that, as long as the only thing using the parameter order are the generated functions it shouldn't matter.

[baggepinnen]

we take care in the example above to use the parameter order defined and returned by generate_control_function, psym. I hope that will be the final and correct order