TCLab State Space Model
![]() | ![]() | ![]() | ![]() | ![]() |
---|
Objective: Derive and simulate a state space model for two heaters and compare the model predictions with the Arduino temperature control lab measurements.
A TCLab state space model is a linear representation of the physics-based model. Putting a model into state space form is the basis for many methods in process dynamics and control analysis and will be used for model predictive control in a later exercise.
![](https://apmonitor.com/pdc/uploads/Main/tclab_sac.png)
A physics-based model is developed in prior exercises. A 2nd-order physics-based model consists of four differential equations from energy balances that include convection, conduction, and thermal radiation.
$$m\,c_p\frac{dT_{H1}}{dt} = U\,A\,\left(T_\infty-T_{H1}\right) + \epsilon\,\sigma\,A\,\left(T_\infty^4-T_{H1}^4\right) + Q_{C12} + Q_{R12} + \alpha_1 Q_1$$
$$m\,c_p\frac{dT_{H2}}{dt} = U\,A\,\left(T_\infty-T_{H2}\right) + \epsilon\,\sigma\,A\,\left(T_\infty^4-T_{H2}^4\right) - Q_{C12} - Q_{R12} + \alpha_2 Q_2$$
$$\tau \frac{dT_{C1}}{dt} = T_{H1} - T_{C1}$$
$$\tau \frac{dT_{C2}}{dt} = T_{H2} - T_{C2}$$
with the following heat transfer terms for simplifying the expressions.
$$Q_{C12} = U_s \, A_s \, \left(T_{H2}-T_{H1}\right)$$
$$Q_{R12} = \epsilon\,\sigma\,A_s\,\left(T_{H2}^4-T_{H1}^4\right)$$
$$\tau = \frac{m_s \, c_{p,s} \, \Delta x}{k\,A_c}$$
Use these equations to develop a state space model of the temperature response to heater changes.
$$\begin{bmatrix}\dot T'_{H1}\\\dot T'_{H2}\\\dot T'_{C1}\\\dot T'_{C2}\end{bmatrix} = \begin{bmatrix}a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4}\\a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4}\\a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4}\\a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4}\end{bmatrix} \begin{bmatrix}T'_{H1}\\T'_{H2}\\T'_{C1}\\T'_{C2}\end{bmatrix} + \begin{bmatrix}b_{1,1} & b_{1,2}\\b_{2,1} & b_{2,2}\\b_{3,1} & b_{3,2}\\b_{4,1} & b_{4,2}\end{bmatrix} \begin{bmatrix}Q'_1\\Q'_2\end{bmatrix}$$
$$\begin{bmatrix}T'_{C1}\\T'_{C2}\end{bmatrix} = \begin{bmatrix}0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{bmatrix} \begin{bmatrix}T'_{H1}\\T'_{H2}\\T'_{C1}\\T'_{C2}\end{bmatrix} + \begin{bmatrix}0&0\\0&0\end{bmatrix} \begin{bmatrix}Q'_1\\Q'_2\end{bmatrix}$$
This state space model:
$$\dot x = A x + B u$$
$$y = C x + D u$$
has four matrices with `A` and `B` that need to be obtained from partial derivatives of the right hand side of the differential equations.
$$A \in \mathbb{R}^{4 \, \mathrm{x} \, 4} = \begin{bmatrix}a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4}\\a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4}\\a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4}\\a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4}\end{bmatrix}$$ $$B \in \mathbb{R}^{4 \, \mathrm{x} \, 2} = \begin{bmatrix}b_{1,1} & b_{1,2}\\b_{2,1} & b_{2,2}\\b_{3,1} & b_{3,2}\\b_{4,1} & b_{4,2}\end{bmatrix} $$ $$C \in \mathbb{R}^{2 \, \mathrm{x} \, 4} = \begin{bmatrix}0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{bmatrix}$$ $$D \in \mathbb{R}^{2 \, \mathrm{x} \, 2} = \begin{bmatrix}0&0\\0&0\end{bmatrix}$$
The individual elements are obtained by evaluating partial derivatives of the differential equations at steady state conditions.
$$a_{1,1} = \frac{\partial f_1}{\partial T_{H1}}\bigg|_{\bar T} = -\frac{U\,A + 4 \epsilon \sigma A T_0^3 + U_s \, A_s + 4 \epsilon \sigma A_s T_0^3}{m\,C_p}$$
$$a_{1,2} = \frac{\partial f_1}{\partial T_{H2}}\bigg|_{\bar T} = \frac{U_s \, A_s + 4 \epsilon \sigma A_s T_0^3}{m\,C_p}$$
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'
data = pd.read_csv(filename)
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)
Q2d[120:320] = 100
Q2d[320:520] = 10
Q2d[520:] = 80
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')
# close connection to Arduino
a.close()
fid.close()
except:
filename = 'https://apmonitor.com/pdc/uploads/Main/tclab_data2.txt'
# read either local file or use link if no TCLab
data = pd.read_csv(filename)
# Fit Parameters of Energy Balance
m = GEKKO() # Create GEKKO Model
# Parameters to Estimate
U = m.FV(value=10,lb=1,ub=20)
Us = m.FV(value=20,lb=5,ub=40)
alpha1 = m.FV(value=0.01,lb=0.001,ub=0.03) # W / % heater
alpha2 = m.FV(value=0.005,lb=0.001,ub=0.02) # W / % heater
tau = m.FV(value=10.0,lb=5.0,ub=60.0)
# Measured inputs
Q1 = m.Param()
Q2 = m.Param()
Ta =23.0+273.15 # K
mass = 4.0/1000.0 # kg
Cp = 0.5*1000.0 # J/kg-K
A = 10.0/100.0**2 # Area not between heaters in m^2
As = 2.0/100.0**2 # Area between heaters in m^2
eps = 0.9 # Emissivity
sigma = 5.67e-8 # Stefan-Boltzmann
TH1 = m.SV()
TH2 = m.SV()
TC1 = m.CV()
TC2 = m.CV()
# Heater Temperatures in Kelvin
T1 = m.Intermediate(TH1+273.15)
T2 = m.Intermediate(TH2+273.15)
# Heat transfer between two heaters
Q_C12 = m.Intermediate(Us*As*(T2-T1)) # Convective
Q_R12 = m.Intermediate(eps*sigma*As*(T2**4-T1**4)) # Radiative
# Energy balances
m.Equation(TH1.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ Q_C12 + Q_R12 \
+ alpha1*Q1))
m.Equation(TH2.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- Q_C12 - Q_R12 \
+ alpha2*Q2))
# Conduction to temperature sensors
m.Equation(tau*TC1.dt() == TH1-TC1)
m.Equation(tau*TC2.dt() == TH2-TC2)
# Options
# STATUS=1 allows solver to adjust parameter
U.STATUS = 1
Us.STATUS = 1
alpha1.STATUS = 1
alpha2.STATUS = 1
tau.STATUS = 1
Q1.value=data['Q1'].values
Q2.value=data['Q2'].values
TH1.value=data['T1'].values[0]
TH2.value=data['T2'].values[0]
TC1.value=data['T1'].values
TC1.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
TC2.value=data['T2'].values
TC2.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
m.time = data['Time'].values
m.options.IMODE = 5 # MHE
m.options.EV_TYPE = 2 # Objective type
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.solve(disp=False) # Solve physics-based model
# Parameter values
print('Estimated Parameters')
print('U : ' + str(U.value[0]))
print('Us : ' + str(Us.value[0]))
print('alpha1: ' + str(alpha1.value[0]))
print('alpha2: ' + str(alpha2.value[-1]))
print('tau: ' + str(tau.value[0]))
print('Constants')
print('Ta: ' + str(Ta))
print('m: ' + str(mass))
print('Cp: ' + str(Cp))
print('A: ' + str(A))
print('As: ' + str(As))
print('eps: ' + str(eps))
print('sigma: ' + str(sigma))
sae = 0.0
for i in range(len(data)):
sae += np.abs(data['T1'][i]-TC1.value[i])
sae += np.abs(data['T2'][i]-TC2.value[i])
print('SAE Energy Balance: ' + str(sae))
# Create plot
plt.figure(figsize=(10,7))
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(data['Time'],data['T1'],'r.',label=r'$T_1$ measured')
plt.plot(m.time,TC1.value,color='black',linestyle='--',\
linewidth=2,label=r'$T_1$ energy balance')
plt.plot(data['Time'],data['T2'],'b.',label=r'$T_2$ measured')
plt.plot(m.time,TC2.value,color='orange',linestyle='--',\
linewidth=2,label=r'$T_2$ energy balance')
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.plot(data['Time'],data['Q2'],'b:',\
linewidth=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.savefig('tclab_2nd_order_physics.png')
plt.show()
Fill in the values of the state space matrices after the parameters (U, Us, `\tau`, `\alpha_1`, and `\alpha_2`) are estimated.
Bm = np.zeros((4,2))
Cm = np.zeros((2,4))
Dm = np.zeros((2,2))
T0 = Ta
c1 = U.value[0]*A
c2 = 4*eps*sigma*A*T0**3
c3 = Us.value[0]*As
c4 = 4*eps*sigma*As*T0**3
c5 = mass*Cp
c6 = 1/tau.value[0]
Am[0,0] = -(c1+c2+c3+c4)/c5
Am[0,1] = (c3+c4)/c5
Am[1,0] =
Am[1,1] =
Am[2,0] = c6
Am[2,2] =
Am[3,1] =
Am[3,3] = -c6
Bm[0,0] = alpha1.value[0]/c5
Bm[1,1] =
Cm[0,2] = 1
Cm[1,3] =
Add a simulation of the linear model with either Gekko or Scipy with the same heater inputs as used during the test.
Simulate Linear Model Am, Bm, Cm, Dm with Gekko
ss = GEKKO()
x,y,u = ss.state_space(Am,Bm,Cm,D=None)
u[0].value = data['Q1'].values
u[1].value = data['Q2'].values
ss.time = data['Time'].values
ss.options.IMODE = 7
ss.solve(disp=False)
Simulate Linear Model Am, Bm, Cm, Dm with Scipy
sys = signal.StateSpace(Am,Bm,Cm,Dm)
tsys = data['Time'].values
Qsys = np.vstack((data['Q1'].values,data['Q2'].values))
Qsys = Qsys.T
tsys,ysys,xsys = signal.lsim(sys,Qsys,tsys)
Add the predictions to the plot by adding the nominal temperature (T0=21oC) to the output to return to oC from deviation variable form.
Solution
![](https://apmonitor.com/pdc/uploads/Main/tclab_state_space.png)
import matplotlib.pyplot as plt
import pandas as pd
from gekko import GEKKO
import tclab
import time
from scipy import signal
# Import data
try:
# try to read local data file first
filename = 'data.csv'
data = pd.read_csv(filename)
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)
Q2d[120:320] = 100
Q2d[320:520] = 10
Q2d[520:] = 80
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')
# close connection to Arduino
a.close()
fid.close()
except:
filename = 'https://apmonitor.com/pdc/uploads/Main/tclab_data2.txt'
# read either local file or use link if no TCLab
data = pd.read_csv(filename)
# Fit Parameters of Energy Balance
m = GEKKO() # Create GEKKO Model
# Parameters to Estimate
U = m.FV(value=10,lb=1,ub=20)
Us = m.FV(value=20,lb=5,ub=40)
alpha1 = m.FV(value=0.01,lb=0.001,ub=0.03) # W / % heater
alpha2 = m.FV(value=0.005,lb=0.001,ub=0.02) # W / % heater
tau = m.FV(value=10.0,lb=5.0,ub=60.0)
# Measured inputs
Q1 = m.Param()
Q2 = m.Param()
Ta =21.0+273.15 # K
mass = 4.0/1000.0 # kg
Cp = 0.5*1000.0 # J/kg-K
A = 10.0/100.0**2 # Area not between heaters in m^2
As = 2.0/100.0**2 # Area between heaters in m^2
eps = 0.9 # Emissivity
sigma = 5.67e-8 # Stefan-Boltzmann
TH1 = m.SV()
TH2 = m.SV()
TC1 = m.CV()
TC2 = m.CV()
# Heater Temperatures in Kelvin
T1 = m.Intermediate(TH1+273.15)
T2 = m.Intermediate(TH2+273.15)
# Heat transfer between two heaters
Q_C12 = m.Intermediate(Us*As*(T2-T1)) # Convective
Q_R12 = m.Intermediate(eps*sigma*As*(T2**4-T1**4)) # Radiative
# Energy balances
m.Equation(TH1.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ Q_C12 + Q_R12 \
+ alpha1*Q1))
m.Equation(TH2.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- Q_C12 - Q_R12 \
+ alpha2*Q2))
# Conduction to temperature sensors
m.Equation(tau*TC1.dt() == TH1-TC1)
m.Equation(tau*TC2.dt() == TH2-TC2)
# Options
# STATUS=1 allows solver to adjust parameter
U.STATUS = 1
Us.STATUS = 1
alpha1.STATUS = 1
alpha2.STATUS = 1
tau.STATUS = 1
Q1.value=data['Q1'].values
Q2.value=data['Q2'].values
TH1.value=data['T1'].values[0]
TH2.value=data['T2'].values[0]
TC1.value=data['T1'].values
TC1.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
TC2.value=data['T2'].values
TC2.FSTATUS = 1 # minimize fstatus * (meas-pred)^2
m.time = data['Time'].values
m.options.IMODE = 5 # MHE
m.options.EV_TYPE = 2 # Objective type
m.options.NODES = 2 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.solve(disp=False) # Solve physics-based model
# Parameter values
print('Estimated Parameters')
print('U : ' + str(U.value[0]))
print('Us : ' + str(Us.value[0]))
print('alpha1: ' + str(alpha1.value[0]))
print('alpha2: ' + str(alpha2.value[-1]))
print('tau: ' + str(tau.value[0]))
print('Constants')
print('Ta: ' + str(Ta))
print('m: ' + str(mass))
print('Cp: ' + str(Cp))
print('A: ' + str(A))
print('As: ' + str(As))
print('eps: ' + str(eps))
print('sigma: ' + str(sigma))
Am = np.zeros((4,4))
Bm = np.zeros((4,2))
Cm = np.zeros((2,4))
Dm = np.zeros((2,2))
T0 = Ta
c1 = U.value[0]*A
c2 = 4*eps*sigma*A*T0**3
c3 = Us.value[0]*As
c4 = 4*eps*sigma*As*T0**3
c5 = mass*Cp
c6 = 1/tau.value[0]
Am[0,0] = -(c1+c2+c3+c4)/c5
Am[0,1] = (c3+c4)/c5
Am[1,0] = (c3+c4)/c5
Am[1,1] = -(c1+c2+c3+c4)/c5
Am[2,0] = c6
Am[2,2] = -c6
Am[3,1] = c6
Am[3,3] = -c6
Bm[0,0] = alpha1.value[0]/c5
Bm[1,1] = alpha2.value[0]/c5
Cm[0,2] = 1
Cm[1,3] = 1
# state space simulation
ss = GEKKO()
x,y,u = ss.state_space(Am,Bm,Cm,D=None)
u[0].value = data['Q1'].values
u[1].value = data['Q2'].values
ss.time = data['Time'].values
ss.options.IMODE = 7
ss.solve(disp=False)
# state space simulation with scipy
sys = signal.StateSpace(Am,Bm,Cm,Dm)
tsys = data['Time'].values
Qsys = np.vstack((data['Q1'].values,data['Q2'].values))
Qsys = Qsys.T
tsys,ysys,xsys = signal.lsim(sys,Qsys,tsys)
sae = 0.0
for i in range(len(data)):
sae += np.abs(data['T1'][i]-TC1.value[i])
sae += np.abs(data['T2'][i]-TC2.value[i])
print('SAE Energy Balance: ' + str(sae))
# Create plot
plt.figure(figsize=(10,7))
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(data['Time'],data['T1'],'r.',label=r'$T_1$ measured')
plt.plot(m.time,TC1.value,color='black',linestyle='--',\
linewidth=2,label=r'$T_1$ energy balance')
plt.plot(ss.time,np.array(y[0].value)+21.0,color='gray',linestyle=':',\
linewidth=2,label=r'$T_1$ state space')
plt.plot(data['Time'],data['T2'],'b.',label=r'$T_2$ measured')
plt.plot(m.time,TC2.value,color='orange',linestyle='--',\
linewidth=2,label=r'$T_2$ energy balance')
plt.plot(ss.time,np.array(y[1].value)+21.0,color='purple',linestyle=':',\
linewidth=2,label=r'$T_2$ state space')
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.plot(data['Time'],data['Q2'],'b:',\
linewidth=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.savefig('tclab_state_space.png')
plt.show()