Computing Sensitivity Matrix in Pyomo.DoE

Recall that Pyomo.DoE utilizes the sensitivity matrix defined as the derivative of the estimates of the measured outputs $\hat{y}$ with respect to the parameter to be estimated $\hat{\theta}$ :

Q = \nabla_{\hat{\theta}} \hat{y}

(1)

to compute the FIM as:

M = Q^T \Sigma_y^{-1}Q

(2)

where $\Sigma_y$ is the measurement covariance matrix

Pyomo.DoE offers two different approaches to computing the sensitivity matrix:

Finite difference method: With Central finite difference method that utilizes parameter perturbations for estimating the sensitivity matrix as:
$Q \approx \frac{\hat{y}(\hat{\theta} + \Delta \hat{\theta})- \hat{y}(\hat{\theta} - \Delta \hat{\theta})}{2 \Delta \hat{\theta}}$
(3)
where $\Delta \hat{\theta}$ is the step size. This method requires $2n_p+1$ copies of the base model to solve. Therefore, its computational complexity increases with the number of parameters
Automatic differentiation method: Pyomo.DoE uses automatic differentiation in PyNumero to compute sensitivities directly from the model equations, avoiding repeated parameter perturbation solves:
$\nabla_{\hat{\theta}}F = \begin{bmatrix}\nabla_{\hat{y}}F\\\nabla_{\hat{x}}F\end{bmatrix}\begin{bmatrix}\nabla_{\hat{\theta}}\hat{y}\\\nabla_{\hat{\theta}}\hat{x}\end{bmatrix}$
(4)
where $F(\hat{x},\hat{y}, \hat{\theta}, u) = 0$ is the model of the system. This method does not require multiple copies of the original model, and therefore, is expected to reduce the total computation time.

Now that we have introduced the two sensitivity calculation methods, we compare their computational cost. Both approaches produce the sensitivity matrix needed to assemble the FIM, but they differ in how many model evaluations they require.

The finite-difference approach estimates sensitivities by repeatedly perturbing parameters and resolving the model. This can become expensive as the number of parameters grows. In contrast, the automatic differentiation approach computes sensitivities directly from the model equations using PyNumero, avoiding repeated finite-difference solves.

Next, we time both methods on the TC Lab example and compare the resulting FIM calculations.

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"

# 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

--2026-06-15 20:42:58--  https://raw.githubusercontent.com/dowlinglab/pyomo-doe/main/notebooks/tclab_pyomo.py
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 46074 (45K) [text/plain]
Saving to: ‘tclab_pyomo.py.4’

tclab_pyomo.py.4    100%[===================>]  44.99K  --.-KB/s    in 0.002s  

2026-06-15 20:42:58 (26.5 MB/s) - ‘tclab_pyomo.py.4’ saved [46074/46074]

idaes was found! No need to install.
Correct git branch of Pyomo.DoE is installed.
Finished installing software

Load Experimental Data (sine test)¶

We will load the sine test data to serve as an initial point. Recall our create model function will use supplied data to initialize the Pyomo model. Carefully initialization is often required for optimization of large-scale dynamic systems.

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)

We will create the data object for the design of experiment problems that follow

# Here, we will induce a step size of 15 seconds, as to not give too many
# degrees of freedom for experimental design.
skip = 8

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

Calculate FIM at Initial Point (Values from $L_2$ Regularization)¶

We will start by re-computing the prior FIM from $L_2$ regularization

# Load Pyomo.DoE class
from pyomo.contrib.doe import DesignOfExperiments, rescale_FIM
import numpy as np
from pyomo.environ import SolverFactory

# Copied from previous notebook
theta_values = {
    "Ua": 0.041705,
    "Ub": 0.012009,
    "inv_CpH": 0.167457,
    "inv_CpS": 4.545432,
}

theta_ref = np.array(
    [
        theta_values["Ua"],
        theta_values["Ub"],
        theta_values["inv_CpH"],
        theta_values["inv_CpS"],
    ]
)

# Defining the prior

import numpy as np

cov = np.array(
    [
        [1.857017e-10, -2.576198e-10, 1.402148e-09, -2.242347e-12],
        [-2.576198e-10, 1.624383e-07, 9.109870e-08, -6.325555e-05],
        [1.402148e-09, 9.109870e-08, 1.031454e-07, -3.890789e-05],
        [-2.242347e-12, -6.325555e-05, -3.890789e-05, 2.499914e-02],
    ]
)
FIM = np.linalg.inv(cov)


SCALE_BY_NOMINAL_PARAM = True
SCALE_CONSTANT_VALUE = 1e-2

if SCALE_BY_NOMINAL_PARAM:
    PRIOR_FIM_SCALED = rescale_FIM(FIM, theta_ref) * SCALE_CONSTANT_VALUE**2
else:
    PRIOR_FIM_SCALED = FIM * SCALE_CONSTANT_VALUE**2

# Create experiment object for design of experiments
doe_experiment = TC_Lab_experiment(
    data=tc_data, theta_initial=theta_values, number_of_states=number_tclab_states
)

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"
solver.options["acceptable_tol"] = 1e-3
grey_box_solver = SolverFactory("cyipopt")
grey_box_solver.options["max_iter"] = 3000
grey_box_solver.options["tol"] = 1e-5
grey_box_solver.options["linear_solver"] = "ma57"
grey_box_solver.options["nlp_scaling_method"] = "gradient-based"
grey_box_solver.options["acceptable_tol"] = 1e-3
# Create the design of experiments object using our experiment instance from above
TC_Lab_DoE = DesignOfExperiments(
    experiment=doe_experiment,
    step=1e-2,
    scale_constant_value=SCALE_CONSTANT_VALUE,
    scale_nominal_param_value=SCALE_BY_NOMINAL_PARAM,
    prior_FIM=PRIOR_FIM_SCALED,
    tee=True,
    solver=solver,
)

# Get the FIM and the jacobian at the initial solve
try:
    init_fim = TC_Lab_DoE.compute_FIM()
    init_jac = TC_Lab_DoE.seq_jac
except:
    init_fim = None
    init_jac = None

Fetching long content....

Optimize Next Experiment (E-Optimality) using Central Finite Difference Method¶

We are now ready to solve the next experiment design problem. Here, we create a new DesignOfExperiments object and pass in the prior_FIM, which represents the information already collected from the previous experiment. By optimizing with this prior information included, Pyomo.DoE identifies the next best experiment to run according to the selected design criterion. In this section, we use E-optimality by setting objective_option="minimum_eigenvalue" and grey box using use_greybox_objective = True.

# Create experiment object for design of experiments
doe_experiment = TC_Lab_experiment(
    data=tc_data, theta_initial=theta_values, number_of_states=number_tclab_states
)

# Create the design of experiments object using our experiment instance from above
TC_Lab_DoE_E = DesignOfExperiments(
    experiment=doe_experiment,
    step=1e-2,
    scale_constant_value=SCALE_CONSTANT_VALUE,
    scale_nominal_param_value=SCALE_BY_NOMINAL_PARAM,
    prior_FIM=PRIOR_FIM_SCALED,
    objective_option="minimum_eigenvalue",
    tee=True,
    grey_box_tee = True,
    solver=solver,
    grey_box_solver=grey_box_solver,
    jac_initial=init_jac,
    fim_initial=init_fim,
    # Use the GreyBox objective implementation
    use_grey_box_objective=True,
)

TC_Lab_DoE_E.run_doe()

Fetching long content....

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/tmp/ipykernel_5012/3459065894.py in <cell line: 0>()
     22 )
     23 
---> 24 TC_Lab_DoE_E.run_doe()

/usr/local/lib/python3.12/dist-packages/pyomo/contrib/doe/doe.py in run_doe(self, model, results_file)
    584         # Solve the full model, which has now been initialized with the square solve
    585         if self.use_grey_box:
--> 586             res = self.grey_box_solver.solve(model, tee=self.grey_box_tee)
    587         else:
    588             res = self.solver.solve(model, tee=self.tee)

/usr/local/lib/python3.12/dist-packages/pyomo/opt/base/solvers.py in solve(self, *args, **kwds)
    625             initial_time = time.time()
    626 
--> 627             self._presolve(*args, **kwds)
    628 
    629             presolve_completion_time = time.time()

/usr/local/lib/python3.12/dist-packages/pyomo/opt/solver/shellcmd.py in _presolve(self, *args, **kwds)
    219         self._define_signal_handlers = kwds.pop('use_signal_handling', None)
    220 
--> 221         OptSolver._presolve(self, *args, **kwds)
    222 
    223         #

/usr/local/lib/python3.12/dist-packages/pyomo/opt/base/solvers.py in _presolve(self, *args, **kwds)
    731             write_start_time = time.time()
    732             self._problem_files, self._problem_format, self._smap_id = (
--> 733                 self._convert_problem(
    734                     args, self._problem_format, self._valid_problem_formats, **kwds
    735                 )

/usr/local/lib/python3.12/dist-packages/pyomo/opt/base/solvers.py in _convert_problem(self, args, problem_format, valid_problem_formats, **kwds)
    783 
    784     def _convert_problem(self, args, problem_format, valid_problem_formats, **kwds):
--> 785         return convert_problem(
    786             args, problem_format, valid_problem_formats, self.has_capability, **kwds
    787         )

/usr/local/lib/python3.12/dist-packages/pyomo/opt/base/convert.py in convert_problem(args, target_problem_type, valid_problem_types, has_capability, **kwds)
     92                     tmpkw = kwds
     93                     tmpkw['capabilities'] = has_capability
---> 94                     problem_files, symbol_map = converter.apply(*tmp, **tmpkw)
     95                     return problem_files, ptype, symbol_map
     96 

/usr/local/lib/python3.12/dist-packages/pyomo/solvers/plugins/converter/model.py in apply(self, *args, **kwds)
    180                     )
    181                 else:
--> 182                     problem_filename, symbol_map_id = instance.write(
    183                         filename=problem_filename,
    184                         format=args[1],

/usr/local/lib/python3.12/dist-packages/pyomo/core/base/block.py in write(self, filename, format, solver_capability, io_options, int_marker)
   1993                 return True
   1994 
-> 1995         filename, smap = problem_writer(self, filename, solver_capability, io_options)
   1996         smap_id = id(smap)
   1997         if not hasattr(self, 'solutions'):

/usr/local/lib/python3.12/dist-packages/pyomo/repn/plugins/nl_writer.py in __call__(self, model, filename, solver_capability, io_options)
    313             _open(col_fname) as COLFILE,
    314         ):
--> 315             info = self.write(model, FILE, ROWFILE, COLFILE, config=config)
    316         if not info.variables:
    317             # This exception is included for compatibility with the

/usr/local/lib/python3.12/dist-packages/pyomo/repn/plugins/nl_writer.py in write(self, model, ostream, rowstream, colstream, **options)
    370         # small objects.
    371         with _NLWriter_impl(ostream, rowstream, colstream, config) as impl:
--> 372             return impl.write(model)
    373 
    374     def _generate_symbol_map(self, info):

/usr/local/lib/python3.12/dist-packages/pyomo/repn/plugins/nl_writer.py in write(self, model)
    561         )
    562         if unknown:
--> 563             raise ValueError(
    564                 "The model ('%s') contains the following active components "
    565                 "that the NL writer does not know how to process:\n\t%s"

ValueError: The model ('unknown') contains the following active components that the NL writer does not know how to process:
	<class 'pyomo.contrib.pynumero.interfaces.external_grey_box.ExternalGreyBoxBlock'>:
		obj_cons.egb_fim_block

Run E-Optimality Experiment using Automatic Differentiation¶

We now use the automatic differentiation method to compute sensitivities by setitng gradient_method = 'pynumero

#
# Create experiment object for design of experiments
doe_experiment = TC_Lab_experiment(
    data=tc_data, theta_initial=theta_values, number_of_states=number_tclab_states
)

# Create the design of experiments object using our experiment instance from above
TC_Lab_DoE_E_AD = DesignOfExperiments(
    experiment=doe_experiment,
    step=1e-2,
    scale_constant_value=SCALE_CONSTANT_VALUE,
    scale_nominal_param_value=SCALE_BY_NOMINAL_PARAM,
    gradient_method="pynumero",
    objective_option="minimum_eigenvalue",
    tee=True,
    grey_box_tee = True,
    solver=solver,
    jac_initial=init_jac,
    fim_initial=init_fim,
    prior_FIM=PRIOR_FIM_SCALED,
    # Use the GreyBox objective implementation
    use_grey_box_objective=True,
    grey_box_solver=grey_box_solver,
)

TC_Lab_DoE_E_AD.run_doe()

Comparison of the Central Finite Difference and Automatic Differentiation Methods’ Computational Time¶

import matplotlib.pyplot as plt
import numpy as np

central_time_metrics = {
    "Build Time": TC_Lab_DoE_E.results.get("Build Time"),
    "Initialization Time": TC_Lab_DoE_E.results.get("Initialization Time"),
    "Solve Time": TC_Lab_DoE_E.results.get("Solve Time"),
    "Wall-clock Time": TC_Lab_DoE_E.results.get("Wall-clock Time"),
}

AD_time_metrics = {
    "Build Time": TC_Lab_DoE_E_AD.results.get("Build Time"),
    "Initialization Time": TC_Lab_DoE_E_AD.results.get("Initialization Time"),
    "Solve Time": TC_Lab_DoE_E_AD.results.get("Solve Time"),
    "Wall-clock Time": TC_Lab_DoE_E_AD.results.get("Wall-clock Time"),
}

labels = []
central_vals = []
AD_vals = []

for key in central_time_metrics:
    c = central_time_metrics.get(key)
    s = AD_time_metrics.get(key)
    if c is None or s is None:
        continue
    try:
        c = float(c)
        s = float(s)
    except (TypeError, ValueError):
        continue
    labels.append(key)
    central_vals.append(c)
    AD_vals.append(s)

x = np.arange(len(labels))
width = 0.36

plt.figure(figsize=(10, 5))
plt.bar(
    x - width / 2,
    central_vals,
    width,
    label="Central finite difference",
    color="#b7cde8",
)
plt.bar(
    x + width / 2, AD_vals, width, label="Automatic Differentiation", color="#4f86c6"
)

plt.xticks(x, labels, rotation=30, ha="right")
plt.ylabel("Time (s)")
plt.title("Time comparison: Central finite difference vs Automatic differentiation")
plt.legend()
plt.tight_layout()
plt.show()

The timing comparison separates the workflow into model build time, initialization time, solve time, and total wall-clock time. In this run, the automatic differentiation formulation takes longer to build, but it solves much faster than the central finite difference formulation.

This tradeoff is expected. Automatic differentiation requires additional setup to construct derivative information directly from the model equations, which increases the build time. However, once that derivative information is available, the optimizer can use it efficiently during the solve.

In contrast, the central finite difference formulation estimates sensitivities by applying parameter perturbations and resolving perturbed model instances. These additional perturbation-based calculations make the solve phase more expensive. For this example, the longer finite-difference solve time dominates the total wall-clock time, so the automatic differentiation workflow is faster overall.

Load Experimental Data (sine test)¶

Calculate FIM at Initial Point (Values from L2L_2L2​ Regularization)¶

Optimize Next Experiment (E-Optimality) using Central Finite Difference Method¶

Run E-Optimality Experiment using Automatic Differentiation¶

Comparison of the Central Finite Difference and Automatic Differentiation Methods’ Computational Time¶

Calculate FIM at Initial Point (Values from $L_2$ Regularization)¶