[Pyomo.DoE] Add simultaneous design of multiple experiments

pyomo/contrib/doe/grey_box_utilities.py

@@ -44,7 +44,14 @@
 class FIMExternalGreyBox(
     ExternalGreyBoxModel if (scipy_available and numpy_available) else object
 ):
-    def __init__(self, doe_object, objective_option="determinant", logger_level=None):
+    def __init__(
+        self,
+        doe_object=None,
+        objective_option="determinant",
+        logger_level=None,
+        parameter_names=None,
+        fim_initial=None,
+    ):
         """
         Grey box model for metrics on the FIM. This methodology reduces
         numerical complexity for the computation of FIM metrics related
@@ -68,17 +75,35 @@ def __init__(self, doe_object, objective_option="determinant", logger_level=None
            default: None, or equivalently, use the logging level of doe_object.
 
            NOTE: Use logging.DEBUG for all messages.
+        parameter_names:
+           Optional ordered iterable of parameter labels. When provided, this
+           lets the grey box object operate on any FIM source with the same
+           ordering instead of assuming the data must come from
+           ``doe_object.model.parameter_names``. This is needed for the
+           multi-experiment grey box path because the linked FIM lives on a
+           scenario block (``scenario.total_fim``), while ``doe_object.model``
+           is the top-level container and does not own ``parameter_names``
+           directly.
+        fim_initial:
+           Optional dense, symmetric FIM used to seed the grey box inputs. This
+           is required when ``doe_object`` is not provided.
         """
 
-        if doe_object is None:
+        if doe_object is None and (parameter_names is None or fim_initial is None):
             raise ValueError(
-                "DoE Object must be provided to build external grey box of the FIM."
+                "Either ``doe_object`` or both ``parameter_names`` and "
+                "``fim_initial`` must be provided to build the FIM grey box."
             )
 
         self.doe_object = doe_object
 
-        # Grab parameter list from the doe_object model
-        self._param_names = [i for i in self.doe_object.model.parameter_names]
+        # Grab parameter ordering from the explicit arguments when available.
+        # Multi-experiment optimization passes the aggregated scenario FIM
+        # directly, so we should not assume the linked FIM always shares the
+        # same location as self.doe_object.model.
+        if parameter_names is None:
+            parameter_names = self.doe_object.model.parameter_names
+        self._param_names = [i for i in parameter_names]
         self._n_params = len(self._param_names)
 
         # Check if the doe_object has model components that are required
@@ -93,15 +118,21 @@ def __init__(self, doe_object, objective_option="determinant", logger_level=None
 
         # If logger level is None, use doe_object's logger level
         if logger_level is None:
-            logger_level = doe_object.logger.level
-
+            if doe_object is not None:
+                logger_level = doe_object.logger.level
+            else:
+                logger_level = logging.WARNING
         self.logger.setLevel(level=logger_level)
 
         # Set initial values for inputs
         # Need a mask structure
-        self._masking_matrix = np.triu(np.ones_like(self.doe_object.fim_initial))
+        if fim_initial is None:
+            fim_initial = self.doe_object.fim_initial
+        fim_initial = np.asarray(fim_initial, dtype=np.float64)
+
+        self._masking_matrix = np.triu(np.ones_like(fim_initial))
         self._input_values = np.asarray(
-            self.doe_object.fim_initial[self._masking_matrix > 0], dtype=np.float64
+            fim_initial[self._masking_matrix > 0], dtype=np.float64
         )
         self._n_inputs = len(self._input_values)
 

pyomo/contrib/doe/tests/experiment_class_example_flags.py

@@ -70,3 +70,80 @@ def get_labeled_model(self):
         m.bad_con_2 = pyo.Constraint(expr=m.hour <= 0.0)
 
         return m
+
+
+class RooneyBieglerMultiExperiment(RooneyBieglerExperiment):
+    """
+    Experiment class based on the multi-experiment Rooney-Biegler prototype.
+
+    This mirrors the implementation in
+    ``examples/multiexperiment-prototype/rooney_biegler_multiexperiment.py``
+    while allowing test-time control over initial hour and bounds.
+    """
+
+    def __init__(
+        self, hour=2.0, y=10.0, theta=None, measure_error=0.1, hour_bounds=(1.0, 10.0)
+    ):
+        data = {'hour': hour, 'y': y}
+        super().__init__(data=data, measure_error=measure_error, theta=theta)
+        self.hour_bounds = hour_bounds
+
+    def get_labeled_model(self):
+        m = super().get_labeled_model()
+        hour_lb, hour_ub = self.hour_bounds
+        m.hour.setlb(hour_lb)
+        m.hour.setub(hour_ub)
+
+        if hasattr(m, "sym_break_cons"):
+            m.sym_break_cons.clear()
+        else:
+            m.sym_break_cons = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.sym_break_cons[m.hour] = None
+        return m
+
+
+class RooneyBieglerMultiInputExperimentFlag(RooneyBieglerExperiment):
+    """
+    Two-input Rooney-Biegler style experiment for symmetry-breaking tests.
+
+    Parameters
+    ----------
+    sym_break_flag : int
+        0 -> do not add ``sym_break_cons`` suffix
+        1 -> add single marker (hour)
+        2 -> add multiple markers (hour and temp)
+    """
+
+    def __init__(self, hour=2.0, temp=300.0, y=10.0, sym_break_flag=1):
+        data = {'hour': hour, 'y': y}
+        super().__init__(
+            data=data, measure_error=0.1, theta={'asymptote': 15, 'rate_constant': 0.5}
+        )
+        self.hour = hour
+        self.temp = temp
+        self.sym_break_flag = sym_break_flag
+
+    def get_labeled_model(self):
+        m = super().get_labeled_model()
+
+        m.hour.setlb(1.0)
+        m.hour.setub(10.0)
+        m.temp = pyo.Var(initialize=self.temp, bounds=(280.0, 340.0))
+        m.temp.fix()
+
+        # Replace base Rooney-Biegler response with two-input synthetic variant
+        # used only for symmetry-breaking tests.
+        m.del_component(m.response_function)
+        m.response_function = pyo.Constraint(
+            expr=m.y
+            == m.asymptote * (1 - pyo.exp(-m.rate_constant * m.hour)) + 0.01 * m.temp
+        )
+        m.experiment_inputs[m.temp] = self.temp
+
+        if self.sym_break_flag in (1, 2):
+            m.sym_break_cons = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+            m.sym_break_cons[m.hour] = None
+            if self.sym_break_flag == 2:
+                m.sym_break_cons[m.temp] = None
+
+        return m