Exercise: ParmEst - Pyomo.DoE Workshop and Tutorials

In the main workshop tutorial, the four heat transfer parameters, $U_a$ , $U_b$ , $C_p^H$ , $C_p^S$ , were estimated from the sine test data using the following TC Lab model:

C_p^H \frac{dT_H}{dt} = U_a (T_{\text{amb}} - T_H) + U_b (T_S - T_H) + \alpha P u(t),

(1)

C_p^S \frac{dT_S}{dt} = U_b (T_H - T_S)

(2)

where $T_{H}$ and $T_{S}$ are the heater and sensor temperatures in $^{\circ} \text{C}$ , respectively, $T_{amb}$ is the ambient temperature in $^{\circ} \text{C}$ , $C_p^H$ and $C_p^S$ are the heat capacity of the heater and sensor in $\text{J}/^{\circ} \text{C}$ , respectively, $U_{a}$ and $U_{b}$ are the heat transfer coefficients from the heater to the sensor and the sensor to ambient (in $\text{W}/^{\circ} \text{C}$ ), respectively, $P$ is the maximum power limit in $\text{bit}$ , $u$ is the heater power in $\%$ , and $\alpha$ is constant that converts the unit of $P u(t)$ into $\text{W}$ .

As discussed in the uncertainty quantification section of the main workshop, the above model is structurally non-identifiable (i.e., the $U_a$ , $U_b$ , $C_p^H$ , and $C_p^S$ parameters cannot be reliably estimated from the mathematical structure of the model). To ensure accurate estimates of the model parameters, we need to explore other structures of the model that simplify parameter correlation.

In this exercise notebook, we present a reformulated, linearized (linear with respect to model parameters) version of the TC Lab model to reliably estimate the original parameters, $U_a$ , $U_b$ , $C_p^H$ , and $C_p^S$ from the sine test data. The structure of the reformulated model is:

\frac{dT_H}{dt} = \beta_1 (T_{\text{amb}} - T_H) + \beta_2 (T_S - T_H) + \beta_4 u(t),

(3)

\frac{dT_S}{dt} = \beta_3 (T_H - T_S)

(4)

where

\beta_1 = \frac{U_a}{C_p^H}, \quad \beta_2 = \frac{U_b}{C_p^H}, \quad \beta_3 = \frac{U_b}{C_p^S}, \quad \beta_4 = \frac{\alpha P}{C_p^H}

(5)

In this exercise, you will practice using ParmEst to estimate the reformulated model and then transform the estimates into the four original parameters ( $U_a$ , $U_b$ , $C_p^H$ , and $C_p^S$ ).

Implementing Reformulated Model in Pyomo¶

Take a few minutes to study the tclab_pyomo.py. Lines 420 to 520 (approximately) in file includes an implementation of the reformulated model. This file also supports regressing either the heat transfer coefficients, $U_a$ and $U_b$ , or their inverses, $1/U_a$ and $1/U_b$ , as model parameters.

Import the Necessary Packages for This Exercise¶

import sys
import numpy as np
import pandas as pd
import logging

# If running on Google Colab, install Pyomo and Ipopt via IDAES
on_colab = "google.colab" in sys.modules
if on_colab:
    !wget "https://raw.githubusercontent.com/dowlinglab/pyomo-doe/main/notebooks/tclab_pyomo.py"
else:
    import os

    if "exercise_solutions" in os.getcwd():
        # Add the "notebooks" folder to the path
        # This is needed for running the solutions from a separate folder
        # You only need this if you run locally
        sys.path.append("../notebooks")

# import TCLab model, simulation, and data analysis functions
from tclab_pyomo import (
    TC_Lab_data,
    TC_Lab_experiment,
    extract_results,
    extract_plot_results,
    plot_profile_likelihood,
    reformulate_parameters,
    recover_original_parameters,
    recover_original_covariance,
)

# set default number of states in the TCLab model
number_tclab_states = 2

Load and Explore Experimental Data (sine test)¶

if on_colab:
    file = "https://raw.githubusercontent.com/dowlinglab/pyomo-doe/main/data/tclab_sine_test_5min_period.csv"
else:
    file = "../data/tclab_sine_test_5min_period.csv"
df = pd.read_csv(file)
df[["Time", "T1", "Q1", "Q2"]].head()  # the Q2 column entries are 0 because
# we are considering the two-state model

Make two plots to visualize the temperature and heat power data as a function of time.

# Add your solution here

# Add your solution here

We’ll now store the data in this custom data class objective. This is a nice trick to help keep data organized, but it is NOT required to use ParmEst or Pyomo data. Alternatively, we could just use a pandas DataFrame.

tc_data = TC_Lab_data(
    name="Sine Wave Test for Heater 1",
    time=df["Time"].values,
    T1=df["T1"].values,
    u1=df["Q1"].values,
    P1=200,
    TS1_data=None,
    T2=df["T2"].values,
    u2=df["Q2"].values,
    P2=200,
    TS2_data=None,
    Tamb=df["T1"].values[0],
)

Our custom data class has a method to export the data as a Pandas Data Frame.

tc_data.to_data_frame().iloc[:, :4].head()

Parameter Estimation with ParmEst¶

Now for the main event: performing nonlinear least squares with ParmEst.

We seek to estimate the original parameters, $C_p^H$ , $C_p^S$ , $U_a$ , and $U_b$ from the reformulated TC Lab model (with $\beta_1$ , $\beta_2$ , $\beta_3$ , and $\beta_4$ as parameters) and the sensor data presented above.

The reformulated TC Lab model predicts the sensor data as follows:

\frac{dT_H}{dt} = \beta_1 (T_{\text{amb}} - T_H) + \beta_2 (T_S - T_H) + \beta_4 u(t),

(6)

\frac{dT_S}{dt} = \beta_3 (T_H - T_S)

(7)

\begin{align*} \text{control input data}\qquad u(t_i) & = \bar{u}_{i}, \forall i \in \mathcal{T} \\ \text{initial condition}\qquad T_H(t_0) & = T_{amb} \\ \text{initial condition}\qquad T_S(t_0) & = T_{amb} \end{align*}

(8)

Remember that

\beta_1 = \frac{U_a}{C_p^H}, \quad \beta_2 = \frac{U_b}{C_p^H}, \quad \beta_3 = \frac{U_b}{C_p^S}, \quad \beta_4 = \frac{\alpha P}{C_p^H}

(9)

In the tclab_pyomo.py model, we defined several helper functions:

extract_results takes a Pyomo model and returns the results stored in an instance of the TC_Lab_data dataclass.
extract_plot_results takes experimental data (stored in a TC_Lab_data instance) and a Pyomo model. The function then generates plots showing the data and model predictions.
results_summary summarizes the Pyomo.DoE results. We’ll use this later in the workshop.
reformulate_parameters reformulates the original parameters.
recover_original_parameters calculates the original parameters from the reformulated ones.
recover_original_covariance computes the covariance matrix of the original parameters from that of the reformulated parameters

import pyomo.contrib.parmest.parmest as parmest

# Set logging level to ERROR to suppress solver output and warnings from parmest
logging.basicConfig(level=logging.ERROR, force=True)

parmest_logger = logging.getLogger("pyomo.contrib.parmest.parmest")
parmest_logger.setLevel(logging.ERROR)

# Solver options used for all parmest estimation problems in this notebook
solver_options = {"linear_solver": "ma57", "max_iter": 1000, "max_cpu_time": 30}

# First, we define an Experiment object within parmest
#
# Hint: when calling TC_lab_experiment, set reparam=True to
# automatically reformulate the problem for better estimation performance
#
# Add your solution here

# Since everything has been labeled properly in the Experiment object, we simply invoke
# parmest's Estimator function to estimate the parameters.
# Add your solution here

parmest_regression_results = extract_plot_results(
    tc_data, pest.ef_instance.exp_scenarios[0], reparam=True
)

# Now that we have estimated the reformulated parameters, we
# need to transform them to the original parameters
theta_orig = recover_original_parameters(
    theta,
    alpha=pest.ef_instance.exp_scenarios[0].alpha,
    P1=pest.ef_instance.exp_scenarios[0].P1,
)

print("The original model parameters are:")
print("Ua =", round(theta_orig["Ua"], 4), "Watts/°C")
print("Ub =", round(theta_orig["Ub"], 4), "Watts/°C")
print("inv_CpH =", round(theta_orig["inv_CpH"], 4), "°C/Joules")
print("inv_CpS =", round(theta_orig["inv_CpS"], 4), "°C/Joules")

Discussion: How do these results compare to our previous analysis? Discuss this in a few sentences.

Quantify the Uncertainty in the Parameter Estimates¶

As mentioned in the uncertainty quantification section of the main workshop, covariance matrix can be used to measure the accuracy of parameter estimates. This is needed to check how close the parameter estimates are to their true values. The leading diagonal of the covariance matrix contains the variance of the estimated parameters.

# Since everything has been labeled properly in the Experiment object, we
# simply use the `cov_est` function to compute the covariance matrix.
# Add your solution here

# Lets check the covariance matrix of the reformulated parameters
# Add a print statement to check the covariance matrix
# Add your solution here

# since the covariance matrix is that of the reformulated parameters, we
# need to transform it to the original parameters
cov_orig = recover_original_covariance(
    theta,
    cov,
    alpha=pest.ef_instance.exp_scenarios[0].alpha,
    P1=pest.ef_instance.exp_scenarios[0].P1,
)

print("The covariance matrix of the original parameters is:\n", cov_orig)
print("\nThe trace of the covariance matrix is:", format(np.trace(cov_orig), ".3e"))

Discussion: How do these results compare to our previous analysis? Discuss this in a few sentences.

Multistart Optimization with ParmEst¶

Use multistart optimization with sobol sampling on the reformulated model.

pest_sobol = parmest.Estimator(
    [TC_Lab_sine_exp], obj_function="SSE", tee=True, solver_options=solver_options
)


# Set some common options for all multistart estimation runs
common_multistart_options = {"n_restarts": 15, "seed": 532, "save_results": False}

# Add your solution here

# Analyze results
print("Best parameter estimates from Sobol sampling multistart:")
print(best_theta_sobol)
print(f"Best objective value from Sobol sampling multistart: {best_obj_sobol}")

# Round the objective values to a reasonable number of decimal places for counting unique minima
results_df_sobol["final objective"] = results_df_sobol["final objective"].round(5)

num_unique_minima_sobol = len(results_df_sobol["final objective"].unique())
print(f"Number of unique minima found with Sobol sampling: {num_unique_minima_sobol}")

# Recalculate the original parameters from the best Sobol sampling result
orig_theta_sobol = recover_original_parameters(
    best_theta_sobol,
    alpha=pest.ef_instance.exp_scenarios[0].alpha,
    P1=pest.ef_instance.exp_scenarios[0].P1,
)

print("Estimated parameters:")
print(f"Ua: {orig_theta_sobol['Ua']:.6f} W/K")
print(f"Ub: {orig_theta_sobol['Ub']:.6f} W/K")
print(f"inv_CpH: {orig_theta_sobol['inv_CpH']:.6f} J/(K*kg)")
print(f"inv_CpS: {orig_theta_sobol['inv_CpS']:.6f} J/(K*kg)")

Profile Likelihood with ParmEst¶

Analyze the profile likeihood of the reformulated model.

# Add your solution here

profiles = profile_results["profiles"]
print("Profile likelihood results:")
profiles.head(5)

# Plot profile curves to visualize profile likelihood for each parameter.

# Add your solution here

How does the estimability of the reformulated model compare to the original model?

Add L2 Regularization to Parameter Estimation Objective¶

Using the same prior from the regularization notebook, converted to the reformulated model parameters, add regularization and compare the optimal solution.

# Original prior:

# ---- Physically intuitive guesses (Cp-space) ----
theta_phys = pd.Series(
    {"Ua": 0.030, "Ub": 0.018, "inv_CpH": 1 / 7.5, "inv_CpS": 1 / 0.22}
)


# Transform to estimator parameterization [beta-space]
theta0_phys_reparam = reformulate_parameters(
    theta_phys,
    alpha=pest.ef_instance.exp_scenarios[0].alpha,
    P1=pest.ef_instance.exp_scenarios[0].P1,
)

# Define diagonal covariance matrix
cov_x = pd.DataFrame(
    np.diag([0.02, 0.01, 0.05, 0.05]),
    index=["beta_1", "beta_2", "beta_3", "beta_4"],
    columns=["beta_1", "beta_2", "beta_3", "beta_4"],
)

# Invert to get the physically informed prior_FIM
prior_FIM_phys = pd.DataFrame(
    np.linalg.inv(cov_x.values), index=cov_x.index, columns=cov_x.columns
)

# Optional scaling factor to tune regularization strength
prior_weight = 2
prior_FIM_phys = prior_weight * prior_FIM_phys


print("theta0_phys_reparam:", theta0_phys_reparam)
print("prior_FIM_phys:\n", prior_FIM_phys)

# Add your solution here

# Compare to unregularized estimation results and multistart results
print("\nUnregularized objective:", obj)
print("Unregularized theta:\n", theta)

print("\nBest Sobol multistart objective:", best_obj_sobol)
print(
    "Best Sobol multistart theta:\n",
    pd.DataFrame(best_theta_sobol, index=["best_theta"]).T,
)

print("\nL2 (physical prior) objective:", obj_phys)
print("L2 (physical prior) theta:\n", theta_phys_est)