TCLab PID Control Tuning
 |  |  |  |
|---|
Objective: Tune a discrete PID controller and test the performance with a series of setpoint changes over 10 min in the sequence from 23oC initially, 50oC at 10 sec, and 40oC at 300 sec. Modify the tuning parameters to achieve a low Integral Absolute Error between the measured temperature and the setpoint.
TCLab PID Control Simulator
A simulator is a useful tool to help evaluate changes in tuning before testing on a physical system. Use the PID simulator to find acceptable control performance that minimizes the integral absolute error between the setpoint and measured temperature.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import ipywidgets as wg
from IPython.display import display
n = 601 # time points to plot
tf = 600.0 # final time
# TCLab FOPDT
Kp = 0.9
taup = 175.0
thetap = 15.0
def process(y,t,u):
dydt = (1.0/taup) * (-(y-23.0) + Kp * u)
return dydt
def pidPlot(Kc,tauI,tauD):
t = np.linspace(0,tf,n) # create time vector
P = np.zeros(n) # initialize proportional term
I = np.zeros(n) # initialize integral term
D = np.zeros(n) # initialize derivative term
e = np.zeros(n) # initialize error
OP = np.zeros(n) # initialize controller output
PV = np.ones(n)*23.0 # initialize process variable
SP = np.ones(n)*23.0 # initialize setpoint
SP[10:300] = 50.0 # step up
SP[300:601] = 40.0 # step down
y0 = 23.0 # initial condition
iae = 0.0
# loop through all time steps
for i in range(1,n):
# simulate process for one time step
ts = [t[i-1],t[i]] # time interval
y = odeint(process,y0,ts,args=(OP[max(0,i-int(thetap))],))
y0 = y[1] # record new initial condition
iae += np.abs(SP[i]-y0[0])
# calculate new OP with PID
PV[i] = y[1] # record PV
e[i] = SP[i] - PV[i] # calculate error = SP - PV
dt = t[i] - t[i-1] # calculate time step
P[i] = Kc * e[i] # calculate proportional term
I[i] = I[i-1] + (Kc/tauI) * e[i] * dt # calculate integral term
D[i] = -Kc * tauD * (PV[i]-PV[i-1])/dt # calculate derivative term
OP[i] = P[i] + I[i] + D[i] # calculate new controller output
if OP[i]>=100:
OP[i] = 100.0
I[i] = I[i-1] # reset integral
if OP[i]<=0:
OP[i] = 0.0
I[i] = I[i-1] # reset integral
# plot PID response
plt.figure(1,figsize=(15,7))
plt.subplot(2,2,1)
plt.plot(t,SP,'k-',linewidth=2,label='Setpoint (SP)')
plt.plot(t,PV,'r:',linewidth=2,label='Temperature (PV)')
plt.ylabel(r'T $(^oC)$')
plt.text(100,30,'Integral Abs Error: ' + str(np.round(iae,2)))
plt.text(400,30,r'$K_c$: ' + str(np.round(Kc,0)))
plt.text(400,27,r'$\tau_I$: ' + str(np.round(tauI,0)) + ' sec')
plt.text(400,24,r'$\tau_D$: ' + str(np.round(tauD,0)) + ' sec')
plt.legend(loc='best')
plt.subplot(2,2,2)
plt.plot(t,P,'g.-',linewidth=2,label=r'Proportional = $K_c \; e(t)$')
plt.plot(t,I,'b-',linewidth=2,label=r'Integral = ' + \
r'$\frac{K_c}{\tau_I} \int_{i=0}^{n_t} e(t) \; dt $')
plt.plot(t,D,'r--',linewidth=2,label=r'Derivative = ' + \
r'$-K_c \tau_D \frac{d(PV)}{dt}$')
plt.legend(loc='best')
plt.subplot(2,2,3)
plt.plot(t,e,'m--',linewidth=2,label='Error (e=SP-PV)')
plt.ylabel(r'$\Delta T$ $(^oC)$')
plt.legend(loc='best')
plt.subplot(2,2,4)
plt.plot(t,OP,'b--',linewidth=2,label='Heater (OP)')
plt.legend(loc='best')
plt.xlabel('time')
Kc_slide = wg.FloatSlider(value=5.0,min=0.0,max=50.0,step=1.0)
tauI_slide = wg.FloatSlider(value=120.0,min=20.0,max=180.0,step=5.0)
tauD_slide = wg.FloatSlider(value=0.0,min=0.0,max=20.0,step=1.0)
wg.interact(pidPlot, Kc=Kc_slide, tauI=tauI_slide, tauD=tauD_slide)
print('PID Simulator: Adjust Kc, tauI, and tauD to achieve lowest Integral Abs Error')
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import ipywidgets as wg
from IPython.display import display
n = 601 # time points to plot
tf = 600.0 # final time
# TCLab FOPDT
Kp = 0.9
taup = 175.0
thetap = 15.0
def process(y,t,u):
dydt = (1.0/taup) * (-(y-23.0) + Kp * u)
return dydt
def pidPlot(Kc,tauI,tauD):
t = np.linspace(0,tf,n) # create time vector
P = np.zeros(n) # initialize proportional term
I = np.zeros(n) # initialize integral term
D = np.zeros(n) # initialize derivative term
e = np.zeros(n) # initialize error
OP = np.zeros(n) # initialize controller output
PV = np.ones(n)*23.0 # initialize process variable
SP = np.ones(n)*23.0 # initialize setpoint
SP[10:300] = 50.0 # step up
SP[300:601] = 40.0 # step down
y0 = 23.0 # initial condition
iae = 0.0
# loop through all time steps
for i in range(1,n):
# simulate process for one time step
ts = [t[i-1],t[i]] # time interval
y = odeint(process,y0,ts,args=(OP[max(0,i-int(thetap))],))
y0 = y[1] # record new initial condition
iae += np.abs(SP[i]-y0[0])
# calculate new OP with PID
PV[i] = y[1] # record PV
e[i] = SP[i] - PV[i] # calculate error = SP - PV
dt = t[i] - t[i-1] # calculate time step
P[i] = Kc * e[i] # calculate proportional term
I[i] = I[i-1] + (Kc/tauI) * e[i] * dt # calculate integral term
D[i] = -Kc * tauD * (PV[i]-PV[i-1])/dt # calculate derivative term
OP[i] = P[i] + I[i] + D[i] # calculate new controller output
if OP[i]>=100:
OP[i] = 100.0
I[i] = I[i-1] # reset integral
if OP[i]<=0:
OP[i] = 0.0
I[i] = I[i-1] # reset integral
# plot PID response
plt.figure(1,figsize=(15,7))
plt.subplot(2,2,1)
plt.plot(t,SP,'k-',linewidth=2,label='Setpoint (SP)')
plt.plot(t,PV,'r:',linewidth=2,label='Temperature (PV)')
plt.ylabel(r'T $(^oC)$')
plt.text(100,30,'Integral Abs Error: ' + str(np.round(iae,2)))
plt.text(400,30,r'$K_c$: ' + str(np.round(Kc,0)))
plt.text(400,27,r'$\tau_I$: ' + str(np.round(tauI,0)) + ' sec')
plt.text(400,24,r'$\tau_D$: ' + str(np.round(tauD,0)) + ' sec')
plt.legend(loc='best')
plt.subplot(2,2,2)
plt.plot(t,P,'g.-',linewidth=2,label=r'Proportional = $K_c \; e(t)$')
plt.plot(t,I,'b-',linewidth=2,label=r'Integral = ' + \
r'$\frac{K_c}{\tau_I} \int_{i=0}^{n_t} e(t) \; dt $')
plt.plot(t,D,'r--',linewidth=2,label=r'Derivative = ' + \
r'$-K_c \tau_D \frac{d(PV)}{dt}$')
plt.legend(loc='best')
plt.subplot(2,2,3)
plt.plot(t,e,'m--',linewidth=2,label='Error (e=SP-PV)')
plt.ylabel(r'$\Delta T$ $(^oC)$')
plt.legend(loc='best')
plt.subplot(2,2,4)
plt.plot(t,OP,'b--',linewidth=2,label='Heater (OP)')
plt.legend(loc='best')
plt.xlabel('time')
Kc_slide = wg.FloatSlider(value=5.0,min=0.0,max=50.0,step=1.0)
tauI_slide = wg.FloatSlider(value=120.0,min=20.0,max=180.0,step=5.0)
tauD_slide = wg.FloatSlider(value=0.0,min=0.0,max=20.0,step=1.0)
wg.interact(pidPlot, Kc=Kc_slide, tauI=tauI_slide, tauD=tauD_slide)
print('PID Simulator: Adjust Kc, tauI, and tauD to achieve lowest Integral Abs Error')
TCLab PID Control Test Script
The TCLab PID control validation script implements the best PID tuning values and calculates the integral absolute error with a sequence of setpoint changes from 23oC initially, 50oC at 10 sec, and 40oC at 300 sec. Update the PID tuning parameters before running the script.
# PID Parameters
Kc = 5.0
tauI = 120.0 # sec
tauD = 2.0 # sec
Kc = 5.0
tauI = 120.0 # sec
tauD = 2.0 # sec
The PID controller is embedded in the script and includes an animated plot that can be generated by running from the command line or from IDLE (Python.org install).
import tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# PID Parameters
Kc = 5.0
tauI = 120.0 # sec
tauD = 2.0 # sec
#-----------------------------------------
# PID Controller
#-----------------------------------------
# inputs ---------------------------------
# sp = setpoint
# pv = current temperature
# pv_last = prior temperature
# ierr = integral error
# dt = time increment between measurements
# outputs --------------------------------
# op = output of the PID controller
# P = proportional contribution
# I = integral contribution
# D = derivative contribution
def pid(sp,pv,pv_last,ierr,dt):
# Parameters in terms of PID coefficients
KP = Kc
KI = Kc/tauI
KD = Kc*tauD
# ubias for controller (initial heater)
op0 = 0
# upper and lower bounds on heater level
ophi = 100
oplo = 0
# calculate the error
error = sp-pv
# calculate the integral error
ierr = ierr + KI * error * dt
# calculate the measurement derivative
dpv = (pv - pv_last) / dt
# calculate the PID output
P = KP * error
I = ierr
D = -KD * dpv
op = op0 + P + I + D
# implement anti-reset windup
if op < oplo or op > ophi:
I = I - KI * error * dt
# clip output
op = max(oplo,min(ophi,op))
# return the controller output and PID terms
return [op,P,I,D]
# save txt file with data and set point
# t = time
# u1,u2 = heaters
# y1,y2 = tempeatures
# sp1,sp2 = setpoints
def save_txt(t, u1, u2, y1, y2, sp1, sp2):
data = np.vstack((t, u1, u2, y1, y2, sp1, sp2)) # vertical stack
data = data.T # transpose data
top = ('Time,Q1,Q2,T1,T2,TSP1,TSP2')
np.savetxt('validate.txt', data, delimiter=',',\
header=top, comments='')
# Connect to Arduino
a = tclab.TCLab()
# Turn LED on
print('LED On')
a.LED(100)
# Run time in minutes
run_time = 10.0
# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)
# Temperature
# set point (degC)
Tsp1 = np.ones(loops) * a.T1
# Heater set point steps
Tsp1[3:] = 50.0
Tsp1[300:] = 40.0
T1 = np.ones(loops) * a.T1 # measured T (degC)
error_sp = np.zeros(loops)
Tsp2 = np.ones(loops) * a.T2 # set point (degC)
T2 = np.ones(loops) * a.T2 # measured T (degC)
# impulse tests (0 - 100%)
Q1 = np.ones(loops) * 0.0
Q2 = np.ones(loops) * 0.0
print('Running Main Loop. Ctrl-C to end.')
print(' Time SP PV Q1 = P + I + D IAE')
print(('{:6.1f} {:6.2f} {:6.2f} ' + \
'{:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}').format( \
tm[0],Tsp1[0],T1[0], \
Q1[0],0.0,0.0,0.0,0.0))
# Main Loop
start_time = time.time()
prev_time = start_time
dt_error = 0.0
# Integral error
ierr = 0.0
# Integral absolute error
iae = 0.0
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
try:
for i in range(1,loops):
# Sleep time
sleep_max = 1.0
sleep = sleep_max - (time.time() - prev_time) - dt_error
if sleep>=1e-4:
time.sleep(sleep-1e-4)
else:
print('exceeded max cycle time by ' + str(abs(sleep)) + ' sec')
time.sleep(1e-4)
# Record time and change in time
t = time.time()
dt = t - prev_time
if (sleep>=1e-4):
dt_error = dt-1.0+0.009
else:
dt_error = 0.0
prev_time = t
tm[i] = t - start_time
# Read temperatures in Kelvin
T1[i] = a.T1
T2[i] = a.T2
# Integral absolute error
iae += np.abs(Tsp1[i]-T1[i])
# Calculate PID output
[Q1[i],P,ierr,D] = pid(Tsp1[i],T1[i],T1[i-1],ierr,dt)
# Write output (0-100)
a.Q1(Q1[i])
a.Q2(0.0)
# Print line of data
print(('{:6.1f} {:6.2f} {:6.2f} ' + \
'{:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}').format( \
tm[i],Tsp1[i],T1[i], \
Q1[i],P,ierr,D,iae))
# Update plot
plt.clf()
# Plot
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm[0:i],Tsp1[0:i],'k--',label=r'$T_1$ set point')
plt.plot(tm[0:i],T1[0:i],'r.',label=r'$T_1$ measured')
plt.ylabel(r'Temperature ($^oC$)')
plt.legend(loc=4)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1[0:i],'b-',label=r'$Q_1$')
plt.ylabel('Heater (%)')
plt.legend(loc=1)
plt.xlabel('Time (sec)')
plt.draw()
plt.pause(0.05)
# Turn off heaters
a.Q1(0)
a.Q2(0)
a.close()
# Save text file
save_txt(tm[0:i],Q1[0:i],Q2[0:i],T1[0:i],T2[0:i],Tsp1[0:i],Tsp2[0:i])
# Save figure
plt.savefig('PID_Control.png')
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
save_txt(tm[0:i],Q1[0:i],Q2[0:i],T1[0:i],T2[0:i],Tsp1[0:i],Tsp2[0:i])
plt.savefig('PID_Control.png')
# Make sure serial connection closes with an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Error: Shutting down')
a.close()
save_txt(tm[0:i],Q1[0:i],Q2[0:i],T1[0:i],T2[0:i],Tsp1[0:i],Tsp2[0:i])
plt.savefig('PID_Control.png')
raise
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# PID Parameters
Kc = 5.0
tauI = 120.0 # sec
tauD = 2.0 # sec
#-----------------------------------------
# PID Controller
#-----------------------------------------
# inputs ---------------------------------
# sp = setpoint
# pv = current temperature
# pv_last = prior temperature
# ierr = integral error
# dt = time increment between measurements
# outputs --------------------------------
# op = output of the PID controller
# P = proportional contribution
# I = integral contribution
# D = derivative contribution
def pid(sp,pv,pv_last,ierr,dt):
# Parameters in terms of PID coefficients
KP = Kc
KI = Kc/tauI
KD = Kc*tauD
# ubias for controller (initial heater)
op0 = 0
# upper and lower bounds on heater level
ophi = 100
oplo = 0
# calculate the error
error = sp-pv
# calculate the integral error
ierr = ierr + KI * error * dt
# calculate the measurement derivative
dpv = (pv - pv_last) / dt
# calculate the PID output
P = KP * error
I = ierr
D = -KD * dpv
op = op0 + P + I + D
# implement anti-reset windup
if op < oplo or op > ophi:
I = I - KI * error * dt
# clip output
op = max(oplo,min(ophi,op))
# return the controller output and PID terms
return [op,P,I,D]
# save txt file with data and set point
# t = time
# u1,u2 = heaters
# y1,y2 = tempeatures
# sp1,sp2 = setpoints
def save_txt(t, u1, u2, y1, y2, sp1, sp2):
data = np.vstack((t, u1, u2, y1, y2, sp1, sp2)) # vertical stack
data = data.T # transpose data
top = ('Time,Q1,Q2,T1,T2,TSP1,TSP2')
np.savetxt('validate.txt', data, delimiter=',',\
header=top, comments='')
# Connect to Arduino
a = tclab.TCLab()
# Turn LED on
print('LED On')
a.LED(100)
# Run time in minutes
run_time = 10.0
# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)
# Temperature
# set point (degC)
Tsp1 = np.ones(loops) * a.T1
# Heater set point steps
Tsp1[3:] = 50.0
Tsp1[300:] = 40.0
T1 = np.ones(loops) * a.T1 # measured T (degC)
error_sp = np.zeros(loops)
Tsp2 = np.ones(loops) * a.T2 # set point (degC)
T2 = np.ones(loops) * a.T2 # measured T (degC)
# impulse tests (0 - 100%)
Q1 = np.ones(loops) * 0.0
Q2 = np.ones(loops) * 0.0
print('Running Main Loop. Ctrl-C to end.')
print(' Time SP PV Q1 = P + I + D IAE')
print(('{:6.1f} {:6.2f} {:6.2f} ' + \
'{:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}').format( \
tm[0],Tsp1[0],T1[0], \
Q1[0],0.0,0.0,0.0,0.0))
# Main Loop
start_time = time.time()
prev_time = start_time
dt_error = 0.0
# Integral error
ierr = 0.0
# Integral absolute error
iae = 0.0
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
try:
for i in range(1,loops):
# Sleep time
sleep_max = 1.0
sleep = sleep_max - (time.time() - prev_time) - dt_error
if sleep>=1e-4:
time.sleep(sleep-1e-4)
else:
print('exceeded max cycle time by ' + str(abs(sleep)) + ' sec')
time.sleep(1e-4)
# Record time and change in time
t = time.time()
dt = t - prev_time
if (sleep>=1e-4):
dt_error = dt-1.0+0.009
else:
dt_error = 0.0
prev_time = t
tm[i] = t - start_time
# Read temperatures in Kelvin
T1[i] = a.T1
T2[i] = a.T2
# Integral absolute error
iae += np.abs(Tsp1[i]-T1[i])
# Calculate PID output
[Q1[i],P,ierr,D] = pid(Tsp1[i],T1[i],T1[i-1],ierr,dt)
# Write output (0-100)
a.Q1(Q1[i])
a.Q2(0.0)
# Print line of data
print(('{:6.1f} {:6.2f} {:6.2f} ' + \
'{:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}').format( \
tm[i],Tsp1[i],T1[i], \
Q1[i],P,ierr,D,iae))
# Update plot
plt.clf()
# Plot
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm[0:i],Tsp1[0:i],'k--',label=r'$T_1$ set point')
plt.plot(tm[0:i],T1[0:i],'r.',label=r'$T_1$ measured')
plt.ylabel(r'Temperature ($^oC$)')
plt.legend(loc=4)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1[0:i],'b-',label=r'$Q_1$')
plt.ylabel('Heater (%)')
plt.legend(loc=1)
plt.xlabel('Time (sec)')
plt.draw()
plt.pause(0.05)
# Turn off heaters
a.Q1(0)
a.Q2(0)
a.close()
# Save text file
save_txt(tm[0:i],Q1[0:i],Q2[0:i],T1[0:i],T2[0:i],Tsp1[0:i],Tsp2[0:i])
# Save figure
plt.savefig('PID_Control.png')
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
save_txt(tm[0:i],Q1[0:i],Q2[0:i],T1[0:i],T2[0:i],Tsp1[0:i],Tsp2[0:i])
plt.savefig('PID_Control.png')
# Make sure serial connection closes with an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Error: Shutting down')
a.close()
save_txt(tm[0:i],Q1[0:i],Q2[0:i],T1[0:i],T2[0:i],Tsp1[0:i],Tsp2[0:i])
plt.savefig('PID_Control.png')
raise
Solution