Update parameters in ModelingToolkit using callbacks

[jairorua]

Hello,

I was trying to create a simple control loop with modelingtoolkit. The idea was to have a generic controller whose setpoint can be changed using callbacks during a dynamic simulation. For example:

using ModelingToolkit, DifferentialEquations, Plots


# Define first order system
function firstOrder(;name, K, tau)

    @parameters t
    @variables f(t), y(t)
    Dt = Differential(t)
    
    eqs = [
           Dt(y) ~ (K*f - y)/tau
          ]

    ODESystem(eqs, t, [f, y], []; name=name)
end

@named model = firstOrder(K=2, tau=10)

# Define controller
function controller(;name, Kp, Ki, Kd)

    @parameters t, setpoint
    @variables error(t), integralError(t), manipulatedVariable(t), controlledVariable(t)
    Dt = Differential(t)

    equations = [
                 error ~ setpoint - controlledVariable
                 Dt(integralError) ~ error
                 Dt(manipulatedVariable) ~ expand_derivatives(Dt(setpoint)) - (manipulatedVariable - Kp * error - Ki * integralError) / Kd
                ]

    ODESystem(equations, t, [error, integralError, manipulatedVariable, controlledVariable], [setpoint]; name=name)

end

@named myController = controller(Kp = 10, Ki = 1, Kd = .1)

# Create control loop
@named loop = ODESystem([
                         model.y ~ myController.controlledVariable,
                         model.f ~ myController.manipulatedVariable
                        ], systems=[model, myController])

# Initial conditions and parameters
x0 = [
      model.f => 0,
      model.y => 0,
      myController.error => 0,
      myController.integralError => 0,
      myController.controlledVariable => 0,
      myController.manipulatedVariable => 0,
     ]

par = [myController.setpoint => 1]

# Create problem
sys = structural_simplify(loop)
prob = ODEProblem(sys, x0, (0.0, 20.0), par)

# Define callback
function condition(u, t, integrator)
    t - 10
end

function affect!(integrator) 
    integrator.p[1] += 1  # Setpoint = 2 after 10 seconds
end

cb = ContinuousCallback(condition, affect!)

# Solve
sol = solve(prob, Tsit5(), callback=cb)

# Plot
plot(sol, vars=[model.y, myController.error])

This results in:

While the controller behaves correctly, i.e. the controlled variable goes firstly to 1 and after 10 seconds moves up to 2 (the setpoint changes at this point thanks to the callback function), the error does not. The error is setpoint - controlledVariable, but it seems that because of the callback, which makes it 2 after 10 seconds, it is constantly 2 throughout the simulation.

The error without a callback starts at 1, as it should:

Why does this happen and how can I actually change the value of my setpoint so the error is 0 whenever I reach the desired steady-state?

[contradict]

The problem is that parameters are not recorded, so when error is re-computed for plotting, it uses the final value of setpoint. Moving the setpoint to a state produces the correct plot:

using ModelingToolkit, DifferentialEquations, Plots


# Define first order system
function firstOrder(;name, K, tau)

    @parameters t
    @variables f(t), y(t)
    Dt = Differential(t)

    eqs = [
           Dt(y) ~ (K*f - y)/tau
          ]

    ODESystem(eqs, t, [f, y], []; name=name)
end

@named model = firstOrder(K=2, tau=10)

# Define controller
function controller(;name, Kp, Ki, Kd)

    @parameters t
    @variables setpoint(t), error(t), integralError(t), manipulatedVariable(t), controlledVariable(t)
    Dt = Differential(t)

    equations = [
                 Dt(setpoint) ~ 0
                 error ~ setpoint - controlledVariable
                 Dt(integralError) ~ error
                 Dt(manipulatedVariable) ~ -(manipulatedVariable - Kp * error - Ki * integralError) / Kd
                ]

    ODESystem(equations, t, [setpoint, error, integralError, manipulatedVariable, controlledVariable], []; name=name)

end

@named myController = controller(Kp = 10, Ki = 1, Kd = .1)

# Create control loop
@named loop = ODESystem([
                         model.y ~ myController.controlledVariable,
                         model.f ~ myController.manipulatedVariable
                        ], systems=[model, myController])

# Initial conditions and parameters
x0 = [
      model.f => 0,
      model.y => 0,
      myController.setpoint => 1,
      myController.error => 0,
      myController.integralError => 0,
      myController.controlledVariable => 0,
      myController.manipulatedVariable => 0,
     ]

# Create problem
sys = structural_simplify(loop)
prob = ODEProblem(sys, x0, (0.0, 20.0))

# Define callback
function condition(u, t, integrator)
    t - 10
end

function affect!(integrator)
    integrator.u[2] = 2  # Setpoint = 2 after 10 seconds
end

cb = ContinuousCallback(condition, affect!)

# Solve
sol = solve(prob, Tsit5(), callback=cb)

# Plot
plot(sol, vars=[model.y, myController.error])

[jairorua]

Thanks for your reply.

I did also think about this possibility, but, at least to me, including a setpoint as a variable feels a bit more like a “work-around” than a proper way of modeling. Is there a way to enforce the solver to record the parameters so one can print them afterwards. The fact that just the last value is what shows up in the plot is misleading, specially for examples that are not trivial. Although this may be the nature of a callback.

Now some discussion/constructive criticism. I guess this discussion boils down to how the modelingtoolkit library considers both variables and parameters. For instance, Modelica has a clear distinction among variables, parameters and constants (see \pi). From the modelicatoolkit documentation is not entirely clear. I do not know if their distinction has already been decided or it is subjected to discussion, but it would be good for the user to have some clear definitions. I personally like quite a lot the trinity variable, parameter, constant. I found these two links (Contextual Variable Types · ModelingToolkit.jl, The AbstractSystem Interface · ModelingToolkit.jl) but I could not find a proper definition.

[ChrisRackauckas]

There is the same trinity, except Julia is a programming language so constants are handled as, well, const.

[contradict]

My mental model is variables change during each solve, parameters may change between solves, and constants are the same for all solves.

[zdenek_hurak]

That is why I have problems to view t as a parameter. If anything, I think it should be viewed as a variable. And indeed, the latest official examples already treat t as a variable: Composing Ordinary Differential Equations · ModelingToolkit.jl.

[jairorua]

Yes, but my concern here is with parameters. One can change it with callbacks and it propagates through the solution, but one cannot see that change in the actual parameter. So parameters behave as parameters but look like constants. So, is there a way to record potential changes?

[ChrisRackauckas]

You can manually record it in callbacks, push! to a vector or something. Parameters are not part of the saving mechanisms by default, though soon I hope to take another look at the saving mechanisms. I will say that if they were tracked by default you’d take a pretty big memory hit in a lot of cases.

[jairorua]

Sounds fair. It would be positive to put this somewhere in the documentation so users of the library have a common-ground and a unified practice for developing libraries based on MKT.

For example, as @zdenek_hurak mentioned in this post, t feels like a variable, but I have seen it many times defined as a parameter, and this is the main reason why I define it like that, although I see this has changed.

[jairorua]

Thanks for the info.

With all this I just to grasp the best practice using MKT, I do not want to pressure into any kind of development.

[zdenek_hurak]

Although not immediately related to the topic of this thread, now that there is some PID controller in the game here, I am taking a closer look. With the common feedback loop configuration and (re)labelling of the variables as in the figure below

             +----------------+    +---------+
 r        e  |                | u  |         |   y
------>o---->| PID Controller +--->|  Plant  +----+--->
       ^     |                |    |         |    |
     - |     +----------------+    +---------+    |
       |                                          |
       +------------------------------------------+

your controller is described in complex (Laplace) domain (with the common abuse of the notation for the time-domain variables) as

(Kd*s+1)*u = (Kp + Ki/s)*e,

which leads to

u = 1/(Kd*s+1)*(Kp + Ki/s)*e,

which is not really a prescription for a PID controller, is it? Instead, it is a PI controller with some low-pass filter.

[zdenek_hurak]

And I would argue similarly for the controller in the original post, which would translate to

(Kd*s+1)*u = Kd*s*r + (Kp + Ki/s)*e,

which only differentiates (and filters) the reference. In fact, sometimes it is preferable to do the opposite, that is, to leave the reference out from differentiation in order to save the controller from derivative kicks in response to step inputs, leading to saturation.

[jairorua]

Yes, I agree that this expression leads to the controller you mention.

\frac{\mathrm{d}u(t)}{\mathrm{d}t} = - (u(t) - K_{\mathrm{P}}\, e(t) - K_{\mathrm{I}}\, \int e(t)\, \mathrm{d}t)\, \frac{1}{K_{\mathrm{D}}} \\ s\, u(s) = - (u(s) - K_{\mathrm{P}}\, e(s) - K_{\mathrm{I}}\, \frac{e(s)}{s})\, \frac{1}{K_{\mathrm{D}}} \\ (K_{\mathrm{D}}\, s + 1)\, u(s) = (K_{\mathrm{P}} + \frac{K_{\mathrm{I}}}{s})\, e(s)

which is different from the Laplacre transform of a PID controller:

u(s) = (K_{\mathrm{P}} + \frac{K_{\mathrm{I}}}{s} + s\, K_{\mathrm{D}})\, e(s)

The idea here though was just to check how to include setpoint changes in a controller. It would be great to have a MKT library with controllers with different functionalities, e.g. anti-windup, no derivative kick, gain scheduling, etc. I may take a look to this.