2.2. Temperature Control Lab, step identification#

The temperature control lab or TCLab is a laboratory device used for thermal projects. It can be used to simulate and develop thermal control and regulation systems. The TCLab features a feedback control application with an Arduino, an LED, 2 heaters and 2 temperature sensors:

  • one heater will represent the TCS heating elements used for temperature control.

  • another will represent external disturbances such as solar radiation in the Cubesat application.

During this tutorial we will start using this card and carry out temperature measurements on control steps. The aim is to identify a dynamic model that can be used to synthesise the control. In this second part of the tutorial, we’ll identify a transfer function to represent the step response previously measured.

Follow the steps below in Colab, using the measurement file provided. When revising, you can repeat these steps with your measurements in Anaconda.

2.2.1. Fitting the Step Change Data to a First Order Model#

For a first-order linear system initially at steady-state, the response to a step input change at \(t=0\) is given by

\[y(t) = y(0) + K(1 - e^{-t/\tau}) \Delta U\]

where \(\Delta U\) is the magnitude of the step change. Converting to notation used for the temperature control lab where \(y(t) = T_1(t)\) and \(\Delta U = \Delta Q_1\)

\[T_1(t) = T_1(0) + K_1(1 - e^{-t/\tau_1}) \Delta Q_1\]

Questions:

  • Recall the differential equation and the transfer function of a 1st order SISO (Single Input Single Output) dynamic system. Use notations compatible with the description of the step response above.

  • Identify the steady state gain \(K_1\) and time constant \(\tau_1\)

  • Compare the 1st order response with the measured data. Conclusion.

# With this line : no print 
%%script echo skipping

import pandas as pd

# Read data
df = pd.read_csv('https://github.com/SizingLab/SystemsEngineeringIntroduction/blob/main/notebook/Step_Test_Data.csv') # provided data 
# df = pd.read_csv('./data/Step_Test_Data.csv') # my data 

df = df.set_index('Time')

T1 = df['T1']
Q1 = df['Q1']

# Static gain K1
# complete

# Time constant tau1
# complete


import matplotlib.pyplot as plt
import numpy as np

exp = np.exp
t = df.index

T1_est = T1.min() + K1*(1 - exp(-t/tau1))*DeltaQ1

plt.figure(figsize=(10,5))
ax = plt.subplot(2,1,1)
df['T1'].plot(ax = ax, grid=True)
plt.plot(t,T1_est)
plt.title('Step Test Data Compared to Model')

plt.subplot(2,1,2)
plt.plot(t,T1_est-T1)
plt.grid()
plt.title('Residual Error')
plt.show()

UsageError: Line magic function `%%script` not found.

2.2.2. Nodal thermal model#

We will show here that a physical model made up of thermal resistances and capacities can explain the dynamic responses observed and the transfer functions identified previously.

Questions:

  • Recall the characteristic equations for thermal resistances and capacities relating temperatures \(T_i\) and heat flows (power) \(P\).

  • Select the thermal model to represent the temperature rise if we assume that the heating element and the temperature sensor are at the same temperature.

  • Determine the values of the thermal resistance \(R_{th}\) and the thermal capacity \(C_{th}\) from the 1st order transfer function identified above. Assume that the maximum power of the heating element (4W) and the outside temperature are known.

Modèle thermique 1 corps

# With this line : no print 
%%script echo skipping

%matplotlib inline
import numpy as np                  # basic package for numerical calculations
import matplotlib.pyplot as plt     # plotting package
from scipy.signal import step       # Import from scipy of the step response

# parameter values and units
# complete

# First order step response
num = [Rth]
den = [Rth*Cth, 1]
t, y = step(system=(num, den))
T1_est = T_ambient  + P1*y

plt.figure(figsize=(10,5))
ax = plt.subplot(1,1,1)
df['T1'].plot(ax = ax, grid=True)
plt.plot(t,T1_est)
plt.title('Step Test Data Compared to Model')
plt.show()

#Print data
print("Thermal resistance Rth = %.2f K/W"%Rth)
print("Thermal capcitance Cth = %.2f J/K"%Cth)

UsageError: Line magic function `%%script` not found.

The parameter values in the above plot were chosen to (at least roughly) reproduce the measured response of temperature control laboratory. The specific heat capacity for solids is typically has values in the range of 0.2 to 0.9 J/K/gram. Using a value of 0.9 that is typical of aluminum and plastics used for electronic products, the estimated mass of the heater/sensor pair would be:

print("Equivalent mass : %.2f g"%(Cth/0.9))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 1
----> 1 print("Equivalent mass : %.2f g"%(Cth/0.9))

NameError: name 'Cth' is not defined

2.2.3. Two State Model#

2.2.3.1. From first order to second order thermal model#

Question: Explain why the 1st order model is not completely satisfactory.

Question: Propose a nodal thermal diagram to represent the observed dynamics, taking into account the response of the temperature sensor. For this third model, we consider the possibility that the heater and sensor may not be at the same temperature. In other words, that the heater/sensor assembly is not at a uniform temperature. To account for this possibility, we introduce \(𝑇_𝐻\) to denote the temperature of heater one and \(𝑇_𝑆\) to denote the temperature of the corresponding sensor. We’ll further assume that sensor mainly exchanges heat with the heater, and the dominant heat transfer to the surroundings is through the heat sink attached to the heater.

Exercice or homework: Determine the form of the transfer function corresponding to this thermal diagram. The input is the temperature rise in %, the output is the temperature measurement in °C.

2.2.3.2. Fitting a Second Order Model by Least Squares#

Graphically determining a large number of parameters can be difficult. Second-order parameters can be determined using the least_square function, which minimises an error by optimisation.

Exercice: Propose a Python code to determine the second-order transfer function best suited to the measurement.

# With this line : no print 
%%script echo skipping

from scipy.optimize import least_squares
import numpy as np
Qmax = 50

def f(x):
# complete
    return resid

ic = [0.86,40,130,20]

r = least_squares(f,ic,bounds=(0,np.inf))
r.x

UsageError: Line magic function `%%script` not found.

For the next tutorial, we’ll use the following transfer function: \( G_1(s) = \frac{0.65}{(27s + 1)(160s + 1)} \)

# With this line : no print 
%%script echo skipping

%matplotlib inline
import numpy as np                  # basic package for numerical calculations
import matplotlib.pyplot as plt     # plotting package
from scipy.signal import step       # Import from scipy of the step response

# First order step response
num = [0.657]
den = [160*27,160+27,1]
t, y = step(system=(num, den),T=np.linspace(0,600,100))
T1_est = 20.25  + y*50

plt.figure(figsize=(10,5))
ax = plt.subplot(1,1,1)
df['T1'].plot(ax = ax, grid=True)
plt.plot(t,T1_est)
plt.title('Step Test Data Compared to Model')
plt.show()

UsageError: Line magic function `%%script` not found.

2.2.4. Correction#

An example of correction can be found here

2.2.5. References#

[Hedengren] Hedengren, J., Process Dynamics and Control, Brigham Young University. Link

[Kantor, 2021] Kantor, JC., Chemical Process Control class materials, University of Notre Dame. Link / Github