Add grey box objectives to Pyomo.DoE

pyomo/contrib/doe/__init__.py

@@ -10,6 +10,7 @@
 #  ___________________________________________________________________________
 from .doe import DesignOfExperiments, ObjectiveLib, FiniteDifferenceStep
 from .utils import rescale_FIM
+from .grey_box_utilities import FIMExternalGreyBox
 
 # Deprecation errors for old Pyomo.DoE interface classes and structures
 from pyomo.common.deprecation import deprecated

pyomo/contrib/doe/examples/reactor_experiment.py

@@ -156,6 +156,7 @@ def finalize_model(self):
         for t in m.t:
             if t in control_points:
                 cv = control_points[t]
+                m.T[t].fix()
             m.T[t].setlb(self.data["T_bounds"][0])
             m.T[t].setub(self.data["T_bounds"][1])
             m.T[t] = cv

pyomo/contrib/doe/examples/result.json

@@ -1 +1 @@
-{"CA0": 5.0, "CA_bounds": [1.0, 5.0], "CB0": 0.0, "CC0": 0.0, "t_range": [0, 1], "control_points": {"0": 500, "0.125": 300, "0.25": 300, "0.375": 300, "0.5": 300, "0.625": 300, "0.75": 300, "0.875": 300, "1": 300}, "T_bounds": [300, 700], "A1": 84.79, "A2": 371.72, "E1": 7.78, "E2": 15.05}
+{"CA0": 5.0, "CA_bounds": [1.0, 5.0], "CB0": 0.0, "CC0": 0.0, "t_range": [0, 1], "control_points": {"0": 500, "0.125": 300, "0.25": 300, "0.375": 300, "0.5": 300, "0.625": 300, "0.75": 300, "0.875": 300, "1": 300}, "T_bounds": [300, 700], "A1": 84.79, "A2": 371.72, "E1": 7.78, "E2": 15.05}

pyomo/contrib/doe/examples/rooney_biegler_example.py

@@ -0,0 +1,148 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2025
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+"""
+Rooney Biegler model, based on Rooney, W. C. and Biegler, L. T. (2001). Design for
+model parameter uncertainty using nonlinear confidence regions. AIChE Journal,
+47(8), 1794-1804.
+"""
+from pyomo.common.dependencies import numpy as np, pathlib
+
+from pyomo.contrib.doe.examples.rooney_biegler_experiment import (
+    RooneyBieglerExperimentDoE,
+)
+from pyomo.contrib.doe import DesignOfExperiments
+
+import pyomo.environ as pyo
+
+import json
+import sys
+
+
+# Example comparing Cholesky factorization
+# (standard solve) with grey box objective
+# solve for the log-deteriminant of the FIM
+# (D-optimality)
+def run_rooney_biegler_doe():
+    # Create a RooneyBiegler Experiment
+    experiment = RooneyBieglerExperimentDoE(data={'hour': 2, 'y': 10.3})
+
+    # Use a central difference, with step size 1e-3
+    fd_formula = "central"
+    step_size = 1e-3
+
+    # Use the determinant objective with scaled sensitivity matrix
+    objective_option = "determinant"
+    scale_nominal_param_value = True
+
+    data = [[1, 8.3], [7, 19.8], [2, 10.3], [5, 15.6], [3, 19.0], [4, 16.0]]
+    FIM_prior = np.zeros((2, 2))
+    # Calculate prior using existing experiments
+    for i in range(len(data)):
+        if i > int(sys.argv[1]):
+            break
+        prev_experiment = RooneyBieglerExperimentDoE(
+            data={'hour': data[i][0], 'y': data[i][1]}
+        )
+        doe_obj = DesignOfExperiments(
+            prev_experiment,
+            fd_formula=fd_formula,
+            step=step_size,
+            objective_option=objective_option,
+            scale_nominal_param_value=scale_nominal_param_value,
+            prior_FIM=None,
+            tee=False,
+        )
+
+        FIM_prior += doe_obj.compute_FIM(method='sequential')
+
+    if sys.argv[1] == 0:
+        FIM_prior[0][0] += 1e-6
+        FIM_prior[1][1] += 1e-6
+
+    # Create the DesignOfExperiments object
+    # We will not be passing any prior information in this example
+    # and allow the experiment object and the DesignOfExperiments
+    # call of ``run_doe`` perform model initialization.
+    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=FIM_prior,
+        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,
+    )
+
+    # Begin optimal DoE
+    ####################
+    doe_obj.run_doe()
+
+    # Print out a results summary
+    print("Optimal experiment values: ")
+    print(
+        "\tOptimal measurement time: {:.2f}".format(
+            doe_obj.results["Experiment Design"][0]
+        )
+    )
+    print("FIM at optimal design:\n {}".format(np.array(doe_obj.results["FIM"])))
+    print(
+        "Objective value at optimal design: {:.2f}".format(
+            pyo.value(doe_obj.model.objective)
+        )
+    )
+
+    print(doe_obj.results["Experiment Design Names"])
+
+    ###################
+    # End optimal DoE
+
+    # Begin optimal greybox DoE
+    ############################
+    doe_obj_gb = DesignOfExperiments(
+        experiment,
+        fd_formula=fd_formula,
+        step=step_size,
+        objective_option=objective_option,
+        use_grey_box_objective=True,
+        scale_nominal_param_value=scale_nominal_param_value,
+        prior_FIM=FIM_prior,
+        tee=False,
+    )
+
+    doe_obj_gb.run_doe()
+
+    print("Optimal experiment values: ")
+    print(
+        "\tOptimal measurement time: {:.2f}".format(
+            doe_obj_gb.results["Experiment Design"][0]
+        )
+    )
+    print("FIM at optimal design:\n {}".format(np.array(doe_obj_gb.results["FIM"])))
+    print(
+        "Objective value at optimal design: {:.2f}".format(
+            np.log10(np.exp(pyo.value(doe_obj_gb.model.objective)))
+        )
+    )
+
+    ############################
+    # End optimal greybox DoE
+
+
+if __name__ == "__main__":
+    run_rooney_biegler_doe()

pyomo/contrib/doe/examples/rooney_biegler_experiment.py

@@ -0,0 +1,100 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2025
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+
+"""
+Rooney Biegler model, based on Rooney, W. C. and Biegler, L. T. (2001). Design for
+model parameter uncertainty using nonlinear confidence regions. AIChE Journal,
+47(8), 1794-1804.
+"""
+
+from pyomo.common.dependencies import pandas as pd
+import pyomo.environ as pyo
+from pyomo.contrib.parmest.experiment import Experiment
+
+
+class RooneyBieglerExperimentDoE(Experiment):
+    def __init__(self, data=None, theta=None):
+        if data is None:
+            self.data = {}
+            self.data['hour'] = 1
+            self.data['y'] = 8.3
+        else:
+            self.data = data
+        if theta is None:
+            self.theta = {}
+            self.theta['asymptote'] = 19.143
+            self.theta['rate constant'] = 0.5311
+        else:
+            self.theta = theta
+        self.model = None
+
+    def create_model(self):
+        # Creates Roony-Biegler model for
+        # individual data points as
+        # an experimental decision.
+        m = self.model = pyo.ConcreteModel()
+
+        # Specify the unknown parameters
+        m.asymptote = pyo.Var(initialize=self.theta['asymptote'])
+        m.rate_constant = pyo.Var(initialize=self.theta['rate constant'])
+
+        # Fix the unknown parameters
+        m.asymptote.fix()
+        m.rate_constant.fix()
+
+        # Add the experiment inputs
+        m.hour = pyo.Var(initialize=self.data['hour'], bounds=(0, 10))
+
+        # Fix the experimental design variable
+        m.hour.fix()
+
+        # Add the experiment outputs
+        m.y = pyo.Var(initialize=self.data['y'])
+
+        # Add governing equation
+        m.response_function = pyo.Constraint(
+            expr=m.y - m.asymptote * (1 - pyo.exp(-m.rate_constant * m.hour)) == 0
+        )
+
+    def finalize_model(self):
+        m = self.model
+        pass
+
+    def label_model(self):
+        m = self.model
+
+        # Add y value as experiment output
+        m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.experiment_outputs[m.y] = m.y()
+
+        # Add measurement error associated with y
+        # We are assuming a flat error of 0.3
+        # or about 1-3 percent
+        m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.measurement_error[m.y] = 1
+
+        # Add hour as experiment input
+        # We are deciding when to sample
+        m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.experiment_inputs[m.hour] = m.hour()
+
+        # Adding the unknown parameters
+        m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.unknown_parameters.update(
+            (k, k.value) for k in [m.asymptote, m.rate_constant]
+        )
+
+    def get_labeled_model(self):
+        self.create_model()
+        self.finalize_model()
+        self.label_model()
+
+        return self.model