Lab 3: Model Identification

Lab 3: Model Identification#

Fitting models to data is an engineering skill that links between the real world of engineering systems to the theory you’ve been learning in the classroom. For this laboratory session you will collect data from a step test experiments, then fit the data to models derived from first-principles energy balances.

Procedure#

  1. Download a copy of this notebook to your laptop and complete the the data collection exercises shown below. The results should be embedded in the notebook. Be sure to ‘save-as-you-go’ to avoid losing your work.

  2. Use your step test data to fit two versions of a model for the temperature control lab:

  • One-state model with one input.

  • Two-state model with one input.

  1. Submit your completed lab notebook via Gradescope.

Pre-lab Exercise#

  1. Read through the entirety of the assignment before coming to lab.

  2. Write pseudocode for estimating the model parameters with uncertainty in Exercise 3.

  3. Bring pseudocode to lab on a piece of paper or a tablet. The instructor and TAs will likely ask to see this as a starting point for debugging.

Exercise 1. Step Tests and Data Collection#

Step 1. Verify operation.#

Execute the following cell to verify that you have a working connection to the temperature control lab hardware. This will test for installation of TCLab.py, connection to the Arduino device, and working firmware within the Arduino.

from tclab import TCLab, clock, Historian, Plotter, setup


TCLab = setup(connected=False)

lab = TCLab()
print("TCLab Temperatures:", lab.T1, lab.T2)
lab.close()

Step 2. Check for steady state#

As discussed in class, for step testing the device must be initially at steady state. Run the following code to verify the heaters are off and that the temperatures are at a steady ambient temperature.

# experimental parameters
tfinal = 30

# perform experiment
with TCLab() as lab:
    lab.U1 = 0
    lab.U2 = 0
    h = Historian(lab.sources)
    p = Plotter(h, tfinal)
    for t in clock(tfinal):
        p.update(t)

Step 3. Step Test#

The step test consists of turning one heater one at 50% power and recording temperature data for at least 900 seconds. Copy the code from Step 2 into the following cell. Then modify as needed to accomplish the step test with P1 at 200 and using 50% of maximum power.

# Add your solution here

Step 4. Save data to a .csv file#

Run the following cell to verify and save your data to a ‘.csv’ file. Be sure you can find and locate the data on your laptop before ending your session. You will need access to this data for subsequent exercises.

import matplotlib.pyplot as plt
import os.path

t, T1, T2, Q1, Q2 = h.fields

plt.plot(t, T1, t, T2, t, Q1, t, Q2)
plt.legend(['T1','T2','Q1','Q2'])
plt.xlabel('Time / seconds')
plt.grid()

# Change the filename here
data_file = 'tclab-model-id.csv'

# Set to True to overwrite the file. Default is False
# to safeguard against accidentally rerunning this cell.
overwrite_file = False

if not overwrite_file and os.path.isfile('./'+data_file):
    raise FileExistsError(data_file + ' already exisits. Either choose a new filename or set overwrite_file = True.')
else:
    h.to_csv(data_file)
    print("Successfully saved data to "+data_file)

Exercise 2. Fitting a First-Order Model#

As discussed in class, a simple energy balance model for T1 is given by

\[C_p \frac{dT_1}{dt} = U_a(T_{amb} - T_1) + \alpha P_1 U_1\]

where the parameter \(\alpha\) has, through independent means, been determined as 0.00016 when U1 is given in percent of full power and power is measured in Watts. Following the notes from Section 2.4, use the results of this experiment to estimate values for \(C_p\) and \(U_a\). Write your answers in the following cell.

You may cut and paste code from Section 2.4 for this task. You will need to modify the code to read your data file. Be sure to report

  1. Final estimates for \(U_a\) and \(C_p\)

  2. An estimate of the time constant computed from \(U_a\) and \(C_p\).

  3. An estimate of the steady-state gain computed from \(U_a\).

# Load and plot data
# Add your solution here

# Implement and test model
# Add your solution here

# Fit model
# Add your solution here

# Plot histogram of residuals
# Add your solution here

# Estimate covariance and correlation of fitted parameters
# Add your solution here

print("Covariance of p = [Ua, Cp]")
print(cov_p)

# Convert to correlation matrix

def covariance_to_correlation(cov):
    ''' Convert covariance matrix into correlation matrix

    Argument:
        cov: covariance matrix

    Returns:
        cor: correlation matrix

    '''

    # Copy matrix
    cor = cov.copy()

    # Get number of rows
    n = cor.shape[0]

    # Loop over rows
    for r in range(n):
        # Loop over columns
        for c in range(n):
            # Scale element
            cor[r,c] = cor[r,c] / np.sqrt(cov[r,r]*cov[c,c])

    return cor

print("Correlation matrix")
print(covariance_to_correlation(cov_p))

# Calculate and print model parameters with uncertainty

# Add your solution here

print("The value of Ua is",round(Ua,4), "+/-", round(np.sqrt(cov_p[0,0]),4), "Watts/degC.")
print("The value of Cp is", round(Cp,3), "+/-", round(np.sqrt(cov_p[1,1]),3), "Joules/degC.")
print("The value of the time constant (tau) is", round(tau,3), "+/-", round(np.sqrt(var_tau),3), "seconds.")
print("The value of the steady state gain (K) is",round(K,3),"+/-", round(np.sqrt(var_K),3), "Watts.")

Exercise 3. Fitting a Second-Order Model#

A second order model is for the heater/sensor combination is given by

\[\begin{split} \begin{align} C^H_p\frac{dT_{H,1}}{dt} & = U_a(T_{amb} - T_{H,1}) + U_b(T_{S,1} - T_{H,1}) + P_1u_1\\ C^S_p\frac{dT_{S,1}}{dt} & = U_b(T_{H,1} - T_{S,1}) \end{align} \end{split}\]

where \(T_{H,1}\) is the heater temperature, \(T_{S,1}\) is the sensor temperature, and \(U_b\) is the heat transfer coefficient between the heater and sensor.

Modify the code you developed for Exercise 2 to fit this second order model. The code shown in Section 2.5 of the notes may be helpful.

Do the following:

  1. Report your best fit for \(U_a\), \(U_b\), \(C^H_p\), and \(C^S_p\) with uncertainty.

  2. Calculate the time constants for the heater and sensor and report with uncertainty.

# Load and plot data
# Add your solution here

# Implement and test model

# Add your solution here

# Regress model

# Add your solution here

# Plot fitted model

# Add your solution here

# Plot histogram of residuals
# Add your solution here

# Estimate covariance and correlation of fitted parameters

# Estimate covariance and correlation of fitted parameters
# Add your solution here

print("\nCovariance of p = [Ua, Ub, CpH, CpS]")
print(cov_p)

# Convert to correlation matrix
print("\nCorrelation matrix")
print(covariance_to_correlation(cov_p))

Discussion Question: What are the time constants for the heater and sensor?

# Calculate and print model parameters with uncertainty

# Add your solution here

print("The value of Ua is", round(Ua,4), "+/-", round(np.sqrt(cov_p[0,0]),4), "Watts/degC.")
print("The value of Ub is", round(Ub,4), "+/-", round(np.sqrt(cov_p[1,1]),4), "Watts/degC.")
print("The value of CpH is", round(CpH,3), "+/-", round(np.sqrt(cov_p[2,2]),3), "Joules/degC.")
print("The value of CpS is", round(CpS,3), "+/-", round(np.sqrt(cov_p[3,3]),3), "Joules/degC.")

print("The value of the time constant TH1 is", round(TH1,3), "+/-", round(np.sqrt(cov_tau[0,0]),3), "seconds.")
print("The value of the time constant TS1 is", round(TS1,3), "+/-", round(np.sqrt(cov_tau[1,1]),3), "seconds.")
print("The value of the time constant TH1 + TS1 is", round(addition,3), "+/-", round(np.sqrt(cov_tau[2,2]),3), "seconds.")

Discussion Question: Do you see any relationship between the heater and sensor time constants and what you found for the first-order model?

Answer:

Contents