Pyomo.DoE: adding more verbose output for sensitivity analysis

pyomo/contrib/doe/doe.py

@@ -45,22 +45,14 @@
 from pyomo.contrib.sensitivity_toolbox.sens import get_dsdp
 
 import pyomo.environ as pyo
+from pyomo.contrib.doe.utils import (
+    check_FIM,
+    compute_FIM_metrics,
+    _SMALL_TOLERANCE_DEFINITENESS,
+)
 
 from pyomo.opt import SolverStatus
 
-# This small and positive tolerance is used when checking
-# if the prior is negative definite or approximately
-# indefinite. It is defined as a tolerance here to ensure
-# consistency between the code below and the tests. The
-# user should not need to adjust it.
-_SMALL_TOLERANCE_DEFINITENESS = 1e-6
-
-# This small and positive tolerance is used to check
-# the FIM is approximately symmetric. It is defined as
-# a tolerance here to ensure consistency between the code
-# below and the tests. The user should not need to adjust it.
-_SMALL_TOLERANCE_SYMMETRY = 1e-6
-
 
 class ObjectiveLib(Enum):
     determinant = "determinant"
@@ -1383,24 +1375,8 @@ def check_model_FIM(self, model=None, FIM=None):
                 )
             )
 
-        # Compute the eigenvalues of the FIM
-        evals = np.linalg.eigvals(FIM)
-
-        # Check if the FIM is positive definite
-        if np.min(evals) < -_SMALL_TOLERANCE_DEFINITENESS:
-            raise ValueError(
-                "FIM provided is not positive definite. It has one or more negative eigenvalue(s) less than -{:.1e}".format(
-                    _SMALL_TOLERANCE_DEFINITENESS
-                )
-            )
-
-        # Check if the FIM is symmetric
-        if not np.allclose(FIM, FIM.T, atol=_SMALL_TOLERANCE_SYMMETRY):
-            raise ValueError(
-                "FIM provided is not symmetric using absolute tolerance {}".format(
-                    _SMALL_TOLERANCE_SYMMETRY
-                )
-            )
+        # Check FIM is positive definite and symmetric
+        check_FIM(FIM)
 
         self.logger.info(
             "FIM provided matches expected dimensions from model and is approximately positive (semi) definite."
@@ -1455,7 +1431,7 @@ def update_FIM_prior(self, model=None, FIM=None):
 
         self.logger.info("FIM prior has been updated.")
 
-    # ToDo: Add an update function for the parameter values? --> closed loop parameter estimation?
+    # TODO: Add an update function for the parameter values? --> closed loop parameter estimation?
     # Or leave this to the user?????
     def update_unknown_parameter_values(self, model=None, param_vals=None):
         raise NotImplementedError(
@@ -1474,12 +1450,36 @@ def compute_FIM_full_factorial(
 
         Parameters
         ----------
-        model: model to perform the full factorial exploration on
-        design_ranges: dict of lists, of the form {<var_name>: [start, stop, numsteps]}
-        method: string to specify which method should be used
-                options are ``kaug`` and ``sequential``
+        model: DoE model, optional
+            model to perform the full factorial exploration on
+        design_ranges: dict
+            dictionary of lists, of the form {<var_name>: [start, stop, numsteps]}
+        method: str, optional
+            to specify which method should be used.
+            Options are ``kaug`` and ``sequential``
 
+        Returns
+        -------
+        fim_factorial_results: dict
+            a dictionary of the results with the following keys and their corresponding
+            values as a list.
+            - keys of model's experiment_inputs
+            - "log10 D-opt": list of log10(D-optimality)
+            - "log10 A-opt": list of log10(A-optimality)
+            - "log10 E-opt": list of log10(E-optimality)
+            - "log10 ME-opt": list of log10(ME-optimality)
+            - "eigval_min": list of minimum eigenvalues
+            - "eigval_max": list of maximum eigenvalues
+            - "det_FIM": list of determinants
+            - "trace_FIM": list of traces
+            - "solve_time": list of solve times
+
+        Raises
+        ------
+        ValueError
+            If the design_ranges' keys do not match the model's experiment_inputs' keys.
         """
+
         # Start timer
         sp_timer = TicTocTimer()
         sp_timer.tic(msg=None)
@@ -1514,15 +1514,19 @@ def compute_FIM_full_factorial(
                 "Design ranges keys must be a subset of experimental design names."
             )
 
-        # ToDo: Add more objective types? i.e., modified-E; G-opt; V-opt; etc?
-        # ToDo: Also, make this a result object, or more user friendly.
+        # TODO: Add more objective types? i.e., modified-E; G-opt; V-opt; etc?
+        # TODO: Also, make this a result object, or more user friendly.
         fim_factorial_results = {k.name: [] for k, v in model.experiment_inputs.items()}
         fim_factorial_results.update(
             {
                 "log10 D-opt": [],
                 "log10 A-opt": [],
                 "log10 E-opt": [],
                 "log10 ME-opt": [],
+                "eigval_min": [],
+                "eigval_max": [],
+                "det_FIM": [],
+                "trace_FIM": [],
                 "solve_time": [],
             }
         )
@@ -1584,24 +1588,9 @@ def compute_FIM_full_factorial(
 
             FIM = self._computed_FIM
 
-            # Compute and record metrics on FIM
-            D_opt = np.log10(np.linalg.det(FIM))
-            A_opt = np.log10(np.trace(FIM))
-            E_vals, E_vecs = np.linalg.eig(FIM)  # Grab eigenvalues
-            E_ind = np.argmin(E_vals.real)  # Grab index of minima to check imaginary
-            # Warn the user if there is a ``large`` imaginary component (should not be)
-            if abs(E_vals.imag[E_ind]) > 1e-8:
-                self.logger.warning(
-                    "Eigenvalue has imaginary component greater than 1e-6, contact developers if this issue persists."
-                )
-
-            # If the real value is less than or equal to zero, set the E_opt value to nan
-            if E_vals.real[E_ind] <= 0:
-                E_opt = np.nan
-            else:
-                E_opt = np.log10(E_vals.real[E_ind])
-
-            ME_opt = np.log10(np.linalg.cond(FIM))
+            det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt = (
+                compute_FIM_metrics(FIM)
+            )
 
             # Append the values for each of the experiment inputs
             for k, v in model.experiment_inputs.items():
@@ -1611,6 +1600,10 @@ def compute_FIM_full_factorial(
             fim_factorial_results["log10 A-opt"].append(A_opt)
             fim_factorial_results["log10 E-opt"].append(E_opt)
             fim_factorial_results["log10 ME-opt"].append(ME_opt)
+            fim_factorial_results["eigval_min"].append(E_vals.min())
+            fim_factorial_results["eigval_max"].append(E_vals.max())
+            fim_factorial_results["det_FIM"].append(det_FIM)
+            fim_factorial_results["trace_FIM"].append(trace_FIM)
             fim_factorial_results["solve_time"].append(time_set[-1])
 
         self.fim_factorial_results = fim_factorial_results

pyomo/contrib/doe/examples/reactor_example.py

@@ -14,13 +14,37 @@
 from pyomo.contrib.doe import DesignOfExperiments
 
 import pyomo.environ as pyo
+import pyomo.environ as pyo
+import idaes.core.solvers.get_solver
 
 import json
 
 
 # Example for sensitivity analysis on the reactor experiment
 # After sensitivity analysis is done, we perform optimal DoE
-def run_reactor_doe():
+def run_reactor_doe(
+    n_points_for_design=9,
+    compute_FIM_full_factorial=True,
+    plot_factorial_results=True,
+    save_plots=True,
+    run_optimal_doe=True,
+):
+    """
+    This function demonstrates how to perform sensitivity analysis on the reactor
+
+    Parameters
+    ----------
+    n_points_for_design : int, optional
+        number of points to use for the design ranges, by default 9
+    compute_FIM_full_factorial : bool, optional
+        whether to compute the full factorial design, by default True
+    plot_factorial_results : bool, optional
+        whether to plot the results of the full factorial design, by default True
+    save_plots : bool, optional
+        whether to save draw_factorial_figure plots, by default True
+    run_optimal_doe : bool, optional
+        whether to run the optimal DoE, by default True
+    """
     # Read in file
     DATA_DIR = pathlib.Path(__file__).parent
     file_path = DATA_DIR / "result.json"
@@ -66,63 +90,76 @@ def run_reactor_doe():
         _Cholesky_option=True,
         _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,
-    )
+    if compute_FIM_full_factorial:
+        # Make design ranges to compute the full factorial design
+        design_ranges = {
+            "CA[0]": [1, 5, n_points_for_design],
+            "T[0]": [300, 700, n_points_for_design],
+        }
+
+        # Compute the full factorial design with the sequential FIM calculation
+        doe_obj.compute_FIM_full_factorial(
+            design_ranges=design_ranges, method="sequential"
+        )
+    if plot_factorial_results:
+        if save_plots:
+            figure_file_name = "example_reactor_compute_FIM"
+        else:
+            figure_file_name = None
+
+        # 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=figure_file_name,
+            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]
+    if run_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(
+            ("\tTemperature values: [" + "{:.2f}, " * 8 + "{:.2f}]").format(
+                *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("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["Experiment Design Names"])
+
+        ###################
+        # End optimal DoE
 
-    ###################
-    # End optimal DoE
+    return doe_obj
 
 
 if __name__ == "__main__":