DiscretePIDs

This package implements a discrete-time PID controller as an approximation of the continuous-time PID controller given by $$U(s) = K \left( w_p R(s) - Y(s) + \dfrac{1}{sT_i} \left( R(s) - Y(s) \right) + \dfrac{sT_d}{1 + s T_d / N}(w_d R(s) - Y(s)) \right) + U_\textrm{ff}(s),$$ where

• $u(t) \leftrightarrow U(s)$ is the control signal
• $y(t) \leftrightarrow Y(s)$ is the measurement signal
• $r(t) \leftrightarrow R(s)$ is the reference / set point
• $u_\textrm{ff}(t) \leftrightarrow U_\textrm{ff}(s)$ is the feed-forward contribution
• $K$ is the proportional gain
• $T_i$ is the integral time
• $T_d$ is the derivative time
• $N$ is a parameter that limits the gain of the derivative term at high frequencies, typically ranges from 2 to 20,
• $w_p \in [0, 1]$ is a parameter that gives the proportion of the reference signal that appears in the proportional term.
• $w_d \in [0, 1]$ is a parameter that gives the proportion of the reference signal that appears in the derivative term (default 0).

Saturation of the controller output is parameterized by $u_{\min}$ and $u_{\max}$, and the integrator anti-windup is parameterized by the tracking time $T_\mathrm{t}$.

Usage

Construct a controller by

pid = DiscretePID(; K = 1, Ti = false, Td = false, Tt = √(Ti*Td), N = 10, wp = 1, wd = 0, umin = -Inf, umax = Inf, Ts, I = 0, D = 0, yold = 0)

and compute the control signal at a given time using

u = pid(r, y, uff)

or

u = calculate_control!(pid, r, y, uff)

The derivative term is by default computed by filtering the measurement $y$, but it can also be sourced externally by setting the keyword argument yd:

u = calculate_control!(pid, r, y, uff; yd)

When yd is provided, no filtering is applied by the PID controller, i.e., $N$ is ignored. This is useful when the derivative is computed externally, e.g., from a velocity sensor or an observer.

The parameters $K$, $T_i$, and $T_d$ may be updated using the functions set_K!, set_Ti!, and set_Td!, respectively.

The numeric type used by the controller (the T in DiscretePID{T}) is determined by the types of the parameters. To use a custom number type, e.g., a fixed-point number type, simply pass the parameters as that type, see example below. The controller will automatically convert measurements and references to this type before performing the control calculations.

The internal state of the controller can be reset to zero using the function reset_state!(pid). If repeated simulations using the same controller object are performed, the state should likely be reset between simulations.

Examples

Example using ControlSystems.jl

The following example simulates a feedback control system containing a PID controller using ControlSystems.jl package. We simulate a response of the closed-loop system to the step disturbance $d(t) = 1$ entering at the plant (the system to be controlled) input, while the reference is $r(t) = 0$.

using DiscretePIDs, ControlSystemsBase, Plots
Tf = 30   # Simulation time
K  = 1    # Proportional gain
Ti = 1    # Integral time
Td = 1    # Derivative time
Ts = 0.01 # sample time

P   = c2d(ss(tf(1, [1, 1])), Ts) # Process to be controlled, discretized using zero-order hold
pid = DiscretePID(; K, Ts, Ti, Td)

ctrl = function(x,t)
    y = (P.C*x)[] # measurement
    d = 1         # disturbance
    r = (t >= 15) # reference
    u = pid(r, y) # control signal
    u + d # Plant input is control signal + disturbance
end

res = lsim(P, ctrl, Tf)

plot(res, plotu=true); ylabel!("u + d", sp=2)

Here we simulated a linear plant, in which case we were able to call ControlSystems.lsim specialized for linear systems. Below, we show two methods for simulation that works with a nonlinear plant, but we still use a linear system to make the comparison easier.

For comparison, we also perform the same simulation with a two degree-of-freedom PID controller

using ControlSystemsBase, DiscretePIDs, Plots
t = 0:Ts:Tf
u = [ones(length(t)) t .>= 15]' # Input signal [d; r]
C = pid_2dof(K, Ti, Td; Ts, N=10)
Gcl = feedback(P, C, W1=1, U2=2, W2=1, Z2=1, pos_feedback=true)
simres = lsim(Gcl, u)
plot(simres, plotu=true, lab=["y" "u" "d" "r"], layout=(2,1), sp=[1 2 2 1], ylabel="")

Please note: The result of simulation with a controller computed by pid_2dof will be slightly different due to a difference in discretization. DiscretePIDs uses a forward-Euler approximation for the integrator and a backward Euler approximation for the derivative, while pid_2dof uses a Tustin approximation (default) for both.

Example using DifferentialEquations.jl

This example is identical to the one above except for using DifferentialEquations.jl for the simulation, which makes it possible to consider more complex plants, in particular nonlinear ones.

There are several different ways one could go about including a discrete-time controller in a continuous-time simulation, in particular, we must choose a way to store the computed control variable. Two common approaches are

1. We use a global variable into which we write the control signal at each discrete time step.
2. We add an extra state variable to the system, and use it to store the control variable.

In this example we choose the latter approach, since it has the added benefit of adding the computed control variable to the solution object.

We use DiffEqCallbacks.PeriodicCallback, in which we perform the PID-controller update, and store the computed control signal in the extra state variable.

using DiscretePIDs, ControlSystemsBase, OrdinaryDiffEq, DiffEqCallbacks, Plots

Tf = 30   # Simulation time
K  = 1    # Proportional gain
Ti = 1    # Integral time
Td = 1    # Derivative time
Ts = 0.01 # sample time

P = ss(tf(1, [1, 1]))  # Process to be controlled in continuous time
A, B, C, D = ssdata(P) # Extract the system matrices
pid = DiscretePID(; K, Ts, Ti, Td)

function dynamics!(dxu, xu, p, t)
    A, B, C, r, d = p   # We store the reference and disturbance in the parameter object
    x = xu[1:P.nx]      # Extract the state
    u = xu[P.nx+1:end]  # Extract the control signal
    dxu[1:P.nx] .= A*x .+ B*(u .+ d) # Plant input is control signal + disturbance
    dxu[P.nx+1:end] .= 0             # The control signal has no dynamics, it's updated by the callback
end

cb = PeriodicCallback(Ts) do integrator
    p = integrator.p    # Extract the parameter object from the integrator
    (; C, d) = p     # Extract the reference and disturbance from the parameter object
    x = integrator.u[1:P.nx] # Extract the state (the integrator uses the variable name `u` to refer to the state, in control theory we typically use the variable name `x`)
    r = (integrator.t >= 15) # Reference
    y = (C*x)[]         # Simulated measurement
    u = pid(r, y)       # Compute the control signal
    integrator.u[P.nx+1:end] .= u # Update the control-signal state variable 
end

parameters = (; A, B, C, d=1) # disturbance = 1
xu0 = zeros(P.nx + P.nu) # Initial state of the system + control signals
prob = ODEProblem(dynamics!, xu0, (0, Tf), parameters, callback=cb) # disturbance = 1
sol = solve(prob, Tsit5(), saveat=Ts)

plot(sol, layout=(2, 1), ylabel=["x" "u"], lab="")

The figure should look more or less identical to the one above, except that we plot the control signal $u$ instead of the combined input $u + d$ like we did above. Due to the fast sample rate $T_s$, the control signal looks continuous, however, increase $T_s$ and you'll notice the zero-order-hold nature of $u$.

Example using SeeToDee.jl

SeeToDee.jl is a library of fixed-time-step integrators useful for "manual" (=one integration step at a time) simulation of control systems. The same example as above is simulated using SeeToDee.Rk4 here. The call to

discrete_dynamics = SeeToDee.Rk4(dynamics, Ts)

considers the continuous-time dynamical system modelled by

\dot x(t) = f(x(t), u(t), p(t), t)

and at a given state $x$ and time $t$ and for a given control $u$, it computes an approximation $x^+$ to the state $x(t+T_s)$ at the next time step $t+T_s$

x(t+T_s) \approx x^+ = \phi(x(t), u(t), p(t), t,T_s).

using DiscretePIDs, ControlSystemsBase, SeeToDee, Plots
Tf = 30   # Simulation time
K  = 1    # Proportional gain
Ti = 1    # Integral time
Td = 1    # Derivative time
Ts = 0.01 # sample time
P  = ss(tf(1, [1, 1]))    # Process to be controlled, in continuous time
A,B,C = ssdata(P)         # Extract the system matrices
p = (; A, B, C, d=1) # reference = 0, disturbance = 1

pid = DiscretePID(; K, Ts, Ti, Td)

ctrl = function(x,p,t)
    r = (t >= 15)   # reference
    y = (p.C*x)[]   # measurement
    pid(r, y)
end

function dynamics(x, u, p, t) # This time we define the dynamics as a function of the state and control signal
    A, B, C, d = p   # We store the reference and disturbance in the parameter object
    A*x .+ B*(u .+ d) # Plant input is control signal + disturbance
end
discrete_dynamics = SeeToDee.Rk4(dynamics, Ts) # Create a discrete-time dynamics function

x = zeros(P.nx) # Initial condition
X, U = [], []   # To store the solution 
t = range(0, step=