Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Exercise: Pyomo.DoE

In this notebook, you will use Pyomo.DoE to compute the A- and D-optimal experiments from the TCLab. In our previous notebook, we used the sine test as a starting point. We will do the same here as well.

Recall, we can computing the next best experiment assuming we already completed one prior experiment. Thus it is important to confirm our optimal experiment design does not change if we change the prior experiment of optimization initial point.

import sys

# 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,
    results_summary,
)

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

Load and Explore Experimental Data (sine test)

We will load the sine test experimental data, similar to our previous notebooks.

import pandas as pd

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.head()
ax = df.plot(x="Time", y=["T1", "T2"], xlabel="Time (s)", ylabel="Temperature (°C)")
ax.grid(True)
ax = df.plot(x="Time", y=["Q1", "Q2"], xlabel="Time (s)", ylabel="Heater Power (%)")
ax.grid(True)
# Here, we will induce a step size of 10 seconds, as to not give too many
# degrees of freedom for experimental design.
skip = 10

# Create the data object considering the new control points every 6 seconds
tc_data = TC_Lab_data(
    name="Sine Wave Test for Heater 1",
    time=df["Time"].values[::skip],
    T1=df["T1"].values[::skip],
    u1=df["Q1"].values[::skip],
    P1=200,
    TS1_data=None,
    T2=df["T2"].values[::skip],
    u2=df["Q2"].values[::skip],
    P2=200,
    TS2_data=None,
    Tamb=df["T1"].values[0],
)

Use Prior Parameter Information

You will use the parameter estimate and the Fisher information matrix (FIM) from previous parmest exercise notebook to pass it as the prior FIM in the DesignOfExperiments object.

import numpy as np

# Parameter values estimated from the regularized regression in the previous parmest
# exercise notebook
theta_values = {
    "beta_1": 0.0133,  # Watts/Joules
    "beta_2": 0.0216,  # Watts/Joules
    "beta_3": 0.0287,  # Watts/Joules
    "beta_4": 0.0102,  # °C.Watts/(Joules.%)
}

# Prior FIM from the notebook (FIM = covariance matrix)^(-1))
PRIOR_FIM = np.array(
    [
        [2.5458e09, -8.4425e08, 1.2854e09, 1.5849e09],
        [-8.4425e08, 3.1042e08, -4.0976e08, -5.2518e08],
        [1.2854e09, -4.0976e08, 6.5803e08, 8.0051e08],
        [1.5849e09, -5.2518e08, 8.0051e08, 9.8675e08],
    ]
)
# Call our custom function to summarize the results
# and compute the eigendecomposition of the FIM
# Set `reparam=True` in this function since we are using the reparameterized model

# Add your solution here

Optimize the Next Experiment (Pseudo-A-Optimality)

Now we are ready to compute the pseudo-A-optimal next best experiment. Why are we starting with pseudo-A-optimality? It runs faster so it is better for debugging syntax.

# Add a solver object to pass to DesignOfExperiments
from pyomo.environ import SolverFactory

# Add your solution here
from pyomo.contrib.doe.utils import rescale_FIM

# Scale the FIM using the nominal parameter values as reference to improve numerical conditioning
# Also, we can apply an additional constant scaling factor to further improve conditioning if needed.
# Add your solution here
from pyomo.contrib.doe import DesignOfExperiments

# Create experiment object for design of experiments
# Add your solution here

# Create the design of experiments object using our experiment instance from above
# Add your solution here

# Run DoE analysis
# Add your solution here
# Extract and plot the results using our custom function
# Add your solution here
# Compute the FIM at the optimal solution and don't forget to unscale it to get
# the correct values for the summary and eigendecomposition
# Add your solution here

Discussion: How do these compare to our previous A-optimal results considering the sine test as the prior experiment?

Optimize the Next Experiment (D-Optimality)

Finally, we are ready to solve the D-optimality problem. This may take 2 minutes to run.

# Create experiment object for design of experiments
# Add your solution here

# Create the design of experiments object using our experiment instance from above
# Add your solution here

# Add your solution here
# Extract and plot the results using our custom function
# Add your solution here
# Extract the FIM at the optimal solution
# Don't forget to unscale it to get the correct values for the summary and eigendecomposition
# Add your solution here

Discussion: How do these results compare to our previous analysis?

Optimize the Next Experiment (ME-Optimality)

Now we look at the modified E-optimality

# Create experiment object for design of experiments
# Add your solution here

# Create the design of experiments object using our experiment instance from above
# Add your solution here

# Add your solution here
# Extract and plot the results using our custom function
# Add your solution here
# Extract the FIM at the optimal solution
# As before, don't forget to unscale it to get the correct values for the summary and eigendecomposition
# Add your solution here

Discussion: How do these results compare to our previous analysis?

Exercise: Multi-experiment design

In this section, you will use the reparameterized tc lab model which you used in the previous parmest exercise and design two experiments using Pyomo.DoE. The structure of the reformulated model is:

dTHdt=β1(Tamb−TH)+β2(TS−TH)+β4u(t),\frac{dT_H}{dt} = \beta_1 (T_{\text{amb}} - T_H) + \beta_2 (T_S - T_H) + \beta_4 u(t),
(1)
dTSdt=β3(TH−TS)\frac{dT_S}{dt} = \beta_3 (T_H - T_S)
(2)

where

β1=UaCpH,β2=UbCpH,β3=UbCpS,β4=αPCpH\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}
(3)

You will use D-optimality (determinant of the FIM) to design the two experiments.

Import necessary functions for the two-state TC Lab model

import sys

# 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_plot_results,
    results_summary,
)

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

Load Experimental Data (sine test)

We will load the sine test experimental data to serve as an initial point.

import pandas as pd

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.head()
# Here, we will induce a step size of 15 seconds, as to make the optimization problem
# more tractable. You can change this value to see how it affects the results
skip = 20

# Create the data object considering the new control points
tc_data_1 = TC_Lab_data(
    name="Sine Wave Test for Heater 1",
    time=df["Time"].values[::skip],
    T1=df["T1"].values[::skip],
    u1=df["Q1"].values[::skip],
    P1=200,
    TS1_data=None,
    T2=df["T2"].values[::skip],
    u2=df["Q2"].values[::skip],
    P2=200,
    TS2_data=None,
    Tamb=df["T1"].values[0],
)

Use Prior Parameter Information

You will use the parameter estimate and the Fisher information matrix (FIM) from previous parmest exercise notebook to pass it as the prior FIM in the DesignOfExperiments object.

import numpy as np

# Parameter values estimated from the regularized regression in the previous parmest
# exercise notebook
theta_values = {
    "beta_1": 0.0133,  # Watts/Joules
    "beta_2": 0.0216,  # Watts/Joules
    "beta_3": 0.0287,  # Watts/Joules
    "beta_4": 0.0102,  # °C.Watts/(Joules.%)
}

# Prior FIM from the notebook (FIM = covariance matrix)^(-1))
PRIOR_FIM = np.array(
    [
        [2.5458e09, -8.4425e08, 1.2854e09, 1.5849e09],
        [-8.4425e08, 3.1042e08, -4.0976e08, -5.2518e08],
        [1.2854e09, -4.0976e08, 6.5803e08, 8.0051e08],
        [1.5849e09, -5.2518e08, 8.0051e08, 9.8675e08],
    ]
)
# Call our custom function to summarize the results
# and compute the eigendecomposition of the FIM
# Set `reparam=True` in this function since we are using the reparameterized model

# Add your solution here

Initialize Two Experiments

You need to create two experiment objects. The first one uses the measured sine-test data which you have already created, and the second uses a different initial input profile so that the optimizer can optimize both experiments simultaneously.

from pyomo.contrib.doe.utils import rescale_FIM

# Scale the FIM using the nominal parameter values as reference to improve numerical conditioning
# Also, we can apply an additional constant scaling factor to further improve conditioning if needed.

# Add your solution here
# Experiment object 1
# Create experiment object for design of experiments and set reparam=True to use the
# reparameterized version of the model
# Add your solution here
from dataclasses import replace

# Experiment object 2
# Choose between "normal" or "uniform" distribution for the new design of u1
# Add your solution here
# Build the experiment 2 with the new design variable and set reparam=True
# Add your solution here

Optimize Two Experiments (D-Optimality)

Now, you create a DesignOfExperiments object for the two experiments and choose a D-optimality (“determinant”) objective. You also pass prior_FIM, which represents information already available from previous parameter estimation.

# Load Pyomo.DoE class
from pyomo.contrib.doe import DesignOfExperiments

# Create the `DesignOfExperiments` object using your experiment instances from above
# Add your solution here
# Optimize two experiments simultaneously
# Add your solution here
# Extract and plot the results using our custom function
### BEGIN SOLUTION
multiexp_results = extract_plot_results(None, TC_Lab_DoE_D)
###
# Extract the FIM and compute the results summary at the optimal solution
# Don't forget to unscale it to get the correct values for the summary and eigendecomposition
# Add your solution here
# Create a helper function to compute FIM metrics
# Add your solution here
# Compare the how much the new experiment improved the FIM metrics compared to the 
# prior experiment and the single experiment optimization with D-optimal design. 
# Add your solution here
Pyomo.DoE Workshop and Tutorials
Planning Multiple Experiments with Pyomo.DoE
Pyomo.DoE Workshop and Tutorials
Exercise Solutions