TCLab Cascade Control
 |  | |  |  |  |
|---|
Objective: Design, tune, and test a cascade controller for the temperature control lab.
Design a cascade controller by labeling a block diagram that includes the inner and outer control loops. The inner loop is a T1 controller that has Q1 as the controller output. The outer loop is a T2 controller that has T1,SP as the controller output.
Tune the cascade controller with the TCLab to meet a temperature 2 (T2) target set point by adjusting heater 1 (Q1). Minimize the integral absolute error (IAE) while staying below 85oC for temperature 1 (T1). Compare the cascade control performance with a standard PI controller.
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 = 1201 # time points to plot
tf = 1200.0 # final time
# TCLab Second-Order
Kp = 0.8473
Kd = 0.3
taus = 51.08
zeta = 1.581
thetap = 0.0
def process(z,t,u):
x1,y1,x2,y2 = z
dx1dt = (1.0/(taus**2)) * (-2.0*zeta*taus*x1-(y1-23.0) + Kp*u + Kd*(y2-y1))
dy1dt = x1
dx2dt = (1.0/(taus**2)) * (-2.0*zeta*taus*x2-(y2-23.0) + Kd*(y1-y2))
dy2dt = x2
return [dx1dt,dy1dt,dx2dt,dy2dt]
def pidPlot(Kc,tauI):
t = np.linspace(0,tf,n) # create time vector
P = np.zeros(n) # initialize proportional term
I = np.zeros(n) # initialize integral term
e = np.zeros(n) # initialize error
ie = np.zeros(n) # initialize integral error
OP = np.zeros(n) # initialize controller output
PV1 = np.ones(n)*23.0 # initialize process variable
PV2 = np.ones(n)*23.0 # initialize process variable
SP2 = np.ones(n)*23.0 # initialize setpoint
SP2[10:] = 35.0 # step up
z0 = [0,23.0,0,23.0] # initial condition
# 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
z = odeint(process,z0,ts,args=(OP[max(0,i-1-int(thetap))],))
z0 = z[1] # record new initial condition
# calculate new OP with PID
PV1[i] = z0[1] # record PV 1
PV2[i] = z0[3] # record PV 2
e[i] = SP2[i] - PV2[i] # calculate error = SP - PV
ie[i] = e[i] + ie[i-1]
dt = t[i] - t[i-1] # calculate time step
P[i] = Kc * e[i] # calculate proportional term
I[i] = Kc/tauI * ie[i]
OP[i] = min(100,max(0,P[i]+I[i])) # calculate new controller output
if OP[i]==100 or OP[i]==0:
ie[i] = ie[i-1] # anti-windup
# plot PID response
plt.figure(1,figsize=(15,5))
plt.subplot(1,2,1)
plt.plot(t,PV1,'r-',linewidth=2,label='Temperature 1')
plt.plot(t,SP2,'k-',linewidth=2,label='T2 Setpoint (SP)')
plt.plot(t,PV2,'b--',linewidth=2,label='Temperature 2')
plt.ylabel(r'T $(^oC)$')
plt.text(800,30,r'$K_c$: ' + str(np.round(Kc,1)))
plt.text(800,26,r'$\tau_I$: ' + str(np.round(tauI,1)))
plt.text(800,40,r'IAE: ' + str(np.round(np.sum(np.abs(e)),1)))
plt.legend(loc=1)
plt.xlabel('time (sec)')
plt.subplot(1,2,2)
plt.plot(t,OP,'r-',linewidth=2,label='Heater 1 (OP)')
plt.ylabel('Heater (%)')
plt.legend(loc='best')
plt.xlabel('time (sec)')
Kc_slide = wg.FloatSlider(value=2.0,min=1.0,max=10.0,step=1.0)
tauI_slide = wg.FloatSlider(value=150.0,min=5.0,max=300.0,step=5.0)
wg.interact(pidPlot, Kc=Kc_slide, tauI=tauI_slide)
print('PI Control: Adjust Kc and tauI')
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import ipywidgets as wg
from IPython.display import display
n = 1201 # time points to plot
tf = 1200.0 # final time
# TCLab Second-Order
Kp = 0.8473
Kd = 0.3
taus = 51.08
zeta = 1.581
thetap = 0.0
def process(z,t,u):
x1,y1,x2,y2 = z
dx1dt = (1.0/(taus**2)) * (-2.0*zeta*taus*x1-(y1-23.0) + Kp*u + Kd*(y2-y1))
dy1dt = x1
dx2dt = (1.0/(taus**2)) * (-2.0*zeta*taus*x2-(y2-23.0) + Kd*(y1-y2))
dy2dt = x2
return [dx1dt,dy1dt,dx2dt,dy2dt]
def pidPlot(Kc,tauI):
t = np.linspace(0,tf,n) # create time vector
P = np.zeros(n) # initialize proportional term
I = np.zeros(n) # initialize integral term
e = np.zeros(n) # initialize error
ie = np.zeros(n) # initialize integral error
OP = np.zeros(n) # initialize controller output
PV1 = np.ones(n)*23.0 # initialize process variable
PV2 = np.ones(n)*23.0 # initialize process variable
SP2 = np.ones(n)*23.0 # initialize setpoint
SP2[10:] = 35.0 # step up
z0 = [0,23.0,0,23.0] # initial condition
# 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
z = odeint(process,z0,ts,args=(OP[max(0,i-1-int(thetap))],))
z0 = z[1] # record new initial condition
# calculate new OP with PID
PV1[i] = z0[1] # record PV 1
PV2[i] = z0[3] # record PV 2
e[i] = SP2[i] - PV2[i] # calculate error = SP - PV
ie[i] = e[i] + ie[i-1]
dt = t[i] - t[i-1] # calculate time step
P[i] = Kc * e[i] # calculate proportional term
I[i] = Kc/tauI * ie[i]
OP[i] = min(100,max(0,P[i]+I[i])) # calculate new controller output
if OP[i]==100 or OP[i]==0:
ie[i] = ie[i-1] # anti-windup
# plot PID response
plt.figure(1,figsize=(15,5))
plt.subplot(1,2,1)
plt.plot(t,PV1,'r-',linewidth=2,label='Temperature 1')
plt.plot(t,SP2,'k-',linewidth=2,label='T2 Setpoint (SP)')
plt.plot(t,PV2,'b--',linewidth=2,label='Temperature 2')
plt.ylabel(r'T $(^oC)$')
plt.text(800,30,r'$K_c$: ' + str(np.round(Kc,1)))
plt.text(800,26,r'$\tau_I$: ' + str(np.round(tauI,1)))
plt.text(800,40,r'IAE: ' + str(np.round(np.sum(np.abs(e)),1)))
plt.legend(loc=1)
plt.xlabel('time (sec)')
plt.subplot(1,2,2)
plt.plot(t,OP,'r-',linewidth=2,label='Heater 1 (OP)')
plt.ylabel('Heater (%)')
plt.legend(loc='best')
plt.xlabel('time (sec)')
Kc_slide = wg.FloatSlider(value=2.0,min=1.0,max=10.0,step=1.0)
tauI_slide = wg.FloatSlider(value=150.0,min=5.0,max=300.0,step=5.0)
wg.interact(pidPlot, Kc=Kc_slide, tauI=tauI_slide)
print('PI Control: Adjust Kc and tauI')
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 = 1201 # time points to plot
tf = 1200.0 # final time
# TCLab Second-Order
Kp = 0.8473
Kd = 0.3
taus = 51.08
zeta = 1.581
thetap = 0.0
def process(z,t,u):
x1,y1,x2,y2 = z
dx1dt = (1.0/(taus**2)) * (-2.0*zeta*taus*x1-(y1-23.0) + Kp*u + Kd*(y2-y1))
dy1dt = x1
dx2dt = (1.0/(taus**2)) * (-2.0*zeta*taus*x2-(y2-23.0) + Kd*(y1-y2))
dy2dt = x2
return [dx1dt,dy1dt,dx2dt,dy2dt]
def pidPlot(Kc1,tauI1,Kc2,tauI2):
t = np.linspace(0,tf,n) # create time vector
P1 = np.zeros(n) # initialize proportional term
I1 = np.zeros(n) # initialize integral term
P2 = np.zeros(n) # initialize proportional term
I2 = np.zeros(n) # initialize integral term
e1 = np.zeros(n) # initialize error
ie1 = np.zeros(n) # initialize integral error
e2 = np.zeros(n) # initialize error
ie2 = np.zeros(n) # initialize integral error
OP = np.zeros(n) # initialize controller output
PV1 = np.ones(n)*23.0 # initialize process variable
PV2 = np.ones(n)*23.0 # initialize process variable
SP1 = np.ones(n)*23.0 # initialize setpoint
SP2 = np.ones(n)*23.0 # initialize setpoint
SP2[10:] = 35.0 # step up
z0 = [0,23.0,0,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
z = odeint(process,z0,ts,args=(OP[max(0,i-1-int(thetap))],))
z0 = z[1] # record new initial condition
dt = t[i] - t[i-1] # calculate time step
# outer loop
PV2[i] = z0[3] # record PV 2
e2[i] = SP2[i] - PV2[i] # calculate error = SP - PV
ie2[i] = e2[i]*dt + ie2[i-1]
P2[i] = Kc2 * e2[i] # calculate proportional term
I2[i] = Kc2/tauI2 * ie2[i]
SP1[i] = min(85,max(23,P2[i]+I2[i])) # calculate new controller output
if SP1[i]==85 or SP1[i]==23:
ie2[i] = ie2[i-1] # anti-windup
# inner loop
PV1[i] = z0[1] # record PV 1
e1[i] = SP1[i] - PV1[i] # calculate error = SP - PV
ie1[i] = e1[i]*dt + ie1[i-1]
P1[i] = Kc1 * e1[i] # calculate proportional term
I1[i] = Kc1/tauI1 * ie1[i]
OP[i] = min(100,max(0,P1[i]+I1[i])) # calculate new controller output
if OP[i]==100 or OP[i]==0:
ie1[i] = ie1[i-1] # anti-windup
# plot PID response
plt.figure(1,figsize=(15,5))
plt.subplot(1,2,1)
plt.plot(t,SP1,'k:',linewidth=2,label='T1 Setpoint (SP)')
plt.plot(t,PV1,'r-',linewidth=2,label='Temperature 1')
plt.plot(t,SP2,'k-',linewidth=2,label='T2 Setpoint (SP)')
plt.plot(t,PV2,'b--',linewidth=2,label='Temperature 2')
plt.ylabel(r'T $(^oC)$')
plt.text(500,30,'Inner Loop 1')
plt.text(500,26,r'$K_{c1}$: ' + str(np.round(Kc1,1)))
plt.text(500,22,r'$\tau_{I1}$: ' + str(np.round(tauI1,1)))
plt.text(800,30,'Outer Loop 2')
plt.text(800,26,r'$K_{c2}$: ' + str(np.round(Kc2,1)))
plt.text(800,22,r'$\tau_{I2}$: ' + str(np.round(tauI2,1)))
plt.text(800,40,r'IAE: ' + str(np.round(np.sum(np.abs(e2)),1)))
plt.legend(loc=1)
plt.xlabel('time (sec)')
plt.subplot(1,2,2)
plt.plot(t,OP,'r-',linewidth=2,label='Heater 1 (OP)')
plt.ylabel('Heater (%)')
plt.legend(loc='best')
plt.xlabel('time (sec)')
Kc1_slide = wg.FloatSlider(value=2.0,min=1.0,max=10.0,step=0.5)
tauI1_slide = wg.FloatSlider(value=200.0,min=5.0,max=300.0,step=5.0)
Kc2_slide = wg.FloatSlider(value=3.0,min=2.0,max=10.0,step=0.5)
tauI2_slide = wg.FloatSlider(value=150.0,min=5.0,max=300.0,step=5.0)
wg.interact(pidPlot, Kc1=Kc1_slide, tauI1=tauI1_slide, \
Kc2=Kc2_slide, tauI2=tauI2_slide)
print('PI Cascade Control: Adjust Kc1, Kc2, tauI1, tauI2')
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import ipywidgets as wg
from IPython.display import display
n = 1201 # time points to plot
tf = 1200.0 # final time
# TCLab Second-Order
Kp = 0.8473
Kd = 0.3
taus = 51.08
zeta = 1.581
thetap = 0.0
def process(z,t,u):
x1,y1,x2,y2 = z
dx1dt = (1.0/(taus**2)) * (-2.0*zeta*taus*x1-(y1-23.0) + Kp*u + Kd*(y2-y1))
dy1dt = x1
dx2dt = (1.0/(taus**2)) * (-2.0*zeta*taus*x2-(y2-23.0) + Kd*(y1-y2))
dy2dt = x2
return [dx1dt,dy1dt,dx2dt,dy2dt]
def pidPlot(Kc1,tauI1,Kc2,tauI2):
t = np.linspace(0,tf,n) # create time vector
P1 = np.zeros(n) # initialize proportional term
I1 = np.zeros(n) # initialize integral term
P2 = np.zeros(n) # initialize proportional term
I2 = np.zeros(n) # initialize integral term
e1 = np.zeros(n) # initialize error
ie1 = np.zeros(n) # initialize integral error
e2 = np.zeros(n) # initialize error
ie2 = np.zeros(n) # initialize integral error
OP = np.zeros(n) # initialize controller output
PV1 = np.ones(n)*23.0 # initialize process variable
PV2 = np.ones(n)*23.0 # initialize process variable
SP1 = np.ones(n)*23.0 # initialize setpoint
SP2 = np.ones(n)*23.0 # initialize setpoint
SP2[10:] = 35.0 # step up
z0 = [0,23.0,0,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
z = odeint(process,z0,ts,args=(OP[max(0,i-1-int(thetap))],))
z0 = z[1] # record new initial condition
dt = t[i] - t[i-1] # calculate time step
# outer loop
PV2[i] = z0[3] # record PV 2
e2[i] = SP2[i] - PV2[i] # calculate error = SP - PV
ie2[i] = e2[i]*dt + ie2[i-1]
P2[i] = Kc2 * e2[i] # calculate proportional term
I2[i] = Kc2/tauI2 * ie2[i]
SP1[i] = min(85,max(23,P2[i]+I2[i])) # calculate new controller output
if SP1[i]==85 or SP1[i]==23:
ie2[i] = ie2[i-1] # anti-windup
# inner loop
PV1[i] = z0[1] # record PV 1
e1[i] = SP1[i] - PV1[i] # calculate error = SP - PV
ie1[i] = e1[i]*dt + ie1[i-1]
P1[i] = Kc1 * e1[i] # calculate proportional term
I1[i] = Kc1/tauI1 * ie1[i]
OP[i] = min(100,max(0,P1[i]+I1[i])) # calculate new controller output
if OP[i]==100 or OP[i]==0:
ie1[i] = ie1[i-1] # anti-windup
# plot PID response
plt.figure(1,figsize=(15,5))
plt.subplot(1,2,1)
plt.plot(t,SP1,'k:',linewidth=2,label='T1 Setpoint (SP)')
plt.plot(t,PV1,'r-',linewidth=2,label='Temperature 1')
plt.plot(t,SP2,'k-',linewidth=2,label='T2 Setpoint (SP)')
plt.plot(t,PV2,'b--',linewidth=2,label='Temperature 2')
plt.ylabel(r'T $(^oC)$')
plt.text(500,30,'Inner Loop 1')
plt.text(500,26,r'$K_{c1}$: ' + str(np.round(Kc1,1)))
plt.text(500,22,r'$\tau_{I1}$: ' + str(np.round(tauI1,1)))
plt.text(800,30,'Outer Loop 2')
plt.text(800,26,r'$K_{c2}$: ' + str(np.round(Kc2,1)))
plt.text(800,22,r'$\tau_{I2}$: ' + str(np.round(tauI2,1)))
plt.text(800,40,r'IAE: ' + str(np.round(np.sum(np.abs(e2)),1)))
plt.legend(loc=1)
plt.xlabel('time (sec)')
plt.subplot(1,2,2)
plt.plot(t,OP,'r-',linewidth=2,label='Heater 1 (OP)')
plt.ylabel('Heater (%)')
plt.legend(loc='best')
plt.xlabel('time (sec)')
Kc1_slide = wg.FloatSlider(value=2.0,min=1.0,max=10.0,step=0.5)
tauI1_slide = wg.FloatSlider(value=200.0,min=5.0,max=300.0,step=5.0)
Kc2_slide = wg.FloatSlider(value=3.0,min=2.0,max=10.0,step=0.5)
tauI2_slide = wg.FloatSlider(value=150.0,min=5.0,max=300.0,step=5.0)
wg.interact(pidPlot, Kc1=Kc1_slide, tauI1=tauI1_slide, \
Kc2=Kc2_slide, tauI2=tauI2_slide)
print('PI Cascade Control: Adjust Kc1, Kc2, tauI1, tauI2')
Test the cascade control on the TCLab. Discuss any differences between the simulated and measured values.
import tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#-----------------------------------------
# PID controller performance for the TCLab
#-----------------------------------------
# PID Parameters
Kc1 = 5.0
tauI1 = 50.0 # sec
Kc2 = 3.0
tauI2 = 100.0 # sec
# Animate Plot?
animate = True
if animate:
try:
from IPython import get_ipython
from IPython.display import display,clear_output
get_ipython().run_line_magic('matplotlib', 'inline')
ipython = True
print('IPython Notebook')
except:
ipython = False
print('Not IPython Notebook')
#-----------------------------------------
# Cascade 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,d,cid):
# Parameters in terms of PID coefficients
if cid==1:
# controller 1
KP = Kc1
KI = Kc1/tauI1
ophi = 100
oplo = 0
else:
# controller 2
KP = Kc2
KI = Kc2/tauI2
ophi = 85
oplo = 23
# ubias for controller (initial heater)
op0 = 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
op = op0 + P + I
# 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]
# 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()
# Wait until temperature is below 25 degC
print('Check that temperatures are < 25 degC before starting')
i = 0
while a.T1>=25.0 or a.T2>=25.0:
print(f'Time: {i} T1: {a.T1} T2: {a.T2}')
i += 10
time.sleep(10)
# Turn LED on
print('LED On')
a.LED(100)
# Run time in minutes
run_time = 15.0
# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)
# Heater set point steps
Tsp1 = np.ones(loops) * a.T1
Tsp2 = np.ones(loops) * a.T2
Tsp2[10:] = 35.0 # step up
T1 = np.ones(loops) * a.T1 # measured T (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
if not animate:
print('Running Main Loop. Ctrl-C to end.')
print(' Time SP1 PV1 Q1 SP2 PV2 Q2 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],Tsp2[0],T2[0],Q2[0],0.0))
# Main Loop
start_time = time.time()
prev_time = start_time
dt_error = 0.0
# Integral error
ierr1 = 0.0
ierr2 = 0.0
# Integral absolute error
iae = 0.0
if not ipython:
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-sleep_max+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
# Disturbances
d1 = T1[i] - 23.0
d2 = T2[i] - 23.0
# Integral absolute error
iae += np.abs(np.abs(Tsp2[i]-T2[i]))
# Calculate PID output
[Tsp1[i],P,ierr2] = pid(Tsp2[i],T2[i],T2[i-1],ierr2,dt,d1,2)
[Q1[i],P,ierr1] = pid(Tsp1[i],T1[i],T1[i-1],ierr1,dt,d2,1)
# Write output (0-100)
a.Q1(Q1[i])
a.Q2(0)
if not animate:
# 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],Tsp2[i],T2[i],Q2[i],iae))
else:
if ipython:
plt.figure(figsize=(10,7))
else:
plt.clf()
# Update 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.plot(tm[0:i],Tsp2[0:i],'k-',label=r'$T_2$ set point')
plt.plot(tm[0:i],T2[0:i],'b.',label=r'$T_2$ measured')
plt.ylabel(r'Temperature ($^oC$)')
plt.text(0,65,'IAE: ' + str(np.round(iae,2)))
plt.legend(loc=4)
plt.ylim([15,100])
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1[0:i],'b-',label=r'$Q_1$')
plt.ylim([-10,110])
plt.ylabel('Heater (%)')
plt.legend(loc=1)
plt.xlabel('Time (sec)')
if ipython:
clear_output(wait=True)
display(plt.gcf())
else:
plt.draw()
plt.pause(0.05)
# Turn off heaters
a.Q1(0)
a.Q2(0)
a.close()
# Save text file
save_txt(tm,Q1,Q2,T1,T2,Tsp1,Tsp2)
# Save Plot
if not animate:
plt.figure(figsize=(10,7))
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm,Tsp1,'k--',label=r'$T_1$ set point')
plt.plot(tm,T1,'r.',label=r'$T_1$ measured')
plt.plot(tm,Tsp2,'k-',label=r'$T_2$ set point')
plt.plot(tm,T2,'b.',label=r'$T_2$ measured')
plt.ylabel(r'Temperature ($^oC$)')
plt.text(0,65,'IAE: ' + str(np.round(iae,2)))
plt.legend(loc=4)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm,Q1,'b-',label=r'$Q_1$')
plt.ylabel('Heater (%)')
plt.legend(loc=1)
plt.xlabel('Time (sec)')
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
print('Cascade control test complete')
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#-----------------------------------------
# PID controller performance for the TCLab
#-----------------------------------------
# PID Parameters
Kc1 = 5.0
tauI1 = 50.0 # sec
Kc2 = 3.0
tauI2 = 100.0 # sec
# Animate Plot?
animate = True
if animate:
try:
from IPython import get_ipython
from IPython.display import display,clear_output
get_ipython().run_line_magic('matplotlib', 'inline')
ipython = True
print('IPython Notebook')
except:
ipython = False
print('Not IPython Notebook')
#-----------------------------------------
# Cascade 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,d,cid):
# Parameters in terms of PID coefficients
if cid==1:
# controller 1
KP = Kc1
KI = Kc1/tauI1
ophi = 100
oplo = 0
else:
# controller 2
KP = Kc2
KI = Kc2/tauI2
ophi = 85
oplo = 23
# ubias for controller (initial heater)
op0 = 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
op = op0 + P + I
# 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]
# 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()
# Wait until temperature is below 25 degC
print('Check that temperatures are < 25 degC before starting')
i = 0
while a.T1>=25.0 or a.T2>=25.0:
print(f'Time: {i} T1: {a.T1} T2: {a.T2}')
i += 10
time.sleep(10)
# Turn LED on
print('LED On')
a.LED(100)
# Run time in minutes
run_time = 15.0
# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)
# Heater set point steps
Tsp1 = np.ones(loops) * a.T1
Tsp2 = np.ones(loops) * a.T2
Tsp2[10:] = 35.0 # step up
T1 = np.ones(loops) * a.T1 # measured T (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
if not animate:
print('Running Main Loop. Ctrl-C to end.')
print(' Time SP1 PV1 Q1 SP2 PV2 Q2 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],Tsp2[0],T2[0],Q2[0],0.0))
# Main Loop
start_time = time.time()
prev_time = start_time
dt_error = 0.0
# Integral error
ierr1 = 0.0
ierr2 = 0.0
# Integral absolute error
iae = 0.0
if not ipython:
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-sleep_max+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
# Disturbances
d1 = T1[i] - 23.0
d2 = T2[i] - 23.0
# Integral absolute error
iae += np.abs(np.abs(Tsp2[i]-T2[i]))
# Calculate PID output
[Tsp1[i],P,ierr2] = pid(Tsp2[i],T2[i],T2[i-1],ierr2,dt,d1,2)
[Q1[i],P,ierr1] = pid(Tsp1[i],T1[i],T1[i-1],ierr1,dt,d2,1)
# Write output (0-100)
a.Q1(Q1[i])
a.Q2(0)
if not animate:
# 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],Tsp2[i],T2[i],Q2[i],iae))
else:
if ipython:
plt.figure(figsize=(10,7))
else:
plt.clf()
# Update 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.plot(tm[0:i],Tsp2[0:i],'k-',label=r'$T_2$ set point')
plt.plot(tm[0:i],T2[0:i],'b.',label=r'$T_2$ measured')
plt.ylabel(r'Temperature ($^oC$)')
plt.text(0,65,'IAE: ' + str(np.round(iae,2)))
plt.legend(loc=4)
plt.ylim([15,100])
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1[0:i],'b-',label=r'$Q_1$')
plt.ylim([-10,110])
plt.ylabel('Heater (%)')
plt.legend(loc=1)
plt.xlabel('Time (sec)')
if ipython:
clear_output(wait=True)
display(plt.gcf())
else:
plt.draw()
plt.pause(0.05)
# Turn off heaters
a.Q1(0)
a.Q2(0)
a.close()
# Save text file
save_txt(tm,Q1,Q2,T1,T2,Tsp1,Tsp2)
# Save Plot
if not animate:
plt.figure(figsize=(10,7))
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm,Tsp1,'k--',label=r'$T_1$ set point')
plt.plot(tm,T1,'r.',label=r'$T_1$ measured')
plt.plot(tm,Tsp2,'k-',label=r'$T_2$ set point')
plt.plot(tm,T2,'b.',label=r'$T_2$ measured')
plt.ylabel(r'Temperature ($^oC$)')
plt.text(0,65,'IAE: ' + str(np.round(iae,2)))
plt.legend(loc=4)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm,Q1,'b-',label=r'$Q_1$')
plt.ylabel('Heater (%)')
plt.legend(loc=1)
plt.xlabel('Time (sec)')
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
print('Cascade control test complete')
Solution
Design: A standard feedback control system includes only one measurement (T2) and one actuator (Q1) that is adjusted to meet the setpoint (T2,SP).
By adding the second measurement (T1) and another PID controller with (T1,SP) as the inner loop, cascade control is implemented.
Tune the cascade controller involves minimizing the integral absolute error (IAE) while staying below 85oC for temperature 1 (T1). The controller output (T1,SP) from the outer controller (T2) is set to 85oC to avoid overshooting the upper limit constraint. The anticipated cascade performance is IAE = 2345.2.
Test: The 15 minute test gives a larger IAE by about 8% ((2537.8-2345.2)/2345.2) than predicted because of unmodeled disturbances that are not included in the simulation. This is evident from the Q1 adjustments even after T2 reaches the setpoint.
# PID Parameters
Kc1 = 6.5
tauI1 = 85.0 # sec
Kc2 = 7.5
tauI2 = 110.0 # sec
Kc1 = 6.5
tauI1 = 85.0 # sec
Kc2 = 7.5
tauI2 = 110.0 # sec
