Optimizing Experiments with Pyomo.DoE

5.3. Optimizing Experiments with `Pyomo.DoE`#

Prepared by: Jialu Wang, Prof. Alex Dowling, Hailey Lynch, and Andrew Marquardt (amarquar@nd.edu, 2024) at the University of Notre Dame.

5.3.1. Introduction and Learning Objectives#

This notebook uses design of experiments for a reactor kinetics experiment with Pyomo.DoE. The user will be able to learn concepts involved in Model-Based Design of Experiments (MBDoE) and practice using Pyomo.DoE from methodology in the notebook. Results will be interpreted throughout the notebook to connect the material with the Pyomo implementation.

The general process we will follow throughout this notebook:

Import Modules
- Step 0: Import Pyomo and Pyomo.DoE Module
Problem Statement
- Step 1: Mathematical Model for the Reaction Kinetics Example
Implementation in Pyomo
- Step 2: Implement Mathematical Model
- Step 3: Generate an Experiment
- Step 4: Run the DOE Module to get Optimal Design Parameters
Methodology
- Step 5: Method for Computing FIM
- Step 6: Running a Full Factorial Design Experiment
Visualizing Results
- Step 7: Evaluating the Full Factorial Heat Maps
Key Takeaways

5.3.2. Step 0: Import Pyomo and Pyomo.DoE Module#

# IPOPT installer
import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper
    helper.install_idaes()
    helper.install_ipopt()
else:
    sys.path.insert(0, '../')
    import helper

    # Import solver
    import idaes
helper.set_plotting_style()

#F# Imports
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyomo.environ as pyo
import pyomo.common.unittest as unittest
from pyomo.contrib.doe import DesignOfExperiments
from pyomo.dae import ContinuousSet, DerivativeVar

5.3.3. Step 1: Mathematical Model for the Reaction Kinetics Example#

5.3.3.1. Reaction Kinetics#

Consider two chemical reactions that convert molecule $A$ to desired product $B$ and a less valuable side-product $C$:

\[A \overset{k_1}{\rightarrow} B \overset{k_2}{\rightarrow} C\]

Goal:

Design a large-scale continuous reactor that maximizes the production of $B$.

The rate laws for these two chemical reactions are:

\[r_A = -k_1 C_A\]

\[r_B = k_1 C_A - k_2 C_B\]

\[r_C = k_2 C_B\]

Here, $C_A$, $C_B$, and $C_C$ are the concentrations of each species. \

The rate constants $k_1$ and $k_2$ depend on temperature as follows:

\[k_1 = A_1 \exp \left({\frac{-1000 E_1}{R T}} \right)\]

\[k_2 = A_2 \exp \left({\frac{-1000 E_2}{R T}}\right)\]

where:

$A_1$ [$s^{-1}$], $A_2$ [$s^{-1}$], $E_1$ [$kJ/mol$], and $E_2$ [$kJ/mol$] are fitted model parameters.
$R$ [$J/mol K$] is the ideal-gas constant.
$T$ [$K$] is absolute temperature.

Objective:

Using the Pyomo ecosystem, we would like to perform uncertainty quantification and design of experiments on a small-scale batch reactor to infer parameters $A_1$, $A_2$, $E_1$, and $E_2$.

5.3.3.2. Batch Reactor#

The concentrations in a batch reactor evolve with time and are modeled by the following differential equations:

\[ \frac{d C_A}{dt} = r_A = -k_1 C_A(t) \]

\[ \frac{d C_B}{dt} = r_B = k_1 C_A(t) - k_2 C_B(t) \]

\[ \frac{d C_C}{dt} = r_C = k_2 C_B(t) \]

We have now established a linear system of differential equations. Next we can write our initial conditions where we assume the feed is only species $A$ such that:

\[C_A(t=0) = C_{A0}, \quad C_B(t=0) = 0, \quad C_C(t=0) = 0\]

When $k_1$ and $k_2$ are at constant temperature, we have the following analytic solution:

\[C_A(t) = C_{A0} \exp(-k_1 t)\]

\[C_B(t) = \frac{k_1}{k_2 - k_1} C_{A0} \left[\exp(-k_1 t) - \exp(-k_2 t) \right]\]

\[C_C(t) = C_{A0} - \frac{k_2}{k_2 - k_1} C_{A0} \exp(-k_1 t) + \frac{k_1}{k_2 - k_1} \exp(-k_2 t) C_{A0} = C_{A0} - C_{A}(t) - C_{B}(t)\]

See the following for more information on batch reactors:

Batch Reactor: AIChE / Wikipedia

5.3.4. Step 2: Implement Mathematical Model#

This mathematical model is comprised of a system of dynamic algebraic equations (DAEs). This system will be solved using Pyomo.DAE.

See the following notebooks from CBE 60499 regarding Pyomo.DAE:

The experiment class build design contains 5 functions:

the class init__() containing assignments to the input values the module needs. These can be specified different than the given defaults in the experiment build call; otherwise, the defaults will be saved as class variables.
create_model() containing the initial essential build scaffolding for the model. The time, data collection points, concentrations (both initial and time-dependent), kinetic expressions, and mass balances are given here. Note the concentration differentials are grouped together in a conditional function; this can also be done with separate expressions if that is more intuitive.
finalize_model() containing essential fixes of the arrhenius parameters (and other values), discretization, and time bin behavior (namely that the end points of each data collection period are continuous- analogous to the technique used on splines).
label_experiment() containing the necessary setup for the experiment to recognize inputs, outputs, measurement errors (this can be constant or functional), and the arrhenius parameters meant to be tested. Note that each section is necessary- the experiment will not run if inputs, outputs, errors, and unknowns are not declared.
get_labeled_model() containing calls to the previous 3 functions and a return of the final, completed model. This is what is used by the DOE module- all parts must be properly built previously and then collected by this master function. Practically, this function doesn’t have any real effect on the model except to combine everything done in the discrete sections.

Note all the functions are in one code block below- the notebook struggled to recognize the class components if this was not the case.

class Reactor_Experiment():
    def __init__(self,t_control=[0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1],control_val=[500,300,300,300,300,300,300,300,300],t_range=[0,1],CA_init=5,C_init=[5,0,0],theta_pe={'A1': 85,'A2': 370,'E1': 8,'E2': 15,},NFE=32,ncp=3):
        """
        Arguments:
        ---------
        t_control: time-dependent design (control) variables, a list of control timepoints
        control_val: control design variable values T at corresponding timepoints
        t_range: time range, h
        CA_init: time-independent design (control) variable, an initial value for CA
        C_init: An initial value for C
        theta_pe: optimized parmest parameters for the arrhenius kinetic equations
        NFE: number of finite elements in the discretizer
        ncp: number of collocation points per finite element in the discretizer
        """
        self.t_range = t_range
        self.CA_init = CA_init
        self.C_init = C_init
        self.NFE = NFE
        self.ncp = ncp
        self.t_control = t_control
        self.control_val = control_val
        self.theta_pe = theta_pe
        self.model = None

    def create_model(self):
        # model build
        m = self.model = pyo.ConcreteModel()

        # concentration intialization
        m.CA_init = self.CA_init

        # Time-independent design variable
        m.t0 = pyo.Set(initialize=[0])
        m.CA0 = pyo.Var(m.t0, initialize=self.CA_init, bounds=(0, 5.0), within=pyo.NonNegativeReals)  # mol/L

        # parameters
        m.R = pyo.Param(mutable=False,initialize=8.314)

        # variables
        m.t = ContinuousSet(bounds=(self.t_range[0], self.t_range[1]))

        # Parameter list
        para_list = ['A1', 'A2', 'E1', 'E2']
        m.para_list = para_list
        # Define parameters as Param
        m.A1 = pyo.Var(within=pyo.NonNegativeReals)
        m.A2 = pyo.Var(within=pyo.NonNegativeReals)
        m.E1 = pyo.Var(within=pyo.NonNegativeReals)
        m.E2 = pyo.Var(within=pyo.NonNegativeReals)

        # Concentration variables under perturbation
        m.C_set = pyo.Set(initialize=['CA', 'CB', 'CC'])
        m.CA = pyo.Var(m.t, initialize=self.C_init[0], within=pyo.NonNegativeReals)
        m.CB = pyo.Var(m.t, initialize=self.C_init[1], within=pyo.NonNegativeReals)
        m.CC = pyo.Var(m.t, initialize=self.C_init[2], within=pyo.NonNegativeReals)

        # Time derivatives
        m.dCAdt = DerivativeVar(m.CA, wrt=m.t)
        m.dCBdt = DerivativeVar(m.CB, wrt=m.t)
        m.dCCdt = DerivativeVar(m.CC, wrt=m.t)

        # Time-dependent design variable, initialized with the first control value
        m.t_control = self.t_control
        # Control time points
        m.t_con = pyo.Set(initialize=m.t_control)
        # Controls
        controls = {}
        for i, t in enumerate(self.t_control):
            controls[t] = self.control_val[i]
        def T_initial(m, t):
            if t in m.t_con:
                return controls[t]
            else:
                # count how many control points are before the current t;
                # locate the nearest neighbouring control point before this t
                j = -1
                for t_con in m.t_con:
                    if t > t_con:
                        j += 1
                neighbour_t = m.t_control[j]
                return controls[neighbour_t]
        m.T = pyo.Var(m.t, initialize=T_initial, bounds=(300, 700), within=pyo.NonNegativeReals)

        m.kp1 = pyo.Var(m.t, within=pyo.NonNegativeReals)
        m.kp2 = pyo.Var(m.t, within=pyo.NonNegativeReals)
        @m.Constraint(m.t)
        def kp1_cons(m,t):
            return m.kp1[t] == m.A1 * pyo.exp(-m.E1 * 1000 / (m.R * m.T[t]))
        @m.Constraint(m.t)
        def kp2_cons(m,t):
            return m.kp2[t] == m.A2 * pyo.exp(-m.E2 * 1000 / (m.R * m.T[t]))

        @m.Constraint(m.C_set,m.t)
        def dCdt_rule(m, y, t):
            """
            Calculate CA in Jacobian matrix analytically

            Arguments:
            z: scenario No.
            y: CA, CB, CC
            t: timepoints

            Return:
            m: Pyomo model
            """
            if y == 'CA':
                return m.dCAdt[t] == -m.kp1[t] * m.CA[t]
            elif y == 'CB':
                return (
                    m.dCBdt[t]
                    == m.kp1[t] * m.CA[t] - m.kp2[t] * m.CB[t]
                )
            elif y == 'CC':
                return m.dCCdt[t] == m.kp2[t] * m.CB[t]

        @m.Constraint(m.t)
        def alge(m, t):
            """
            The algebraic equation for mole balance

            Arguments:
            z: m.pert
            t: time

            Return:
            m: Pyomo model
            """
            return m.CA[t] + m.CB[t] + m.CC[t] == m.CA0[0]

    def finalize_model(self):
        m = self.model

        # fix arrhenius parameters (any changes necessary will be done automatically by pyomo- this is necessary to ensure that the model DoF is 0)
        m.A1.fix(self.theta_pe['A1'])
        m.A2.fix(self.theta_pe['A2'])
        m.E1.fix(self.theta_pe['E1'])
        m.E2.fix(self.theta_pe['E2'])

        # fix the starting concentration value (similar to arrhenius parameters, though not strictly necessary in the same way)
        m.CA0.fix(self.CA_init)

        # Boundary Conditions
        m.CB[0.0].fix(0.0)
        m.CC[0.0].fix(0.0)

        @m.Constraint(m.t)
        def T_control(m, t):
            """
            T at interval timepoint equal to the T of the control time point at the beginning of this interval
            Count how many control points are before the current t;
            locate the nearest neighbouring control point before this t

            Arguments:
            m: model
            t: time

            Return:
            m: Pyomo model
            """
            if t in m.t_con:
                return pyo.Constraint.Skip
            else:
                j = -1
                for t_con in m.t_con:
                    if t > t_con:
                        j += 1
                neighbour_t = m.t_control[j]
                return m.T[t] == m.T[neighbour_t]


        # Discretization
        discretizer = pyo.TransformationFactory('dae.collocation')
        discretizer.apply_to(m, nfe=self.NFE, ncp=self.ncp, wrt=m.t)
        for t in m.t:
            m.dCdt_rule["CC", t].deactivate()

    def label_experiment(self):
        m = self.model

        # Concentration measurement labels
        m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
        m.experiment_outputs.update((m.CA[t], None) for t in m.t_control)
        m.experiment_outputs.update((m.CB[t], None) for t in m.t_control)
        m.experiment_outputs.update((m.CC[t], None) for t in m.t_control)

        # measurement values
        m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
        m.measurement_error.update((m.CA[t], 1e-4) for t in m.t_control)
        m.measurement_error.update((m.CB[t], 1e-4) for t in m.t_control)
        m.measurement_error.update((m.CC[t], 1e-4) for t in m.t_control)

        # Design variables
        m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
        m.experiment_inputs.update((m.CA0[t], None) for t in m.t0)
        m.experiment_inputs.update((m.T[t], None) for t in m.t_control)

        # unkown parameter labels
        m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
        m.unknown_parameters.update((k,pyo.value(k)) for k in [m.A1,m.A2,m.E1,m.E2])

    def get_labeled_model(self):
        if self.model is None:
            self.create_model()
            self.finalize_model()
            self.label_experiment()
        return self.model

5.3.5. Step 3: Generate an Experiment#

Recall the unknown model parameters, $\theta = (A_1, A_2, E_1, E_2)$.

Our goal is to maximize the precision of $\theta$ by measuring the model outputs called $y$ for the initial conditions at each time point such that $y(t)=(C_A(t), C_B(t), C_C(t))$.

The code block below generates an experiment and then runs the DOE module on the generated experiment class. The given parameters are the default parameters for DesignOfExperiments, the step size is a generally acceptable 0.001, and the finite difference formula used is the most robust first order basic method (central difference).

experiment = Reactor_Experiment()
fd_formula = "central"
step_size = 10**-3
objective_option = "determinant"
scale_nominal_param_value = True
doe_obj = DesignOfExperiments(
        experiment,
        fd_formula=fd_formula,
        step=step_size,
        objective_option=objective_option,
        scale_constant_value=1,
        scale_nominal_param_value=scale_nominal_param_value,
        prior_FIM=None,
        jac_initial=None,
        fim_initial=None,
        L_diagonal_lower_bound=1e-7,
        solver=None,
        tee=False,
        get_labeled_model_args=None,
        _Cholesky_option=True,
        _only_compute_fim_lower=True,)

5.3.6. Step 4: Run the DOE Module to get Optimal Design Parameters#

The ultimate goal of an experiment of this type is to determine the optimum design parameters. Later, we will see what it looks like to run an ensemble of experiments (which can give more information), but for now, one will suffice.

doe_obj.run_doe()

for i,j in enumerate(doe_obj.results['Experiment Design']):
    if i == 0:
        print("Ideal CA initial concentration: {} M".format(j))
    else:
        print("Ideal Temperature at checkpoint {}: {} K".format(i-1,round(j,2)))

Ideal CA initial concentration: 5.0 M
Ideal Temperature at checkpoint 0: 455.51 K
Ideal Temperature at checkpoint 1: 300.0 K
Ideal Temperature at checkpoint 2: 300.0 K
Ideal Temperature at checkpoint 3: 300.0 K
Ideal Temperature at checkpoint 4: 300.0 K
Ideal Temperature at checkpoint 5: 300.0 K
Ideal Temperature at checkpoint 6: 300.0 K
Ideal Temperature at checkpoint 7: 511.39 K
Ideal Temperature at checkpoint 8: 700.0 K

From here, the optimum parameters are to start the experiment at an A concentration of 5.0 M and an initial temperature of about 500K, with all subsequent temperatures being dropped to 300K.

5.3.6.1. Understanding the Output for Running an Experiment#

Let’s take some time to understand exactly what just happened there from a theory standpoint. When we generate experiments, we will get a class containing for the FIM, determinant, and eigenvalues.

See the following notebooks from CBE 60499 and CBE 60258 for more information on determinants and eigenvalues:

Determinant: CBE 60499 / CBE 60258
Eigenvalue: CBE 60499 / CBE 60258

5.3.6.2. Fisher Information Matrix (FIM):#

Objective:

The FIM measures the information content for the unknown parameters $\theta$ from the model output $y_i$ such that:

\[y_i=f(\psi_i, \theta)\]

In order to quantify the uncertainty of the estimated parameters for parameter estimation, consider the covariance matrix for the parameters:

\[V(\hat{\theta},\psi) = \left[\sum_{r}^{N_{r}}\sum_{r}^{N_{r}} \tilde{\sigma}_{(r,r')}Q_{r}^{T}Q_{r'}+V_{\theta}(\hat{\theta})^{-1}\right]^{-1}\]

where:

$\hat{\theta}$: estimated parameters.
$\tilde{\sigma}$: element in the inverse of the observational covariance matrix.
$r,$ $r'$: measurements.
$Q$: dynamic sensitivity.
$V_{\theta}$: prior information.

The inverse of $V$ estimates the FIM such that:

\[V(\hat{\theta},\psi) = [M(\hat{\theta},\psi)]^{-1}\]

where:

$\psi$: design vector from a DAE system

For sequential design of experiments, consider prior information such that after $N_e$ experiments, the FIM is calculated by: $$M= \sum_{k=1}^{N_e-1}M_k+M_{N_e}(\hat{\theta},\psi_{N_e}) = K+M_{N_e}(\hat{\theta},\psi_{N_e})$$

where:

$N_e - 1$: previous experiments
$K$: constant matrix encoding information from all $N_e - 1$

Key Takeaways: \ In regards to parameter estimation,

A large FIM value denotes more information about $\theta$ is gained from the model.

See the following notebook from CBE 60258 for more information on FIM:

Computing FIM

5.3.7. Step 5: Method for Computing FIM#

This method computes a FIM-based MBDoE optimization problem with no degrees of freedom.

5.3.7.1. Optimality Conditions#

In our results, we will use the four optimality conditions below:

Optimality Condition	Definition	Computation	Geometry
D-Optimality	Maximizes the determinant of $M$ or minimizes the determinant of $V$	Determinant	Minimizes the volume of the confidence ellipsoid
A-Optimality	Maximizes the trace of $M$ or minimizes the trace of $V$	Trace	Minimizes the dimensions of the enclosing box around the confidence ellipsoid
E-Optimality	Minimizes the variance of the most uncertain parameter	Eigenvalue	Minimizes the size of the major axis of the confidence ellipsoid
Modified E-Optimality	Reduces the correlations between parameters	Condition Number	Transforms the confidence ellipsoid into a round sphere

See the following notebooks from CBE 60499 and CBE 60258 for more information on condition numbers:

Condition Number: CBE 60499 / CBE 60258

For additional information on condition number and trace:

Wikipedia: Condition Number / Trace

Now that our analysis is complete, we can take a further look at our results. First, the initial FIM matrix, with only the diagonal populated.

doe_obj.fim_initial

array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

Next, let’s take a look at a single FIM matrix. After the experiment is run, the FIM is populated with values showing that the experiment has yielded information about the 4 arrhenius parameters.

doe_obj.compute_FIM()

array([[  186908.66237351,   158693.0963409 ,  -388214.6707347 ,
         -749813.95132776],
       [  158693.0963409 ,   344092.54462273,  -294525.36474362,
        -1598447.48509604],
       [ -388214.6707347 ,  -294525.36474362,   821096.94785898,
         1399756.9225552 ],
       [ -749813.95132776, -1598447.48509604,  1399756.9225552 ,
         7502382.99492581]])

5.3.8. Step 6: Running a Full Factorial Design Experiment (Exploratory Analysis)#

Next, we will run a full factorial experiment. A full factorial design tests all possible combinations of the design variable inputs to give a full picture of behavior as both are left to vary.

It is possible to run a partial factorial of any fraction, though these designs are not always advisable since they give less information (potentially leading to persisting confounding aspects of the design). However, the full factorial is by far the most computationally expensive evaluation, so it is sometimes not possible to complete if individual experiments are resource-intensive. Here, we can run a full factorial in an acceptable time, so it will be completed.

The first step in a full factorial is to designate design ranges and the number of points for each of the design variables.

# Make design ranges to compute the full factorial design
design_ranges = {"CA0[0]": [0.5,5,10], "T[0]": [300, 700, 9]}

Now, to evaluate the FIM full factorial. Using the design ranges, the object can be modified to contain the full factorial.

# Compute the full factorial design with the sequential FIM calculation
doe_obj.compute_FIM_full_factorial(design_ranges=design_ranges, method="sequential")

{'CA0[0]': [np.float64(0.5),
  np.float64(0.5),
  np.float64(0.5),
  np.float64(0.5),
  np.float64(0.5),
  np.float64(0.5),
  np.float64(0.5),
  np.float64(0.5),
  np.float64(0.5),
  np.float64(1.0),
  np.float64(1.0),
  np.float64(1.0),
  np.float64(1.0),
  np.float64(1.0),
  np.float64(1.0),
  np.float64(1.0),
  np.float64(1.0),
  np.float64(1.0),
  np.float64(1.5),
  np.float64(1.5),
  np.float64(1.5),
  np.float64(1.5),
  np.float64(1.5),
  np.float64(1.5),
  np.float64(1.5),
  np.float64(1.5),
  np.float64(1.5),
  np.float64(2.0),
  np.float64(2.0),
  np.float64(2.0),
  np.float64(2.0),
  np.float64(2.0),
  np.float64(2.0),
  np.float64(2.0),
  np.float64(2.0),
  np.float64(2.0),
  np.float64(2.5),
  np.float64(2.5),
  np.float64(2.5),
  np.float64(2.5),
  np.float64(2.5),
  np.float64(2.5),
  np.float64(2.5),
  np.float64(2.5),
  np.float64(2.5),
  np.float64(3.0),
  np.float64(3.0),
  np.float64(3.0),
  np.float64(3.0),
  np.float64(3.0),
  np.float64(3.0),
  np.float64(3.0),
  np.float64(3.0),
  np.float64(3.0),
  np.float64(3.5),
  np.float64(3.5),
  np.float64(3.5),
  np.float64(3.5),
  np.float64(3.5),
  np.float64(3.5),
  np.float64(3.5),
  np.float64(3.5),
  np.float64(3.5),
  np.float64(4.0),
  np.float64(4.0),
  np.float64(4.0),
  np.float64(4.0),
  np.float64(4.0),
  np.float64(4.0),
  np.float64(4.0),
  np.float64(4.0),
  np.float64(4.0),
  np.float64(4.5),
  np.float64(4.5),
  np.float64(4.5),
  np.float64(4.5),
  np.float64(4.5),
  np.float64(4.5),
  np.float64(4.5),
  np.float64(4.5),
  np.float64(4.5),
  np.float64(5.0),
  np.float64(5.0),
  np.float64(5.0),
  np.float64(5.0),
  np.float64(5.0),
  np.float64(5.0),
  np.float64(5.0),
  np.float64(5.0),
  np.float64(5.0)],
 'T[0]': [np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0),
  np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0),
  np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0),
  np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0),
  np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0),
  np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0),
  np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0),
  np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0),
  np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0),
  np.float64(300.0),
  np.float64(350.0),
  np.float64(400.0),
  np.float64(450.0),
  np.float64(500.0),
  np.float64(550.0),
  np.float64(600.0),
  np.float64(650.0),
  np.float64(700.0)],
 'T[0.125]': [300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300],
 'T[0.25]': [300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300],
 'T[0.375]': [300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300],
 'T[0.5]': [300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300],
 'T[0.625]': [300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300],
 'T[0.75]': [300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300],
 'T[0.875]': [300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300],
 'T[1]': [300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300,
  300],
 'log10 D-opt': [np.float64(2.2255801249152545),
  np.float64(9.549738128233253),
  np.float64(10.879904123755436),
  np.float64(11.384829138648877),
  np.float64(11.429057071876878),
  np.float64(11.174065002637425),
  np.float64(10.717496346471625),
  np.float64(10.126533183901312),
  np.float64(9.449997295889082),
  np.float64(7.1784084049900105),
  np.float64(11.95797809355288),
  np.float64(13.288144089066181),
  np.float64(13.79306910384321),
  np.float64(13.83729703672706),
  np.float64(13.582304966703033),
  np.float64(13.125736309064692),
  np.float64(12.534773144224005),
  np.float64(11.858237252822262),
  np.float64(8.58713895994145),
  np.float64(13.366708166000416),
  np.float64(14.696874161509157),
  np.float64(15.201799176290809),
  np.float64(15.24602710917269),
  np.float64(14.99103503912183),
  np.float64(14.534466381348881),
  np.float64(13.943503216178811),
  np.float64(13.266967324195152),
  np.float64(9.586643868883643),
  np.float64(14.366218058853294),
  np.float64(15.69638405437815),
  np.float64(16.201309069154576),
  np.float64(16.245537002039615),
  np.float64(15.990544931989215),
  np.float64(15.533976274212334),
  np.float64(14.943013108979834),
  np.float64(14.266477216851964),
  np.float64(10.361933591953962),
  np.float64(15.141498162939131),
  np.float64(16.4716641584397),
  np.float64(16.97658917321968),
  np.float64(17.02081710610387),
  np.float64(16.765825036053226),
  np.float64(16.309256378276956),
  np.float64(15.718293213040905),
  np.float64(15.041757320871744),
  np.float64(10.995382262297488),
  np.float64(15.774948131318448),
  np.float64(17.105114126820013),
  np.float64(17.610039141601302),
  np.float64(17.654267074484174),
  np.float64(17.399275004434102),
  np.float64(16.94270634665801),
  np.float64(16.351743181423362),
  np.float64(15.675207289246952),
  np.float64(11.530957975341293),
  np.float64(16.310522448367486),
  np.float64(17.640688443863073),
  np.float64(18.14561345864236),
  np.float64(18.189841391528216),
  np.float64(17.934849321479),
  np.float64(17.47828066370255),
  np.float64(16.88731749846797),
  np.float64(16.210781606292116),
  np.float64(11.994893150027677),
  np.float64(16.77445802417301),
  np.float64(18.104624019691123),
  np.float64(18.60954903446663),
  np.float64(18.653776967351472),
  np.float64(18.39878489730097),
  np.float64(17.94221623952475),
  np.float64(17.35125307428907),
  np.float64(16.674717182113394),
  np.float64(12.404113132908455),
  np.float64(17.183678203764376),
  np.float64(18.513844199269766),
  np.float64(19.018769214044763),
  np.float64(19.062997146928776),
  np.float64(18.80800507687958),
  np.float64(18.351436419102953),
  np.float64(17.760473253868),
  np.float64(17.08393736169353),
  np.float64(12.77016816487489),
  np.float64(17.54973812825557),
  np.float64(18.879904123752354),
  np.float64(19.38482913853232),
  np.float64(19.429057071416),
  np.float64(19.17406500136531),
  np.float64(18.717496343588405),
  np.float64(18.126533178353643),
  np.float64(17.44999728617809)],
 'log10 A-opt': [np.float64(4.960694446283322),
  np.float64(4.989231476927692),
  np.float64(5.007765325324744),
  np.float64(4.9995778942580476),
  np.float64(4.947209778444284),
  np.float64(4.843968417284895),
  np.float64(4.6939113309168485),
  np.float64(4.507021882902733),
  np.float64(4.29519258624265),
  np.float64(5.562785680625736),
  np.float64(5.59129146825561),
  np.float64(5.609825316652699),
  np.float64(5.601637885585729),
  np.float64(5.54926976977038),
  np.float64(5.4460284086056125),
  np.float64(5.295971322224454),
  np.float64(5.109081874181917),
  np.float64(4.897252577465214),
  np.float64(5.914968198701419),
  np.float64(5.94347398636703),
  np.float64(5.962007834764038),
  np.float64(5.953820403697108),
  np.float64(5.901452287881727),
  np.float64(5.798210926716919),
  np.float64(5.64815384033502),
  np.float64(5.461264392289517),
  np.float64(5.249435095564961),
  np.float64(6.164845670747029),
  np.float64(6.19335145958361),
  np.float64(6.211885307980684),
  np.float64(6.203697876913671),
  np.float64(6.1513297610983075),
  np.float64(6.048088399933528),
  np.float64(5.898031313551529),
  np.float64(5.71114186550573),
  np.float64(5.499312568779676),
  np.float64(6.35866569895547),
  np.float64(6.387171485599771),
  np.float64(6.405705333996549),
  np.float64(6.39751790292977),
  np.float64(6.3451497871144715),
  np.float64(6.241908425949612),
  np.float64(6.091851339567684),
  np.float64(5.904961891521902),
  np.float64(5.6931325947954825),
  np.float64(6.517028191447347),
  np.float64(6.54553397769498),
  np.float64(6.564067826092011),
  np.float64(6.555880395025108),
  np.float64(6.503512279209689),
  np.float64(6.40027091804487),
  np.float64(6.250213831662943),
  np.float64(6.063324383617121),
  np.float64(5.851495086890726),
  np.float64(6.650921769866683),
  np.float64(6.679427556956249),
  np.float64(6.697961405353224),
  np.float64(6.689773974286344),
  np.float64(6.637405858470937),
  np.float64(6.534164497306053),
  np.float64(6.384107410924057),
  np.float64(6.197217962878299),
  np.float64(5.985388666151961),
  np.float64(6.766905663059889),
  np.float64(6.79541145091158),
  np.float64(6.813945299308622),
  np.float64(6.805757868241632),
  np.float64(6.753389752426268),
  np.float64(6.6501483912614985),
  np.float64(6.500091304879496),
  np.float64(6.313201856833693),
  np.float64(6.101372560107306),
  np.float64(6.869210709315597),
  np.float64(6.897716495806382),
  np.float64(6.9162503442033945),
  np.float64(6.908062913136452),
  np.float64(6.855694797321057),
  np.float64(6.752453436156159),
  np.float64(6.602396349774256),
  np.float64(6.415506901728434),
  np.float64(6.203677605002089),
  np.float64(6.960725688947834),
  np.float64(6.989231476927729),
  np.float64(7.007765325324571),
  np.float64(6.999577894257735),
  np.float64(6.947209778442441),
  np.float64(6.843968417277565),
  np.float64(6.6939113308956175),
  np.float64(6.507021882849856),
  np.float64(6.295192586123432)],
 'log10 E-opt': [np.float64(-4.106392072778763),
  np.float64(-0.2686138138673309),
  np.float64(0.7368536413016509),
  np.float64(1.133651821121558),
  np.float64(1.1574468599662964),
  np.float64(1.108772492698694),
  np.float64(1.021460905928118),
  np.float64(0.9111173437424053),
  np.float64(0.7866328308820192),
  np.float64(-2.198365443268814),
  np.float64(0.33344617746570865),
  np.float64(1.3389136326308357),
  np.float64(1.7357118122741713),
  np.float64(1.7595068507515264),
  np.float64(1.710832482680246),
  np.float64(1.6235208944544095),
  np.float64(1.5131773301142029),
  np.float64(1.3886928141894206),
  np.float64(-1.8461826994348454),
  np.float64(0.6856286955797929),
  np.float64(1.6910961507415478),
  np.float64(2.0878943303860633),
  np.float64(2.111689368862331),
  np.float64(2.0630150007580013),
  np.float64(1.9757034123808714),
  np.float64(1.8653598476974294),
  np.float64(1.7408753312027407),
  np.float64(-1.596303432431236),
  np.float64(0.9355061687826229),
  np.float64(1.9409736239587831),
  np.float64(2.337771803602034),
  np.float64(2.3615668420797813),
  np.float64(2.3128924739741366),
  np.float64(2.2255808855931916),
  np.float64(2.1152373208362723),
  np.float64(1.9907528041899203),
  np.float64(-1.4024870494805368),
  np.float64(1.1293261948222466),
  np.float64(2.1347936499737457),
  np.float64(2.5315918296177227),
  np.float64(2.5553868680944367),
  np.float64(2.5067124999908628),
  np.float64(2.4194009116093866),
  np.float64(2.309057346849236),
  np.float64(2.1845728301505623),
  np.float64(-1.2441244934382105),
  np.float64(1.2876886869131605),
  np.float64(2.2931561420692566),
  np.float64(2.6899543217137465),
  np.float64(2.713749360190004),
  np.float64(2.665074992085396),
  np.float64(2.577763403705297),
  np.float64(2.4674198389451316),
  np.float64(2.342935322240347),
  np.float64(-1.1102309227677472),
  np.float64(1.4215822661768922),
  np.float64(2.427049721326405),
  np.float64(2.823847900974552),
  np.float64(2.8476429394517817),
  np.float64(2.7989685713479084),
  np.float64(2.711656982965921),
  np.float64(2.601313418206244),
  np.float64(2.4768289015014195),
  np.float64(-0.9942469530278655),
  np.float64(1.5375661601181623),
  np.float64(2.5430336152873254),
  np.float64(2.939831794929547),
  np.float64(2.963626833407538),
  np.float64(2.9149524653025463),
  np.float64(2.8276408769213033),
  np.float64(2.71729731216166),
  np.float64(2.5928127954567075),
  np.float64(-0.8919418607049382),
  np.float64(1.6398712050221116),
  np.float64(2.645338660181363),
  np.float64(3.0421368398249777),
  np.float64(3.0659318783020963),
  np.float64(3.0172575101970223),
  np.float64(2.9299459218155217),
  np.float64(2.8196023570562723),
  np.float64(2.6951178403518052),
  np.float64(-0.8004252448250803),
  np.float64(1.731386186151668),
  np.float64(2.7368536413014937),
  np.float64(3.1336518209462776),
  np.float64(3.157446859423134),
  np.float64(3.1087724913187116),
  np.float64(3.021460902937315),
  np.float64(2.9111173381776103),
  np.float64(2.7866328214729483)],
 'log10 ME-opt': [np.float64(8.915359601025779),
  np.float64(5.149643468356237),
  np.float64(4.200186197797772),
  np.float64(3.8195616500810443),
  np.float64(3.7569090036921984),
  np.float64(3.7092379845969963),
  np.float64(3.6496748495985996),
  np.float64(3.574181707997031),
  np.float64(3.486625228949382),
  np.float64(7.609568995961778),
  np.float64(5.149643468349458),
  np.float64(4.200186197796994),
  np.float64(3.8195616502563188),
  np.float64(3.7569090042329525),
  np.float64(3.7092379859349296),
  np.float64(3.649674852375437),
  np.float64(3.5741817128925955),
  np.float64(3.4866252368382615),
  np.float64(7.609568768963021),
  np.float64(5.149643468352362),
  np.float64(4.200186197797348),
  np.float64(3.819561650255849),
  np.float64(3.7569090042334996),
  np.float64(3.7092379859685445),
  np.float64(3.649674852559465),
  np.float64(3.5741817134160843),
  np.float64(3.4866252379216665),
  np.float64(7.609566975088596),
  np.float64(5.1496434683628705),
  np.float64(4.200186197796811),
  np.float64(3.8195616502564995),
  np.float64(3.7569090042325133),
  np.float64(3.7092379859689477),
  np.float64(3.6496748525636638),
  np.float64(3.57418171349349),
  np.float64(3.4866252381486746),
  np.float64(7.609570621684523),
  np.float64(5.149643468337591),
  np.float64(4.200186197797803),
  np.float64(3.8195616502567202),
  np.float64(3.756909004234074),
  np.float64(3.70923798596826),
  np.float64(3.649674852563642),
  np.float64(3.574181713496656),
  np.float64(3.486625238203747),
  np.float64(7.609570557261331),
  np.float64(5.1496434683415115),
  np.float64(4.200186197797587),
  np.float64(3.819561650256158),
  np.float64(3.7569090042337),
  np.float64(3.7092379859689877),
  np.float64(3.6496748525630527),
  np.float64(3.574181713495834),
  np.float64(3.4866252382092617),
  np.float64(7.609570566046614),
  np.float64(5.149643468342439),
  np.float64(4.200186197801948),
  np.float64(3.819561650256549),
  np.float64(3.7569090042333557),
  np.float64(3.7092379859675986),
  np.float64(3.649674852563387),
  np.float64(3.5741817134959764),
  np.float64(3.486625238209388),
  np.float64(7.609570490503468),
  np.float64(5.149643468354963),
  np.float64(4.200186197796167),
  np.float64(3.8195616502567766),
  np.float64(3.7569090042328583),
  np.float64(3.7092379859684654),
  np.float64(3.6496748525634697),
  np.float64(3.574181713495981),
  np.float64(3.486625238209561),
  np.float64(7.609570441181276),
  np.float64(5.14964346834127),
  np.float64(4.200186197796839),
  np.float64(3.819561650256194),
  np.float64(3.7569090042331506),
  np.float64(3.7092379859686577),
  np.float64(3.649674852564037),
  np.float64(3.5741817134960256),
  np.float64(3.486625238209141),
  np.float64(7.60956880674363),
  np.float64(5.149643468335421),
  np.float64(4.200186197798042),
  np.float64(3.819561650256228),
  np.float64(3.756909004233304),
  np.float64(3.7092379859684486),
  np.float64(3.649674852563485),
  np.float64(3.5741817134961495),
  np.float64(3.4866252382093386)],
 'solve_time': [0.3887155420379713,
  0.31036762497387826,
  0.35544995800592005,
  0.31027691590134054,
  0.3091701250523329,
  0.31066708301659673,
  0.31494841701351106,
  0.31176258402410895,
  0.3090299169998616,
  0.3306788749760017,
  0.32615229196380824,
  0.3225219160085544,
  0.3215128749143332,
  0.32514083303976804,
  0.3152942080050707,
  0.3214532919228077,
  0.3539591250009835,
  0.3200489589944482,
  0.3378319159382954,
  0.3141996250487864,
  0.3126658749533817,
  0.3177407500334084,
  0.3137342919362709,
  0.3141678749816492,
  0.3264159579994157,
  0.32318554201629013,
  0.32300170802045614,
  0.33766987500712276,
  0.32201145798899233,
  0.3064495000289753,
  0.3335783750517294,
  0.3337524589151144,
  0.33812750002834946,
  0.3391392088960856,
  0.31692454195581377,
  0.33777674997691065,
  0.3451724579790607,
  0.3124734580051154,
  0.3133432500762865,
  0.3193974159657955,
  0.31442708300892264,
  0.3109232500428334,
  0.31132012500893325,
  0.31164483400061727,
  0.3096928750164807,
  0.33850658300798386,
  0.3115229170070961,
  0.30144129192922264,
  0.315865374985151,
  0.2999213330913335,
  0.31695312494412065,
  0.30790612497366965,
  0.30875312502030283,
  0.31051137496251613,
  0.3306241250829771,
  0.2994608749868348,
  0.30201999994460493,
  0.30542208394035697,
  0.31475062505342066,
  0.3460419160546735,
  0.3137564589269459,
  0.3105030000442639,
  0.3014305420219898,
  0.3378682079492137,
  0.30100712506100535,
  0.30140591599047184,
  0.29906833299901336,
  0.2995164589956403,
  0.2988877500174567,
  0.2980056250235066,
  0.3007318750023842,
  0.29879241599701345,
  0.3300924169598147,
  0.301105041988194,
  0.3033839170821011,
  0.29826008388772607,
  0.3058716249652207,
  0.3013390420237556,
  0.3060440829722211,
  0.30388787493575364,
  0.3008562079630792,
  0.3370975418947637,
  0.2986111659556627,
  0.30381266702897847,
  0.30961641704197973,
  0.2995606670156121,
  0.341500207898207,
  0.3235001669963822,
  0.31139237503521144,
  0.3073393749073148]}

Why do this? It is useful to get an idea of the optimum values for the design variables, but by evaluating each over a range and filling in a full rectangle of results, we get an idea of the shape of the optimization surface. We also get an idea of what directions of optimization work on each of the optimalities, and whether the overall optimal solution performs badly at any of the metrics.

5.3.9. Step 7: Evaluating the Full Factorial Heatmaps (Exploratory Analysis)#

A heatmap shows the change of the objective function or the experimental information content in the design region.

Heatmaps can be drawn by two design variables while fixing the other design variables.

5.3.9.1. Interpreting Heatmaps#

The horizontal and vertical axes represent the two design variables, while the color of each grid shows the experimental information content.

Using the FIM criteria calculated in the last code block, the optimality conditions can be drawn in heat maps. You can see each plotted below. Note that darker colors signify better performance on that particular optimality condition.

# Plot the results
doe_obj.draw_factorial_figure(
    sensitivity_design_variables=["CA0[0]", "T[0]"],
    fixed_design_variables={
        "T[0.125]": 300,
        "T[0.25]": 300,
        "T[0.375]": 300,
        "T[0.5]": 300,
        "T[0.625]": 300,
        "T[0.75]": 300,
        "T[0.875]": 300,
        "T[1]": 300,
    },
    title_text="Reactor Example",
    xlabel_text="Concentration of A (M)",
    ylabel_text="Initial Temperature (K)",
    figure_file_name="example_reactor_compute_FIM",
    log_scale=False,
)

../../_images/062c1025524a9cd128e8a5346e04a0bd04ef390a15cb76e0a6f6268639c8580a.png

../../_images/f1c8a2362404dd3d744f60bcba1a4ff2e85b9d35c765086d05d978e35e6e2731.png

../../_images/f6541c20fb405a4fd9c4415ff6818854ae796cedfce37629e593e1834ac499f0.png

../../_images/bc4fed4ed639bea21f8abb366c534cf440309f774fa968a83855a7f5c65f5c84.png

As an example, the Figure from the Reactor case - D optimality shows:

The most informative region is around $C_{A0}=5.0 \ \text{M},$ $T=500.0 \ \text{K}$. \ The least informative region is around $C_{A0}=0.5 \ \text{M},$ $T=300.0 \ \text{K}$.

This can be a lot to take in at once, but think of each square within each plot as a unique experiment. To know what the optimality values are for an experiment set at an initial A concentration of 2.5 M and a temperature of 700 K, merely requires finding that location on the plots.

The power of viewing multiple experiments at once is that it gives a more contextualized view of the design space. Viewing one value is also possible, but it tells less of a complete story out of the surrounding space’s context.

5.3.10. Step 8: Compute D-Optimal Experiment Design (Optimization)#

5.3.10.1. Computational Optimization#

# Toggle on printing
doe_obj.tee = True

# Notice we specified the "determinant" objective when creating the DesignOfExperiments object
doe_obj.run_doe()

WARNING: Implicitly replacing the Component attribute obj_cons (type=<class
'pyomo.core.base.block.ScalarBlock'>) on block unknown with a new Component
(type=<class 'pyomo.core.base.block.ScalarBlock'>). This is usually indicative
of a modelling error. To avoid this warning, use block.del_component() and
block.add_component().
WARNING: Implicitly replacing the Component attribute objective (type=<class
'pyomo.core.base.objective.ScalarObjective'>) on block unknown with a new
Component (type=<class 'pyomo.core.base.objective.ScalarObjective'>). This is
usually indicative of a modelling error. To avoid this warning, use
block.del_component() and block.add_component().
WARNING: Implicitly replacing the Component attribute dummy_obj (type=<class
'pyomo.core.base.objective.ScalarObjective'>) on block unknown with a new
Component (type=<class 'pyomo.core.base.objective.ScalarObjective'>). This is
usually indicative of a modelling error. To avoid this warning, use
block.del_component() and block.add_component().
Ipopt 3.13.2: linear_solver=ma57
halt_on_ampl_error=yes
max_iter=3000


******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma57.

Number of nonzeros in equality constraint Jacobian...:    25881
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     2581

Total number of variables............................:     7076
                     variables with only lower bounds:     3864
                variables with lower and upper bounds:      774
                     variables with only upper bounds:        0
Total number of equality constraints.................:     7076
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 4.00e+00 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
Reallocating memory for MA57: lfact (286238)
   1  0.0000000e+00 2.99e+01 8.25e-02  -1.0 1.44e+04    -  9.80e-01 9.90e-01h  1
   2  0.0000000e+00 2.96e+01 1.98e+00  -1.0 1.42e+04    -  1.00e+00 9.90e-01h  1
   3  0.0000000e+00 2.72e-03 3.27e-06  -1.0 1.35e+02    -  1.00e+00 1.00e+00h  1
   4  0.0000000e+00 9.31e-10 2.83e-08  -2.5 2.80e-03    -  1.00e+00 1.00e+00h  1
Cannot recompute multipliers for feasibility problem.  Error in eq_mult_calculator

Number of Iterations....: 4

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   5.6568541627013561e+04    5.6568541627013561e+04
Constraint violation....:   6.1395254091954413e-13    9.3132257461547852e-10
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   6.1395254091954413e-13    5.6568541627013561e+04


Number of objective function evaluations             = 5
Number of objective gradient evaluations             = 5
Number of equality constraint evaluations            = 5
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 5
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 4
Total CPU secs in IPOPT (w/o function evaluations)   =      0.062
Total CPU secs in NLP function evaluations           =      0.001

EXIT: Optimal Solution Found.
Ipopt 3.13.2: linear_solver=ma57
halt_on_ampl_error=yes
max_iter=3000


******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma57.

Number of nonzeros in equality constraint Jacobian...:    26194
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     2610

Reallocating memory for MA57: lfact (271756)
Total number of variables............................:     7096
                     variables with only lower bounds:     3868
                variables with lower and upper bounds:      784
                     variables with only upper bounds:        0
Total number of equality constraints.................:     7086
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -1.9340418e+01 1.03e+00 1.01e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
Reallocating memory for MA57: lfact (300708)
   1 -1.9340276e+01 6.76e+02 5.90e+00  -1.0 2.94e+05    -  9.76e-01 4.21e-01h  1
   2 -1.9216117e+01 1.40e+03 4.51e+00  -1.0 3.63e+04    -  9.38e-01 1.00e+00f  1
   3 -1.8686210e+01 4.32e+04 2.16e+01  -1.0 1.08e+06    -  2.18e-01 4.14e-01f  1
   4 -1.7970661e+01 4.30e+05 3.54e+01  -1.0 3.26e+05    -  1.00e+00 1.00e+00f  1
   5 -1.7495968e+01 2.30e+05 2.63e+02  -1.0 7.73e+05    -  1.00e+00 1.00e+00h  1
   6 -1.7457567e+01 1.74e+05 1.10e+02  -1.0 1.30e+06    -  1.00e+00 1.00e+00h  1
   7 -1.7383351e+01 1.09e+04 3.16e+01  -1.0 2.83e+05    -  1.00e+00 1.00e+00h  1
   8 -1.7410531e+01 3.61e+01 2.76e-01  -1.0 1.47e+04    -  1.00e+00 1.00e+00h  1
   9 -1.7495161e+01 1.01e+02 5.04e-01  -1.7 1.91e+04    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -1.7727475e+01 1.25e+03 1.71e+00  -1.7 6.84e+04    -  1.00e+00 1.00e+00h  1
  11 -1.7731011e+01 1.10e+01 1.80e-02  -1.7 1.24e+04    -  1.00e+00 1.00e+00h  1
  12 -1.7731051e+01 6.89e-03 1.99e-05  -1.7 4.93e+01    -  1.00e+00 1.00e+00h  1
  13 -1.8012138e+01 3.09e+03 2.62e+00  -3.8 7.57e+04    -  8.49e-01 1.00e+00f  1
  14 -1.8673635e+01 3.34e+05 7.04e+01  -3.8 2.92e+05    -  6.15e-01 1.00e+00h  1
  15 -1.8862837e+01 1.02e+05 4.50e+01  -3.8 4.39e+05    -  6.17e-01 1.00e+00h  1
  16 -1.9109565e+01 1.04e+04 4.36e+00  -3.8 6.72e+04    -  1.00e+00 1.00e+00h  1
  17 -1.9287347e+01 1.09e+05 1.25e+01  -3.8 8.40e+05    -  8.47e-01 7.70e-01h  1
Reallocating memory for MA57: lfact (315877)
  18 -1.9242152e+01 6.93e+03 1.30e+00  -3.8 2.55e+05    -  1.00e+00 1.00e+00h  1
  19 -1.9250112e+01 1.47e+02 8.58e-02  -3.8 1.14e+04    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20 -1.9250104e+01 9.04e-01 7.72e-04  -3.8 6.51e+02    -  1.00e+00 1.00e+00h  1
  21 -1.9250104e+01 1.55e-05 4.00e-09  -3.8 4.41e+00    -  1.00e+00 1.00e+00h  1
  22 -1.9335879e+01 8.97e+03 6.95e-01  -5.7 5.96e+04    -  7.99e-01 1.00e+00f  1
  23 -1.9339832e+01 2.77e+03 8.26e-02  -5.7 4.70e+04    -  1.00e+00 1.00e+00h  1
  24 -1.9339328e+01 7.33e+00 4.71e-04  -5.7 1.89e+03    -  1.00e+00 1.00e+00h  1
  25 -1.9339326e+01 4.09e-05 5.03e-09  -5.7 7.86e+00    -  1.00e+00 1.00e+00h  1
Reallocating memory for MA57: lfact (333219)
  26 -1.9340417e+01 7.36e+00 3.42e-04  -8.6 1.75e+03    -  9.97e-01 9.99e-01h  1
Reallocating memory for MA57: lfact (363246)
Reallocating memory for MA57: lfact (448020)
Reallocating memory for MA57: lfact (475112)
  27 -1.9340422e+01 5.93e-04 3.90e-08  -8.6 2.25e+01    -  1.00e+00 1.00e+00h  1
  28 -1.9340422e+01 5.68e-04 5.39e-04  -8.6 2.71e+00  -4.0 4.37e-02 4.30e-02h  1
  29 -1.9340422e+01 5.28e-04 3.36e-01  -8.6 3.62e+00  -4.5 1.00e+00 7.05e-02f  2
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30 -1.9340421e+01 3.38e-04 1.55e-01  -8.6 3.62e+00  -5.0 9.29e-01 3.72e-01h  1
  31 -1.9340419e+01 2.05e-05 6.88e-02  -8.6 2.43e+00  -5.4 5.24e-01 1.00e+00f  1
  32 -1.9340419e+01 5.09e-07 4.58e-02  -8.6 4.96e-02  -5.9 7.09e-01 1.00e+00h  1
  33 -1.9340418e+01 3.93e-07 8.31e-03  -8.6 7.38e-02  -6.4 8.25e-01 1.00e+00h  1
  34 -1.9340418e+01 2.35e-07 1.02e-07  -8.6 4.52e-02  -6.9 1.00e+00 1.00e+00h  1
  35 -1.9340418e+01 1.11e-08 4.87e-09  -8.6 8.17e-03  -7.3 1.00e+00 1.00e+00h  1

Number of Iterations....: 35

                                   (scaled)                 (unscaled)
Objective...............:  -1.9340417993573684e+01   -1.9340417993573684e+01
Dual infeasibility......:   4.8741362352397229e-09    4.8741362352397229e-09
Constraint violation....:   2.2171843028218064e-09    1.1085921514109032e-08
Complementarity.........:   2.6084073511108184e-09    2.6084073511108184e-09
Overall NLP error.......:   4.8741362352397229e-09    1.1085921514109032e-08


Number of objective function evaluations             = 37
Number of objective gradient evaluations             = 36
Number of equality constraint evaluations            = 37
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 36
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 35
Total CPU secs in IPOPT (w/o function evaluations)   =      0.628
Total CPU secs in NLP function evaluations           =      0.011

EXIT: Optimal Solution Found.

5.3.10.2. Extract and Summarize FIM#

Next, we can summarize the FIM.

def results_summary(result):
    ''' Summarize the results of the experiment design
    Arguments:
    ----------
    result: The FIM matrix from the optimized experiment

    Returns:
    --------
    None
    
    Note: Taken from https://github.com/dowlinglab/pyomo-doe/blob/main/notebooks/tclab_pyomo.py

    '''
    eigenvalues, eigenvectors = np.linalg.eig(result)

    min_eig = min(eigenvalues)

    print("======Results Summary======")
    print("Four design criteria log10() value:")
    print("A-optimality:", np.log10(np.trace(result)))
    print("D-optimality:", np.log10(np.linalg.det(result)))
    print("E-optimality:", np.log10(min_eig))
    print("Modified E-optimality:", np.log10(np.linalg.cond(result)))
    print("\nFIM:\n", np.array(result))

    print("\neigenvalues:\n", eigenvalues)

    print("\neigenvectors:\n", eigenvectors)

results_summary(doe_obj.results['FIM'])

======Results Summary======
Four design criteria log10() value:
A-optimality: 6.9316222021521385
D-optimality: 19.34041799357368
E-optimality: 3.138893050718441
Modified E-optimality: 3.741473277791075

FIM:
 [[  177727.96998605    99340.89223057  -413453.83379003  -501967.36430149]
 [   99340.89223057   292556.10957161  -199008.0877333  -1435221.89445854]
 [ -413453.83379003  -199008.0877333    975791.78913822  1008806.63136677]
 [ -501967.36430149 -1435221.89445854  1008806.63136677  7097156.16924074]]

eigenvalues:
 [7.59217707e+06 9.47467067e+05 1.37687036e+03 2.21102927e+03]

eigenvectors:
 [[-0.07674202  0.37300899 -0.91723726 -0.11683652]
 [-0.19507385 -0.04024196 -0.12386907  0.97210248]
 [ 0.15779223 -0.91113555 -0.37684189 -0.05407228]
 [ 0.96496553  0.17051948 -0.03636558  0.19606679]]

5.3.10.3. Extract and Visualize Optimal Dynamic Experiment#

Finally, let’s extract and plot the results of the dynamic experiment.

def plot_experiment_results(model):
    """
    Plot the results of the experiment design

    Arguments:
    ----------
    model: The Pyomo model

    Returns:
    --------
    None
    """
    # Extract the timepoints
    time = np.array([pyo.value(t) for t in model.t])

    # Extract the concentrations
    CA = np.array([pyo.value(model.CA[t]) for t in model.t])
    CB = np.array([pyo.value(model.CB[t]) for t in model.t])
    CC = np.array([pyo.value(model.CC[t]) for t in model.t])
    T = np.array([pyo.value(model.T[t]) for t in model.t])

        # Create the plot
    fig, ax1 = plt.subplots()

    # Plot concentrations on the left y-axis
    ax1.plot(time, CA, label="$C_A$", color="blue")
    ax1.plot(time, CB, label="$C_B$", color="green")
    ax1.plot(time, CC, label="$C_C$", color="purple")
    ax1.set_xlabel("Time (h)")
    ax1.set_ylabel("Concentration (mol/L)")
    ax1.legend(loc="upper left")

    # Create a second y-axis for temperature
    ax2 = ax1.twinx()
    ax2.plot(time, T, label="$T$", color="red", linestyle="--")
    ax2.set_ylabel("Temperature (K)", color="red")
    ax2.tick_params(axis='y', labelcolor="red")

    # Add legend for the temperature line
    ax2.legend(loc="upper right")

    plt.title("Experiment Results")
    plt.show()

# Plot the results of the experiment design
plot_experiment_results(doe_obj.model.scenario_blocks[0])

../../_images/1e7776db187c9f795d1cf90d1f47effc6304436f067e2d72a54b0814304b090b.png

5.3.11. Key Takeaways#

DOE is helpful for guiding decision-making by maximizing information yield in experimental design.
FIM allows us to gain information about the data from a mathematical model in an experiment.
Optimality conditions tells us the most and least informative regions in regards to the measurements in experimental design.
Heatmaps enable us to visualize the most informative parameters using the optimality conditions.

Optimality Condition	Definition	Computation	Geometry
D-Optimality	Maximizes the determinant of \(M\) or minimizes the determinant of \(V\)	Determinant	Minimizes the volume of the confidence ellipsoid
A-Optimality	Maximizes the trace of \(M\) or minimizes the trace of \(V\)	Trace	Minimizes the dimensions of the enclosing box around the confidence ellipsoid
E-Optimality	Minimizes the variance of the most uncertain parameter	Eigenvalue	Minimizes the size of the major axis of the confidence ellipsoid
Modified E-Optimality	Reduces the correlations between parameters	Condition Number	Transforms the confidence ellipsoid into a round sphere

Optimizing Experiments with Pyomo.DoE

Contents

5.3. Optimizing Experiments with Pyomo.DoE#

5.3.1. Introduction and Learning Objectives#

5.3.2. Step 0: Import Pyomo and Pyomo.DoE Module#

5.3.3. Step 1: Mathematical Model for the Reaction Kinetics Example#

5.3.3.1. Reaction Kinetics#

5.3.3.2. Batch Reactor#

5.3.4. Step 2: Implement Mathematical Model#

5.3.5. Step 3: Generate an Experiment#

5.3.6. Step 4: Run the DOE Module to get Optimal Design Parameters#

5.3.6.1. Understanding the Output for Running an Experiment#

5.3.6.2. Fisher Information Matrix (FIM):#

5.3.7. Step 5: Method for Computing FIM#

5.3.7.1. Optimality Conditions#

5.3.8. Step 6: Running a Full Factorial Design Experiment (Exploratory Analysis)#

5.3.9. Step 7: Evaluating the Full Factorial Heatmaps (Exploratory Analysis)#

5.3.9.1. Interpreting Heatmaps#

5.3.10. Step 8: Compute D-Optimal Experiment Design (Optimization)#

5.3.10.1. Computational Optimization#

5.3.10.2. Extract and Summarize FIM#

5.3.10.3. Extract and Visualize Optimal Dynamic Experiment#

5.3.11. Key Takeaways#

5.3. Optimizing Experiments with `Pyomo.DoE`#