## TCLab Stability Analysis

Objective: Determine closed-loop behavior of a P-only controller with the TCLab including oscillatory behavior and stability. Background information on stability analysis may be useful for completing this exercise.

For the temperature control lab, determine a second-order model between heater 1 (Q1) and temperature sensor 1 (TC1).

$$G_p(s) = \frac{T_{C1}(s)}{Q_1(s)} = \frac{K_p}{\tau_s^2\,s^2+2\,\zeta\,\tau_s\,s+1}$$

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from gekko import GEKKO
import tclab
import time

# Import data
try:
# try to read local data file first
filename = 'data.csv'
except:
# heater steps
Q1d = np.zeros(601)
Q1d[10:200] = 80
Q1d[200:280] = 20
Q1d[280:400] = 70
Q1d[400:] = 50

Q2d = np.zeros(601)

try:
# Connect to Arduino
a = tclab.TCLab()
fid = open(filename,'w')
fid.write('Time,Q1,Q2,T1,T2\n')
fid.close()

# run step test (10 min)
for i in range(601):
# set heater values
a.Q1(Q1d[i])
a.Q2(Q2d[i])
print('Time: ' + str(i) + \
' Q1: ' + str(Q1d[i]) + \
' Q2: ' + str(Q2d[i]) + \
' T1: ' + str(a.T1)   + \
' T2: ' + str(a.T2))
# wait 1 second
time.sleep(1)
fid = open(filename,'a')
fid.write(str(i)+','+str(Q1d[i])+','+str(Q2d[i])+',' \
+str(a.T1)+','+str(a.T2)+'\n')
fid.close()
# close connection to Arduino
a.close()
except:

# Second order model of TCLab
m = GEKKO()
m.time = data['Time'].values
Kp   = m.FV(1.0,lb=0.5,ub=2.0)
taus = m.FV(50,lb=10,ub=200)
zeta =  m.FV(1.2,lb=1.1,ub=5)
T0 = data['T1'][0]
Q1 = m.Param(0)
x = m.Var(0); TC1 = m.CV(T0)
m.Equation(x==TC1.dt())
m.Equation((taus**2)*x.dt()+2*zeta*taus*TC1.dt()+(TC1-T0) == Kp*Q1)
m.options.IMODE = 5
m.options.NODES = 2
m.options.EV_TYPE = 2 # Objective type
Kp.STATUS = 1
taus.STATUS = 1
zeta.STATUS = 1
TC1.FSTATUS = 1
Q1.value=data['Q1'].values
TC1.value=data['T1'].values
m.solve(disp=True)

# Parameter values
print('Estimated Parameters')
print('Kp  : ' + str(Kp.value[0]))
print('taus: ' + str(taus.value[0]))
print('zeta: ' + str(zeta.value[0]))

# Create plot
plt.figure(figsize=(10,7))
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(data['Time'],data['T1'],'b.',label=r'$T_1$ measured')
plt.plot(m.time,TC1.value,color='orange',linestyle='--',\
linewidth=2,label=r'$T_1$ second order')
plt.ylabel(r'T ($^oC$)')
plt.legend(loc=2)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(data['Time'],data['Q1'],'r-',\
linewidth=3,label=r'$Q_1$')
plt.ylabel('Heater (%)')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.savefig('tclab_2nd_order.png')
plt.show()

With the second-order TCLab model, determine the range of controller gains (Kc) where a P-only controller oscillates and is stable (does not diverge). Traditional stability analysis does not include information about controller output saturation (0-100% heater).

P-Only Simulator

Use the P-Only simulator to test the stability and oscillation predictions. The simulator shows the stability and oscillation of a P-only controller with a TCLab second-order model with K_p=0.8473 oC/%, \tau_s=51.08 sec, \zeta=1.581 sec, and \theta_p=0.0 sec. Use the second-order model parameters from your own TCLab device for a more accurate simulation. The controller gain K_c is adjusted with a slider to compute the updated temperature response.

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 Second-Order
Kp = 0.8473
taus = 51.08
zeta = 1.581
thetap = 0.0

def process(z,t,u):
x,y = z
dxdt = (1.0/(taus**2)) * (-2.0*zeta*taus*x-(y-23.0) + Kp * u)
dydt = x
return [dxdt,dydt]

def pidPlot(Kc):
t = np.linspace(0,tf,n) # create time vector
P = np.zeros(n)         # initialize proportional 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:] = 60.0          # step up
z0 = [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
PV[i] = z0[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
OP[i] = min(100,max(0,P[i])) # calculate new controller output

P = np.zeros(n)         # initialize proportional term
e = np.zeros(n)         # initialize error
OPu = np.zeros(n)       # initialize controller output
PVu = np.ones(n)*23.0   # initialize process variable
SP = np.ones(n)*23.0    # initialize setpoint
SP[10:] = 60.0          # step up
z0 = [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=(OPu[max(0,i-1-int(thetap))],))
z0 = z[1]                  # record new initial condition
# calculate new OP with PID
PVu[i] = z0[1]             # record PV
e[i] = SP[i] - PVu[i]       # calculate error = SP - PV
dt = t[i] - t[i-1]         # calculate time step
P[i] = Kc * e[i]           # calculate proportional term
OPu[i] = P[i]               # calculate new controller output

# plot PID response
plt.figure(1,figsize=(15,5))
plt.subplot(1,2,1)
plt.plot(t,SP,'k-',linewidth=2,label='Setpoint (SP)')
plt.plot(t,PV,'r-',linewidth=2,label='Temperature - OP Limits (PV)')
plt.plot(t,PVu,'b--',linewidth=2,label='Temperature - No OP Limits (PV)')
plt.ylabel(r'T $(^oC)$')
plt.text(100,30,'OP Limit Offset: ' + str(np.round(SP[-1]-PVu[-1],2)))
M = SP[-1]-SP[0]
pred_offset = M*(1-Kp*Kc/(1+Kp*Kc))
plt.text(100,25,'No OP Limit Offset: ' + str(np.round(pred_offset,2)))
plt.text(400,30,r'$K_c$: ' + str(np.round(Kc,1)))
plt.legend(loc=1)
plt.xlabel('time (sec)')
plt.subplot(1,2,2)
plt.plot(t,OP,'r-',linewidth=2,label='Heater - OP Limits (OP)')
plt.plot(t,OPu,'b--',linewidth=2,label='Heater - No OP Limits (OP)')
plt.ylabel('Heater (%)')
plt.legend(loc='best')
plt.xlabel('time (sec)')

Kc_slide = wg.FloatSlider(value=2.0,min=-2.0,max=10.0,step=0.1)
wg.interact(pidPlot, Kc=Kc_slide)
print('P-only Simulator with and without OP Limits: Adjust Kc')

Solution

The second-order parameter regression gives the following parameters:

$$K_p= 0.8473$$

$$\tau_s= 51.08$$

$$\zeta = 1.581$$

$$G_p(s)=\frac{K_p}{\tau_s^2\,s^2+2\,\zeta\,\tau_s\,s+1}$$

The controller stability range is analyzed with a Routh Array and Root Locus plot. The Root Locus plot also shows where the controller begins to oscillate.

Routh Array

The closed-loop transfer function with a P-only controller gain K_c is:

$$G_{cl}(s)=\frac{T_{C1}(s)}{T_{SP}(s)}=\frac{K_c\,G_p(s)}{1+K_c\,G_p(s)}=\frac{\frac{K_c\,K_p}{\tau_s^2\,s^2+2\,\zeta\,\tau_s\,s+1}}{1+\frac{K_c\,K_p}{\tau_s^2\,s^2+2\,\zeta\,\tau_s\,s+1}}=\frac{K_c\,K_p}{\tau_s^2\,s^2+2\,\zeta\,\tau_s\,s+1+K_c\,K_p}$$

The coefficients of the denominator are analyzed to determine the stability range.

 a_n=\tau_s^2 a_{n-2}=1+K_c\,K_p a_{n-1}=2 \zeta \tau_s a_{n-3}=0 b_{1}=\frac{(2 \zeta \tau_s)(1+K_c K_p)}{2 \zeta \tau_s}=1+K_c K_p b_{2}=0 c_{1}=0 0

The leading edge cannot change signs for the system to be stable. Therefore, the following conditions must be met:

$$a_n=\tau_s^2=2609 > 0$$

$$a_{n-1}=2\,\zeta\,\tau_s=161.4 > 0$$

$$b_{1}=1+K_c\,K_p > 0$$

$$c_{1}=0 > 0$$

The positive constraint on b_1 leads to K_c>-\frac{1}{K_p}. Therefore the following range is acceptable for the controller stability.

$$K_c > -1.1803$$

Root Locus Plot

One of the real portion of the roots becomes positive when K_c<-1.19. There is no positive value of the controller gain that leads to an unstable system. This implies that the only way to make a second-order system unstable is to incorrectly specify the controller as a Direct acting controller with a value at least equal to the inverse of the process gain. The unstable controller does not have imaginary roots so the response is expect to go to +\infty or -\infty without oscillation.

There are imaginary roots and oscillation when K_c>1.78. The real portion of the roots are negative so the system is stable and converges but not necessarily to the new set point.

# Stability Analysis
import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button, RadioButtons

Kp   = 0.84726
taus = 51.0800
zeta = 1.58068

# open loop
num = [Kp]
den = [taus**2,2*zeta*taus,1]
# demoninator of the closed loop
# Gcl = direct/(1+loop)
def dcl(K):
return [taus**2,2*zeta*taus,1.0+K*Kp]

# root locus plot
n = 10000 # number of points to plot
nr = len(den)-1 # number of roots
rs = np.zeros((n,2*nr))   # store results
Kc1 = -5.0
Kc2 = 5.0
Kc = np.linspace(Kc1,Kc2,n)  # Kc values
for i in range(n):        # cycle through n times
roots = np.roots(dcl(Kc[i]))
for j in range(nr):   # store roots
rs[i,j] = roots[j].real # store real
rs[i,j+nr] = roots[j].imag # store imaginary

# create the image
fig, ax = plt.subplots()
ls = []
for i in range(nr):
plt.plot(rs[:,i],rs[:,i+nr],'r.',markersize=2)
# this handle is required to update the plot
if math.isclose(rs[0,i+nr],0.0):
lbl = f'{rs[0,i]:0.4f}'
else:
lbl = f'{rs[0,i]:0.4f}, {rs[0,i+nr]:0.4f}i'
l, = plt.plot(rs[0,i],rs[0,i+nr], 'ks', markersize=5,label=lbl)
ls.append(l)
leg = plt.legend(loc='best')
plt.xlabel('Root (real)')
plt.ylabel('Root (imag)')
plt.grid(b=True, which='major', color='b', linestyle='-',alpha=0.5)
plt.grid(b=True, which='minor', color='r', linestyle='--',alpha=0.5)

# slider creation
axcolor = 'lightgoldenrodyellow'
axKc = plt.axes([0.10, 0.1, 0.80, 0.03], facecolor=axcolor)
sKc = Slider(axKc, 'Kc', Kc1, Kc2, valinit=0, valstep=0.01)

def update(val):
Kc_val= sKc.val
indx = (np.abs(Kc-Kc_val)).argmin()
for i in range(nr):
ls[i].set_ydata(rs[indx,i+nr])
ls[i].set_xdata(rs[indx,i])
if math.isclose(rs[indx,i+nr],0.0):
lbl = f'{rs[indx,i]:0.4f}'
else:
lbl = f'{rs[indx,i]:0.4f}, {rs[indx,i+nr]:0.4f}i'
leg.texts[i].set_text(lbl)
fig.canvas.draw_idle()
sKc.on_changed(update)

resetax = plt.axes([0.8, 0.025, 0.1, 0.04])
button = Button(resetax, 'Reset', color=axcolor, hovercolor='0.975')

def reset(event):
sKc.reset()
button.on_clicked(reset)

plt.show()

Extra: Calculate Offset

If there were a set point change of +M with T_{SP}(s)=\frac{M}{s} then the temperature response is predicted to increase:

$$T_{C1,\infty} = \lim_{s \to 0} s \, T_{C1}(s) = \lim_{s \to 0} s \, G_{cl}(s)\frac{M}{s} = \frac{M\,K_p\,K_c}{1+K_p\,K_c}$$

This gives an offset of:

$$\mathrm{Offset} = T_{SP} - T_{C1,\infty} = M \left(1 - \frac{K_p\,K_c}{1+K_p\,K_c} \right)$$