Making a PD Controller Component in ModelingToolkit.jl

If I were to make a PD controller component, I’d probably do something like this:

using ModelingToolkit
function PD_controller(Pgain, Dgain, num_inputs; name=:pd)
    @parameters t
    D = Differential(t)
    @variables u[1:num_inputs](t)
    @variables o[1:num_inputs](t)

    eqs = o .~ Pgain.*(u) + Dgain.* D.(u)
    return ODESystem(eqs, t, vcat(u,o), []; name = name)
end

When I hook up an actual control system involving a PD controller, I cant structural_simplify the system. (It’s fine without the derivative component). The stacktrace is

… Internal error in Tearing.jl: vs[1] = 0.

I guess this is because the component needs an explicit equation for D(u), which it doesn’t have. Now in my case, D(u) is defined, since u ~ another_component.x , where another_component.x has its dynamics defined. But it’s defined in the equations for a different component. Is there a way for structural_simplify to be extended to deal with such cases, or is there a particular coding pattern that allows for it without requiring the PD controller to receive D(u) as a separate input?

Thanks in advance! It really feels like ModelingToolkit.jl could become a more flexible, fully differentiable version of Simulink in the near future.

1 Like

Normally, PD controllers are implemented with a filter to make them proper. Does that help? In other words:

u(s) = K_\mathrm{p}\frac{1+T_\mathrm{d}s}{1+\gamma T_\mathrm{d}s}\cdot e(s)

with typical value \gamma = 0.1, instead of the “pure” PD controller:

u(s) = K_\mathrm{p}(1+T_\mathrm{d}s)\cdot e(s)
3 Likes

Indeed.

I will also add that while the filtering part makes the controller proper, which is plausible from a mathematical/fundamental viewpoint, it also does the very practical job of filtering, which is certainly plausible from an engineering viewpoint.

And while discussing these practical aspects, it is perhaps worth mentioning that in practical implementations it is often wise to apply the (filtered) derivative part of the controller to the measured plant outputs y only. The proportional (and possibly the integration part, if present) of the controller are applied to the regulation error e = r-y (sometimes just the integral part is applied to the regulation error). The reason is that if there are abrupt changes in the reference signal r (Heaviside steps), the derivative part gets immediately saturated. Some discussion of this is in chapter 11.5 of http://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-public_24Jul2020.pdf (paragraph on Setpoint weighting) or elsewhere.

2 Likes

That would be great. Let’s not forget, however, what the key reason for the popularity of Simulink in the control systems community is (and I really mean not just the academic community but also at least some industries, certainly the automotive one). It is certainly not (just) the convenience and/or performance of its numerical solvers for this and that. There are other software tools – commercial and FOSS – that offer competitive features and performance for the mere control design related computations. But what matters is that from a controller designed in Matlab and implemented in Simulink an engineer can generate a code optimized for running in real time on a wide selection of target platforms.

Thanks, yes, that is a good point! I guess the PD controller I made was just a minimal example highlighting the issue I was asking about:

Component A has input u. The dynamics of component A depend upon du/dt.

du/dt is defined in the equations for Component X, but not component A. One could explicitly feed du/dt as an extra input to component X if necessary. But it would be nice if one could just take the derivative of u inside component A, and when connecting the systems and using structural_simplify, the simplification algorithm realises that du/dt is defined elsewhere, and hence the differential equations are complete and can compile.

Thanks for the insight and the references

Yes definitely, I’m speaking with the bias of somebody who isn’t implementing code on a target platform in the wild!

This is my naive version of a classic PID:

function controlPID(;name)
    @parameters Kc τI τD
    @variables ym(t) ε(t) c(t) x(t)
    eqs = [
        D(ym) ~ expand_derivatives(D(sp)) - (c - Kc*ε - Kc/τI*x)/(Kc*τD),
        D(x) ~ ε
    ]
    ODESystem(eqs, name=name)
end

The error is defined as \varepsilon = y_{sp} -y_m. To avoid calculating \frac{\mathrm{d} \varepsilon}{\mathrm{d}t}, I differentiated the error \frac{\mathrm{d} \varepsilon}{\mathrm{d}t} = \frac{\mathrm{d} y_{sp}}{\mathrm{d}t} -\frac{\mathrm{d} y_m}{\mathrm{d}t}.

c is the output of the controller.

I let ModelingToolkit to differentiate y_{sp}.

I hope that everything is cleat, if not I can share a more complete example.

1 Like

I haven’t started looking too much into the ModelingToolkit yet, so I’m curious as to how you connect it to a simple systems in a feedback configuration?

  • Where is the setpoint (sp) provided?
  • Where is the control error (ε) defined?

A simple example with, say, controlling a second order system, would be very useful for me and others to try to set aside time to absorb the possibilities of ModelingToolkit :-).

This is a long step-by-step tutorial to simulate a PID loop. Keep in mind that my knowledge of Julia and ModelingToolkit is very limited. I’m pretty sure that there are a lot of things to be improved.

You can use Pluto to run this tutorial.

First, we will load the libraries:

using ModelingToolkit, DifferentialEquations, Plots

Then, we will define time t, as a parameter of the simulation and a couple of functions that we will need: error \varepsilon(t) and setpoint y_{sp}(t) as variables of the simulation:

## Removed to define t as a variable
## @parameters t
@variables t, ε(t), ysp(t)
D = Differential(t)

As set point we can use any function, like Heaviside, sin, or something more complex:

# Set-point function:

# sp = t*(1-0.5*(1+tanh((t-5)/.1)))

sp = 1

The control loop will have the following elements:

  • Summing point, to calculate the error: \varepsilon = y_{sp} - y_m.
  • PID controller: K_c = 10. \tau_I = 2., and \tau_D = 0.1
  • Plant: First order (K_p = 2 and \tau_p = 10)
  • Sensor: First order (K_m = 1 and \tau_m = 1)

We set up all these parameters:

# Simulation parameters

begin
	Kp = 2
	τp = 10
	Km = 1
	τm = 1
	Kc = 10.
	τI = 2.
	τD = 0.1

	tend = 20.
end

Now, we will define a generic function of a first order system:

# First order system: G(s) = y(t)/f(s) = K/(τ*s +1)
#                     τ dy(t)/dt +y(t) = K f(t)

function firstOrder(;name)
    @parameters K τ
    @variables f(t) y(t)
    eqs = [
        D(y) ~ (K*f - y)/τ
        ]
    ODESystem(eqs, name=name)
end

And we will define the PID controller:

# PID controller: Gc(s) = c(s)/ε(s) = Kc*(1 + 1/(τI*s) + τD*s)
#                 c(t) = Kc*(ε(t) + 1/τI*∫ε(t)dt + τD*dε(t)/dt
#                 x(t) = ∫ε(t)dt
#                 ε(t) = ysp(t) - ym(t)

function controlPID(;name)
    @parameters Kc τI τD
    @variables ym(t) ε(t) c(t) x(t)
    eqs = [
        D(ym) ~ expand_derivatives(D(sp)) - (c - Kc*ε - Kc/τI*x)/(Kc*τD),
        D(x) ~ ε
    ]
    ODESystem(eqs, name=name)
end

The next step is setting up the plant, sensor and controller:

# Process (Plant): First order system
# Input: c(t)
# Output: y(t)

@named proc = firstOrder()
# Sensor: First order system
# Input: y(t)
# Output: ym(t) 

@named sensor = firstOrder()
# PID controller:
# Inputs: ysp(t), ym(t)
# Output: c(t)

@named control_PID = controlPID()

After defining all the elements of our control loop, we will define the loop itself:

@named loopPID = ODESystem([
        sensor.f ~ proc.y,
        proc.f ~ control_PID.c,
        control_PID.ε ~ ε,
        control_PID.ym ~ sensor.y,
        ε ~ ysp - sensor.y,
        ysp ~ sp
        ], t, systems=[proc, sensor, control_PID]);

I think that the name of the variables is clear. If not, I can prepare a drawing to make them clearer.

Now we have defined our system, but is a DAE and it has more equations than necessary. We can use ModelingToolkit to make all the simplifications:

sys = structural_simplify(dae_index_lowering(loopPID))

We have almost finished. We need to define the initial conditions and set up the parameters of the loop:

# Initial conditions

u0_PID = [
    proc.y => 0.0,
    proc.f => .0,
    sensor.y => 0.0,
    sensor.f => 0.0,
    control_PID.c => 0.0,
    control_PID.ε => 0,
    control_PID.x => 0,
    control_PID.ym => 0,
    ε => 0
    ]
p_PID = [
    proc.τ => τp,
    proc.K => Kp,
    sensor.τ => τm,
    sensor.K => Km,
    control_PID.Kc => Kc,
    control_PID.τI => τI,
    control_PID.τD => τD
    ]

We have everything set, so we can ask ModelingToolkit to create the problem:

probPID = ODEProblem(sys, u0_PID, (0.0, tend), p_PID)

Finally, we solve the problem:

solPID = solve(probPID, ImplicitEuler())

Done!

We can plot the solution easily:

plot(solPID)

If we want to choose which variable to be plotted:

plot(solPID, vars=[proc.y,sensor.y])

We can test if the solution is right. To do it, we will use ControlSystems.jl:

using ControlSystems
s = tf("s")
Gp = Kp/(τp*s+1)
Gm = Km/(τm*s+1)
Gc = Kc*(1 + 1/(τI*s) + τD*s)
G = Gc*Gp/(1+Gc*Gp*Gm)

begin
	stepplot(G, tend)
	plot!(solPID, vars=[proc.y])
end

Please let me know which points are not clear or which ones could be improved and I will edit this post.

5 Likes

I find it confusing to regard (and declare) the time as a parameter. The parameters are constant throughout the simulation, aren’t they? Shouldn’t the time t be better viewed as a variable? The example at Composing Ordinary Differential Equations · ModelingToolkit.jl suggests so. Note also that the issue has been filed Comment on terminology: parameter vs. (independent) variable t · Issue #564 · SciML/ModelingToolkit.jl · GitHub.

I think you’re right. Probably I learnt that in an old version of ModelingToolkit documentation. I’m going to edit my post to reflect your suggestion .

1 Like

The meaning of parameter varies among scientific fields. In Modelica, there are two types of constants: constant and parameter, where the former is a physics/mathematics constant such as \pi, Planck’s constant, etc., while the latter are constant that vary from one simulation to the next (e.g., geometric quantities, etc.). Thus, in Modelica (and some mathematics fields), a “parameter” is a constant.

In some physics fields, “parameter” means a different thing. As an example, in “distributed parameter system”.

I haven’t really seen a formal definition of parameter in ModelingToolkit, but the word seems to include constants as well as independent variables (including time). Possibly also functions of independent variables, but I’m not sure. (That is, according to some older documentation of ModelingToolkit.jl. It is possible that the use of concepts has changed over time.)

The concept of parameter doesn’t exists in Symbolics.jl (as far as I know), and is an extension of variables from Symbolics.jl — I think Chris R. mentioned that parameters are variables with some added meta information, but this is not entirely clear to me.

1 Like

Is there a typo in this sequence? (Should there be a keyword here at the end of the last line?)

When I remove the word here, your code runs fine in IJulia.

Thank you for noticing that mistake. I have removed it.

1 Like

One question: when you create models such as firstOrder

# First order system: G(s) = y(t)/f(s) = K/(τ*s +1)
#                     τ dy(t)/dt +y(t) = K f(t)

function firstOrder(;name)
    @parameters K τ
    @variables f(t) y(t)
    eqs = [
        D(y) ~ (K*f - y)/τ
        ]
    ODESystem(eqs, name=name)
end

… have you found a way to give default values to your parameters K and τ that kick in if you don’t specify the value in your proc/sensor parameter specification…

p_PID = [
    proc.τ => τp,
    proc.K => Kp,
    sensor.τ => τm,
    sensor.K => Km,
    control_PID.Kc => Kc,
    control_PID.τI => τI,
    control_PID.τD => τD
    ]

?

Just after Saving my answer, I noticed that I was answering a different question… Sorry! Next time I’ll read the question more carefully.

I have found a way, I’m not sure if it’s the best way but it works.

You can define the parameters for each block of the simulation. First, the plant (process):

p_PROC = [proc.K => Kp, proc.τ => τp]

Then, the sensor:

p_SENSOR = [sensor.K => Km, sensor.τ => τm]

And, finally, the controller:

p_PID = [
    control_PID.Kc => Kc,
    control_PID.τI => τI,
    control_PID.τD => τD
    ]

To define the simulation problem, we need to put them together:

p_LOOP = union(p_PROC, p_SENSOR, p_PID)

So, now the problem is:

probPID = ODEProblem(sys, u0_PID, (0.0, tend), p_LOOP)

We could do the same for the initial conditions.

1 Like

I’ll reply more fully too all the comments on this thread a bit later.

For now, what I would say for defaults is:

for each component ODESystem, when you generate the ODESystem, you can add a keyword argument: defaults = Dict( tau => 1., K => 2.), and so on. The entries in this dict comprise either/both initial conditions and parameters. When you connect up the system, the default values are inherited

1 Like

Does this mean:

# First order system: G(s) = y(t)/f(s) = K/(τ*s +1)
#                     τ dy(t)/dt +y(t) = K f(t)

function firstOrder(;name)
    @parameters K τ
    @variables f(t) y(t)
    eqs = [
        D(y) ~ (K*f - y)/τ
        ]
    ODESystem(eqs, name=name; defaults = Dict(tau => 1., K => 2.))
end

?

Yes, see ‘Defaults’ section here: Composing Ordinary Differential Equations · ModelingToolkit.jl

1 Like