## TCLab Second Order Regression     Objective: Develop an adaptive regression that recalculates parameters of a second order model as new data is measured. Discuss how adaptive parameter regression could be used to adapt a PID controller with gain scheduling.

Gain scheduling is practice of changing the PID control gain based on many factors including time, distance from setpoint, operating regime such as startup or shutdown, based on a process variable zone. Gain scheduling is may be used for nonlinear processes where the process gain or time constant is different at high or low values. Gain scheduling is one example of adaptive control that includes other strategies that sense and modify control behavior based on observed system response. The purpose of this exercise is to adaptively estimate a second-order model as new data arrives, not as a batch regression when all of the data is collected.

To solve the adaptive parameter regression with either Scipy.opitimize.minimize or Python Gekko, modify the function to include the analytic solution for an overdamped 2nd order response to a step input of 80%.

# 2nd order step response
def model(T0,t,M,Kp,taus,zeta):
# T0 = initial T
# t  = time
# M  = magnitude of the step
# Kp = gain
# taus = second order time constant
# zeta = damping factor (zeta>1 for overdamped)
T = ?  # fill in the analytic solution
return T

import tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import random

# Second order model of TCLab
# initial parameter guesses
Kp = 1.0
taus = 50.0
zeta = 1.2

# Detect session is IPython
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')

# magnitude of step
M = 80

# 2nd order step response
def model(T0,t,M,Kp,taus,zeta):
# T0 = initial T
# t  = time
# M  = magnitude of the step
# Kp = gain
# taus = second order time constant
# zeta = damping factor (zeta>1 for overdamped)
T = ?
return T

# define objective for optimizer
def objective(p,tm,ymeas):
# p = optimization parameters
Kp = p
taus = p
zeta = p
# tm = time points
# ymeas = measurements
# ypred = predicted values
n = np.size(tm)
ypred = np.ones(n)*ymeas
for i in range(1,n):
ypred[i] = model(ymeas,tm[i],M,Kp,taus,zeta)
sse = sum((ymeas-ypred)**2)
# penalize bound violation
if taus<10.0:
sse = sse + 100.0 * (10.0-taus)**2
if taus>200.0:
sse = sse + 100.0 * (200.0-taus)**2
if zeta<=1.1:
sse = sse + 1e6 * (1.0-zeta)**2
if zeta>=5.0:
sse = sse + 1e6 * (5.0-zeta)**2
if Kp>=2.0:
sse = sse + 1e6 * (2.0-Kp)**2
if Kp<=0.5:
sse = sse + 1e6 * (0.5-Kp)**2
return sse

# Connect to Arduino
a = tclab.TCLab()

# Get Version
print(a.version)

# Turn LED on
print('LED On')
a.LED(100)

# Run time in minutes
run_time = 5

# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)
z = np.zeros(loops)

# Temperature (K)
T1 = np.ones(loops) * a.T1 # measured T (degC)
T1p = np.ones(loops) * a.T1 # predicted T (degC)

# step test (0 - 100%)
Q1 = np.ones(loops) * 0.0
Q1[1:] = M # magnitude of the step

print('Running Main Loop. Ctrl-C to end.')
print('  Time   Kp    taus    zeta')
print('{:6.1f} {:6.2f} {:6.2f} {:6.2f}'.format(tm,Kp,taus,zeta))

# Create plot
if not ipython:
plt.figure(figsize=(10,7))
plt.ion()
plt.show()

# Main Loop
start_time = time.time()
prev_time = start_time
try:
for i in range(1,loops):
# Sleep time
sleep_max = 1.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep)
else:
time.sleep(0.01)

# Record time and change in time
t = time.time()
dt = t - prev_time
prev_time = t
tm[i] = t - start_time

T1[i] = a.T1

###############################
### CONTROLLER or ESTIMATOR ###
###############################
# Estimate parameters after 15 cycles and every 3 steps
if i>=15 and (np.mod(i,3)==0):
# randomize guess values
r = random.random()-0.5  # random number -0.5 to 0.5
Kp = Kp + r*0.05
taus = taus + r*1.0
zeta = zeta + r*0.01
p0=[Kp,taus,zeta]  # initial parameters
solution = minimize(objective,p0,args=(tm[0:i+1],T1[0:i+1]))
p = solution.x
Kp = p
taus = max(10.0,min(200.0,p))  # clip to >10, <=200
zeta = max(1.1,min(5.0,p)) # clip to >=1.1, <=5

# Update 2nd order prediction
for j in range(1,i+1):
T1p[j] = model(T1p,tm[j],M,Kp,taus,zeta)

# Write output (0-100)
a.Q1(Q1[i])

# Print line of data
if np.mod(i,15)==0:
print('  Time   Kp    taus    zeta')
print('{:6.1f} {:6.2f} {:6.2f} {:6.2f}'.format(tm[i],Kp,taus,zeta))

# Plot
if ipython:
plt.figure(figsize=(10,7))
else:
plt.clf()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm[0:i],T1[0:i],'ro',label=r'$T_1 \, Meas$')
plt.plot(tm[0:i],T1p[0:i],'k-',label=r'$T_1 \, Pred$')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1[0:i],'b-',label=r'$Q_1$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
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 figure
plt.savefig('test_Second_Order.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()
plt.savefig('test_Second_Order.png')

# Make sure serial connection still closes when there's an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Error: Shutting down')
a.close()
plt.savefig('test_Second_Order.png')
raise

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

# Detect session is IPython
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')

# magnitude of step
M = 80

# 2nd order step response
def model(T0,t,M,Kp,taus,zeta):
# T0 = initial T
# t  = time
# M  = magnitude of the step
# Kp = gain
# taus = second order time constant
# zeta = damping factor (zeta>1 for overdamped)
T = ?
return T

# Connect to Arduino
a = tclab.TCLab()

# Second order model of TCLab
m = GEKKO(remote=False)
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)
y0 = a.T1
u = m.MV(0)
x = m.Var(y0); y = m.CV(y0)
m.Equation(x==y.dt())
m.Equation((taus**2)*x.dt()+2*zeta*taus*y.dt()+(y-y0) == Kp*u)
m.options.IMODE = 5
m.options.NODES = 2
m.time = np.linspace(0,200,101)
m.solve(disp=False)
y.FSTATUS = 1

# Turn LED on
print('LED On')
a.LED(100)

# Run time in minutes
run_time = 5

# Number of cycles
loops = int(30.0*run_time)
tm = np.zeros(loops)
z = np.zeros(loops)

# Temperature (K)
T1 = np.ones(loops) * a.T1 # measured T (degC)
T1p = np.ones(loops) * a.T1 # predicted T (degC)

# step test (0 - 100%)
Q1 = np.ones(loops) * 0.0
Q1[1:] = M # magnitude of the step

print('Running Main Loop. Ctrl-C to end.')
print('  Time   Kp    taus    zeta')

# Create plot
if not ipython:
plt.figure(figsize=(10,7))
plt.ion()
plt.show()

# Main Loop
start_time = time.time()
prev_time = start_time
try:
for i in range(1,loops):
# Sleep time
sleep_max = 2.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep)
else:
time.sleep(0.01)

# Record time and change in time
t = time.time()
dt = t - prev_time
prev_time = t
tm[i] = t - start_time

T1[i] = a.T1

# Regression with Gekko
if i==5:
Kp.STATUS = 1
taus.STATUS = 1
zeta.STATUS = 1
u.meas = Q1[i]
y.meas = T1[i]
m.solve(disp=False)

# Update 2nd order prediction
Kpm = Kp.value; tausm = taus.value; zetam = zeta.value
for j in range(1,i+1):
T1p[j] = model(T1p,tm[j],M,Kpm,tausm,zetam)

# Write output (0-100)
a.Q1(Q1[i])

# Print line of data
if np.mod(i,15)==0:
print('  Time   Kp    taus    zeta')
print('{:6.1f} {:6.2f} {:6.2f} {:6.2f}'\
.format(tm[i],Kpm,tausm,zetam))

# Plot
if ipython:
plt.figure(figsize=(10,7))
else:
plt.clf()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm[0:i],T1[0:i],'ro',label=r'$T_1 \, Meas$')
plt.plot(tm[0:i],T1p[0:i],'k-',label=r'$T_1 \, Pred$')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1[0:i],'b-',label=r'$Q_1$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
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 figure
plt.savefig('test_Second_Order.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()
plt.savefig('test_Second_Order.png')

# Make sure serial connection still closes when there's an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Error: Shutting down')
a.close()
plt.savefig('test_Second_Order.png')
raise

Use a regression method to fit a 2nd order model to the closed-loop response by finding K_p, \zeta, and \tau_s. Dead-time \theta_p is not needed for this model.

$$\tau_s^2 \frac{d^2T_1}{dt^2} + 2 \zeta \tau_s \frac{dT_1}{dt} + T_1 = K_p \, Q_1$$

Insert the analytic solution of an overdamped zeta>1 second order step response with a heat step of 80% and run the code to observe the estimation of parameters. See Second Order Systems for additional information on analytic solutions.

Solution

The analytic solution of an overdamped zeta>1 second order step (80%) response is:

$$T_1(t) = 80 \, K_p \left( 1-e^{-\zeta\,t/\tau_s} \left[ \cosh\left( \frac{t}{\tau_s}\sqrt{\zeta^2 - 1} \right) + \frac{\zeta}{\sqrt{\zeta^2-1}} \sinh\left( \frac{t}{\tau_s}\sqrt{\zeta^2 - 1} \right) \right] \right) + T_{1,0}$$ # 2nd order step response
def model(y0,t,M,Kp,taus,zeta):
# y0 = initial y
# t  = time
# M  = magnitude of the step
# Kp = gain
# taus = second order time constant
# zeta = damping factor (zeta>1 for overdamped)
a = np.exp(-zeta*t/taus)
b = np.sqrt(zeta**2-1.0)
c = (t/taus)*b
y = Kp * M * (1.0 - a * (np.cosh(c)+(zeta/b)*np.sinh(c))) + y0
return y