Description
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)
endwhere 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)
endand 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)
endThis 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