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 = 2Load 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()Loading...
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
### BEGIN SOLUTION
results_summary(PRIOR_FIM, reparam=True)
### END SOLUTION======Results Summary======
Five design criteria log10() value:
Pseudo A-optimality: 9.653309012938479
A-optimality: -3.7180317015494597
D-optimality: 25.766913035762006
E-optimality: 3.759129284988475
Modified E-optimality: 5.890357174341899
FIM:
[[ 2.5458e+09 -8.4425e+08 1.2854e+09 1.5849e+09]
[-8.4425e+08 3.1042e+08 -4.0976e+08 -5.2518e+08]
[ 1.2854e+09 -4.0976e+08 6.5803e+08 8.0051e+08]
[ 1.5849e+09 -5.2518e+08 8.0051e+08 9.8675e+08]]
eigenvalues:
[4.46155714e+09 3.93791714e+07 5.79469681e+04 5.74287396e+03]
Eigenvector matrix:
eigvec_1 eigvec_2 eigvec_3 eigvec_4
beta_1 0.7554 0.0245 -0.6547 0.0137
beta_2 -0.2508 0.8710 -0.2496 0.3409
beta_3 0.3813 0.4899 0.4447 -0.6456
beta_4 0.4703 0.0279 0.5580 0.6832
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
### BEGIN SOLUTION
solver = SolverFactory("ipopt")
solver.options["max_iter"] = 3000
solver.options["tol"] = 1e-5
solver.options["linear_solver"] = "ma57"
solver.options["nlp_scaling_method"] = "gradient-based"
### END SOLUTIONfrom 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.
### BEGIN SOLUTION
theta_ref = np.array(
[
theta_values["beta_1"],
theta_values["beta_2"],
theta_values["beta_3"],
theta_values["beta_4"],
]
)
SCALE_BY_NOMINAL_PARAM = True
SCALE_CONSTANT_VALUE = 1e-2
if SCALE_BY_NOMINAL_PARAM:
PRIOR_FIM_SCALED = rescale_FIM(PRIOR_FIM, theta_ref) * SCALE_CONSTANT_VALUE**2
else:
PRIOR_FIM_SCALED = PRIOR_FIM * SCALE_CONSTANT_VALUE**2
### END SOLUTIONfrom pyomo.contrib.doe import DesignOfExperiments
# Create experiment object for design of experiments
### BEGIN SOLUTION
doe_experiment = TC_Lab_experiment(
data=tc_data,
theta_initial=theta_values,
number_of_states=number_tclab_states,
reparam=True,
)
### END SOLUTION
# Create the design of experiments object using our experiment instance from above
### BEGIN SOLUTION
TC_Lab_DoE_pA = DesignOfExperiments(
experiment=doe_experiment,
step=1e-2,
scale_constant_value=SCALE_CONSTANT_VALUE,
scale_nominal_param_value=SCALE_BY_NOMINAL_PARAM,
objective_option="pseudo_trace", # Now we specify a type of objective, A-opt = "trace"
prior_FIM=PRIOR_FIM_SCALED, # We use the prior information from the same existing experiment as in the D-optimal case!
tee=True,
solver=solver,
)
### END SOLUTION
# Run DoE analysis
### BEGIN SOLUTION
TC_Lab_DoE_pA.run_doe()
### END SOLUTIONFetching long content....
# Extract and plot the results using our custom function
### BEGIN SOLUTION
pAopt_pyomo_doe_results = extract_plot_results(
None, TC_Lab_DoE_pA.model.fd_scenario_blocks[0], reparam=True
)
### END SOLUTION
Model parameters:
Beta_1 = 0.0134 Watts/Joules
Beta_2 = 0.0216 Watts/Joules
Beta_3 = 0.0287 Watts/Joules
Beta_4 = 0.0102 °C.Watts/(Joules.%)
# Compute the FIM at the optimal solution and don't forget to unscale it to get
# the correct values for the summary and eigendecomposition
### BEGIN SOLUTION
FIM_pA = np.asarray(TC_Lab_DoE_pA.results["FIM"])
if SCALE_BY_NOMINAL_PARAM:
D = np.diag(theta_ref)
FIM_unscaled = D @ FIM_pA @ D / SCALE_CONSTANT_VALUE**2
else:
FIM_unscaled = FIM_pA / SCALE_CONSTANT_VALUE**2
results_summary(np.asarray(FIM_unscaled))
### END SOLUTION======Results Summary======
Five design criteria log10() value:
Pseudo A-optimality: 9.653309035847789
A-optimality: -3.719818051423565
D-optimality: 25.770219332375724
E-optimality: 3.760945246419933
Modified E-optimality: 5.888541214259059
FIM:
[[ 2.54580012e+09 -8.44249988e+08 1.28539998e+09 1.58489989e+09]
[-8.44249988e+08 3.10420003e+08 -4.09760007e+08 -5.25180013e+08]
[ 1.28539998e+09 -4.09760007e+08 6.58030017e+08 8.00510027e+08]
[ 1.58489989e+09 -5.25180013e+08 8.00510027e+08 9.86750095e+08]]
eigenvalues:
[4.46155715e+09 3.93791717e+07 5.81461613e+04 5.76693752e+03]
Eigenvector matrix:
eigvec_1 eigvec_2 eigvec_3 eigvec_4
Ua 0.7554 0.0245 -0.6547 0.0145
Ub -0.2508 0.8710 -0.2492 0.3412
inv_CpH 0.3813 0.4899 0.4440 -0.6461
inv_CpS 0.4703 0.0279 0.5588 0.6825
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
### BEGIN SOLUTION
doe_experiment = TC_Lab_experiment(
data=tc_data,
theta_initial=theta_values,
number_of_states=number_tclab_states,
reparam=True,
)
### END SOLUTION
# Create the design of experiments object using our experiment instance from above
### BEGIN SOLUTION
TC_Lab_DoE_D = DesignOfExperiments(
experiment=doe_experiment,
step=1e-2,
scale_constant_value=SCALE_CONSTANT_VALUE,
scale_nominal_param_value=SCALE_BY_NOMINAL_PARAM,
objective_option="determinant", # D-opt = "determinant"
prior_FIM=PRIOR_FIM_SCALED, # prior information from the existing experiment!
tee=True,
solver=solver,
)
### END SOLUTION
### BEGIN SOLUTION
TC_Lab_DoE_D.run_doe()
### END SOLUTIONFetching long content....
# Extract and plot the results using our custom function
### BEGIN SOLUTION
dopt_pyomo_doe_results = extract_plot_results(
None, TC_Lab_DoE_D.model.fd_scenario_blocks[0], reparam=True
)
### END SOLUTION
Model parameters:
Beta_1 = 0.0134 Watts/Joules
Beta_2 = 0.0216 Watts/Joules
Beta_3 = 0.0287 Watts/Joules
Beta_4 = 0.0102 °C.Watts/(Joules.%)
# Extract the FIM at the optimal solution
# Don't forget to unscale it to get the correct values for the summary and eigendecomposition
### BEGIN SOLUTION
FIM_D = np.asarray(TC_Lab_DoE_D.results["FIM"])
if SCALE_BY_NOMINAL_PARAM:
D = np.diag(theta_ref)
FIM_unscaled_D = D @ FIM_D @ D / SCALE_CONSTANT_VALUE**2
else:
FIM_unscaled_D = FIM_D / SCALE_CONSTANT_VALUE**2
results_summary(np.asarray(FIM_unscaled_D))
### END SOLUTION======Results Summary======
Five design criteria log10() value:
Pseudo A-optimality: 9.653309111948408
A-optimality: -3.727219189374622
D-optimality: 25.782395202854836
E-optimality: 3.7686379681597773
Modified E-optimality: 5.8808485003811155
FIM:
[[ 2.54580053e+09 -8.44249964e+08 1.28539994e+09 1.58489956e+09]
[-8.44249964e+08 3.10420019e+08 -4.09760041e+08 -5.25180049e+08]
[ 1.28539994e+09 -4.09760041e+08 6.58030093e+08 8.00510092e+08]
[ 1.58489956e+09 -5.25180049e+08 8.00510092e+08 9.86750386e+08]]
eigenvalues:
[4.46155723e+09 3.93791732e+07 5.87494997e+04 5.86999820e+03]
Eigenvector matrix:
eigvec_1 eigvec_2 eigvec_3 eigvec_4
Ua 0.7554 0.0245 -0.6546 0.0169
Ub -0.2508 0.8710 -0.2479 0.3421
inv_CpH 0.3813 0.4899 0.4416 -0.6478
inv_CpS 0.4703 0.0279 0.5613 0.6805
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
### BEGIN SOLUTION
doe_experiment = TC_Lab_experiment(
data=tc_data,
theta_initial=theta_values,
number_of_states=number_tclab_states,
reparam=True,
)
### END SOLUTION
# Create the design of experiments object using our experiment instance from above
### BEGIN SOLUTION
TC_Lab_DoE_ME = DesignOfExperiments(
experiment=doe_experiment,
step=1e-2,
scale_constant_value=1,
scale_nominal_param_value=True,
objective_option="condition_number", # Now we specify a type of objective, ME-opt = "condition_number"
prior_FIM=PRIOR_FIM, # We use the prior information from the existing experiment!
tee=True,
)
### END SOLUTION
### BEGIN SOLUTION
TC_Lab_DoE_ME.use_grey_box = True # Use greybox for ME
TC_Lab_DoE_ME.run_doe()
### END SOLUTIONFetching long content....
# Extract and plot the results using our custom function
### BEGIN SOLUTION
meopt_pyomo_doe_results = extract_plot_results(
None, TC_Lab_DoE_ME.model.fd_scenario_blocks[0], reparam=True
)
### END SOLUTION
Model parameters:
Beta_1 = 0.0134 Watts/Joules
Beta_2 = 0.0216 Watts/Joules
Beta_3 = 0.0287 Watts/Joules
Beta_4 = 0.0102 °C.Watts/(Joules.%)
# 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
### BEGIN SOLUTION
SCALE_CONSTANT_VALUE = 1
FIM_ME = np.asarray(TC_Lab_DoE_ME.results['FIM'])
if SCALE_BY_NOMINAL_PARAM:
D = np.diag(theta_ref)
FIM_unscaled_ME = D @ FIM_ME @ D / SCALE_CONSTANT_VALUE**2
else:
FIM_unscaled_ME = FIM_ME / SCALE_CONSTANT_VALUE**2
results_summary(np.asarray(FIM_unscaled_ME))
### END SOLUTION======Results Summary======
Five design criteria log10() value:
Pseudo A-optimality: 6.093661307557583
A-optimality: -1.3294323156659942
D-optimality: 14.586640923832137
E-optimality: 1.3438521070122227
Modified E-optimality: 4.742053255059707
FIM:
[[ 450725.85179762 -242518.7309079 490633.91261288 214689.74486591]
[-242518.7309079 144857.76130071 -254079.92230317 -115752.36413873]
[ 490633.91261288 -254079.92230317 542151.27130027 234423.60777105]
[ 214689.74486591 -115752.36413873 234423.60777105 102949.47556696]]
eigenvalues:
[1.21872400e+06 2.12633735e+04 6.74918430e+02 2.20725296e+01]
Eigenvector matrix:
eigvec_1 eigvec_2 eigvec_3 eigvec_4
Ua 0.6078 -0.0782 0.6684 0.4215
Ub -0.3256 0.8584 0.1709 0.3578
inv_CpH 0.6636 0.5061 -0.2280 -0.5015
inv_CpS 0.2902 -0.0305 -0.6870 0.6655
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:
where
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 = 2Load 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()Loading...
# 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
### BEGIN SOLUTION
results_summary(PRIOR_FIM, reparam=True)
### END SOLUTION======Results Summary======
Five design criteria log10() value:
Pseudo A-optimality: 9.653309012938479
A-optimality: -3.7180317015494597
D-optimality: 25.766913035762006
E-optimality: 3.759129284988475
Modified E-optimality: 5.890357174341899
FIM:
[[ 2.5458e+09 -8.4425e+08 1.2854e+09 1.5849e+09]
[-8.4425e+08 3.1042e+08 -4.0976e+08 -5.2518e+08]
[ 1.2854e+09 -4.0976e+08 6.5803e+08 8.0051e+08]
[ 1.5849e+09 -5.2518e+08 8.0051e+08 9.8675e+08]]
eigenvalues:
[4.46155714e+09 3.93791714e+07 5.79469681e+04 5.74287396e+03]
Eigenvector matrix:
eigvec_1 eigvec_2 eigvec_3 eigvec_4
beta_1 0.7554 0.0245 -0.6547 0.0137
beta_2 -0.2508 0.8710 -0.2496 0.3409
beta_3 0.3813 0.4899 0.4447 -0.6456
beta_4 0.4703 0.0279 0.5580 0.6832
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.
### BEGIN SOLUTION
theta_ref = np.array(
[
theta_values["beta_1"],
theta_values["beta_2"],
theta_values["beta_3"],
theta_values["beta_4"],
]
)
# Because the parameter values differ substantially in magnitude, we can use parameter
# scaling in DesignOfExperiments to improve numerical conditioning.
SCALE_BY_NOMINAL_PARAM = False
SCALE_CONSTANT_VALUE = 1e-2
if SCALE_BY_NOMINAL_PARAM:
PRIOR_FIM_SCALED = rescale_FIM(PRIOR_FIM, theta_ref) * SCALE_CONSTANT_VALUE**2
else:
PRIOR_FIM_SCALED = PRIOR_FIM * SCALE_CONSTANT_VALUE**2
### END SOLUTION# Experiment object 1
# Create experiment object for design of experiments and set reparam=True to use the
# reparameterized version of the model
### BEGIN SOLUTION
doe_experiment_1 = TC_Lab_experiment(
data=tc_data_1,
theta_initial=theta_values,
number_of_states=number_tclab_states,
reparam=True,
)
### END SOLUTIONfrom dataclasses import replace
# Experiment object 2
# Choose between "normal" or "uniform" distribution for the new design of u1
### BEGIN SOLUTION
SEED: int = 11 # Choose any integer seed you like
PDF = "normal"
def get_initial_u1_design(PDF="uniform", seed=11, sample_data=tc_data_1):
rng = np.random.default_rng(seed=seed)
if PDF == "normal":
u1_design = rng.normal(loc=50.0, scale=10.0, size=len(sample_data.time))
# Normal distribution is unbounded, so we need to clip the values to be within the
# bounds of the control input, which is [0, 100] in this case.
u1_design = np.clip(u1_design, 0.0, 100.0)
elif PDF == "uniform":
u1_design = rng.uniform(low=0.0, high=100.0, size=len(sample_data.time))
# To break permutation symmetry, we can enforce that the first value of u1_design
# is greater than the first value of the original u1 data.
if u1_design[0] <= sample_data.u1[0] and sample_data.u1[0] < 90.0:
u1_design[0] = sample_data.u1[0] + 10.0
else:
u1_design[0] = 100.0
# Create a new data object with only u1 replaced
tc_data_new = replace(sample_data, u1=u1_design)
return tc_data_new
tc_data_2 = get_initial_u1_design(PDF=PDF, seed=SEED, sample_data=tc_data_1)
### END SOLUTION
# Build the experiment 2 with the new design variable and set reparam=True
### BEGIN SOLUTION
doe_experiment_2 = TC_Lab_experiment(
data=tc_data_2,
theta_initial=theta_values,
number_of_states=number_tclab_states,
reparam=True,
)
### END SOLUTIONOptimize 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
### BEGIN SOLUTION
TC_Lab_DoE_D = DesignOfExperiments(
experiment=[doe_experiment_1, doe_experiment_2],
# We are optimizing two experiments simultaneously!
step=1e-2,
scale_constant_value=SCALE_CONSTANT_VALUE,
scale_nominal_param_value=SCALE_BY_NOMINAL_PARAM,
objective_option="determinant",
# Now we specify a type of objective, D-opt = "determinant"
prior_FIM=PRIOR_FIM,
tee=True,
solver=solver,
)
### END SOLUTION# Optimize two experiments simultaneously
### BEGIN SOLUTION
TC_Lab_DoE_D.optimize_experiments()
### END SOLUTIONFetching long content....
# 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
### BEGIN SOLUTION
FIM_D = np.asarray(TC_Lab_DoE_D.results["solution"]["param_scenarios"][0]["total_fim"])
if SCALE_BY_NOMINAL_PARAM:
D = np.diag(theta_ref)
FIM_unscaled_mexp = D @ FIM_D @ D / SCALE_CONSTANT_VALUE**2
else:
FIM_unscaled_mexp = FIM_D / SCALE_CONSTANT_VALUE**2
results_summary(np.asarray(FIM_unscaled_mexp))
### END SOLUTION======Results Summary======
Five design criteria log10() value:
Pseudo A-optimality: 13.653933167405627
A-optimality: -8.31302292779217
D-optimality: 44.36934452682199
E-optimality: 8.314637798889118
Modified E-optimality: 5.334854546024938
FIM:
[[ 2.54780124e+13 -8.44179783e+12 1.28532234e+13 1.58196622e+13]
[-8.44179783e+12 3.10428128e+12 -4.09769983e+12 -5.25311912e+12]
[ 1.28532234e+13 -4.09769983e+12 6.58042604e+12 8.00661788e+12]
[ 1.58196622e+13 -5.25311912e+12 8.00661788e+12 9.91201373e+12]]
eigenvalues:
[4.46161760e+13 2.06365835e+08 6.45573303e+10 3.93793752e+11]
Eigenvector matrix:
eigvec_1 eigvec_2 eigvec_3 eigvec_4
Ua 0.7554 -0.3603 -0.5468 0.0242
Ub -0.2508 -0.4215 -0.0301 0.8710
inv_CpH 0.3813 0.7833 0.0323 0.4899
inv_CpS 0.4703 -0.2811 0.8361 0.0282
# Create a helper function to compute FIM metrics
### BEGIN SOLUTION
def log10_FIM_metric(FIM):
pseudo_A_opt = np.log10(np.trace(FIM))
A_opt = np.log10(np.trace(np.linalg.inv(FIM)))
D_opt = np.log10(np.linalg.det(FIM))
E_opt = np.log10(np.min(np.linalg.eigvals(FIM)))
ME_opt = np.log10(np.max(np.linalg.eigvals(FIM)) / np.min(np.linalg.eigvals(FIM)))
return pseudo_A_opt, A_opt, D_opt, E_opt, ME_opt
### END SOLUTION# 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.
### BEGIN SOLUTION
import pandas as pd
prior_metrics = log10_FIM_metric(PRIOR_FIM)
optimized_metrics_mexp = log10_FIM_metric(FIM_unscaled_mexp)
optimized_metrics_single_exp = log10_FIM_metric(FIM_unscaled_D)
metric_names = ["pseudo_A_opt", "A_opt", "D_opt", "E_opt", "ME_opt"]
metrics_df = pd.DataFrame(
{
"prior": list(prior_metrics),
"optimized_mexp": list(optimized_metrics_mexp),
"optimized_single_exp": list(optimized_metrics_single_exp),
},
index=metric_names,
)
metrics_df.index.name = "metric"
metrics_df.round(2)
### END SOLUTIONLoading...