PID controller implemented digitally

I've been reading articles about digital PID implementation. So, in the discrete domain, we can write the parallel form PID law, as:

$$U(z) = E(z)*\left[K_p + \frac{K_iT_s}{1 - z^{-1}} + K_d\frac{1 - z^{-1}}{T_s}\right]$$

Where $T_s$ is the sampling time. It's also said that $K_p$, $K_i$ and $K_d$ are related as:

$$K_i=\frac{K_p}{T_i}\\ K_d = K_p\cdot T_d$$

where $T_i$ is the integral time and $T_d$ is the derivative time.

So, in the above equation I would get $$K_i=K_p\frac{T_s}{T_i} \\ K_d=K_p\frac{T_d}{T_s}$$

Let's say I implement my control law, in discrete time (k), as:

$$u[k] = K_p\cdot e[k] + K_i\big(e[k] + e[k-1]\big) + K_d\big(e[k] - e[k-1]\big)$$

My questions:

• is the integral time and derivative time ($T_s$ and $T_d$) defined by the clock period of the circuit that perform these computations?

• if my error signal is sampled by the same clock that computes the integral and derivative parts of my law, $T_i=T_d=T_s$?

• If $T_i=T_d=T_s$, $K_p = K_i = K_d$? It sounds really odd.

• There is no relationship between the sampling increment, Ts, and the PID gains, Kp, Ki, and Kd. The PID gains are independent of the platform (would be the same numerical values if implemented using op-amps, for example). Ts is used to calculate the instantaneous values of the derivative and integral of the error signal. Your expression for the integral term in the difference equation is not correct, btw. – Chu Oct 31 '15 at 22:10
• Thank you. What about Ti and Ts? How they relate to system clock? Look at eq.5 of this paper: scilab.ninja/doc/b4/pid.pdf. If Td=Ti, Kp=Ki=Kd – user3845287 Oct 31 '15 at 23:28
• What do you mean by 'system clock'? – Chu Nov 1 '15 at 0:41
• The signal that clocks all the blocks of my closed loop-system., supposing it's the same. – user3845287 Nov 1 '15 at 16:46
• So, is the system clock the device that determines Ts? In any case, neither of these are directly related to Tp, Ti, or Td. – Chu Nov 1 '15 at 17:17

is the integral time and derivative time (Ts and Td) defined by the clock period of the circuit that perform these computations?

The integral time and the derivative time are related to any associated strobes used around those block.

Digital electronics take a finite time to perform the required operations. Take a trapezoid discrete integrator who's difference equation looks like this:

$y_n = y_{n-1} + K_i*T_s*\frac{x_n + x_{n-1}}{2}$

There is two additions, two multiplications and one division (although a cheap shift right division) for one integration step. This must be completed before the next sample $x_n$ arrives. Now to save on a multiplication that $K_i * T_s$ could be combined into one scalar, which is standard practice & this results on the gain variable being related to the sample time $T_s$

What sets your $T_s$ ? . Those multiplications and additions take a finite time & they MUST be completed by the time the next sample is available $x_n$

if my error signal is sampled by the same clock that computes the integral and derivative parts of my law, Ti=Td=Ts?

Yes it would be, if Ts is the data acquisition or interpolator rate, of one is used)

If Ti=Td=Ts, Kp = Ki = Kd? It sounds really odd.

No because Ts is only part of the discrete integral equation as there is then the gain factor. As previously mentioned the gain factor AND the sample period can be pre-calculated to save on a multiplication step

--edit--

To help clarify. Below is a piece of python code that implements a discrete integrator (choice of 3, but Trap more often than not is the one you want). The Sample time, $T_s$ here is 1us (as well as 100us to prove a point). In this instance $T_i = T_s = 1us (or 100us)$. Likewise the "clock" of the processor, the python virtual machine clock speed, for simplicity, can be considered to be 2.9GHz. This is a reasonable example for a FPGA/uP running at a higher speed and the ADC's are sampling at a slower rate BUT a rate that is slow enough for the main processor to complete its computational steps.

#!/usr/bin/env python
import numpy as np
from matplotlib import pyplot as plt

def fwdEuler(X=0, dt=0):
'''y(n) = y(n-1) + [t(n)-t(n-1)]*x(n-1)'''
x[0] = X*Ki # apply gain prior to integration
y[0] = self.y[-1] + dt*x[-1]
x[-1] = float(x[0])
y[-1] = float(y[0])
return y[0]

def revEuler(X=0, dt=0):
'''y(n) = y(n-1) + [t(n)-t(n-1)]*x(n)'''
x[0] = X*Ki # apply gain prior to integration
y[0] = y[-1] + dt*x[0]
x[-1] = float(x[0])
y[-1] = float(y[0])
return y[0]

def Trap(X=0, dt=0):
'''y(n) = y(n-1) + K*[t(n)-t(n-1)]*[x(n)+x(n-1)]/2'''
x[0] = X*Ki # apply gain prior to integration
y[0] = y[-1] + dt*(x[0]+x[-1])/2
x[-1] = float(x[0])
y[-1] = float(y[0])
return y[0]

plt.hold(True)
plt.grid(True)
f,p = plt.subplots(4)
x = [0,0]
y = [0,0]
Ki = 1

t = np.arange(0,1,1e-6)  # 1million samples at 1us step
data = np.ones(t.size)*1 #
p[0].plot(t,data)
p[0].set_title('Simple straight line to integrate')

int_1 = np.zeros(t.size)
for i in range(t.size):
int_1[i] = Trap(data[i],t[1]-t[0])
p[1].set_title('Trap integration with Ki=1 and Ti=1u')
p[1].plot(t,int_1)

x = [0,0]
y = [0,0]
Ki = 1
t = np.arange(0,1,100e-6)  # 1million samples at 1us step
data = np.ones(t.size)*1 #
int_2 = np.zeros(t.size)
for i in range(t.size):
int_2[i] = Trap(data[i],t[1]-t[0])
p[2].set_title('Trap integration with Ki=1 and Ti=100u')
p[2].plot(t,int_2)

x = [0,0]
y = [0,0]
Ki = 2
t = np.arange(0,1,1e-6)  # 1million samples at 1us step
data = np.ones(t.size)*1 #
int_3 = np.zeros(t.size)
for i in range(t.size):
int_3[i] = Trap(data[i],t[1]-t[0])
p[3].set_title('Trap integration with Ki=2 and Ti=1u')
p[3].plot(t,int_3)

plt.show()


As you can see, for a valid integral difference algorith, the actual integral gain is unity and is thus can be ... tuned, via an additional gain, $K_i$. This is equally true for differential terms.

• JonRB, isn't the difference equation of the differential therm written as: Kd*(e[k] - e[k-1])/Ts ?? – user3845287 Nov 1 '15 at 16:56
• For a differentiator, not quite. I think you might be confusing differentiator with difference equation. I posted the difference equation for an integrator. The difference equation for a differentiator would be like: y[n] = Kd*(x[n] -x[n-1])/Ts – JonRB Nov 1 '15 at 17:52
• Yes thats what i mean... The difference equation for a differentiator also takes Ts into account. – user3845287 Nov 1 '15 at 17:57
• it is a difference equation. It isn't the greatest (just like a really poor discrete integrator/accumulator that is referenced alot...) . A better one would be: y[n] = -x[n]/16 + x[n-2] - x[n-4] +x[n-6]/16 . This has Kd and Td wrapped into the main sample period. – JonRB Nov 1 '15 at 20:13

Yes, the coefficients in integrating and derivating are subject to sampling time, but not the way you suggest. Ti, Td remain always the same as they are defined in countinous time - like an opamp PID regulator works identicaly the same compared to the discrete PID regulator with same Kp,Ti,Td. What realy cahnges are the weigts of an equaition that change with respect of sampling time.
A PID regualtor is looking very simple to implement, but it has also traps, one known is the saturation of integrator, therefore special algorithms like anti-windup are called. Integrator can be done with simpson formula, while diferentiator is also tricky since it has to be damped, therefore an aproximative diferentiator is used.
Regulators with parameter Kp, Ti, Td are European type, while the US uses Kp, ki, kd. They are different because in EU type the Kp is multiplied to integrator and differentiator as you noted ki=Kp/Ti, kd=Kp*Td and Ti,Td rapresent time constants [s] and can be calculated with methods like Ziegler-Nichols, while ki,kd have no units.

Here you can find source code of implemented controllers as well to the book.

//TO...sample time
// at initial time all: ek1, ek2, ek3, e4, uk1 are zero

q0 := Kp*(1 + T0/Ti + Td/(6*T0));
q1 := Kp*(-1 + Td/(3*T0));
q2 := Kp*(-Td/T0);
q3 := Kp*(Td/(3*T0));
q4 := Kp*(Td/(6*T0));

ek4 := ek3
ek3 := ek2;
ek2 := ek1;
ek1 := ek;
ek := setpoint - actual;

u := q0*ek + q1*ek1 + q2*ek2 + q3*ek3 + q4*ek4 + uk1;

if u > u_max then
u := u_max;
elsif u < u_min then
u := u_min;
end_if;

uk1 := u;


This is Ziegler-Nichols regulaltor (Zn3fpd.m) Controller is based on forward rectangular method of discretization replacing derivation by a four-point difference from link, rearanged by me in different language. I don't give any guarantee.

This a code of zn2fr.m, similarly with your link, but with sampling time taken into account, while the algorithm in the link you posted does not provide the same bahaviour whit different sampling times. The suitable sampling time refered by the author is to be 5 to 20 times the sytem rise time.

q0 = Kp*(1+Td/T0);
q1 = -Kp*(1-T0/Ti+2*Td/T0);
q2 = Kp*(Td/T0);

ek2 := ek1;
ek1 := ek;
ek := setpoint - actual;

u = q0*ek + q1*ek1 + q2*ek2 + uk1;

if u > u_max then
u := u_max;
elsif u < u_min then
u := u_min;
end_if;

uk1 := u;


% zn3fd % Controller for 3rd order processes with filtration of D-component using % Tustin approximation. Time constant of the filter Tf = Td/alfa.

This is controller is widely used in industry, it has aproximated D-component (like real opamp leaking), alfa is in between 3...30, at 30 it will behave almost non filtered D, at 3 the D-component is almost useless (too filtered)

Tf = Td/alfa;
cf = Tf/T0;
ci = T0/Ti;
cd = Td/T0;
p1 = -4*cf/(1+2*cf);
p2 = (2*cf-1)/(1+2*cf);
q0 = Kp * (1 + 2*(cf+cd) + (ci/2)*(1+2*cf))/(1+2*cf);
q1 = Kp * (ci/2-4*(cf+cd))/(1+2*cf);
q2 = Kp * (cf*(2-ci) + 2*cd + ci/2 - 1)/(1+2*cf);

ek4 := ek3
ek3 := ek2;
ek2 := ek1;
ek1 := ek;
ek := setpoint - actual;

u = q0*ek + q1*ek1 + q2*ek2 - p1*uk1 - p2*uk2;

if u > u_max then
u := u_max;
elsif u < u_min then
u := u_min;
end_if;

uk2 := uk1;
uk1 := u;

• Incorrect. Ti, Td are not "EU" types. Ti is the integrator time and Td is the differentiator time. Equally the gains Kp, Ki and Kd are not initless... If the input is amps and the output is volts: Kp units=V/A. Likewise Ki units equals Vs/A and Kd units are V/As. – JonRB Nov 1 '15 at 17:50
• Ok, let they are parallel or serial or whatever, while at university they told us that more common type in US is parallel with kp,ki,kd, while we never adopted it since we used different topology kp,ti,td. en.wikipedia.org/wiki/PID_controller – Marko Buršič Nov 1 '15 at 18:01
• I'm just wondering if the relationship Ki=Kp.Ts/Ti and Kd=Kp.Td/Ts are true. If so, for a closed loop digital system, which all blocks are clocked by the same clock at Frequency = 1/Ts, so Ti = Td =Ts, and Kp=Ki=Kd. I wonder if that's true. If not, why? – user3845287 Nov 1 '15 at 18:35
• If Ti=1s, Kd=1s and Ts=1s then, and only then Kp=Ki=Kd, but I don't see any significance in that. – Marko Buršič Nov 1 '15 at 18:39
• Ok lets say my error is sampled by an ADC with a CLK signal 1/Ts. – user3845287 Nov 1 '15 at 19:12