Digital PID Control - Part 1

Digital PID Control - Part 1#

Required imports#

from IPython.core.display import HTML
from sympy import *
from mathprint import *

Define symbols#

Kp, Ki, Kd = symbols('K_p K_i K_d', positive=True)

Ts, N = symbols('T_s N', positive=True)
s = symbols('s', complex=True)
z = symbols('z')

e, u = symbols('e u')
e0, e1, e2 = symbols('e_k e_{k-1} e_{k-2}')
u0, u1, u2 = symbols('u_k u_{k-1} u_{k-2}')

Continuous-time PID#

Let us start with defining our PID control in a continuous form. Here, we are going to follow MATLAB Simulink implementation. In Simulink, the derivative term is a first-order low-pass filter.

\[ U(s) = \Bigg( K_p + K_i \underbrace{\boxed{\frac{1}{s}}}_{\text{integrator}} + \underbrace{\boxed{s}}_{\text{differentiator}} K_d \Bigg) E(s) \]

where:

\(K_p\) is the proportional gain
\(K_i\) is the integral gain
\(K_d\) is the derivative gain
\(e(t)\) or \(E(s)\) is the input to the controller, which is the error between the target and the actual output of the system
\(u(t)\) or \(U(s)\) is the computed control signal that will be sent to the system

Discrete-time PID#

Let us use the buliear transform method for both the integral term and the derivative term.

\[ \frac{1}{s} \longleftarrow \frac{T_s}{2} \frac{z+1}{z-1} \]

Gpid = Kp + Ki/s + Kd*s     # PID-control
mprint("G_{PID}=", latex(Gpid))

Gpid = Gpid.subs(1/s, Ts/2 * (z+1)/(z-1) )
mprint("G_{PID}=", latex(Gpid))

Gpid = factor(Gpid)

Gpid = collect(numer(Gpid), [z**2, z]) /  collect(expand(denom(Gpid)), [2*Ts])
mprintb("G_{PID}=", latex(Gpid))

\[\displaystyle G_{PID}=K_{d} s + \frac{K_{i}}{s} + K_{p}\]

\[\displaystyle G_{PID}=\frac{2 K_{d} \left(z - 1\right)}{T_{s} \left(z + 1\right)} + \frac{K_{i} T_{s} \left(z + 1\right)}{2 \left(z - 1\right)} + K_{p}\]

\[\displaystyle \boxed{G_{PID}=\frac{4 K_{d} + K_{i} T_{s}^{2} - 2 K_{p} T_{s} + z^{2} \left(4 K_{d} + K_{i} T_{s}^{2} + 2 K_{p} T_{s}\right) + z \left(- 8 K_{d} + 2 K_{i} T_{s}^{2}\right)}{2 T_{s} \left(z^{2} - 1\right)}}\]

eq = Eq(u, Gpid * e)
mprint("\\small ", latex(eq))

eq = Eq(numer(eq.rhs), eq.lhs * denom(eq.rhs))
mprint("\\small ", latex(eq))

eq = expand(Eq(eq.lhs/z**2/Ts/2, eq.rhs/z**2/Ts/2))
mprint("\\small ", latex(eq))

eq = eq.subs(e/z**2, e2).subs(e/z, e1).subs(e, e0).subs(u/z**2, u2).subs(u/z, u1).subs(u, u0)
mprint("\\small ", latex(eq))

eq_ = Eq(collect(eq.rhs, [u0, u1, u2]), collect(expand(eq.lhs), [e0, e1, e2]))
mprintb("\\small ", latex(eq_))

\[\displaystyle \small u = \frac{e \left(4 K_{d} + K_{i} T_{s}^{2} - 2 K_{p} T_{s} + z^{2} \left(4 K_{d} + K_{i} T_{s}^{2} + 2 K_{p} T_{s}\right) + z \left(- 8 K_{d} + 2 K_{i} T_{s}^{2}\right)\right)}{2 T_{s} \left(z^{2} - 1\right)}\]

\[\displaystyle \small e \left(4 K_{d} + K_{i} T_{s}^{2} - 2 K_{p} T_{s} + z^{2} \left(4 K_{d} + K_{i} T_{s}^{2} + 2 K_{p} T_{s}\right) + z \left(- 8 K_{d} + 2 K_{i} T_{s}^{2}\right)\right) = 2 T_{s} u \left(z^{2} - 1\right)\]

\[\displaystyle \small \frac{2 K_{d} e}{T_{s}} - \frac{4 K_{d} e}{T_{s} z} + \frac{2 K_{d} e}{T_{s} z^{2}} + \frac{K_{i} T_{s} e}{2} + \frac{K_{i} T_{s} e}{z} + \frac{K_{i} T_{s} e}{2 z^{2}} + K_{p} e - \frac{K_{p} e}{z^{2}} = u - \frac{u}{z^{2}}\]

\[\displaystyle \small \frac{2 K_{d} e_{k}}{T_{s}} - \frac{4 K_{d} e_{k-1}}{T_{s}} + \frac{2 K_{d} e_{k-2}}{T_{s}} + \frac{K_{i} T_{s} e_{k}}{2} + K_{i} T_{s} e_{k-1} + \frac{K_{i} T_{s} e_{k-2}}{2} + K_{p} e_{k} - K_{p} e_{k-2} = u_{k} - u_{k-2}\]

\[\displaystyle \boxed{\small u_{k} - u_{k-2} = e_{k} \left(\frac{2 K_{d}}{T_{s}} + \frac{K_{i} T_{s}}{2} + K_{p}\right) + e_{k-1} \left(- \frac{4 K_{d}}{T_{s}} + K_{i} T_{s}\right) + e_{k-2} \left(\frac{2 K_{d}}{T_{s}} + \frac{K_{i} T_{s}}{2} - K_{p}\right)}\]

Continuous-time “PID with derivative filter”#

Let us start with defining our PID control in a continuous form. Here, we are going to follow MATLAB Simulink implementation. In Simulink, the derivative term is a first-order low-pass filter.

\[ U(s) = \Bigg( K_p + K_i \underbrace{\boxed{\frac{1}{s}}}_{\text{integrator}} + K_d \frac{N}{1+N \underbrace{\boxed{\frac{1}{s}}}_{\text{filter}}} \Bigg) E(s) \]

where:
\(K_p\) is the proportional gain
\(K_i\) is the integral gain
\(K_d\) is the derivative gain
\(N\) is the filter coefficient
\(e(t)\) or \(E(s)\) is the input to the controller, which is the error between the target and the actual output of the system
\(u(t)\) or \(U(s)\) is the computed control signal that will be sent to the system

As for the discretization methods, we will use the tarpezoidal (bilinear transform or Tustin) metod. From this point onwards, let us introduce \(T_s\) as the sampling period of the discrete-time PID control.

Discrete-time “PID with derivative filter”#

Let us use the buliear transform method for both the integral term and the filter.

\[ \frac{1}{s} \longleftarrow \frac{T_s}{2} \frac{z+1}{z-1} \]

Gpid = Kp + Ki/s + Kd*N/(1+N/s)     # PID-control
mprint("G_{PID}=", latex(Gpid))

Gpid = Gpid.subs(1/s, Ts/2 * (z+1)/(z-1) )
mprint("G_{PID}=", latex(Gpid))

Gpid = factor(Gpid)

Gpid = collect(numer(Gpid), [z**2, z]) /  collect(expand(denom(Gpid)), [z**2, z])
mprintb("\\small ", "G_{PID}=", latex(Gpid))

\[\displaystyle G_{PID}=\frac{K_{d} N}{\frac{N}{s} + 1} + \frac{K_{i}}{s} + K_{p}\]

\[\displaystyle G_{PID}=\frac{K_{d} N}{\frac{N T_{s} \left(z + 1\right)}{2 \left(z - 1\right)} + 1} + \frac{K_{i} T_{s} \left(z + 1\right)}{2 \left(z - 1\right)} + K_{p}\]

\[\displaystyle \boxed{\small G_{PID}=\frac{4 K_{d} N + K_{i} N T_{s}^{2} - 2 K_{i} T_{s} - 2 K_{p} N T_{s} + 4 K_{p} + z^{2} \left(4 K_{d} N + K_{i} N T_{s}^{2} + 2 K_{i} T_{s} + 2 K_{p} N T_{s} + 4 K_{p}\right) + z \left(- 8 K_{d} N + 2 K_{i} N T_{s}^{2} - 8 K_{p}\right)}{- 2 N T_{s} + z^{2} \left(2 N T_{s} + 4\right) - 8 z + 4}}\]

Now, let us apply the input to the output and turn the transfer function into algorithm (or difference equation).

eq = Eq(u, Gpid * e)
mprint("\\small ", latex(eq))

eq = Eq(numer(eq.rhs), eq.lhs * denom(eq.rhs))
mprint("\\small ", latex(eq))

eq = expand(Eq(eq.lhs/z**2/Ts/N/8, eq.rhs/z**2/Ts/N/8))
mprint("\\small ", latex(eq))

eq = eq.subs(e/z**2, e2).subs(e/z, e1).subs(e, e0).subs(u/z**2, u2).subs(u/z, u1).subs(u, u0)
mprint("\\small ", latex(eq))

eq = Eq(collect(eq.rhs, [u0, u1, u2]), collect(expand(eq.lhs), [e0, e1, e2]))
mprintb("\\small ", latex(eq))

\[\displaystyle \small u = \frac{e \left(4 K_{d} N + K_{i} N T_{s}^{2} - 2 K_{i} T_{s} - 2 K_{p} N T_{s} + 4 K_{p} + z^{2} \left(4 K_{d} N + K_{i} N T_{s}^{2} + 2 K_{i} T_{s} + 2 K_{p} N T_{s} + 4 K_{p}\right) + z \left(- 8 K_{d} N + 2 K_{i} N T_{s}^{2} - 8 K_{p}\right)\right)}{- 2 N T_{s} + z^{2} \left(2 N T_{s} + 4\right) - 8 z + 4}\]

\[\displaystyle \small e \left(4 K_{d} N + K_{i} N T_{s}^{2} - 2 K_{i} T_{s} - 2 K_{p} N T_{s} + 4 K_{p} + z^{2} \left(4 K_{d} N + K_{i} N T_{s}^{2} + 2 K_{i} T_{s} + 2 K_{p} N T_{s} + 4 K_{p}\right) + z \left(- 8 K_{d} N + 2 K_{i} N T_{s}^{2} - 8 K_{p}\right)\right) = u \left(- 2 N T_{s} + z^{2} \left(2 N T_{s} + 4\right) - 8 z + 4\right)\]

\[\displaystyle \small \frac{K_{d} e}{2 T_{s}} - \frac{K_{d} e}{T_{s} z} + \frac{K_{d} e}{2 T_{s} z^{2}} + \frac{K_{i} T_{s} e}{8} + \frac{K_{i} T_{s} e}{4 z} + \frac{K_{i} T_{s} e}{8 z^{2}} + \frac{K_{i} e}{4 N} - \frac{K_{i} e}{4 N z^{2}} + \frac{K_{p} e}{4} - \frac{K_{p} e}{4 z^{2}} + \frac{K_{p} e}{2 N T_{s}} - \frac{K_{p} e}{N T_{s} z} + \frac{K_{p} e}{2 N T_{s} z^{2}} = \frac{u}{4} - \frac{u}{4 z^{2}} + \frac{u}{2 N T_{s}} - \frac{u}{N T_{s} z} + \frac{u}{2 N T_{s} z^{2}}\]

\[\displaystyle \small \frac{K_{d} e_{k}}{2 T_{s}} - \frac{K_{d} e_{k-1}}{T_{s}} + \frac{K_{d} e_{k-2}}{2 T_{s}} + \frac{K_{i} T_{s} e_{k}}{8} + \frac{K_{i} T_{s} e_{k-1}}{4} + \frac{K_{i} T_{s} e_{k-2}}{8} + \frac{K_{i} e_{k}}{4 N} - \frac{K_{i} e_{k-2}}{4 N} + \frac{K_{p} e_{k}}{4} - \frac{K_{p} e_{k-2}}{4} + \frac{K_{p} e_{k}}{2 N T_{s}} - \frac{K_{p} e_{k-1}}{N T_{s}} + \frac{K_{p} e_{k-2}}{2 N T_{s}} = \frac{u_{k}}{4} - \frac{u_{k-2}}{4} + \frac{u_{k}}{2 N T_{s}} - \frac{u_{k-1}}{N T_{s}} + \frac{u_{k-2}}{2 N T_{s}}\]

\[\displaystyle \boxed{\small u_{k} \left(\frac{1}{4} + \frac{1}{2 N T_{s}}\right) + u_{k-2} \left(- \frac{1}{4} + \frac{1}{2 N T_{s}}\right) - \frac{u_{k-1}}{N T_{s}} = e_{k} \left(\frac{K_{d}}{2 T_{s}} + \frac{K_{i} T_{s}}{8} + \frac{K_{i}}{4 N} + \frac{K_{p}}{4} + \frac{K_{p}}{2 N T_{s}}\right) + e_{k-1} \left(- \frac{K_{d}}{T_{s}} + \frac{K_{i} T_{s}}{4} - \frac{K_{p}}{N T_{s}}\right) + e_{k-2} \left(\frac{K_{d}}{2 T_{s}} + \frac{K_{i} T_{s}}{8} - \frac{K_{i}}{4 N} - \frac{K_{p}}{4} + \frac{K_{p}}{2 N T_{s}}\right)}\]

Separating the control terms#

The proportional term#

By setting \(Ki = 0\) and \(K_d=0\), we obtain:

eq.subs(([Ki,0], [Kd,0]))

\[\displaystyle u_{k} \left(\frac{1}{4} + \frac{1}{2 N T_{s}}\right) + u_{k-2} \left(- \frac{1}{4} + \frac{1}{2 N T_{s}}\right) - \frac{u_{k-1}}{N T_{s}} = - \frac{K_{p} e_{k-1}}{N T_{s}} + e_{k} \left(\frac{K_{p}}{4} + \frac{K_{p}}{2 N T_{s}}\right) + e_{k-2} \left(- \frac{K_{p}}{4} + \frac{K_{p}}{2 N T_{s}}\right)\]

The integral term#

By setting \(Kp = 0\) and \(K_d=0\), we obtain:

eq.subs(([Kp,0], [Kd,0]))

\[\displaystyle u_{k} \left(\frac{1}{4} + \frac{1}{2 N T_{s}}\right) + u_{k-2} \left(- \frac{1}{4} + \frac{1}{2 N T_{s}}\right) - \frac{u_{k-1}}{N T_{s}} = \frac{K_{i} T_{s} e_{k-1}}{4} + e_{k} \left(\frac{K_{i} T_{s}}{8} + \frac{K_{i}}{4 N}\right) + e_{k-2} \left(\frac{K_{i} T_{s}}{8} - \frac{K_{i}}{4 N}\right)\]

The derivative term#

By setting \(Kp = 0\) and \(K_i=0\), we obtain:

eq.subs(([Kp,0], [Ki,0]))

\[\displaystyle u_{k} \left(\frac{1}{4} + \frac{1}{2 N T_{s}}\right) + u_{k-2} \left(- \frac{1}{4} + \frac{1}{2 N T_{s}}\right) - \frac{u_{k-1}}{N T_{s}} = \frac{K_{d} e_{k}}{2 T_{s}} - \frac{K_{d} e_{k-1}}{T_{s}} + \frac{K_{d} e_{k-2}}{2 T_{s}}\]

Comparison with Simulink#

The following figure will be our general PID subsytem, where we can implement the discrete PID algorithm inside the MATLAB function block.

The MATLAB function below has already been equipped with an integral windup prevention (clamping method):

function [u0, u] = PID(u1, u2, e0, e1, e2, params)
% u1 is u(k-1)
% u2 is u(k-2)
% e0 is e(k)
% e1 is e(k-1)
% e2 is e(k-2)

Kp = params(1);
Ki = params(2);  
Kd = params(3); 
N  = params(4); 
UB = params(5); % lower-bound value for the saturated output
LB = params(6); % upper-bound value for the saturated output

Ts  = params(7);
NTs = N * Ts;

% PID
E0 = Kd/(2*Ts) + Ki*Ts/8 + Ki/(4*N) + Kp/4 + Kp/(2*NTs);
E1 = -Kd/Ts + Ki*Ts/4 - Kp/NTs;
E2 = Kd/(2*Ts) + Ki*Ts/8 - Ki/(4*N) - Kp/4 + Kp/(2*NTs);
lhs = (-1/NTs) * u1 + (-1/4 + 1/(2*NTs))*u2;
rhs = E0 * e0 + E1 * e1 + E2 * e2;
u_pid = (rhs - lhs) / (1/4 + 1/(2*NTs));

% With disabled integrator, Ki = 0
E0 = Kd/(2*Ts) + Kp/4 + Kp/(2*NTs);
E1 = -Kd/Ts - Kp/NTs;
E2 = Kd/(2*Ts) - Kp/4 + Kp/(2*NTs);
rhs = E0 * e0 + E1 * e1 + E2 * e2;
u_pd = (rhs - lhs) / (1/4 + 1/(2*NTs));

if u_pid > UB
    u  = UB;
    u0 = u_pd;
elseif u_pid < LB
    u  = LB;
    u0 = u_pd;
else
    u  = u_pid;
    u0 = u_pid;
end

For testing, we will use the following Simulink program:

As shown in the figure, as our plant, we use a first order system: \(\frac{5/255}{2s+1}\). Identical control parameters are used for the testing.

The result is as follows:

Downloads#

Simulink file for bilinear PID

MATLAB version: R2024b