Better integration with ModellingToolkit

[ufechner7]

baggepinnen wrote on Discourse what you need to do if you want to use an MTK model with this package:

• Build a dynamics function that includes inputs, docs here.
• Write a wrapper that changes the order of the input arguments to match the order assumed by ModelPredictiveControl.jl
• Possibly build an output function, docs on the same page
• Manually keep track of the order of state variables in the generated function

As a first step it would be nice to have a simple, working example for this.

[baggepinnen]

Please consider working through and contributing this example @ufechner7

[franckgaga]

FYI, I'm working on some (final) breaking change for the v1.0.0 release that should simplify these steps. I will change the order of the input arguments to match the order of MTK and I will officially support a time $t$ and a parameter (instead of a closure, like in the manual) arguments.

[baggepinnen]

A problem I see with including an MTK example is that it is rather involved to cover all corner cases that MTK exposes. Committing to an interface to MTK is a substantial amount of nontrivial work, and we do not want to give users false indications of a well supported interface.

[franckgaga]

I don't plan on making an official interface to MTK, but a section in the manual that exposes the steps for that would be clearly useful (with a disclaimer that it's not an official interface). I personally don't have much experience with MTK so it would not be easy task for me (and my little experience with MTK was rather bad in term of clear error messages)

Anyway, there is already JuliaSimControl that is tightly integrated with MTK by design. No need to overlap on this specific feature.

[1-Bart-1]

I can try making a working example. And yes, JuliaSimControl is tightly integrated with MTK, but it is not open source. So I think more people would be interested in documentation for using MTK with ModelPredictiveControl.jl.

[1-Bart-1]

@franckgaga do you know how soon you will make those breaking changes?

[baggepinnen]

You can use this example as a guide https://baggepinnen.github.io/ControlSystemIdentification.jl/dev/nonlinear/#Identifying-parameters-in-a-ModelingToolkit-model

[franckgaga]

@franckgaga do you know how soon you will make those breaking changes?

I will work on these feature this week and in the next week. I would say in max 2 weeks these changes will be released.

[1-Bart-1]

The following example works until the last 4 lines. In these last lines, ForwardDiff tries to run f! with ForwardDiff.Dual numbers, but this does not work with the generated f! function. Any ideas on how to fix this?

using ModelPredictiveControl
using ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t

@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

function generate_f_h(model, par)
    @named mtkmodel = model()
    mtkmodel = complete(mtkmodel)
    inputs = [mtkmodel.τ]
    outputs = [mtkmodel.y]
    (f_, f_ip), dvs, psym, io_sys = ModelingToolkit.generate_control_function(mtkmodel, inputs, split=false; outputs=outputs)
    h_ = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs = inputs)
    function f!(dx, x, u, _)
        @show dx
        dx .= f_(x, u, par, 1)
    end
    function h!(y, x, _)
        y .= h_(x, 1, par, 1)
    end
    return f!, h!
end

const par = (9.8, 0.4, 1.2, 0.3)
f!, h! = generate_f_h(Pendulum, par)
nu, nx, ny, Ts = 1, 2, 1, 0.1
vu, vx, vy = ["\$τ\$ (Nm)"], ["\$θ\$ (rad)", "\$ω\$ (rad/s)"], ["\$θ\$ (°)"]
model = setname!(NonLinModel(f!, h!, Ts, nu, nx, ny); u=vu, x=vx, y=vy)

using Plots
u = [0.5]
N = 35
res = sim!(model, N, u)
display(plot(res, plotu=false))

α=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)

const par_plant = (par[1], par[2], 1.25*par[3], par[4])
f_plant, h_plant = generate_f_h(Pendulum, par_plant)
plant = setname!(NonLinModel(f_plant, h_plant, Ts, nu, nx, ny); u=vu, x=vx, y=vy)
@time res = sim!(estim, N, [0.5], plant=plant, y_noise=[0.5])
display(plot(res, plotu=false, plotxwithx̂=true))

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)

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

res_yd = sim!(nmpc, N, [180.0], plant=plant, x_0=[π, 0], x̂_0=[π, 0, 0], y_step=[10])
display(plot(res_yd))

[baggepinnen]

In

    function f!(dx, x, u, _)
        @show dx
        dx .= f_(x, u, par, 1)
    end

you could use f_ip instead, this is "in place" and accepts dx as the first argument. What type does dx have?

This works

function generate_f_h(model, par)
    @named mtkmodel = model()
    mtkmodel = complete(mtkmodel)
    inputs = [mtkmodel.τ]
    outputs = [mtkmodel.y]
    (f_, f_ip), dvs, psym, io_sys = ModelingToolkit.generate_control_function(mtkmodel, inputs, split=false; outputs=outputs)
    any(ModelingToolkit.is_alg_equation, equations(io_sys)) && error("Systems with algebraic equations are not supported")
    h_ = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs = inputs)
    function f!(dx, x, u, _)
        f_ip(dx, x, u, par, 1)
        dx
    end
    function h!(y, x, _)
        y .= h_(x, 1, par, 1)
    end
    return f!, h!
end

[baggepinnen]

Another note, MTK does not guarantee that the state realization remains constant between versions (even patch versions are allowed to change this), the way you assign vx is thus fragile and subject to breakage at any time. Use something like vx = string.(dvs) instead. nx = 2 is also not a valid assumption, nx = length(dvs) is more robust. For safety, you should also check whether or not the simplified system contains algebraic equations, if it does, an error should be thrown since the RK4 integrator used by NonLinModel does not support algebraic equations. If someone tries to use this example on a system with algebraic equations they will be silently treated as differential equations and you get a very confusing situation :) You check this by

any(ModelingToolkit.is_alg_equation, equations(io_sys)) && error("Systems with algebraic equations are not supported")

[1-Bart-1]

The updated example with your comments works well:

using ModelPredictiveControl
using ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t

@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, par, inputs, outputs)
    (_, f_ip), dvs, _, io_sys = ModelingToolkit.generate_control_function(model, inputs, split=false; outputs=outputs)
    any(ModelingToolkit.is_alg_equation, equations(io_sys)) && error("Systems with algebraic equations are not supported")
    h_ = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs = inputs)
    nx = length(dvs)
    vx = string.(dvs)
    function f!(dx, x, u, _)
        f_ip(dx, x, u, par, 1)
        nothing
    end
    function h!(y, x, _)
        y .= h_(x, 1, par, 1)
        nothing
    end
    return f!, h!, nx, vx
end

inputs, outputs = [mtk_model.τ], [mtk_model.y]

const par = (9.8, 0.4, 1.2, 0.3)
f!, h!, nx, vx = generate_f_h(mtk_model, par, inputs, outputs)
nu, ny, Ts = 1, 1, 0.1
vu, vy = ["\$τ\$ (Nm)"], ["\$θ\$ (°)"]
model = setname!(NonLinModel(f!, h!, Ts, nu, nx, ny); u=vu, x=vx, y=vy)

using Plots
u = [0.5]
N = 35
res = sim!(model, N, u)
display(plot(res, plotu=false))

α=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)

const par_plant = (par[1], par[2], 1.25*par[3], par[4])
f_plant, h_plant, _, _ = generate_f_h(mtk_model, par_plant, inputs, outputs)
plant = setname!(NonLinModel(f_plant, h_plant, Ts, nu, nx, ny); u=vu, x=vx, y=vy)
res = sim!(estim, N, [0.5], plant=plant, y_noise=[0.5])
display(plot(res, plotu=false, plotxwithx̂=true))

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)

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

res_yd = sim!(nmpc, N, [180.0], plant=plant, x_0=[π, 0], x̂_0=[π, 0, 0], y_step=[10])
display(plot(res_yd))

[baggepinnen]

Where does this particular order of parameters come from?

const par = (9.8, 0.4, 1.2, 0.3)

it is most likely incorrect, the generated dynamics have parameters in this order

julia> psym
4-element Vector{Num}:
 L
 g
 K
 m

This is a safer way to get the parameters in the correct order

julia> getindex.(Ref(defaults(io_sys)), psym)
4-element Vector{Float64}:
 0.4
 9.8
 1.2
 0.3

It is generally unsafe to assume the order of elements of any variable in MTK. The one exception is when the order is explicitly indicated by the user, such as the order of inputs and outputs to generate_control_function and build_explicit_observed_function. This page contains some utilities and explanations as to why
https://docs.sciml.ai/ModelingToolkit/stable/basics/FAQ/

The final system that has all the required information about unknowns and variable order is io_sys, and the unknowns dvs and parameters psym are returned from generate_control_function precisely because the layout of those variables will be required for correct usage of the generated function.

[1-Bart-1]

But in which order does f_ip want its parameters? In the Ref(defaults(io_sys)) order or in the psym order? And somehow I get a different order than you:
getindex.(Ref(defaults(io_sys)), psym) = [9.8, 0.4, 1.2, 0.3]
Ref(defaults(io_sys)) = Base.RefValue{Dict{Any, Any}}(Dict{Any, Any}(L => 0.4, g => 9.8, K => 1.2, m => 0.3))

[baggepinnen]

defaults(io_sys) is a Dict so it does not have an order defined.

And somehow I get a different order than you

This is exactly why it is tricky to get things like this right with MTK ;)

But in which order does f_ip want its parameters? In the Ref(defaults(io_sys)) order or in the psym order?

in the psym order

[1-Bart-1]

Should be correct now:

using ModelPredictiveControl
using 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) = 0.0 # 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=outputs)
    any(ModelingToolkit.is_alg_equation, equations(io_sys)) && error("Systems with algebraic equations are not supported")
    h_ = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs = inputs)
    nx = length(dvs)
    vx = string.(dvs)
    @show par = varmap_to_vars(defaults(io_sys), psym)
    function f!(dx, x, u, _)
        f_ip(dx, x, u, par, 1)
        nothing
    end
    function h!(y, x, _)
        y .= h_(x, 1, par, 1)
        nothing
    end
    return f!, h!, nx, vx
end

inputs, outputs = [mtk_model.τ], [mtk_model.y]

f!, h!, nx, vx = generate_f_h(mtk_model, inputs, outputs)
nu, ny, Ts = 1, 1, 0.1
vu, vy = ["\$τ\$ (Nm)"], ["\$θ\$ (°)"]
model = setname!(NonLinModel(f!, h!, Ts, nu, nx, ny); u=vu, x=vx, y=vy)

using Plots
u = [0.5]
N = 35
res = sim!(model, N, u)
display(plot(res, plotu=false))

α=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)

mtk_model.K = defaults(mtk_model)[mtk_model.K] * 1.25
f_plant, h_plant, _, _ = generate_f_h(mtk_model, inputs, outputs)
plant = setname!(NonLinModel(f_plant, h_plant, Ts, nu, nx, ny); u=vu, x=vx, y=vy)
res = sim!(estim, N, [0.5], plant=plant, y_noise=[0.5])
display(plot(res, plotu=false, plotxwithx̂=true))

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)

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

res_yd = sim!(nmpc, N, [180.0], plant=plant, x_0=[π, 0], x̂_0=[π, 0, 0], y_step=[10])
display(plot(res_yd))

[1-Bart-1]

generate_control_function works fine on the small model in this example, but when trying it on a bigger model (1757 equations) it gets extremely slow. Any idea why and what can I do about it?

[baggepinnen]

Is it slower than structural_simplify on the same model? If so, could you use @profview (if in vscode) on the call to see what takes the most time?

[1-Bart-1]

The problem is that the call never completes. I will try to use @profview and let it run overnight if you want. When I interrupt the function with ctrl+c it is stuck in line 223 in generate_control_function: eqs = [eq for eq in full_equations(sys)]. And yes, it is a lot slower than structural_simplify.

[baggepinnen]

1757 equations isn't terribly big, but the MTK compiler is known to scale quite poorly. If you're an academic user, JuliaSimCompiler could potentially be an option to make it significantly faster. Unfortunately, generate_control_function for JSCompiler is still not merged (up for review) so you can't try it right away. If your model is public, I could potentially try it for you, I have it implemented locally.

[1-Bart-1]

The model is public, it is part of KiteModels.jl. It is not part of the main branch yet, but you can add the branch:

julia> ] add https://github.com/ufechner7/KiteModels.jl/tree/feat/xfoil-polars

And the code I am trying to run:
src/kite.jl
https://github.com/Albatross-Kite-Transport/KitePredictiveControl

But don't worry if it is too much work to run the model, I can also wait for the function to be merged.

[baggepinnen]

It was just merged, but we need a new release as well for it to become available. I tried loading your model but couldn't get past kite::KPS4_3L = KPS4_3L(KCU(se("system_3l.yaml"))), despite giving the absolute path to the file.

Once a new release of JSCompiler is made, you can use it like so

using JuliaSimCompiler
(_, f_ip), dvs, psym, io_sys = ModelingToolkit.generate_control_function(IRSystem(model), inputs, split=false)

[ufechner7]

@baggepinnen You need to call:

set_data_path(joinpath(pwd(), "data"))

which is commented in the code. This command assumes you execute the script kite.jl with

include("src/kite.jl")

Giving the absolute path will not work.

[baggepinnen]

I did try that, but I think these functions assume something about the local paths that does not hold the way I installed things. Anyway, you can already now try structural_simplify from JSCompiler by calling

using JuliaSimCompiler
structural_simplify(IRSystem(model), (inputs, []))

this is the internal processing that takes place inside of generate_control_function

[1-Bart-1]

Here is the result of the profile on the current generate_control_function. The left part is structural_simplify, while the right part is symbolics_tearing. So around 60% of the runtime is structural_simplify.

Somehow, when running stuctural_simplify inside KiteModels.jl, it is way faster:

@show length(equations(sys))
@time (s.simple_sys, _) = structural_simplify(sys, (inputs, []); simplify=true)


length(equations(sys)) = 1757
 44.193423 seconds (59.91 M allocations: 3.580 GiB, 2.88% gc time, 85.55% compilation time: 7% of which was recompilation)

[baggepinnen]

The simplify=true part might be adding a lot of time, can you try without that? It typically makes rather small difference for large models. If that does not speed things up sufficiently, JSCompiler is probably the best bet.

[1-Bart-1]

kite_model, inputs = KiteModels.model!(kite, pos, vel)

inputs, outputs = inputs, [kite_model.pos[2, kite.num_A]]
@time structural_simplify(kite_model, (inputs, []))

results:

5.424258 seconds (12.56 M allocations: 717.108 MiB, 2.76% gc time)

Somehow, structural_simplify is way slower inside generate_control_function.

@time (_, f_ip), dvs, psym, io_sys = ModelingToolkit.generate_control_function(model, inputs, split=false; outputs=outputs)
8306.809924 seconds (5.55 G allocations: 260.320 GiB, 80.75% gc time)

How is this possible? Structural_simplify is being called with the exact same arguments and model inside generate_control_function.

[1-Bart-1]

@baggepinnen I submitted an issue for this here: SciML/ModelingToolkit.jl#3077

[1-Bart-1]

My solution for now: splitting up into two functions and using Serialization

function get_control_function(model, inputs; filename="control_function.bin")
    if isfile(filename)
        println("deserializing control function")
        return deserialize(filename)
    else
        println("generating control function: sit back, relax and enjoy...")
        @time (_, f_ip), dvs, psym, io_sys = ModelingToolkit.generate_control_function(model, inputs; split=false)
        any(ModelingToolkit.is_alg_equation, equations(io_sys)) && error("Systems with algebraic equations are not supported")
        println("serializing control function")
        @time serialize(filename, (f_ip, dvs, psym, io_sys))
        return (f_ip, dvs, psym, io_sys)
    end
end

function generate_f_h(inputs, outputs, f_ip, dvs, psym, io_sys)
    @time h_ = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs = inputs)
    nx = length(dvs)
    vx = string.(dvs)
    @show par = varmap_to_vars(defaults(io_sys), psym)
    function f!(dx, x, u, _)
        f_ip(dx, x, u, par, 1)
        nothing
    end
    function h!(y, x, _)
        y .= h_(x, 1, par, 1)
        nothing
    end
    return f!, h!, nx, vx
end

[franckgaga]

@1-Bart-1 @baggepinnen so what would be the most reliable chunk of code to add in the manual as an exemple of MTK integration on a simple nonlinear model ?

P.S. I almost finished the breaking changes for v1.0.0. The NonLinModel function signature won't be identical to MTK since I want to keep a clear distinction between the manipulated inputs $$\mathbf{u}$$ and the measured disturbances $\mathbf{d}$. They are truly different in practice, $$\mathbf{u}$$ will typically include a DAC with zero-order hold, and $\mathbf{d}$ an ADC. Since there is no measured disturbances in MTK, it's impossible to perfectly match the signature.

The breaking change will add a model parameter (any type) in the signature:

f(x, u, d, p) -> xdot   # or f!(xdot, x, u, d, p) -> nothing
h(x, d, p) -> y         # or h!(y, x, d, p) -> nothing

A similar breaking change is introduced in the custom economic function $J_E$ of NonLinMPC objects. I chose to exclude the time $t$ in the signature since a direct dependence is not common for control systems and introducing a new measured disturbance defined as $d(t)=t$ is an easy workaround.

[1-Bart-1]

This one should be the most reliable / simple. Uses the same pendulum problem as the nonlinear example in the documentation.

Should be correct now:

using ModelPredictiveControl
using 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) = 0.0 # 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=outputs)
    any(ModelingToolkit.is_alg_equation, equations(io_sys)) && error("Systems with algebraic equations are not supported")
    h_ = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs = inputs)
    nx = length(dvs)
    vx = string.(dvs)
    par = varmap_to_vars(defaults(io_sys), psym)
    function f!(dx, x, u, _, _)
        f_ip(dx, x, u, par, 1)
        nothing
    end
    function h!(y, x, _, _)
        y .= h_(x, 1, par, 1)
        nothing
    end
    return f!, h!, nx, vx
end

inputs, outputs = [mtk_model.τ], [mtk_model.y]

f!, h!, nx, vx = generate_f_h(mtk_model, inputs, outputs)
nu, ny, Ts = 1, 1, 0.1
vu, vy = ["\$τ\$ (Nm)"], ["\$θ\$ (°)"]
model = setname!(NonLinModel(f!, h!, Ts, nu, nx, ny); u=vu, x=vx, y=vy)

using Plots
u = [0.5]
N = 35
res = sim!(model, N, u)
display(plot(res, plotu=false))

α=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)

mtk_model.K = defaults(mtk_model)[mtk_model.K] * 1.25
f_plant, h_plant, _, _ = generate_f_h(mtk_model, inputs, outputs)
plant = setname!(NonLinModel(f_plant, h_plant, Ts, nu, nx, ny); u=vu, x=vx, y=vy)
res = sim!(estim, N, [0.5], plant=plant, y_noise=[0.5])
display(plot(res, plotu=false, plotxwithx̂=true))

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)

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

res_yd = sim!(nmpc, N, [180.0], plant=plant, x_0=[π, 0], x̂_0=[π, 0, 0], y_step=[10])
display(plot(res_yd))

[franckgaga]

Thanks a lot @1-Bart-1 !

[1-Bart-1]

No problem :) I updated the code above to reflect the breaking change (just added an extra underscore to the arguments in f and h).

[1-Bart-1]

I found a mistake in the example code...

      function h!(y, x, _, _)
          y .= h_(x, 1, par, 1)
          nothing
      end

This should be changed to account for different input lengths:

    function h!(y, x, _, _)
        y .= h_(x, ones(Int, length(inputs)), par, 1)
        nothing
    end

And do I understand correctly that the output/observed variable cannot be directly dependent on the input? Because the input variables are not passed as an argument to h!.

[baggepinnen]

What are these ones doing there in the first place? The second argument to h from MTK is supposed to be u

[1-Bart-1]

The ModelPredictiveControl.jl h! function does not have a u argument. As I understand it, MPC.jl assumes the output only depends on state variables, and not on input variables. This should hold true for most control functions. In the pendulum example, the output y only depends on the state x, and not on input u. It does not matter which numbers you pass to the second argument, because y does not depend on u.

[baggepinnen]

Ok, in that case it would probably be good to mention this in the example and also possibly error if the output does depend on inputs. Maybe passing nothing instead of ones would be more clear, and give the desired error if the output would depend on inputs

[1-Bart-1]

How can you check if outputs depend on inputs? And @franckgaga are there any situations where this would happen? If so, shouldn't h! and h have an argument u?

[baggepinnen]

High-pass systems and systems with acceleration as output often have outputs that depend on inputs

[1-Bart-1]

But does it help to make an error message if MPC.jl doesn't support outputs that depend on inputs anyways?

On a different note: par = ModelingToolkit.varmap_to_vars(defaults(sys), psym; check=false) doesn't seem to work with the JuliaSimCompiler version of generate_control_function.

julia> par = ModelingToolkit.varmap_to_vars(defaults(io_sys), psym; check=false)
ERROR: KeyError: key g not found
julia> typeof(io_sys)
JuliaSimCompiler.ScheduledSystem{Nothing, JuliaSimCompiler.DiscreteInfo{Nothing, Nothing}}

julia> typeof(psym)
Vector{IRElement} (alias for Array{JuliaSimCompiler.ADT.IRElement, 1})

julia> defaults(io_sys)
Dict{Any, Any} with 8 entries:
  m    => 0.3
  L    => 0.4
  θ(t) => 0.0
  τ(t) => 0.0
  g    => 9.8
  K    => 1.2
  ω(t) => 0.0
  y(t) => 0.0

julia> psym
4-element Vector{JuliaSimCompiler.ADT.IRElement}:
 g
 L
 m
 K
 
 julia> defaults(io_sys)[psym[1]]
ERROR: KeyError: key g not found

[baggepinnen]

But does it help to make an error message if MPC.jl doesn't support outputs that depend on inputs anyways?

Yes, because the MTK model has no issues with having outputs depending on inputs. If the user uses this script with such a model, they will now silently get the wrong answer where an arrays of ones is used as input.

For the varmap issue, use

defaults(mtk_model) 

Instead

[1-Bart-1]

I get the same error when using defaults(mtk_model) and defaults(ir_model).

[franckgaga]

How can you check if outputs depend on inputs? And @franckgaga are there any situations where this would happen? If so, shouldn't h! and h have an argument u?

No, this package does not support model outputs $\mathbf{y}$ that are a function of the manipulated inputs $\mathbf{u}$ (or, a direct transmission from the inputs). I'm not keen to add this feature because:

1. This is a very unlikely assumption in practice, from a control perspective. It means that you could instantaneously affect a measured output by changing a manipulated input. We can expect that an input change will impact the measured outputs at least the next period, not right now. This is even weirder when we remember that all discrete controllers (1) samples a measurement (2) computes a new action and (3) applies this action on the real plant. A direct transmission would violate causality (the input would affect the outputs even before computing it).
2. There is practically always a zero-order hold downstream of the controller. It can be explicit (in the DAC) or implicit (the computed action is kept constant between the current discrete time and the next one).
3. A direct transmission from $\mathbf{u}$ can be easily removed by adding a new state. The new state can be a pure unit delay or fast first-order dynamics, and the transformed plant model will be presumably more realistic.

High-pass systems and systems with acceleration as output often have outputs that depend on inputs

If we do not neglect the aspects above in the modelling, the discretized version of these systems will have zero direct transmission from the inputs $\mathbf{u}$

[franckgaga]

But does it help to make an error message if MPC.jl doesn't support outputs that depend on inputs anyways?

Yes, because the MTK model has no issues with having outputs depending on inputs. If the user uses this script with such a model, they will now silently get the wrong answer where an arrays of ones is used as input.

From that perspective, sending nothing to h_ function seems to be the best solution. That is, an error will be thrown if this argument is used inside h_, else everything will run as normal. Do you agree ?

[baggepinnen]

You are of course correct in what you say, but I still wouldn't be so quick in dismissing the direct feedthrough use case. You are not considering the following

Not all systems that are controlled by mpc are purely physical systems. A very common case that isn't is when an mpc controller is used in a slower outer loop, and a much faster inner controller is tracking references provided by mpc. In this case, you cannot form any reasonable model of the inner closed loop without direct feedthrough. The inner loop in this case will interpolate between the samples in the reference trajectory in a non causal way, and the inner loop is in this case best discretized using a Tustin or FOH discretization which introduces nonzero D matrices.

You can of course add a one sample delay, but if you're using a cascade controller because it's too slow to run mpc at the faster rate, an additional sample delay at the slow rate is a very poor solution to an artificial problem introduced by the MPC tool.

[1-Bart-1]

I get the same error when using defaults(mtk_model) and defaults(ir_model).

Any suggestions for how to solve this @baggepinnen ?

[baggepinnen]

I can check carefully tomorrow, but for now, check if parameters(mtk_model) are the same as psym and if they are, use the former. I have to double check the validity of this, even of they happen to agree in this case.

[1-Bart-1]

They are not the same.

julia> parameters(mtk_model)
4-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 L
 g
 K
 m

julia> psym
4-element Vector{JuliaSimCompiler.ADT.IRElement}:
 g
 L
 m
 K
 
 julia> parameters(ir_model)
4-element Vector{JuliaSimCompiler.ADT.IRElement}:
 L
 g
 K
 m
 
 julia> parameters(io_sys)
5-element Vector{JuliaSimCompiler.ADT.IRElement}:
 g
 L
 m
 K
 τ

[baggepinnen]

In that case, you can maybe use JSC.initial_conditions (or something like that, I'm not in office to check) to get the numeric parameter vector.

[1-Bart-1]

Could you check as well how to run build_explicit_observed_function on an IRSystem?

[baggepinnen]

Could you check as well how to run build_explicit_observed_function on an IRSystem?

The same as for an ODESystem? Possibly with keyword arg target = JSC.JuliaTarget().

[1-Bart-1]

Thanks, this is the updated code using IRSystem JSC:

function get_control_function(model, inputs)
    println("generating control function")
    f_ip, dvs, psym, io_sys = ModelingToolkit.generate_control_function(IRSystem(model), inputs)
    # any(ModelingToolkit.is_alg_equation, equations(io_sys)) && error("Systems with algebraic equations are not supported")
    return (f_ip, dvs, psym, io_sys)
end

function generate_f_h(sys, inputs, outputs, f_ip, dvs, psym)
    h_ = JuliaSimCompiler.build_explicit_observed_function(sys, outputs; inputs = inputs, target = JuliaSimCompiler.JuliaTarget())
    nu = length(inputs)
    ny = length(outputs)
    nx = length(dvs)
    vu = string.(inputs)
    vy = string.(outputs)
    vx = string.(dvs)
    par = JuliaSimCompiler.initial_conditions(io_sys, defaults(io_sys), psym)[2]
    function f!(dx, x, u, _, _)
        f_ip(dx, x, u, par, 1)
        nothing
    end
    function h!(y, x, _, _)
        y .= h_(x, fill(nothing, length(inputs)), par, 1)
        nothing
    end
    return f!, (h!, nu, ny, nx, vu, vy, vx)
end

One last comment on the example code: the order of the states is determined by io_sys. So to set an initial state, it is better to use something like this:

x_0 = zeros(nx)
x̂_0 = zeros(nx + ny)
x_0[ModelingToolkit.variable_index(io_sys, :θ)] = π
x̂_0[ModelingToolkit.variable_index(io_sys, :θ)] = π
res_yd = sim!(nmpc, N, [180.0], plant=plant, x_0=x_0, x̂_0=x̂_0, y_step=[10])

[franckgaga]

Thanks, this is the updated code using IRSystem JSC:

function get_control_function(model, inputs)
    println("generating control function")
    f_ip, dvs, psym, io_sys = ModelingToolkit.generate_control_function(IRSystem(model), inputs)
    # any(ModelingToolkit.is_alg_equation, equations(io_sys)) && error("Systems with algebraic equations are not supported")
    return (f_ip, dvs, psym, io_sys)
end

function generate_f_h(sys, inputs, outputs, f_ip, dvs, psym)
    h_ = JuliaSimCompiler.build_explicit_observed_function(sys, outputs; inputs = inputs, target = JuliaSimCompiler.JuliaTarget())
    nu = length(inputs)
    ny = length(outputs)
    nx = length(dvs)
    vu = string.(inputs)
    vy = string.(outputs)
    vx = string.(dvs)
    par = JuliaSimCompiler.initial_conditions(io_sys, defaults(io_sys), psym)[2]
    function f!(dx, x, u, _, _)
        f_ip(dx, x, u, par, 1)
        nothing
    end
    function h!(y, x, _, _)
        y .= h_(x, fill(nothing, length(inputs)), par, 1)
        nothing
    end
    return f!, (h!, nu, ny, nx, vu, vy, vx)
end

One last comment on the example code: the order of the states is determined by io_sys. So to set an initial state, it is better to use something like this:

x_0 = zeros(nx)
x̂_0 = zeros(nx + ny)
x_0[ModelingToolkit.variable_index(io_sys, :θ)] = π
x̂_0[ModelingToolkit.variable_index(io_sys, :θ)] = π
res_yd = sim!(nmpc, N, [180.0], plant=plant, x_0=x_0, x̂_0=x̂_0, y_step=[10])

I'm a little bit lost here, since I don't know what is JuliaSimCompiler (is it a proprietary tool?). Should I need to update the example in the manual with this code ?

[baggepinnen]

Yeah JSC is a proprietary tool, no need to update the example with JSC, we discussed it here only since JSC is better able to handle very large systems that @1-Bart-1 tried to use this functionality on.