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. In this notebook, we will use the step test as the starting point.
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 (step test)#
We will load the step 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_step_test.csv"
else:
file = '../data/tclab_step_test.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],
)
Analyze FIM with Pyomo.DoE at initial point (step test)#
To get started, compute and analyze the FIM of the step test experiment.
# Load Pyomo.DoE class
from pyomo.contrib.doe import DesignOfExperiments
from pyomo.environ import SolverFactory
# Copied from previous notebook
theta_values = {
'Ua': 0.0417051733576387,
'Ub': 0.009440714239773074,
'inv_CpH': 0.1659093525658045,
'inv_CpS': 5.8357556063605465,
}
# 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
# Calculate the FIM
# Add your solution here
# Call our custom function to summarize the results
# and compute the eigendecomposition of the FIM
# Add your solution here
Discussion: How does this FIM compare to the sine test experiment we previously analyzed?
Optimize the next experiment (A-optimality)#
Now we are ready to compute the A-optimal next best experiment. Why are we starting with A-optimality? It runs faster so it is better for debugging syntax.
# 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
# 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
# Compute the FIM at the optimal solution
# Add your solution here
Discussion: How do these results compare to our previous analysis?