Coupling MTK with time discrete control

[ValentinKaisermayer]

Let's say I want to simulate an MPC controller where the plant is an MTK system. I would need to access the current measurements and would have to modify control inputs via a callback in every time step according to some optimization problem.

For example a modified RC circuit:

using ModelingToolkit, Plots, DifferentialEquations

@parameters t
@connector function Pin(;name)
    sts = @variables v(t) = 1.0 i(t) = 1.0
    ODESystem(Equation[], t, sts, []; name=name)
end

function ModelingToolkit.connect(::Type{Pin}, ps...)
    eqs = [
           0 ~ sum(p -> p.i, ps) # KCL
          ]
    # KVL
    for i in 1:length(ps) - 1
        push!(eqs, ps[i].v ~ ps[i + 1].v)
    end

    return eqs
end

function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name), g)
end

function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t) = 1.0 i(t) = 1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name), p, n)
end

function Resistor(;name, R=1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R = R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name), oneport)
end

function Capacitor(;name, C=1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C = C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name), oneport)
end

function ControlledVoltageSource(;name, V=1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    vars = @variables V(t) = V
    D = Differential(t)
    eqs = [
           V ~ v
           D(V) ~ 0 # is constant 
          ]
    extend(ODESystem(eqs, t, vars, []; name), oneport)
end

@named resistor = Resistor(R=1)
@named capacitor = Capacitor(C=1)
@named source = ControlledVoltageSource(V=1)
@named ground = Ground()

rc_eqs = [
    connect(source.p, resistor.p)
    connect(resistor.n, capacitor.p)
    connect(capacitor.n, source.n, ground.g)
]

@named rc_model = compose(ODESystem(rc_eqs, t), [resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
    capacitor.v => 0.0
]
prob = ODAEProblem(sys, u0, (0, 10.0))

@inline indexof(sym, syms) = findfirst(isequal(sym), syms)
function cb!(int) # the callback function
    int.u[indexof(source.V, states(sys))] = sin(int.t) # this is where the actual MPC problem would be solved
    return
end

sol = solve(prob, Tsit5(); callback=PeriodicCallback(cb!, 1.))
plot(sol, vars=[capacitor.v, resistor.v, source.v])

Notice the ControlledVoltageSource implementation:

function ControlledVoltageSource(;name, V=1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    vars = @variables V(t) = V
    D = Differential(t)
    eqs = [
           V ~ v
           D(V) ~ 0 # is constant 
          ]
    extend(ODESystem(eqs, t, vars, []; name), oneport)
end

where I had to make V a state in order to be able to access its value in the solution. For this to work the equation D(V) ~ 0 has to be added. This is actually the solution from: https://discourse.julialang.org/t/update-parameters-in-modelingtoolkit-using-callbacks/63770 An unexperienced user might not be able to do this.

The control input can not be a parameter or else I can not access it correctly in the solution but at the same time the controls kwarg requires it to be one. For a linear MPC I would need the system matrix A, input matrix B and output matrix C (and very rarely feedthrough matrix D), i.e. the jacobi matrices with respect to the states, inputs and outputs of the sate and output equations at some operating point. MTK should be able to provide those.

If I where to implement a time discrete PID controller or any other time discrete control strategy it would be pretty much the same.

A more Modelicaesk solution would be to define a connector for the inputs and outputs like:

@connector function RealInput(;name, u0=0)
    sts = @variables u(t) = u0
    ODESystem([Differential(t)(u) ~ 0], t, sts, []; name)
end

and use it in the component like:

function ControlledVoltageSource(;name, V=1.0)
    @named oneport = OnePort()
    @named input = RealInput(u0=V)
    @unpack v = oneport
    eqs = [
        v ~ input.u
    ]
    extend(compose(ODESystem(eqs, t, [], []; name), input), oneport)
end

This would at least solve the issue of the additional equation, that one has to add for this to work. However, maybe I am doing something wrong but I have to access the variable in a funny way in the callback:

function cb!(int) # the callback function
    int.u[indexof(source.input₊u, states(sys))] = sin(int.t)
    return
end

[ChrisRackauckas]

Did you try using the Difference operator? (Didn't read through all of this)

[ValentinKaisermayer]

Did you try using the Difference operator? (Didn't read through all of this)

I did not

[ChrisRackauckas]

There's not much documented on it right now (anything?) but that is the approach that was developed for this.

[ValentinKaisermayer]

I will have a look at it. I found https://github.com/SciML/ModelingToolkit.jl/blob/5a47dd69a6c401f3ba53a8d980b4ad860ebdba7f/test/odesystem.jl#L403 and https://github.com/SciML/ModelingToolkit.jl/blob/5a47dd69a6c401f3ba53a8d980b4ad860ebdba7f/test/discretesystem.jl#L14

[ValentinKaisermayer]

I have to say the doc string:

Represents a difference operator.

at https://github.com/JuliaSymbolics/Symbolics.jl/blob/c2a4e89cc919cd3a6c9e0b19300b9b8d8d03770a/src/difference.jl#L4 Is not very helpful. 😉

[ValentinKaisermayer]

What exactly does Difference(t;dt=0.01)(x) ~ y mean? Is it x[t+0.01] = y[t]?

[ChrisRackauckas]

x[t+0.01]-x[t] = y[t]

[ValentinKaisermayer]

So i would have to write Difference(t;dt=0.01)(x) ~ f(t) - x if I want to set x to the value of f(t) at some discrete time steps? In other words, what is the intended way of implementing a first-order hold (FOH)?

[ValentinKaisermayer]

How should I implement the Step block such that the discrete time event is handled correctly? My approach would be to put the RHS (i.e. offset + IfElse.ifelse(t < start_time, 0, height)) in a callback and set the value of y in the integrator to it. However, this approach needs y to be a state, which I can enforce via Differential(t)(y) ~ 0. In generate_difference_cb() this is done very similar.: https://github.com/SciML/ModelingToolkit.jl/blob/6d31f4bba00a350b0d1ae1d7fff9483ae04bde39/src/systems/diffeqs/abstractodesystem.jl#L149

using ModelingToolkit, IfElse, DifferentialEquations, Plots

@parameters t

function Plant(;name, x0=zeros(2))
    D = Differential(t)
    sts = @variables x1(t)=x0[1] x2(t)=x0[2] y(t) u(t)
    eqs= [
        D(x1) ~ x2
        D(x2) ~ -x1 - 0.5 * x2 + u
        y ~ 0.9 * x1 + x2
    ]
    ODESystem(eqs, t, sts, []; name)
end

function Step(;name, height=1.0, offset=0.0, start_time=0.0)
    sts = @variables y(t)=0.0
    eqs = [
        y ~ offset + IfElse.ifelse(t < start_time, 0, height) # needs event to work properly
    ]
    ODESystem(eqs, t, sts, []; name)
end

@named plant = Plant()
@named step_block = Step(start_time=1)
connections = [
    plant.u ~ step_block.y
]
@named sys = compose(ODESystem(connections, t), [plant, step_block])
sys = structural_simplify(model)
prob = ODAEProblem(sys, [], (0, 10.0))

sol = solve(prob, Tsit5())
plot(sol, vars=[step_block.y, plant.y])

I do not see how the Difference operator can help me in this case. Not all discrete control schemes need a difference equation. For instance a time discrete P-controller. Nevertheless, I tried but got an error:

function DiscreteTimeStep(;name, height=1.0, offset=0.0, start_time=0.0, dt=0.1)
    sts = @variables y(t)=0.0
    D = Difference(t; dt)
    eqs = [
        D(y) ~ -y + offset + IfElse.ifelse(t < start_time, 0, height)
    ]
    ODESystem(eqs, t, sts, []; name)
end

@named plant = Plant()
@named step_block = DiscreteTimeStep(start_time=1, dt=0.1)
connections = [
    plant.u ~ step_block.y
]
@named sys = compose(ODESystem(connections, t), [plant, step_block])
sys = initialize_system_structure(sys)
prob = ODAEProblem(sys, [], (0, 10.0))

ERROR: LoadError: KeyError: key 6 not found
Stacktrace:
 [1] getindex
   @ .\dict.jl:482 [inlined]
 [2] torn_system_jacobian_sparsity(sys::ODESystem)
   @ ModelingToolkit.StructuralTransformations .julia\packages\ModelingToolkit\45oYI\src\structural_transformation\codegen.jl:83
 [3] build_torn_function(sys::ODESystem; expression::Bool, jacobian_sparsity::Bool, checkbounds::Bool, kw::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit.StructuralTransformations .julia\packages\ModelingToolkit\45oYI\src\structural_transformation\codegen.jl:232
 [4] build_torn_function
   @ .julia\packages\ModelingToolkit\45oYI\src\structural_transformation\codegen.jl:189 [inlined]
 [5] ODAEProblem{true}(sys::ODESystem, u0map::Vector{Any}, tspan::Tuple{Int64, Float64}, parammap::SciMLBase.NullParameters; kw::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit.StructuralTransformations .julia\packages\ModelingToolkit\45oYI\src\structural_transformation\codegen.jl:350
 [6] ODAEProblem (repeats 2 times)
   @ .julia\packages\ModelingToolkit\45oYI\src\structural_transformation\codegen.jl:341 [inlined]
 [7] ODAEProblem(::ODESystem, ::Vararg{Any, N} where N; kw::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit.StructuralTransformations .julia\packages\ModelingToolkit\45oYI\src\structural_transformation\codegen.jl:333
 [8] ODAEProblem(::ODESystem, ::Vararg{Any, N} where N)
   @ ModelingToolkit.StructuralTransformations .julia\packages\ModelingToolkit\45oYI\src\structural_transformation\codegen.jl:333

I think the actual solution could be that IfElse.ifelse generates the callback, as you hinted in #38 and #893.

[ChrisRackauckas]

So i would have to write Difference(t;dt=0.01)(x) ~ f(t) - x if I want to set x to the value of f(t) at some discrete time steps? In other words, what is the intended way of implementing a first-order hold (FOH)?

Yes, it's a little odd but it also a bit more naturally matches the differential forms, and it simplifies anyways.

That IfElse issue doesn't look related. It looks like a structural_simplify bug. It's all not documented yet for a reason 😅, but if you're willing to keep building test cases that would be helpful.

[ValentinKaisermayer]

but if you're willing to keep building test cases that would be helpful.

Sure! 😄

[ValentinKaisermayer]

@variables t
function DiscreteSin(dt; name)
    @variables y(t)=0.0
    U = DiscreteUpdate(t; dt=dt)
    eqs = [
        U(y) ~ sin(t)
    ]
    ODESystem(eqs, t, name=name)
end
function Plant(;name, x0=zeros(2))
    D = Differential(t)
    sts = @variables x1(t)=x0[1] x2(t)=x0[2] y(t) u(t)
    eqs= [
        D(x1) ~ x2
        D(x2) ~ -x1 - 0.5 * x2 + u
        y ~ 0.9 * x1 + x2
    ]
    ODESystem(eqs, t, sts, []; name)
end

dt = 1
@named plant_block = Plant()
@named sin_block = DiscreteSin(dt)
@named model = ODESystem([plant_block.u ~ sin_block.y], systems=[sin_block, plant_block])
sys = structural_simplify(model)
prob = ODAEProblem(sys, [], (0.0, 10.0))
sol = solve(prob, Tsit5())

errors:

ERROR: LoadError: MethodError: no method matching extract_derivative(::Type{ForwardDiff.Tag{Base.Fix2{NonlinearFunction{false, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#347"), Symbol("##arg#348")), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x16a3eeb1, 0x218cb4b5, 0x9fd455e0, 0x71ae8127, 0x109ddc8c)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Tuple{Float64}}, Float64}}, ::SymbolicUtils.Mul{ForwardDiff.Dual{ForwardDiff.Tag{Base.Fix2{NonlinearFunction{false, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#347"), Symbol("##arg#348")), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x16a3eeb1, 0x218cb4b5, 0x9fd455e0, 0x71ae8127, 0x109ddc8c)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Tuple{Float64}}, Float64}, Float64, 1}, Int64, Dict{Any, Number}, Nothing})

(test) pkg> st
      Status `****\Project.toml`
  [0c46a032] DifferentialEquations v6.20.0
  [961ee093] ModelingToolkit v7.1.3
  [91a5bcdd] Plots v1.23.6

[ValentinKaisermayer]

I hacked this simple control loop together, maybe it is of help for someone.

using ModelingToolkit, DifferentialEquations, Plots

@variables t
D = Differential(t)

function Plant(;name)
    sts = @variables x(t)=0 y(t)=0 u(t)=0
    eqs= [
        D(x) ~ u
        y ~ x
    ]
    ODESystem(eqs, t, sts, []; name)
end

function PController(;name, kp=1)
    sts = @variables e(t)=0 u(t)=0
    eqs= [
        u ~ kp * e 
    ]
    ODESystem(eqs, t, sts, []; name)
end

function Feedback(;name)
    sts = @variables y(t)=0 r(t)=0 e(t)=0
    eqs= [
        e ~ r - y
    ]
    ODESystem(eqs, t, sts, []; name)
end

The time continuous model is easy:

@named plant = Plant()
@named controller = PController(kp=1)
@named feedback = Feedback()
@named model = ODESystem([plant.u ~ controller.u, feedback.y ~ plant.y, controller.e ~ feedback.e, feedback.r ~ 1], systems=[plant, controller, feedback])
sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0.0, 10.0))
sol = solve(prob, Tsit5())

plot(sol, vars=[plant.u, plant.y])

For the time discrete part I had to do some tricks. I tried the DiscreteUpdate operator but it seems to be not fully fleshed out.

function DiscretePController(;name, kp=1, dt=1)
    sts = @variables u(t)=0 e(t)=0
    pars = @parameters kp=kp
    U = DiscreteUpdate(t; dt)
    eqs= [
        U(u) ~ kp * e
    ]
    ODESystem(eqs, t, sts, pars; name)
end

@named plant = Plant()
@named controller = DiscretePController(kp=1, dt=0.25)
@named feedback = Feedback()
@named model = ODESystem([plant.u ~ controller.u, feedback.y ~ plant.y, controller.e ~ feedback.e, feedback.r ~ 1], systems=[plant, controller, feedback])
sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0.0, 10.0))
sol = solve(prob; callback=prob.kwargs[:callback])

ERROR: LoadError: type NamedTuple has no field callback

So I came up with this approach:

function DiscretePController(;name, kp=1, dt=1)
    sts = @variables u(t)=0 e(t)=0
    pars = @parameters kp=kp
    eqs= [
        D(u) ~ 0 # force u to be a state, apply control law in callback
    ]
    ODESystem(eqs, t, sts, pars; name)
end

@named plant = Plant()
@named controller = DiscretePController(kp=1)
@named feedback = Feedback()
@named model = ODESystem([plant.u ~ controller.u, feedback.y ~ plant.y, controller.e ~ feedback.e, feedback.r ~ 1], systems=[plant, controller, feedback])
sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0.0, 10.0)) 

indexof(sym,syms) = findfirst(isequal(sym),syms)
function cb!(int)
    kp = 1
    int.u[indexof(controller.u, states(sys))] = kp * prob.f.observed(feedback.e, int.u, int.p, int.t)
end

sol = solve(prob, Tsit5(); callback=PeriodicCallback(cb!, 0.25; initial_affect=true))

plot(sol, vars=[plant.u, plant.y])

The D(u) ~ 0 forces the controller output to be a state and be constant. In the callback only states can be accessed. But via the observable function from the problem, I can access the control error e using the values provided by the integrator.

I know this is a hack but maybe it is useful for someone.

[ChrisRackauckas]

See @baggepinnen 's many issues on this. The old version was deprecated