adding E-opt reformulation to doe

Added a reformulation suggested in a textbook that would appear to be really, really bad at first glance.

Author: djlaky

pyomo/contrib/doe/doe.py

@@ -1135,6 +1135,7 @@ def create_objective_function(self, model=None):
         if self.objective_option not in [
             ObjectiveLib.determinant,
             ObjectiveLib.trace,
+            ObjectiveLib.minimum_eigenvalue,
             ObjectiveLib.zero,
         ]:
             raise AttributeError(
@@ -1237,6 +1238,15 @@ def determinant_general(b):
                 for d in range(len(list_p))
             )
             return m.determinant == det_perm
+        
+        def min_eig_diagonal_constraints(b, p):
+            """
+            Set bounding constraints for minimum diagonal
+            value in the FIM.
+            """
+            m = b.model()
+            return m.fim[p, p] >= m.minimum_diagonal
+
 
         if self.Cholesky_option and self.objective_option == ObjectiveLib.determinant:
             model.obj_cons.cholesky_cons = pyo.Constraint(
@@ -1265,6 +1275,26 @@ def determinant_general(b):
                 expr=pyo.log10(model.trace), sense=pyo.maximize
             )
 
+        elif self.objective_option == ObjectiveLib.minimum_eigenvalue:
+            # Use the fact that the minimum eigenvalue is
+            # smaller than or equal to the minimum diagonal
+            # element. Inspired by formulation 7.27 in 
+            # Convex Optimization by Boyd and Vandenberghe
+            #
+            # Minimum eigenvalue bounded by minimum diagonal element
+            # To-Do: initialize minimum diagonal better
+            model.minimum_diagonal = pyo.Var(
+                initialize=1, bounds=(small_number, None)
+            )
+            # Add constraints on each diagonal element
+            model.obj_cons.minimum_diagonal_constraints = pyo.Constraint(
+                model.parameter_names, rule=min_eig_diagonal_constraints
+            )
+            model.objective = pyo.Objective(
+                expr=pyo.log10(model.minimum_diagonal), sense=pyo.maximize
+            )
+
+
         elif self.objective_option == ObjectiveLib.zero:
             # add dummy objective function
             model.objective = pyo.Objective(expr=0)

pyomo/contrib/doe/examples/reactor_example_E-opt.py

@@ -0,0 +1,137 @@
+#  ___________________________________________________________________________
+#
+#  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.
+#  ___________________________________________________________________________
+from pyomo.common.dependencies import numpy as np
+
+from pyomo.contrib.doe.examples.reactor_experiment import ReactorExperiment
+from pyomo.contrib.doe import DesignOfExperiments
+
+import pyomo.environ as pyo
+
+import json
+from pathlib import Path
+
+
+# Example for sensitivity analysis on the reactor experiment
+# After sensitivity analysis is done, we perform optimal DoE
+def run_reactor_doe():
+    # Read in file
+    DATA_DIR = Path(__file__).parent
+    file_path = DATA_DIR / "result.json"
+
+    with open(file_path) as f:
+        data_ex = json.load(f)
+
+    # Put temperature control time points into correct format for reactor experiment
+    data_ex["control_points"] = {
+        float(k): v for k, v in data_ex["control_points"].items()
+    }
+
+    # Create a ReactorExperiment object; data and discretization information are part
+    # of the constructor of this object
+    experiment = ReactorExperiment(data=data_ex, nfe=10, ncp=3)
+
+    # Use a central difference, with step size 1e-3
+    fd_formula = "central"
+    step_size = 1e-3
+
+    # Use the minimum_eigenvalue objective with scaled sensitivity matrix
+    objective_option = "minimum_eigenvalue"
+    scale_nominal_param_value = True
+
+    # 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=None,
+        jac_initial=None,
+        fim_initial=None,
+        L_diagonal_lower_bound=1e-7,
+        solver=None,
+        tee=True,
+        get_labeled_model_args=None,
+        _Cholesky_option=False,
+        _only_compute_fim_lower=True,
+    )
+
+    # Make design ranges to compute the full factorial design
+    design_ranges = {"CA[0]": [1, 5, 9], "T[0]": [300, 700, 9]}
+
+    # Compute the full factorial design with the sequential FIM calculation
+    doe_obj.compute_FIM_full_factorial(design_ranges=design_ranges, method="sequential")
+
+    # Plot the results
+    # doe_obj.draw_factorial_figure(
+    #     sensitivity_design_variables=["CA[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,
+    # )
+
+    ###########################
+    # End sensitivity analysis
+
+    # Begin optimal DoE
+    ####################
+    doe_obj.run_doe()
+
+    # Print out a results summary
+    print("Optimal experiment values: ")
+    print(
+        "\tInitial concentration: {:.2f}".format(
+            doe_obj.results["Experiment Design"][0]
+        )
+    )
+    print(
+        ("\tTemperature values: [" + "{:.2f}, " * 8 + "{:.2f}]").format(
+            *doe_obj.results["Experiment Design"][1:]
+        )
+    )
+    #print(doe_obj.results["Experiment Design"][1:])
+    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"])
+    print(doe_obj.results["log10 E-opt"])
+
+    doe_obj.model.obj_cons.pprint()
+    doe_obj.model.fim.pprint()
+    doe_obj.model.minimum_diagonal.pprint()
+    doe_obj.model.objective.pprint()
+
+    ###################
+    # End optimal DoE
+
+
+if __name__ == "__main__":
+    run_reactor_doe()

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