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.