TCLab C - Parameter Estimation
The TCLab is a hands-on application of machine learning and advanced temperature control with two heaters and two temperature sensors. The labs reinforce principles of model development, estimation, and advanced control methods. This is the third exercise and it involves estimating parameters in a multi-variate energy balance model. The predictions are aligned to the measured values through an optimizer that adjusts the parameters to minimize a sum of squared error objective. This lab builds upon the TCLab B by using the model of two heaters and two temperature sensors.
Lab Problem Statement
Data and Solutions
import numpy as np
import pandas as pd
import tclab
import time
import matplotlib.pyplot as plt
# generate step test data on Arduino
filename = 'tclab_dyn_data2.csv'
# 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)
Q2d[120:320] = 100
Q2d[320:520] = 10
Q2d[520:] = 80
# Connect to Arduino
a = tclab.TCLab()
fid = open(filename,'w')
fid.write('Time,H1,H2,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) + \
' H1: ' + str(Q1d[i]) + \
' H2: ' + 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')
# close connection to Arduino
a.close()
# read data file
data = pd.read_csv(filename)
# plot measurements
plt.figure()
plt.subplot(2,1,1)
plt.plot(data['Time'],data['H1'],'r-',label='Heater 1')
plt.plot(data['Time'],data['H2'],'b--',label='Heater 2')
plt.ylabel('Heater (%)')
plt.legend(loc='best')
plt.subplot(2,1,2)
plt.plot(data['Time'],data['T1'],'r.',label='Temperature 1')
plt.plot(data['Time'],data['T2'],'b.',label='Temperature 2')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
plt.xlabel('Time (sec)')
plt.savefig('tclab_dyn_meas2.png')
plt.show()
import pandas as pd
import tclab
import time
import matplotlib.pyplot as plt
# generate step test data on Arduino
filename = 'tclab_dyn_data2.csv'
# 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)
Q2d[120:320] = 100
Q2d[320:520] = 10
Q2d[520:] = 80
# Connect to Arduino
a = tclab.TCLab()
fid = open(filename,'w')
fid.write('Time,H1,H2,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) + \
' H1: ' + str(Q1d[i]) + \
' H2: ' + 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')
# close connection to Arduino
a.close()
# read data file
data = pd.read_csv(filename)
# plot measurements
plt.figure()
plt.subplot(2,1,1)
plt.plot(data['Time'],data['H1'],'r-',label='Heater 1')
plt.plot(data['Time'],data['H2'],'b--',label='Heater 2')
plt.ylabel('Heater (%)')
plt.legend(loc='best')
plt.subplot(2,1,2)
plt.plot(data['Time'],data['T1'],'r.',label='Temperature 1')
plt.plot(data['Time'],data['T2'],'b.',label='Temperature 2')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
plt.xlabel('Time (sec)')
plt.savefig('tclab_dyn_meas2.png')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from gekko import GEKKO
# Import or generate data
filename = 'tclab_dyn_data2.csv'
try:
data = pd.read_csv(filename)
except:
url = 'https://apmonitor.com/do/uploads/Main/tclab_dyn_data2.txt'
data = pd.read_csv(url)
# Create GEKKO Model
m = GEKKO()
m.time = data['Time'].values
# Parameters to Estimate
U = m.FV(value=10,lb=1,ub=20)
alpha1 = m.FV(value=0.01,lb=0.003,ub=0.03) # W / % heater
alpha2 = m.FV(value=0.005,lb=0.002,ub=0.02) # W / % heater
# STATUS=1 allows solver to adjust parameter
U.STATUS = 1
alpha1.STATUS = 1
alpha2.STATUS = 1
# Measured inputs
Q1 = m.MV(value=data['H1'].values)
Q2 = m.MV(value=data['H2'].values)
# State variables
TC1 = m.CV(value=data['T1'].values)
TC1.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
TC2 = m.CV(value=data['T2'].values)
TC2.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
Ta = m.Param(value=19.0+273.15) # K
mass = m.Param(value=4.0/1000.0) # kg
Cp = m.Param(value=0.5*1000.0) # J/kg-K
A = m.Param(value=10.0/100.0**2) # Area not between heaters in m^2
As = m.Param(value=2.0/100.0**2) # Area between heaters in m^2
eps = m.Param(value=0.9) # Emissivity
sigma = m.Const(5.67e-8) # Stefan-Boltzmann
# Heater temperatures in Kelvin
T1 = m.Intermediate(TC1+273.15)
T2 = m.Intermediate(TC2+273.15)
# Heat transfer between two heaters
Q_C12 = m.Intermediate(U*As*(T2-T1)) # Convective
Q_R12 = m.Intermediate(eps*sigma*As*(T2**4-T1**4)) # Radiative
# Semi-fundamental correlations (energy balances)
m.Equation(TC1.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ Q_C12 + Q_R12 \
+ alpha1*Q1))
m.Equation(TC2.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- Q_C12 - Q_R12 \
+ alpha2*Q2))
# Options
m.options.IMODE = 5 # MHE
m.options.EV_TYPE = 2 # Objective type
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
# Solve
m.solve(disp=True)
# Parameter values
print('U : ' + str(U.value[0]))
print('alpha1: ' + str(alpha1.value[0]))
print('alpha2: ' + str(alpha2.value[0]))
# Create plot
plt.figure()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(data['Time'],data['T1'],'ro',label=r'$T_1$ measured')
plt.plot(m.time,TC1.value,color='purple',linestyle='--',\
lw=3,label=r'$T_1$ predicted')
plt.plot(data['Time'],data['T2'],'bx',label=r'$T_2$ measured')
plt.plot(m.time,TC2.value,color='orange',linestyle='--',\
lw=3,label=r'$T_2$ predicted')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(data['Time'],data['H1'],'r-',\
lw=3,label=r'$Q_1$')
plt.plot(data['Time'],data['H2'],'b:',\
lw=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.show()
import matplotlib.pyplot as plt
import pandas as pd
from gekko import GEKKO
# Import or generate data
filename = 'tclab_dyn_data2.csv'
try:
data = pd.read_csv(filename)
except:
url = 'https://apmonitor.com/do/uploads/Main/tclab_dyn_data2.txt'
data = pd.read_csv(url)
# Create GEKKO Model
m = GEKKO()
m.time = data['Time'].values
# Parameters to Estimate
U = m.FV(value=10,lb=1,ub=20)
alpha1 = m.FV(value=0.01,lb=0.003,ub=0.03) # W / % heater
alpha2 = m.FV(value=0.005,lb=0.002,ub=0.02) # W / % heater
# STATUS=1 allows solver to adjust parameter
U.STATUS = 1
alpha1.STATUS = 1
alpha2.STATUS = 1
# Measured inputs
Q1 = m.MV(value=data['H1'].values)
Q2 = m.MV(value=data['H2'].values)
# State variables
TC1 = m.CV(value=data['T1'].values)
TC1.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
TC2 = m.CV(value=data['T2'].values)
TC2.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
Ta = m.Param(value=19.0+273.15) # K
mass = m.Param(value=4.0/1000.0) # kg
Cp = m.Param(value=0.5*1000.0) # J/kg-K
A = m.Param(value=10.0/100.0**2) # Area not between heaters in m^2
As = m.Param(value=2.0/100.0**2) # Area between heaters in m^2
eps = m.Param(value=0.9) # Emissivity
sigma = m.Const(5.67e-8) # Stefan-Boltzmann
# Heater temperatures in Kelvin
T1 = m.Intermediate(TC1+273.15)
T2 = m.Intermediate(TC2+273.15)
# Heat transfer between two heaters
Q_C12 = m.Intermediate(U*As*(T2-T1)) # Convective
Q_R12 = m.Intermediate(eps*sigma*As*(T2**4-T1**4)) # Radiative
# Semi-fundamental correlations (energy balances)
m.Equation(TC1.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ Q_C12 + Q_R12 \
+ alpha1*Q1))
m.Equation(TC2.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- Q_C12 - Q_R12 \
+ alpha2*Q2))
# Options
m.options.IMODE = 5 # MHE
m.options.EV_TYPE = 2 # Objective type
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
# Solve
m.solve(disp=True)
# Parameter values
print('U : ' + str(U.value[0]))
print('alpha1: ' + str(alpha1.value[0]))
print('alpha2: ' + str(alpha2.value[0]))
# Create plot
plt.figure()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(data['Time'],data['T1'],'ro',label=r'$T_1$ measured')
plt.plot(m.time,TC1.value,color='purple',linestyle='--',\
lw=3,label=r'$T_1$ predicted')
plt.plot(data['Time'],data['T2'],'bx',label=r'$T_2$ measured')
plt.plot(m.time,TC2.value,color='orange',linestyle='--',\
lw=3,label=r'$T_2$ predicted')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(data['Time'],data['H1'],'r-',\
lw=3,label=r'$Q_1$')
plt.plot(data['Time'],data['H2'],'b:',\
lw=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.show()
Note: Switch to make_mp4 = True to make an MP4 movie animation. This requires imageio and ffmpeg (install available through Python). It creates a folder named figures in your run directory. You can delete this folder after the run is complete.
import numpy as np
import time
import matplotlib.pyplot as plt
import random
# get gekko package with:
# pip install gekko
from gekko import GEKKO
# get tclab package with:
# pip install tclab
from tclab import TCLab
# Connect to Arduino
a = TCLab()
# Make an MP4 animation?
make_mp4 = False
if make_mp4:
import imageio # required to make animation
import os
try:
os.mkdir('./figures')
except:
pass
# Final time
tf = 10 # min
# number of data points (1 pt every 3 seconds)
n = tf * 20 + 1
# Configure heater levels
# Percent Heater (0-100%)
Q1s = np.zeros(n)
Q2s = np.zeros(n)
# Heater random steps every 50 sec
# Alternate steps by Q1 and Q2
Q1s[3:] = 100.0
Q1s[50:] = 0.0
Q1s[100:] = 80.0
Q2s[25:] = 60.0
Q2s[75:] = 100.0
Q2s[125:] = 25.0
# rapid, random changes every 5 cycles between 50 and 100
for i in range(130,180):
if i%10==0:
Q1s[i:i+10] = random.random() * 100.0
if (i+5)%10==0:
Q2s[i:i+10] = random.random() * 100.0
# Record initial temperatures (degC)
T1m = a.T1 * np.ones(n)
T2m = a.T2 * np.ones(n)
# Store MHE values for plots
Tmhe1 = T1m[0] * np.ones(n)
Tmhe2 = T2m[0] * np.ones(n)
Umhe = 10.0 * np.ones(n)
amhe1 = 0.01 * np.ones(n)
amhe2 = 0.0075 * np.ones(n)
#########################################################
# Initialize Model as Estimator
#########################################################
# Use remote=False for local solve (Windows, Linux, ARM)
# remote=True for remote solve (All platforms)
m = GEKKO(name='tclab-mhe',remote=False)
# 60 second time horizon, 20 steps
m.time = np.linspace(0,60,21)
# Parameters to Estimate
U = m.FV(value=10,name='u')
U.STATUS = 0 # don't estimate initially
U.FSTATUS = 0 # no measurements
U.DMAX = 1
U.LOWER = 5
U.UPPER = 15
alpha1 = m.FV(value=0.01,name='a1') # W / % heater
alpha1.STATUS = 0 # don't estimate initially
alpha1.FSTATUS = 0 # no measurements
alpha1.DMAX = 0.001
alpha1.LOWER = 0.003
alpha1.UPPER = 0.03
alpha2 = m.FV(value=0.0075,name='a2') # W / % heater
alpha2.STATUS = 0 # don't estimate initially
alpha2.FSTATUS = 0 # no measurements
alpha2.DMAX = 0.001
alpha2.LOWER = 0.002
alpha2.UPPER = 0.02
# Measured inputs
Q1 = m.MV(value=0,name='q1')
Q1.STATUS = 0 # don't estimate
Q1.FSTATUS = 1 # receive measurement
Q2 = m.MV(value=0,name='q2')
Q2.STATUS = 0 # don't estimate
Q2.FSTATUS = 1 # receive measurement
# Measurements for model alignment
TC1 = m.CV(value=T1m[0],name='tc1')
TC1.STATUS = 1 # minimize error between simulation and measurement
TC1.FSTATUS = 1 # receive measurement
TC1.MEAS_GAP = 0.1 # measurement deadband gap
TC1.LOWER = 0
TC1.UPPER = 200
TC2 = m.CV(value=T2m[0],name='tc2')
TC2.STATUS = 1 # minimize error between simulation and measurement
TC2.FSTATUS = 1 # receive measurement
TC2.MEAS_GAP = 0.1 # measurement deadband gap
TC2.LOWER = 0
TC2.UPPER = 200
Ta = m.Param(value=23.0+273.15) # K
mass = m.Param(value=4.0/1000.0) # kg
Cp = m.Param(value=0.5*1000.0) # J/kg-K
A = m.Param(value=10.0/100.0**2) # Area not between heaters in m^2
As = m.Param(value=2.0/100.0**2) # Area between heaters in m^2
eps = m.Param(value=0.9) # Emissivity
sigma = m.Const(5.67e-8) # Stefan-Boltzmann
# Heater temperatures
T1 = m.Intermediate(TC1+273.15)
T2 = m.Intermediate(TC2+273.15)
# Heat transfer between two heaters
Q_C12 = m.Intermediate(U*As*(T2-T1)) # Convective
Q_R12 = m.Intermediate(eps*sigma*As*(T2**4-T1**4)) # Radiative
# Semi-fundamental correlations (energy balances)
m.Equation(TC1.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ Q_C12 + Q_R12 \
+ alpha1*Q1))
m.Equation(TC2.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- Q_C12 - Q_R12 \
+ alpha2*Q2))
# Global Options
m.options.IMODE = 5 # MHE
m.options.EV_TYPE = 2 # Objective type
m.options.NODES = 3 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.options.COLDSTART = 1 # COLDSTART on first cycle
##################################################################
# Create plot
plt.figure(figsize=(12,7))
plt.ion()
plt.show()
# Main Loop
start_time = time.time()
prev_time = start_time
tm = np.zeros(n)
try:
for i in range(1,n):
# Sleep time
sleep_max = 3.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep-0.01)
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
# Read temperatures in Celsius
T1m[i] = a.T1
T2m[i] = a.T2
# Insert measurements
TC1.MEAS = T1m[i]
TC2.MEAS = T2m[i]
Q1.MEAS = Q1s[i-1]
Q2.MEAS = Q2s[i-1]
# Start estimating U after 10 cycles (20 sec)
if i==10:
U.STATUS = 1
alpha1.STATUS = 1
alpha2.STATUS = 1
# Predict Parameters and Temperatures with MHE
m.solve()
if m.options.APPSTATUS == 1:
# Retrieve new values
Tmhe1[i] = TC1.MODEL
Tmhe2[i] = TC2.MODEL
Umhe[i] = U.NEWVAL
amhe1[i] = alpha1.NEWVAL
amhe2[i] = alpha2.NEWVAL
else:
# Solution failed, copy prior solution
Tmhe1[i] = Tmhe1[i-1]
Tmhe2[i] = Tmhe1[i-1]
Umhe[i] = Umhe[i-1]
amhe1[i] = amhe1[i-1]
amhe2[i] = amhe2[i-1]
# Write new heater values (0-100)
a.Q1(Q1s[i])
a.Q2(Q2s[i])
# Plot
plt.clf()
ax=plt.subplot(3,1,1)
ax.grid()
plt.plot(tm[0:i],T1m[0:i],'ro',label=r'$T_1$ measured')
plt.plot(tm[0:i],Tmhe1[0:i],'k-',label=r'$T_1$ MHE')
plt.plot(tm[0:i],T2m[0:i],'bx',label=r'$T_2$ measured')
plt.plot(tm[0:i],Tmhe2[0:i],'k--',label=r'$T_2$ MHE')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(3,1,2)
ax.grid()
plt.plot(tm[0:i],Umhe[0:i],'k-',label='Heat Transfer Coeff')
plt.plot(tm[0:i],amhe1[0:i]*1000,'r--',label=r'$\alpha_1$x1000')
plt.plot(tm[0:i],amhe2[0:i]*1000,'b--',label=r'$\alpha_2$x1000')
plt.ylabel('Parameters')
plt.legend(loc='best')
ax=plt.subplot(3,1,3)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.draw()
plt.pause(0.05)
if make_mp4:
filename='./figures/plot_'+str(i+10000)+'.png'
plt.savefig(filename)
# Turn off heaters
a.Q1(0)
a.Q2(0)
# Save figure
plt.savefig('tclab_mhe.png')
# generate mp4 from png figures in batches of 350
if make_mp4:
images = []
iset = 0
for i in range(1,n):
filename='./figures/plot_'+str(i+10000)+'.png'
images.append(imageio.imread(filename))
if ((i+1)%350)==0:
imageio.mimsave('results_'+str(iset)+'.mp4', images)
iset += 1
images = []
if images!=[]:
imageio.mimsave('results_'+str(iset)+'.mp4', images)
# 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('tclab_mhe.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('tclab_mhe.png')
raise
import time
import matplotlib.pyplot as plt
import random
# get gekko package with:
# pip install gekko
from gekko import GEKKO
# get tclab package with:
# pip install tclab
from tclab import TCLab
# Connect to Arduino
a = TCLab()
# Make an MP4 animation?
make_mp4 = False
if make_mp4:
import imageio # required to make animation
import os
try:
os.mkdir('./figures')
except:
pass
# Final time
tf = 10 # min
# number of data points (1 pt every 3 seconds)
n = tf * 20 + 1
# Configure heater levels
# Percent Heater (0-100%)
Q1s = np.zeros(n)
Q2s = np.zeros(n)
# Heater random steps every 50 sec
# Alternate steps by Q1 and Q2
Q1s[3:] = 100.0
Q1s[50:] = 0.0
Q1s[100:] = 80.0
Q2s[25:] = 60.0
Q2s[75:] = 100.0
Q2s[125:] = 25.0
# rapid, random changes every 5 cycles between 50 and 100
for i in range(130,180):
if i%10==0:
Q1s[i:i+10] = random.random() * 100.0
if (i+5)%10==0:
Q2s[i:i+10] = random.random() * 100.0
# Record initial temperatures (degC)
T1m = a.T1 * np.ones(n)
T2m = a.T2 * np.ones(n)
# Store MHE values for plots
Tmhe1 = T1m[0] * np.ones(n)
Tmhe2 = T2m[0] * np.ones(n)
Umhe = 10.0 * np.ones(n)
amhe1 = 0.01 * np.ones(n)
amhe2 = 0.0075 * np.ones(n)
#########################################################
# Initialize Model as Estimator
#########################################################
# Use remote=False for local solve (Windows, Linux, ARM)
# remote=True for remote solve (All platforms)
m = GEKKO(name='tclab-mhe',remote=False)
# 60 second time horizon, 20 steps
m.time = np.linspace(0,60,21)
# Parameters to Estimate
U = m.FV(value=10,name='u')
U.STATUS = 0 # don't estimate initially
U.FSTATUS = 0 # no measurements
U.DMAX = 1
U.LOWER = 5
U.UPPER = 15
alpha1 = m.FV(value=0.01,name='a1') # W / % heater
alpha1.STATUS = 0 # don't estimate initially
alpha1.FSTATUS = 0 # no measurements
alpha1.DMAX = 0.001
alpha1.LOWER = 0.003
alpha1.UPPER = 0.03
alpha2 = m.FV(value=0.0075,name='a2') # W / % heater
alpha2.STATUS = 0 # don't estimate initially
alpha2.FSTATUS = 0 # no measurements
alpha2.DMAX = 0.001
alpha2.LOWER = 0.002
alpha2.UPPER = 0.02
# Measured inputs
Q1 = m.MV(value=0,name='q1')
Q1.STATUS = 0 # don't estimate
Q1.FSTATUS = 1 # receive measurement
Q2 = m.MV(value=0,name='q2')
Q2.STATUS = 0 # don't estimate
Q2.FSTATUS = 1 # receive measurement
# Measurements for model alignment
TC1 = m.CV(value=T1m[0],name='tc1')
TC1.STATUS = 1 # minimize error between simulation and measurement
TC1.FSTATUS = 1 # receive measurement
TC1.MEAS_GAP = 0.1 # measurement deadband gap
TC1.LOWER = 0
TC1.UPPER = 200
TC2 = m.CV(value=T2m[0],name='tc2')
TC2.STATUS = 1 # minimize error between simulation and measurement
TC2.FSTATUS = 1 # receive measurement
TC2.MEAS_GAP = 0.1 # measurement deadband gap
TC2.LOWER = 0
TC2.UPPER = 200
Ta = m.Param(value=23.0+273.15) # K
mass = m.Param(value=4.0/1000.0) # kg
Cp = m.Param(value=0.5*1000.0) # J/kg-K
A = m.Param(value=10.0/100.0**2) # Area not between heaters in m^2
As = m.Param(value=2.0/100.0**2) # Area between heaters in m^2
eps = m.Param(value=0.9) # Emissivity
sigma = m.Const(5.67e-8) # Stefan-Boltzmann
# Heater temperatures
T1 = m.Intermediate(TC1+273.15)
T2 = m.Intermediate(TC2+273.15)
# Heat transfer between two heaters
Q_C12 = m.Intermediate(U*As*(T2-T1)) # Convective
Q_R12 = m.Intermediate(eps*sigma*As*(T2**4-T1**4)) # Radiative
# Semi-fundamental correlations (energy balances)
m.Equation(TC1.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ Q_C12 + Q_R12 \
+ alpha1*Q1))
m.Equation(TC2.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- Q_C12 - Q_R12 \
+ alpha2*Q2))
# Global Options
m.options.IMODE = 5 # MHE
m.options.EV_TYPE = 2 # Objective type
m.options.NODES = 3 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.options.COLDSTART = 1 # COLDSTART on first cycle
##################################################################
# Create plot
plt.figure(figsize=(12,7))
plt.ion()
plt.show()
# Main Loop
start_time = time.time()
prev_time = start_time
tm = np.zeros(n)
try:
for i in range(1,n):
# Sleep time
sleep_max = 3.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep-0.01)
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
# Read temperatures in Celsius
T1m[i] = a.T1
T2m[i] = a.T2
# Insert measurements
TC1.MEAS = T1m[i]
TC2.MEAS = T2m[i]
Q1.MEAS = Q1s[i-1]
Q2.MEAS = Q2s[i-1]
# Start estimating U after 10 cycles (20 sec)
if i==10:
U.STATUS = 1
alpha1.STATUS = 1
alpha2.STATUS = 1
# Predict Parameters and Temperatures with MHE
m.solve()
if m.options.APPSTATUS == 1:
# Retrieve new values
Tmhe1[i] = TC1.MODEL
Tmhe2[i] = TC2.MODEL
Umhe[i] = U.NEWVAL
amhe1[i] = alpha1.NEWVAL
amhe2[i] = alpha2.NEWVAL
else:
# Solution failed, copy prior solution
Tmhe1[i] = Tmhe1[i-1]
Tmhe2[i] = Tmhe1[i-1]
Umhe[i] = Umhe[i-1]
amhe1[i] = amhe1[i-1]
amhe2[i] = amhe2[i-1]
# Write new heater values (0-100)
a.Q1(Q1s[i])
a.Q2(Q2s[i])
# Plot
plt.clf()
ax=plt.subplot(3,1,1)
ax.grid()
plt.plot(tm[0:i],T1m[0:i],'ro',label=r'$T_1$ measured')
plt.plot(tm[0:i],Tmhe1[0:i],'k-',label=r'$T_1$ MHE')
plt.plot(tm[0:i],T2m[0:i],'bx',label=r'$T_2$ measured')
plt.plot(tm[0:i],Tmhe2[0:i],'k--',label=r'$T_2$ MHE')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(3,1,2)
ax.grid()
plt.plot(tm[0:i],Umhe[0:i],'k-',label='Heat Transfer Coeff')
plt.plot(tm[0:i],amhe1[0:i]*1000,'r--',label=r'$\alpha_1$x1000')
plt.plot(tm[0:i],amhe2[0:i]*1000,'b--',label=r'$\alpha_2$x1000')
plt.ylabel('Parameters')
plt.legend(loc='best')
ax=plt.subplot(3,1,3)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',label=r'$Q_1$')
plt.plot(tm[0:i],Q2s[0:i],'b:',label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.draw()
plt.pause(0.05)
if make_mp4:
filename='./figures/plot_'+str(i+10000)+'.png'
plt.savefig(filename)
# Turn off heaters
a.Q1(0)
a.Q2(0)
# Save figure
plt.savefig('tclab_mhe.png')
# generate mp4 from png figures in batches of 350
if make_mp4:
images = []
iset = 0
for i in range(1,n):
filename='./figures/plot_'+str(i+10000)+'.png'
images.append(imageio.imread(filename))
if ((i+1)%350)==0:
imageio.mimsave('results_'+str(iset)+'.mp4', images)
iset += 1
images = []
if images!=[]:
imageio.mimsave('results_'+str(iset)+'.mp4', images)
# 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('tclab_mhe.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('tclab_mhe.png')
raise
See also:
Generative AI Learning
Use these prompts around the parameter estimation lab. The estimation setup and the identifiability argument are yours to own.
"Quiz me with 5 questions, one at a time, on estimating U, alpha1, and alpha2 in the TCLab MIMO energy balance: why estimating all three from a single simultaneous step of both heaters is ill-posed (which parameter pairs are correlated), what input excitation breaks the correlation, how the objective (sum of squared errors on T1 and T2) weights the two sensors, why the estimated U absorbs unmodeled conduction, and how moving horizon estimation differs from batch estimation on the same data. Grade my answers and list my misconceptions."
"Here is my estimation result: U = {value}, alpha1 = {value}, alpha2 = {value}, objective {value}, from the step sequence {describe your Q1/Q2 schedule}. As an adversarial reviewer: tell me whether my input schedule can actually distinguish U from alpha (explain using the steady-state balance), propose the smallest change to the experiment that improves identifiability, and ask me to predict how each parameter estimate would shift if the ambient temperature drifted 2 degC during the run. Grade my predictions."
App: Hand-fit the same parameters in the TCLab Simulation Studio before running GEKKO - if the optimizer's U and alpha differ from your hand fit by more than the app's objective difference justifies, one of the two setups has a bug worth finding.
What to Turn In
Submit a short report (PDF, 1-2 pages) that curates your results into a demonstration of what you learned. You may use Generative AI to help write the report, but you must guide it to the correct visualizations, justifications, and assumptions. Answer these questions:
- Show the fit: measured vs predicted T1 and T2 with your estimated U, alpha1, alpha2 (units), the input schedule, and whether the data is hardware or simulated.
- Make the identifiability argument: what in your input design lets the optimizer separate U from alpha1/alpha2? Support it with a second estimation run on a worse input schedule (e.g., both heaters stepped together) and compare the estimates.
- Compare batch estimation vs MHE parameter trajectories on the same dataset. Where do they disagree and why is that expected?
- The largest residual in your fit: defend it physically or concede a model defect - and say which model term you would add first.