Add core symbolic gradient support to Pyomo.DoE

pyomo/contrib/doe/__init__.py

@@ -6,8 +6,8 @@
 # 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 .doe import DesignOfExperiments, ObjectiveLib, FiniteDifferenceStep
-from .utils import rescale_FIM
+from .doe import DesignOfExperiments, ObjectiveLib, FiniteDifferenceStep, GradientMethod
+from .utils import rescale_FIM, ExperimentGradients
 from .grey_box_utilities import FIMExternalGreyBox
 
 # Deprecation errors for old Pyomo.DoE interface classes and structures

pyomo/contrib/doe/examples/polynomial.py

@@ -0,0 +1,99 @@
+# ____________________________________________________________________________________
+#
+# Pyomo: Python Optimization Modeling Objects
+# Copyright (c) 2008-2026 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.
+# ____________________________________________________________________________________
+
+import pyomo.environ as pyo
+
+from pyomo.contrib.doe import DesignOfExperiments
+from pyomo.contrib.parmest.experiment import Experiment
+
+
+class PolynomialExperiment(Experiment):
+    """A small algebraic experiment used for symbolic-gradient testing."""
+
+    def __init__(self):
+        self.model = None
+
+    def get_labeled_model(self):
+        """Build and label the experiment model on first access."""
+        if self.model is None:
+            self.create_model()
+            self.label_experiment()
+        return self.model
+
+    def create_model(self):
+        """Define a polynomial model for testing symbolic sensitivities.
+
+        y = a*x1 + b*x2 + c*x1*x2 + d
+        """
+
+        m = self.model = pyo.ConcreteModel()
+
+        m.x1 = pyo.Var(bounds=(-5, 5), initialize=2.0)
+        m.x2 = pyo.Var(bounds=(-5, 5), initialize=3.0)
+
+        m.a = pyo.Var(bounds=(-5, 5), initialize=2)
+        m.a.fix()
+        m.b = pyo.Var(bounds=(-5, 5), initialize=-1)
+        m.b.fix()
+        m.c = pyo.Var(bounds=(-5, 5), initialize=0.5)
+        m.c.fix()
+        m.d = pyo.Var(bounds=(-5, 5), initialize=-1)
+        m.d.fix()
+
+        m.y = pyo.Var(initialize=0)
+
+        @m.Constraint()
+        def output_equation(m):
+            return m.y == m.a * m.x1 + m.b * m.x2 + m.c * m.x1 * m.x2 + m.d
+
+    def label_experiment(self):
+        """Attach the standard DoE suffixes to the polynomial model."""
+        m = self.model
+
+        m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.experiment_outputs[m.y] = None
+
+        m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.measurement_error[m.y] = 1
+
+        m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.experiment_inputs[m.x1] = None
+        m.experiment_inputs[m.x2] = None
+
+        m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.unknown_parameters.update((k, pyo.value(k)) for k in [m.a, m.b, m.c, m.d])
+
+
+def run_polynomial_doe():
+    """Run a small symbolic DoE FIM calculation for the polynomial model."""
+    experiment = PolynomialExperiment()
+
+    doe_obj = DesignOfExperiments(
+        experiment,
+        gradient_method="pynumero",
+        step=1e-3,
+        objective_option="determinant",
+        scale_constant_value=1,
+        scale_nominal_param_value=False,
+        prior_FIM=None,
+        jac_initial=None,
+        fim_initial=None,
+        L_diagonal_lower_bound=1e-7,
+        solver=pyo.SolverFactory("ipopt"),
+        tee=False,
+        get_labeled_model_args=None,
+        _Cholesky_option=True,
+        _only_compute_fim_lower=True,
+    )
+
+    return doe_obj.compute_FIM()
+
+
+if __name__ == "__main__":
+    run_polynomial_doe()

pyomo/contrib/doe/tests/test_doe_build.py

@@ -6,7 +6,6 @@
 # Solutions of Sandia, LLC, the U.S. Government retains certain rights in this
 # software.  This software is distributed under the 3-clause BSD License.
 # ____________________________________________________________________________________
-import json
 import os.path
 
 from pyomo.common.dependencies import (
@@ -17,20 +16,16 @@
     scipy_available,
 )
 
-from pyomo.common.fileutils import this_file_dir
 import pyomo.common.unittest as unittest
 
 if not (numpy_available and scipy_available):
     raise unittest.SkipTest("Pyomo.DoE needs scipy and numpy to run tests")
 
-if scipy_available:
-    from pyomo.contrib.doe import DesignOfExperiments
-    from pyomo.contrib.doe.examples.reactor_example import (
-        ReactorExperiment as FullReactorExperiment,
-    )
-    from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
-        RooneyBieglerExperiment,
-    )
+from pyomo.contrib.doe import DesignOfExperiments
+from pyomo.contrib.doe.examples.polynomial import PolynomialExperiment
+from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
+    RooneyBieglerExperiment,
+)
 
 from pyomo.contrib.doe.examples.rooney_biegler_doe_example import run_rooney_biegler_doe
 import pyomo.environ as pyo
@@ -39,14 +34,6 @@
 
 ipopt_available = SolverFactory("ipopt").available()
 
-currdir = this_file_dir()
-file_path = os.path.join(currdir, "..", "examples", "result.json")
-
-with open(file_path) as f:
-    data_ex = json.load(f)
-
-data_ex["control_points"] = {float(k): v for k, v in data_ex["control_points"].items()}
-
 
 def get_rooney_biegler_experiment():
     """Get a fresh RooneyBieglerExperiment instance for testing.
@@ -148,8 +135,7 @@ def get_standard_args(experiment, fd_method, obj_used):
 
 
 @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
-@unittest.skipIf(not numpy_available, "Numpy is not available")
-@unittest.skipIf(not scipy_available, "scipy is not available")
+# Tests require NumPy/SciPy, but availability is checked by the file-level SkipTest.
 @unittest.skipIf(not pandas_available, "pandas is not available")
 class TestDoeBuild(unittest.TestCase):
     def test_rooney_biegler_fd_central_check_fd_eqns(self):
@@ -267,6 +253,35 @@ def test_rooney_biegler_fd_forward_check_fd_eqns(self):
                 other_param_val = pyo.value(k)
                 self.assertAlmostEqual(other_param_val, v)
 
+    def test_polynomial_example_labels(self):
+        experiment = PolynomialExperiment()
+        model = experiment.get_labeled_model()
+
+        self.assertEqual(len(model.experiment_outputs), 1)
+        self.assertEqual(len(model.measurement_error), 1)
+        self.assertEqual(len(model.experiment_inputs), 2)
+        self.assertEqual(len(model.unknown_parameters), 4)
+
+        self.assertIn(model.y, model.experiment_outputs)
+        self.assertIn(model.y, model.measurement_error)
+        self.assertIn(model.x1, model.experiment_inputs)
+        self.assertIn(model.x2, model.experiment_inputs)
+
+    def test_polynomial_example_create_doe_model_pynumero(self):
+        experiment = PolynomialExperiment()
+
+        DoE_args = get_standard_args(
+            experiment, fd_method="central", obj_used="determinant"
+        )
+        DoE_args["gradient_method"] = "pynumero"
+
+        doe_obj = DesignOfExperiments(**DoE_args)
+        doe_obj.create_doe_model()
+
+        model = doe_obj.model
+        self.assertEqual(len(model.scenarios), 1)
+        self.assertTrue(hasattr(model.scenario_blocks[0], "jac_variables_wrt_param"))
+
     def test_rooney_biegler_fd_central_design_fixing(self):
         fd_method = "central"
         obj_used = "pseudo_trace"
@@ -510,21 +525,26 @@ def test_generate_blocks_without_model(self):
             )
 
 
-class TestReactorExample(unittest.TestCase):
-    def test_reactor_update_suffix_items(self):
-        """Test the reactor example with updating suffix items."""
-        from pyomo.contrib.doe.examples.update_suffix_doe_example import main
+class TestRooneyBieglerExample(unittest.TestCase):
+    @unittest.skipIf(not pandas_available, "pandas is not available")
+    def test_rooney_biegler_update_suffix_items(self):
+        """Test updating suffix items on the lightweight Rooney-Biegler model."""
+        from pyomo.contrib.parmest.utils.model_utils import update_model_from_suffix
+
+        experiment = get_rooney_biegler_experiment()
+        model = experiment.get_labeled_model()
+        suffix_obj = model.measurement_error
+        orig_vals = np.array(list(suffix_obj.values()))
+        new_vals = orig_vals + 1
 
-        # Run the reactor update suffix items example
-        suffix_obj, _, new_vals = main()
+        update_model_from_suffix(suffix_obj, new_vals)
 
-        # Check that the suffix object has been updated correctly
         for i, v in enumerate(suffix_obj.values()):
             self.assertAlmostEqual(v, new_vals[i], places=6)
 
 
 @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
-@unittest.skipIf(not numpy_available, "Numpy is not available")
+# Test requires NumPy, but availability is checked by the file-level SkipTest.
 @unittest.skipIf(not pandas_available, "pandas is not available")
 class TestDoEObjectiveOptions(unittest.TestCase):
     def test_trace_constraints(self):