From bf27dc5ad4fa2209b18d7251f95e7fdd658a2905 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Tue, 25 Nov 2025 12:22:48 -0500
Subject: [PATCH 01/17] delete the rooney_biegler_with_constraint.py and
 updated the test_parmest and test_examples.py files.

---
 .../rooney_biegler_with_constraint.py         | 109 ------------------
 pyomo/contrib/parmest/tests/test_examples.py  |   7 --
 pyomo/contrib/parmest/tests/test_parmest.py   |   2 +-
 scenarios.csv                                 |  11 ++
 4 files changed, 12 insertions(+), 117 deletions(-)
 delete mode 100644 pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler_with_constraint.py
 create mode 100644 scenarios.csv

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler_with_constraint.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler_with_constraint.py
deleted file mode 100644
index 3416682dca3..00000000000
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler_with_constraint.py
+++ /dev/null
@@ -1,109 +0,0 @@
-#  ___________________________________________________________________________
-#
-#  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
-
-
-def rooney_biegler_model_with_constraint(data):
-    model = pyo.ConcreteModel()
-
-    model.asymptote = pyo.Var(initialize=15)
-    model.rate_constant = pyo.Var(initialize=0.5)
-
-    model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True)
-    model.y = pyo.Param(within=pyo.PositiveReals, mutable=True)
-
-    model.response_function = pyo.Var(data.hour, initialize=0.0)
-
-    # changed from expression to constraint
-    def response_rule(m, h):
-        return m.response_function[h] == m.asymptote * (
-            1 - pyo.exp(-m.rate_constant * h)
-        )
-
-    model.response_function_constraint = pyo.Constraint(data.hour, rule=response_rule)
-
-    def SSE_rule(m):
-        return sum(
-            (data.y[i] - m.response_function[data.hour[i]]) ** 2 for i in data.index
-        )
-
-    model.SSE = pyo.Objective(rule=SSE_rule, sense=pyo.minimize)
-
-    return model
-
-
-class RooneyBieglerExperiment(Experiment):
-
-    def __init__(self, data):
-        self.data = data
-        self.model = None
-
-    def create_model(self):
-        # rooney_biegler_model_with_constraint expects a dataframe
-        data_df = self.data.to_frame().transpose()
-        self.model = rooney_biegler_model_with_constraint(data_df)
-
-    def label_model(self):
-
-        m = self.model
-
-        m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
-        m.experiment_outputs.update(
-            [(m.hour, self.data['hour']), (m.y, self.data['y'])]
-        )
-
-        m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
-        m.unknown_parameters.update(
-            (k, pyo.ComponentUID(k)) for k in [m.asymptote, m.rate_constant]
-        )
-
-    def finalize_model(self):
-
-        m = self.model
-
-        # Experiment output values
-        m.hour = self.data['hour']
-        m.y = self.data['y']
-
-    def get_labeled_model(self):
-        self.create_model()
-        self.label_model()
-        self.finalize_model()
-
-        return self.model
-
-
-def main():
-    # These were taken from Table A1.4 in Bates and Watts (1988).
-    data = pd.DataFrame(
-        data=[[1, 8.3], [2, 10.3], [3, 19.0], [4, 16.0], [5, 15.6], [7, 19.8]],
-        columns=['hour', 'y'],
-    )
-
-    model = rooney_biegler_model_with_constraint(data)
-    solver = pyo.SolverFactory('ipopt')
-    solver.solve(model)
-
-    print('asymptote = ', model.asymptote())
-    print('rate constant = ', model.rate_constant())
-
-
-if __name__ == '__main__':
-    main()
diff --git a/pyomo/contrib/parmest/tests/test_examples.py b/pyomo/contrib/parmest/tests/test_examples.py
index 8d9331a6166..cf052f76bc0 100644
--- a/pyomo/contrib/parmest/tests/test_examples.py
+++ b/pyomo/contrib/parmest/tests/test_examples.py
@@ -38,13 +38,6 @@ def test_model(self):
 
         rooney_biegler.main()
 
-    def test_model_with_constraint(self):
-        from pyomo.contrib.parmest.examples.rooney_biegler import (
-            rooney_biegler_with_constraint,
-        )
-
-        rooney_biegler_with_constraint.main()
-
     @unittest.skipUnless(pynumero_ASL_available, "test requires libpynumero_ASL")
     @unittest.skipUnless(seaborn_available, "test requires seaborn")
     def test_parameter_estimation_example(self):
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index db71d280f7c..e7de72fe59c 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -1349,7 +1349,7 @@ def test_covariance(self):
 @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
 class TestSquareInitialization_RooneyBiegler(unittest.TestCase):
     def setUp(self):
-        from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler_with_constraint import (
+        from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
             RooneyBieglerExperiment,
         )
 
diff --git a/scenarios.csv b/scenarios.csv
new file mode 100644
index 00000000000..af286781a20
--- /dev/null
+++ b/scenarios.csv
@@ -0,0 +1,11 @@
+Name,Probability,k1,k2,E1,E2
+ExpScen0,0.1,25.800350784967552,14.144215235968407,31505.74904933868,35000.0
+ExpScen1,0.1,25.1283730831486,149.99999951481198,31452.3366518825,41938.78130161935
+ExpScen2,0.1,22.225574065242643,130.92739780149637,30948.66911165926,41260.15420926035
+ExpScen3,0.1,100.0,149.9999996987801,35182.7313074055,41444.52600370866
+ExpScen4,0.1,82.99114366257251,45.95424665356903,34810.857217160396,38300.63334950135
+ExpScen5,0.1,100.0,150.0,35142.202191502525,41495.411057950805
+ExpScen6,0.1,2.8743643265327625,149.99999474412596,25000.0,41431.61195917287
+ExpScen7,0.1,2.754580914039618,14.381786098093363,25000.0,35000.0
+ExpScen8,0.1,2.8743643265327625,149.99999474412596,25000.0,41431.61195917287
+ExpScen9,0.1,2.6697808222410906,150.0,25000.0,41514.74761132933

From 5cae869215fe191bbf679dc6916292f15d6b3a48 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Tue, 25 Nov 2025 13:52:08 -0500
Subject: [PATCH 02/17] fix the `parmest/rooney_biegler.py`. Previously, it was
 not calling the experiment class to compute the parameters. it is now fixed.

---
 .../examples/rooney_biegler/rooney_biegler.py     | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index 9f77d75a501..22cfec24c47 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -50,15 +50,17 @@ def __init__(self, data, measure_error=None):
 
     def create_model(self):
         # rooney_biegler_model expects a dataframe
-        data_df = self.data.to_frame().transpose()
-        self.model = rooney_biegler_model(data_df)
+        # data_df = self.data.to_frame().transpose()
+        self.model = rooney_biegler_model(data=self.data)
 
     def label_model(self):
 
         m = self.model
 
         m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
-        m.experiment_outputs.update([(m.y[self.data['hour']], self.data['y'])])
+        m.experiment_outputs.update(
+            [(m.y[t], self.data.loc[i, 'y']) for i, t in enumerate(self.data['hour'])]
+        )
 
         m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.unknown_parameters.update(
@@ -66,7 +68,9 @@ def label_model(self):
         )
 
         m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
-        m.measurement_error.update([(m.y[self.data['hour']], self.measure_error)])
+        m.measurement_error.update(
+            [(m.y[t], self.measure_error) for t in self.data['hour']]
+        )
 
     def finalize_model(self):
 
@@ -90,7 +94,8 @@ def main():
         columns=['hour', 'y'],
     )
 
-    model = rooney_biegler_model(data)
+    # model = rooney_biegler_model(data)
+    model = RooneyBieglerExperiment(data).get_labeled_model()
     solver = pyo.SolverFactory('ipopt')
     solver.solve(model)
 

From fb38feed24cfb16c19bbc318f80daec9e599b058 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Tue, 25 Nov 2025 17:15:52 -0500
Subject: [PATCH 03/17] 1. delete name==__main__ block in rooney_biegler.py,
 and corresponding test in test_examples.py

2. Fix the index issue in test_parmest. New column mismathc issue arises.
---
 .../parameter_estimation_example.py           | 31 ++++++++------
 .../examples/rooney_biegler/rooney_biegler.py | 42 ++++---------------
 pyomo/contrib/parmest/tests/test_examples.py  |  5 ---
 pyomo/contrib/parmest/tests/test_parmest.py   | 17 ++++++--
 4 files changed, 40 insertions(+), 55 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py b/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py
index 8f56a9d849c..dec044e2106 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py
@@ -9,6 +9,7 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
+from numpy import cov
 from pyomo.common.dependencies import pandas as pd
 import pyomo.contrib.parmest.parmest as parmest
 from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
@@ -26,9 +27,7 @@ def main():
 
     # Sum of squared error function
     def SSE(model):
-        expr = (
-            model.experiment_outputs[model.y[model.hour]] - model.y[model.hour]
-        ) ** 2
+        expr = (model.experiment_outputs[model.y] - model.y) ** 2
         return expr
 
     # Create an experiment list
@@ -45,16 +44,17 @@ def SSE(model):
 
     # Parameter estimation and covariance
     n = 6  # total number of data points used in the objective (y in 6 scenarios)
-    obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=n)
+    obj, theta = pest.theta_est()
+    # cov = pest.cov_est()
 
-    # Plot theta estimates using a multivariate Gaussian distribution
-    parmest.graphics.pairwise_plot(
-        (theta, cov, 100),
-        theta_star=theta,
-        alpha=0.8,
-        distributions=['MVN'],
-        title='Theta estimates within 80% confidence region',
-    )
+    if parmest.graphics.seaborn_available:
+        parmest.graphics.pairwise_plot(
+            (theta, cov, 100),
+            theta_star=theta,
+            alpha=0.8,
+            distributions=['MVN'],
+            title='Theta estimates within 80% confidence region',
+        )
 
     # Assert statements compare parameter estimation (theta) to an expected value
     relative_error = abs(theta['asymptote'] - 19.1426) / 19.1426
@@ -62,6 +62,11 @@ def SSE(model):
     relative_error = abs(theta['rate_constant'] - 0.5311) / 0.5311
     assert relative_error < 0.01
 
+    return obj, theta
+
 
 if __name__ == "__main__":
-    main()
+    obj, theta = main()
+    print("Estimated parameters (theta):", theta)
+    print("Objective function value at theta:", obj)
+    # print("Covariance of parameter estimates:", cov)
diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index 22cfec24c47..3f7799b9b60 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -26,17 +26,17 @@ def rooney_biegler_model(data):
     model.asymptote = pyo.Var(initialize=15)
     model.rate_constant = pyo.Var(initialize=0.5)
 
-    model.y = pyo.Var(data.hour, within=pyo.PositiveReals, initialize=5)
+    model.y = pyo.Var(within=pyo.PositiveReals, initialize=5)
 
     def response_rule(m, h):
-        return m.y[h] == m.asymptote * (1 - pyo.exp(-m.rate_constant * h))
+        return m.y == m.asymptote * (1 - pyo.exp(-m.rate_constant * h))
 
     model.response_function = pyo.Constraint(data.hour, rule=response_rule)
 
-    def SSE_rule(m):
-        return sum((data.y[i] - m.y[data.hour[i]]) ** 2 for i in data.index)
+    # def SSE_rule(m):
+    #     return (data.y - m.y) ** 2
 
-    model.SSE = pyo.Objective(rule=SSE_rule, sense=pyo.minimize)
+    # model.SSE = pyo.Objective(rule=SSE_rule, sense=pyo.minimize)
 
     return model
 
@@ -50,17 +50,15 @@ def __init__(self, data, measure_error=None):
 
     def create_model(self):
         # rooney_biegler_model expects a dataframe
-        # data_df = self.data.to_frame().transpose()
-        self.model = rooney_biegler_model(data=self.data)
+        data_df = self.data.to_frame().transpose()
+        self.model = rooney_biegler_model(data_df)
 
     def label_model(self):
 
         m = self.model
 
         m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
-        m.experiment_outputs.update(
-            [(m.y[t], self.data.loc[i, 'y']) for i, t in enumerate(self.data['hour'])]
-        )
+        m.experiment_outputs.update([(m.y, self.data.loc['y'])])
 
         m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.unknown_parameters.update(
@@ -68,9 +66,7 @@ def label_model(self):
         )
 
         m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
-        m.measurement_error.update(
-            [(m.y[t], self.measure_error) for t in self.data['hour']]
-        )
+        m.measurement_error.update([(m.y, self.measure_error)])
 
     def finalize_model(self):
 
@@ -85,23 +81,3 @@ def get_labeled_model(self):
         self.finalize_model()
 
         return self.model
-
-
-def main():
-    # These were taken from Table A1.4 in Bates and Watts (1988).
-    data = pd.DataFrame(
-        data=[[1, 8.3], [2, 10.3], [3, 19.0], [4, 16.0], [5, 15.6], [7, 19.8]],
-        columns=['hour', 'y'],
-    )
-
-    # model = rooney_biegler_model(data)
-    model = RooneyBieglerExperiment(data).get_labeled_model()
-    solver = pyo.SolverFactory('ipopt')
-    solver.solve(model)
-
-    print('asymptote = ', model.asymptote())
-    print('rate constant = ', model.rate_constant())
-
-
-if __name__ == '__main__':
-    main()
diff --git a/pyomo/contrib/parmest/tests/test_examples.py b/pyomo/contrib/parmest/tests/test_examples.py
index cf052f76bc0..ce790b7ddb7 100644
--- a/pyomo/contrib/parmest/tests/test_examples.py
+++ b/pyomo/contrib/parmest/tests/test_examples.py
@@ -33,11 +33,6 @@ def setUpClass(self):
     def tearDownClass(self):
         pass
 
-    def test_model(self):
-        from pyomo.contrib.parmest.examples.rooney_biegler import rooney_biegler
-
-        rooney_biegler.main()
-
     @unittest.skipUnless(pynumero_ASL_available, "test requires libpynumero_ASL")
     @unittest.skipUnless(seaborn_available, "test requires seaborn")
     def test_parameter_estimation_example(self):
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index e7de72fe59c..e4291930034 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -9,6 +9,7 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
+from pyexpat import model
 import sys
 import os
 import subprocess
@@ -310,9 +311,7 @@ def setUp(self):
 
         # Sum of squared error function
         def SSE(model):
-            expr = (
-                model.experiment_outputs[model.y[model.hour]] - model.y[model.hour]
-            ) ** 2
+            expr = (model.experiment_outputs[model.y] - model.y) ** 2
             return expr
 
         # Create an experiment list
@@ -1363,8 +1362,9 @@ def setUp(self):
         def SSE(model):
             expr = (
                 model.experiment_outputs[model.y]
-                - model.response_function[model.experiment_outputs[model.hour]]
+                - pyo.value(model.response_function[model.hour])
             ) ** 2
+
             return expr
 
         exp_list = []
@@ -1378,9 +1378,18 @@ def SSE(model):
             exp_list, obj_function=SSE, solver_options=solver_options, tee=True
         )
 
+        # debugging print
+        o, j = self.pest.theta_est()
+        print("Initial run obj:", o)
+        print("Initial run theta:", j)
+
     def test_theta_est_with_square_initialization(self):
         obj_init = self.pest.objective_at_theta(initialize_parmest_model=True)
         objval, thetavals = self.pest.theta_est()
+        print("*" * 30)
+        print("objval:", objval)
+        print("thetavals:", thetavals)
+        print("*" * 30)
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
         self.assertAlmostEqual(

From 859137192526f247a82e836c3e5ab2cc10334bce Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Tue, 25 Nov 2025 17:41:05 -0500
Subject: [PATCH 04/17] change the SSE function in `test_parmest.py`.
 Everything works now and the test is successful except "FileNotFoundError:
 [Errno 2] No such file or directory: 'mpiexec'"

---
 .../parmest/examples/rooney_biegler/rooney_biegler.py  |  5 -----
 pyomo/contrib/parmest/tests/test_parmest.py            | 10 +---------
 2 files changed, 1 insertion(+), 14 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index 3f7799b9b60..f29d98d9388 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -33,11 +33,6 @@ def response_rule(m, h):
 
     model.response_function = pyo.Constraint(data.hour, rule=response_rule)
 
-    # def SSE_rule(m):
-    #     return (data.y - m.y) ** 2
-
-    # model.SSE = pyo.Objective(rule=SSE_rule, sense=pyo.minimize)
-
     return model
 
 
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index e4291930034..3cce39f56f1 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -1360,10 +1360,7 @@ def setUp(self):
 
         # Sum of squared error function
         def SSE(model):
-            expr = (
-                model.experiment_outputs[model.y]
-                - pyo.value(model.response_function[model.hour])
-            ) ** 2
+            expr = (model.experiment_outputs[model.y] - model.y) ** 2
 
             return expr
 
@@ -1378,11 +1375,6 @@ def SSE(model):
             exp_list, obj_function=SSE, solver_options=solver_options, tee=True
         )
 
-        # debugging print
-        o, j = self.pest.theta_est()
-        print("Initial run obj:", o)
-        print("Initial run theta:", j)
-
     def test_theta_est_with_square_initialization(self):
         obj_init = self.pest.objective_at_theta(initialize_parmest_model=True)
         objval, thetavals = self.pest.theta_est()

From 108b9b7be490a16b9ed1662018472f1ac6a2a439 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Wed, 26 Nov 2025 13:06:04 -0500
Subject: [PATCH 05/17] update parmest/ rooney_biegler, added theta in the
 constructor

---
 .../examples/rooney_biegler/rooney_biegler.py | 28 +++++++++++++++----
 1 file changed, 22 insertions(+), 6 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index f29d98d9388..f1d6449260b 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -20,11 +20,15 @@
 from pyomo.contrib.parmest.experiment import Experiment
 
 
-def rooney_biegler_model(data):
+def rooney_biegler_model(data, theta=None):
+    # if theta is  None:
+    #     theta = {}
+    #     theta['asymptote'] = 15
+    #     theta['rate constant'] = 0.5
     model = pyo.ConcreteModel()
 
-    model.asymptote = pyo.Var(initialize=15)
-    model.rate_constant = pyo.Var(initialize=0.5)
+    model.asymptote = pyo.Var(initialize=theta['asymptote'])
+    model.rate_constant = pyo.Var(initialize=theta['rate constant'])
 
     model.y = pyo.Var(within=pyo.PositiveReals, initialize=5)
 
@@ -38,15 +42,27 @@ def response_rule(m, h):
 
 class RooneyBieglerExperiment(Experiment):
 
-    def __init__(self, data, measure_error=None):
+    def __init__(self, data=None, measure_error=None, theta=None):
         self.data = data
+
+        if measure_error is None:
+            self.measure_error = 1
+        else:
+            self.measure_error = measure_error
+
+        if theta is None:
+            self.theta = {}
+            self.theta['asymptote'] = 15.0
+            self.theta['rate constant'] = 0.5
+        else:
+            self.theta = theta
+
         self.model = None
-        self.measure_error = measure_error
 
     def create_model(self):
         # rooney_biegler_model expects a dataframe
         data_df = self.data.to_frame().transpose()
-        self.model = rooney_biegler_model(data_df)
+        self.model = rooney_biegler_model(data_df, theta=self.theta)
 
     def label_model(self):
 

From 2afd81938892ca8e00d41a5393f6c6b9662f2b94 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Wed, 26 Nov 2025 13:06:23 -0500
Subject: [PATCH 06/17] Revert "update parmest/ rooney_biegler, added theta in
 the constructor"

This reverts commit 108b9b7be490a16b9ed1662018472f1ac6a2a439.
---
 .../examples/rooney_biegler/rooney_biegler.py | 28 ++++---------------
 1 file changed, 6 insertions(+), 22 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index f1d6449260b..f29d98d9388 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -20,15 +20,11 @@
 from pyomo.contrib.parmest.experiment import Experiment
 
 
-def rooney_biegler_model(data, theta=None):
-    # if theta is  None:
-    #     theta = {}
-    #     theta['asymptote'] = 15
-    #     theta['rate constant'] = 0.5
+def rooney_biegler_model(data):
     model = pyo.ConcreteModel()
 
-    model.asymptote = pyo.Var(initialize=theta['asymptote'])
-    model.rate_constant = pyo.Var(initialize=theta['rate constant'])
+    model.asymptote = pyo.Var(initialize=15)
+    model.rate_constant = pyo.Var(initialize=0.5)
 
     model.y = pyo.Var(within=pyo.PositiveReals, initialize=5)
 
@@ -42,27 +38,15 @@ def response_rule(m, h):
 
 class RooneyBieglerExperiment(Experiment):
 
-    def __init__(self, data=None, measure_error=None, theta=None):
+    def __init__(self, data, measure_error=None):
         self.data = data
-
-        if measure_error is None:
-            self.measure_error = 1
-        else:
-            self.measure_error = measure_error
-
-        if theta is None:
-            self.theta = {}
-            self.theta['asymptote'] = 15.0
-            self.theta['rate constant'] = 0.5
-        else:
-            self.theta = theta
-
         self.model = None
+        self.measure_error = measure_error
 
     def create_model(self):
         # rooney_biegler_model expects a dataframe
         data_df = self.data.to_frame().transpose()
-        self.model = rooney_biegler_model(data_df, theta=self.theta)
+        self.model = rooney_biegler_model(data_df)
 
     def label_model(self):
 

From 2684599cd70a548c0ed497befddac4a01ab4b4fc Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Wed, 26 Nov 2025 15:58:53 -0500
Subject: [PATCH 07/17] imported the parmest/rooney_biegler in  the
 doe/greybox_test.py file, and changed the file accordingly.

---
 pyomo/contrib/doe/tests/test_greybox.py       | 14 +++---
 .../examples/rooney_biegler/rooney_biegler.py | 44 ++++++++++++++-----
 2 files changed, 39 insertions(+), 19 deletions(-)

diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py
index 1336b266685..2da72ba41ef 100644
--- a/pyomo/contrib/doe/tests/test_greybox.py
+++ b/pyomo/contrib/doe/tests/test_greybox.py
@@ -32,8 +32,8 @@
     from pyomo.contrib.doe.examples.reactor_example import (
         ReactorExperiment as FullReactorExperiment,
     )
-    from pyomo.contrib.doe.examples.rooney_biegler_example import (
-        RooneyBieglerExperimentDoE,
+    from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
+        RooneyBieglerExperiment,
     )
 
 import pyomo.environ as pyo
@@ -304,7 +304,9 @@ def make_greybox_and_doe_objects_rooney_biegler(objective_option):
     fd_method = "central"
     obj_used = objective_option
 
-    experiment = RooneyBieglerExperimentDoE(data={'hour': 2, 'y': 10.3})
+    data = pd.DataFrame(data=[[2, 10.3]], columns=['hour', 'y'])
+
+    experiment = RooneyBieglerExperiment(data=data.loc[0, :])
 
     DoE_args = get_standard_args(experiment, fd_method, obj_used)
     DoE_args["use_grey_box_objective"] = True
@@ -318,13 +320,11 @@ def make_greybox_and_doe_objects_rooney_biegler(objective_option):
     # Add the grey box solver to DoE_args
     DoE_args["grey_box_solver"] = grey_box_solver
 
-    data = [[1, 8.3], [7, 19.8]]
+    data = pd.DataFrame(data=[[1, 8.3], [7, 19.8]], columns=['hour', 'y'])
     FIM_prior = np.zeros((2, 2))
     # Calculate prior using existing experiments
     for i in range(len(data)):
-        prev_experiment = RooneyBieglerExperimentDoE(
-            data={'hour': data[i][0], 'y': data[i][1]}
-        )
+        prev_experiment = RooneyBieglerExperiment(data=data.loc[i, :])
         doe_obj = DesignOfExperiments(
             **get_standard_args(prev_experiment, fd_method, obj_used)
         )
diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index f29d98d9388..443205a3197 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -20,33 +20,48 @@
 from pyomo.contrib.parmest.experiment import Experiment
 
 
-def rooney_biegler_model(data):
+def rooney_biegler_model(data, theta=None):
     model = pyo.ConcreteModel()
 
-    model.asymptote = pyo.Var(initialize=15)
-    model.rate_constant = pyo.Var(initialize=0.5)
+    if theta is None:
+        theta = {'asymptote': 15, 'rate_constant': 0.5}
 
-    model.y = pyo.Var(within=pyo.PositiveReals, initialize=5)
+    model.asymptote = pyo.Var(initialize=theta['asymptote'])
+    model.rate_constant = pyo.Var(initialize=theta['rate_constant'])
 
-    def response_rule(m, h):
-        return m.y == m.asymptote * (1 - pyo.exp(-m.rate_constant * h))
+    # Fix the unknown parameters
+    model.asymptote.fix()
+    model.rate_constant.fix()
 
-    model.response_function = pyo.Constraint(data.hour, rule=response_rule)
+    # Add the experiment inputs
+    model.hour = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10))
+
+    # Fix the experiment inputs
+    model.hour.fix()
+
+    # Add the response variable
+    model.y = pyo.Var(within=pyo.PositiveReals, initialize=data["y"].iloc[0])
+
+    def response_rule(m):
+        return m.y == m.asymptote * (1 - pyo.exp(-m.rate_constant * m.hour))
+
+    model.response_function = pyo.Constraint(rule=response_rule)
 
     return model
 
 
 class RooneyBieglerExperiment(Experiment):
 
-    def __init__(self, data, measure_error=None):
+    def __init__(self, data, measure_error=None, theta=None):
         self.data = data
         self.model = None
         self.measure_error = measure_error
+        self.theta = theta
 
     def create_model(self):
         # rooney_biegler_model expects a dataframe
         data_df = self.data.to_frame().transpose()
-        self.model = rooney_biegler_model(data_df)
+        self.model = rooney_biegler_model(data_df, theta=self.theta)
 
     def label_model(self):
 
@@ -63,12 +78,17 @@ def label_model(self):
         m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.measurement_error.update([(m.y, self.measure_error)])
 
+        # Add hour as an experiment input
+        m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.experiment_inputs.update([(m.hour, self.data.loc['hour'])])
+
     def finalize_model(self):
+        pass
 
-        m = self.model
+        # m = self.model
 
-        # Experiment input values
-        m.hour = self.data['hour']
+        # # Experiment input values
+        # m.hour = self.data['hour']
 
     def get_labeled_model(self):
         self.create_model()

From 0226681598feba2f386ee614c72c51758b70c210 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Wed, 26 Nov 2025 16:35:24 -0500
Subject: [PATCH 08/17] add comments to clarify suffix purposes, and solved the
 cov issue for pairplot. Now the parmest tests runs fine.

---
 .../parameter_estimation_example.py           | 20 +++++++++----------
 .../examples/rooney_biegler/rooney_biegler.py | 13 +++++++-----
 2 files changed, 18 insertions(+), 15 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py b/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py
index dec044e2106..3da7e4c47ca 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py
@@ -9,7 +9,6 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
-from numpy import cov
 from pyomo.common.dependencies import pandas as pd
 import pyomo.contrib.parmest.parmest as parmest
 from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
@@ -25,10 +24,12 @@ def main():
         columns=['hour', 'y'],
     )
 
+    # If desired, define a custom objective function, and then pass its name
+    # as a string to the Estimator constructor. e.g., obj_function= SSE
     # Sum of squared error function
-    def SSE(model):
-        expr = (model.experiment_outputs[model.y] - model.y) ** 2
-        return expr
+    # def SSE(model):
+    #     expr = (model.experiment_outputs[model.y] - model.y) ** 2
+    #     return expr
 
     # Create an experiment list
     exp_list = []
@@ -40,12 +41,11 @@ def SSE(model):
     # exp0_model.pprint()
 
     # Create an instance of the parmest estimator
-    pest = parmest.Estimator(exp_list, obj_function=SSE)
+    pest = parmest.Estimator(exp_list, obj_function="SSE")
 
     # Parameter estimation and covariance
-    n = 6  # total number of data points used in the objective (y in 6 scenarios)
     obj, theta = pest.theta_est()
-    # cov = pest.cov_est()
+    cov = pest.cov_est()
 
     if parmest.graphics.seaborn_available:
         parmest.graphics.pairwise_plot(
@@ -62,11 +62,11 @@ def SSE(model):
     relative_error = abs(theta['rate_constant'] - 0.5311) / 0.5311
     assert relative_error < 0.01
 
-    return obj, theta
+    return obj, theta, cov
 
 
 if __name__ == "__main__":
-    obj, theta = main()
+    obj, theta, cov = main()
     print("Estimated parameters (theta):", theta)
     print("Objective function value at theta:", obj)
-    # print("Covariance of parameter estimates:", cov)
+    print("Covariance of parameter estimates:", cov)
diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index 443205a3197..31244a2a558 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -67,29 +67,32 @@ def label_model(self):
 
         m = self.model
 
+        # Add experiment outputs as a suffix
+        # Experiment outputs suffix is required for parest
         m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.experiment_outputs.update([(m.y, self.data.loc['y'])])
 
+        # Add unknown parameters as a suffix
+        # Unknown parameters suffix is required for both Pyomo.DoE and parmest
         m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.unknown_parameters.update(
             (k, pyo.ComponentUID(k)) for k in [m.asymptote, m.rate_constant]
         )
 
+        # Add measurement error as a suffix
+        # Measurement error suffix is required for Pyomo.DoE and
+        #  `cov` estimation in parmest
         m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.measurement_error.update([(m.y, self.measure_error)])
 
         # Add hour as an experiment input
+        # Experiment inputs suffix is required for Pyomo.DoE
         m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.experiment_inputs.update([(m.hour, self.data.loc['hour'])])
 
     def finalize_model(self):
         pass
 
-        # m = self.model
-
-        # # Experiment input values
-        # m.hour = self.data['hour']
-
     def get_labeled_model(self):
         self.create_model()
         self.label_model()

From a0242ca998a71cb963c3ed6671ad3d0de3c95e7f Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Wed, 26 Nov 2025 17:20:57 -0500
Subject: [PATCH 09/17] add nominal theta values in doe/test_greybox.py to pass
 inside the experiment constructor.

---
 pyomo/contrib/doe/tests/test_greybox.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py
index 2da72ba41ef..bb175c517a7 100644
--- a/pyomo/contrib/doe/tests/test_greybox.py
+++ b/pyomo/contrib/doe/tests/test_greybox.py
@@ -305,8 +305,9 @@ def make_greybox_and_doe_objects_rooney_biegler(objective_option):
     obj_used = objective_option
 
     data = pd.DataFrame(data=[[2, 10.3]], columns=['hour', 'y'])
+    theta = {'asymptote': 19.143, 'rate_constant': 0.5311}
 
-    experiment = RooneyBieglerExperiment(data=data.loc[0, :])
+    experiment = RooneyBieglerExperiment(data=data.loc[0, :], theta=theta)
 
     DoE_args = get_standard_args(experiment, fd_method, obj_used)
     DoE_args["use_grey_box_objective"] = True
@@ -321,10 +322,11 @@ def make_greybox_and_doe_objects_rooney_biegler(objective_option):
     DoE_args["grey_box_solver"] = grey_box_solver
 
     data = pd.DataFrame(data=[[1, 8.3], [7, 19.8]], columns=['hour', 'y'])
+    theta = {'asymptote': 19.143, 'rate_constant': 0.5311}
     FIM_prior = np.zeros((2, 2))
     # Calculate prior using existing experiments
     for i in range(len(data)):
-        prev_experiment = RooneyBieglerExperiment(data=data.loc[i, :])
+        prev_experiment = RooneyBieglerExperiment(data=data.loc[i, :], theta=theta)
         doe_obj = DesignOfExperiments(
             **get_standard_args(prev_experiment, fd_method, obj_used)
         )

From ce24565adb34c069719a9f276e97faf806da1e13 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Wed, 3 Dec 2025 14:49:27 -0500
Subject: [PATCH 10/17] add measure_error  argument in rooney_biegler  in
 test_greybox. solve the skipping issue

---
 pyomo/contrib/doe/tests/test_greybox.py        | 18 +++++++++++++++---
 .../examples/rooney_biegler/rooney_biegler.py  |  2 +-
 2 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py
index bb175c517a7..beedcc39dfa 100644
--- a/pyomo/contrib/doe/tests/test_greybox.py
+++ b/pyomo/contrib/doe/tests/test_greybox.py
@@ -307,7 +307,9 @@ def make_greybox_and_doe_objects_rooney_biegler(objective_option):
     data = pd.DataFrame(data=[[2, 10.3]], columns=['hour', 'y'])
     theta = {'asymptote': 19.143, 'rate_constant': 0.5311}
 
-    experiment = RooneyBieglerExperiment(data=data.loc[0, :], theta=theta)
+    experiment = RooneyBieglerExperiment(
+        data=data.loc[0, :], theta=theta, measure_error=1
+    )
 
     DoE_args = get_standard_args(experiment, fd_method, obj_used)
     DoE_args["use_grey_box_objective"] = True
@@ -326,7 +328,9 @@ def make_greybox_and_doe_objects_rooney_biegler(objective_option):
     FIM_prior = np.zeros((2, 2))
     # Calculate prior using existing experiments
     for i in range(len(data)):
-        prev_experiment = RooneyBieglerExperiment(data=data.loc[i, :], theta=theta)
+        prev_experiment = RooneyBieglerExperiment(
+            data=data.loc[i, :], theta=theta, measure_error=1
+        )
         doe_obj = DesignOfExperiments(
             **get_standard_args(prev_experiment, fd_method, obj_used)
         )
@@ -349,7 +353,13 @@ def make_greybox_and_doe_objects_rooney_biegler(objective_option):
 # linear solvers.
 bad_message = "Invalid option encountered."
 cyipopt_call_working = True
-if numpy_available and scipy_available and ipopt_available and cyipopt_available:
+if (
+    numpy_available
+    and scipy_available
+    and ipopt_available
+    and cyipopt_available
+    and pandas_available
+):
     try:
         objective_option = "determinant"
         doe_object, _ = make_greybox_and_doe_objects_rooney_biegler(
@@ -1112,6 +1122,7 @@ def test_solve_A_optimality_trace_of_inverse(self):
     @unittest.skipIf(
         not cyipopt_call_working, "cyipopt is not properly accessing linear solvers"
     )
+    @unittest.skipIf(not pandas_available, "pandas is not available")
     def test_solve_E_optimality_minimum_eigenvalue(self):
         # Two locally optimal design points exist
         # (time, optimal objective value)
@@ -1153,6 +1164,7 @@ def test_solve_E_optimality_minimum_eigenvalue(self):
     @unittest.skipIf(
         not cyipopt_call_working, "cyipopt is not properly accessing linear solvers"
     )
+    @unittest.skipIf(not pandas_available, "pandas is not available")
     def test_solve_ME_optimality_condition_number(self):
         # Two locally optimal design points exist
         # (time, optimal objective value)
diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index 31244a2a558..bcc82d86b12 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -76,7 +76,7 @@ def label_model(self):
         # Unknown parameters suffix is required for both Pyomo.DoE and parmest
         m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.unknown_parameters.update(
-            (k, pyo.ComponentUID(k)) for k in [m.asymptote, m.rate_constant]
+            (k, pyo.value(k)) for k in [m.asymptote, m.rate_constant]
         )
 
         # Add measurement error as a suffix

From 99f36c591babac3b61af33945c20055603d13f52 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Wed, 3 Dec 2025 15:22:21 -0500
Subject: [PATCH 11/17] delete copies of `rooney_biegler` in `doe/examples`

---
 .../doe/examples/rooney_biegler_example.py    | 148 ------------------
 .../doe/examples/rooney_biegler_experiment.py | 100 ------------
 2 files changed, 248 deletions(-)
 delete mode 100644 pyomo/contrib/doe/examples/rooney_biegler_example.py
 delete mode 100644 pyomo/contrib/doe/examples/rooney_biegler_experiment.py

diff --git a/pyomo/contrib/doe/examples/rooney_biegler_example.py b/pyomo/contrib/doe/examples/rooney_biegler_example.py
deleted file mode 100644
index 954f0eed435..00000000000
--- a/pyomo/contrib/doe/examples/rooney_biegler_example.py
+++ /dev/null
@@ -1,148 +0,0 @@
-#  ___________________________________________________________________________
-#
-#  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()
diff --git a/pyomo/contrib/doe/examples/rooney_biegler_experiment.py b/pyomo/contrib/doe/examples/rooney_biegler_experiment.py
deleted file mode 100644
index 4918100b78d..00000000000
--- a/pyomo/contrib/doe/examples/rooney_biegler_experiment.py
+++ /dev/null
@@ -1,100 +0,0 @@
-#  ___________________________________________________________________________
-#
-#  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

From 95b235486365f74f53fce80b08d584db5d814dbd Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Wed, 3 Dec 2025 16:53:12 -0500
Subject: [PATCH 12/17] remove indexing problem of y data

---
 .../parmest/examples/rooney_biegler/bootstrap_example.py      | 4 +---
 .../examples/rooney_biegler/likelihood_ratio_example.py       | 4 +---
 .../contrib/parmest/examples/rooney_biegler/rooney_biegler.py | 4 ++--
 pyomo/contrib/parmest/tests/test_examples.py                  | 1 +
 4 files changed, 5 insertions(+), 8 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/bootstrap_example.py b/pyomo/contrib/parmest/examples/rooney_biegler/bootstrap_example.py
index c74ce32d5b4..2d7c4dd18c0 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/bootstrap_example.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/bootstrap_example.py
@@ -26,9 +26,7 @@ def main():
 
     # Sum of squared error function
     def SSE(model):
-        expr = (
-            model.experiment_outputs[model.y[model.hour]] - model.y[model.hour]
-        ) ** 2
+        expr = (model.experiment_outputs[model.y] - model.y) ** 2
         return expr
 
     # Create an experiment list
diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/likelihood_ratio_example.py b/pyomo/contrib/parmest/examples/rooney_biegler/likelihood_ratio_example.py
index 3b9209b366e..4d5814fd09c 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/likelihood_ratio_example.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/likelihood_ratio_example.py
@@ -27,9 +27,7 @@ def main():
 
     # Sum of squared error function
     def SSE(model):
-        expr = (
-            model.experiment_outputs[model.y[model.hour]] - model.y[model.hour]
-        ) ** 2
+        expr = (model.experiment_outputs[model.y] - model.y) ** 2
         return expr
 
     # Create an experiment list
diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index bcc82d86b12..6ce70a56dc0 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -70,7 +70,7 @@ def label_model(self):
         # Add experiment outputs as a suffix
         # Experiment outputs suffix is required for parest
         m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
-        m.experiment_outputs.update([(m.y, self.data.loc['y'])])
+        m.experiment_outputs.update([(m.y, self.data['y'])])
 
         # Add unknown parameters as a suffix
         # Unknown parameters suffix is required for both Pyomo.DoE and parmest
@@ -88,7 +88,7 @@ def label_model(self):
         # Add hour as an experiment input
         # Experiment inputs suffix is required for Pyomo.DoE
         m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
-        m.experiment_inputs.update([(m.hour, self.data.loc['hour'])])
+        m.experiment_inputs.update([(m.hour, self.data['hour'])])
 
     def finalize_model(self):
         pass
diff --git a/pyomo/contrib/parmest/tests/test_examples.py b/pyomo/contrib/parmest/tests/test_examples.py
index ce790b7ddb7..987638de60b 100644
--- a/pyomo/contrib/parmest/tests/test_examples.py
+++ b/pyomo/contrib/parmest/tests/test_examples.py
@@ -15,6 +15,7 @@
 from pyomo.contrib.pynumero.asl import AmplInterface
 from pyomo.opt import SolverFactory
 
+print("seaborn available", seaborn_available)
 ipopt_available = SolverFactory("ipopt").available()
 pynumero_ASL_available = AmplInterface.available()
 

From 92b8a25fe1266cbe8ce44f62d5f86eb6148546fd Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Thu, 4 Dec 2025 12:35:11 -0500
Subject: [PATCH 13/17] Solve the parmest/datarec.rst issue.

---
 doc/OnlineDocs/explanation/analysis/parmest/datarec.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/OnlineDocs/explanation/analysis/parmest/datarec.rst b/doc/OnlineDocs/explanation/analysis/parmest/datarec.rst
index 380bd7cb876..847b90022e9 100644
--- a/doc/OnlineDocs/explanation/analysis/parmest/datarec.rst
+++ b/doc/OnlineDocs/explanation/analysis/parmest/datarec.rst
@@ -44,8 +44,8 @@ The following example returns model values from a Pyomo Expression.
 
    >>> # Define objective
    >>> def SSE(model):
-   ...     expr = (model.experiment_outputs[model.y[model.hour]]
-   ...             - model.y[model.hour]
+   ...     expr = (model.experiment_outputs[model.y]
+   ...             - model.y
    ...            ) ** 2
    ...     return expr
 

From 8c2709fffb664eb7ff5bcaa1f4b7f1df5d1f608a Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Tue, 9 Dec 2025 15:50:18 -0500
Subject: [PATCH 14/17] Delete debugging print statement

---
 pyomo/contrib/parmest/tests/test_examples.py |  1 -
 pyomo/contrib/parmest/tests/test_parmest.py  |  4 ----
 scenarios.csv                                | 11 -----------
 3 files changed, 16 deletions(-)
 delete mode 100644 scenarios.csv

diff --git a/pyomo/contrib/parmest/tests/test_examples.py b/pyomo/contrib/parmest/tests/test_examples.py
index 987638de60b..ce790b7ddb7 100644
--- a/pyomo/contrib/parmest/tests/test_examples.py
+++ b/pyomo/contrib/parmest/tests/test_examples.py
@@ -15,7 +15,6 @@
 from pyomo.contrib.pynumero.asl import AmplInterface
 from pyomo.opt import SolverFactory
 
-print("seaborn available", seaborn_available)
 ipopt_available = SolverFactory("ipopt").available()
 pynumero_ASL_available = AmplInterface.available()
 
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index 3cce39f56f1..7f71a56e1d2 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -1378,10 +1378,6 @@ def SSE(model):
     def test_theta_est_with_square_initialization(self):
         obj_init = self.pest.objective_at_theta(initialize_parmest_model=True)
         objval, thetavals = self.pest.theta_est()
-        print("*" * 30)
-        print("objval:", objval)
-        print("thetavals:", thetavals)
-        print("*" * 30)
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
         self.assertAlmostEqual(
diff --git a/scenarios.csv b/scenarios.csv
deleted file mode 100644
index af286781a20..00000000000
--- a/scenarios.csv
+++ /dev/null
@@ -1,11 +0,0 @@
-Name,Probability,k1,k2,E1,E2
-ExpScen0,0.1,25.800350784967552,14.144215235968407,31505.74904933868,35000.0
-ExpScen1,0.1,25.1283730831486,149.99999951481198,31452.3366518825,41938.78130161935
-ExpScen2,0.1,22.225574065242643,130.92739780149637,30948.66911165926,41260.15420926035
-ExpScen3,0.1,100.0,149.9999996987801,35182.7313074055,41444.52600370866
-ExpScen4,0.1,82.99114366257251,45.95424665356903,34810.857217160396,38300.63334950135
-ExpScen5,0.1,100.0,150.0,35142.202191502525,41495.411057950805
-ExpScen6,0.1,2.8743643265327625,149.99999474412596,25000.0,41431.61195917287
-ExpScen7,0.1,2.754580914039618,14.381786098093363,25000.0,35000.0
-ExpScen8,0.1,2.8743643265327625,149.99999474412596,25000.0,41431.61195917287
-ExpScen9,0.1,2.6697808222410906,150.0,25000.0,41514.74761132933

From a48eeb45ed03aa518ce741a2403643da5a5f59d8 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Wed, 10 Dec 2025 16:23:12 -0500
Subject: [PATCH 15/17] Remove the commented custom SSE objective function.

---
 .../rooney_biegler/parameter_estimation_example.py    | 11 -----------
 1 file changed, 11 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py b/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py
index 3da7e4c47ca..ece018923d4 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/parameter_estimation_example.py
@@ -24,22 +24,11 @@ def main():
         columns=['hour', 'y'],
     )
 
-    # If desired, define a custom objective function, and then pass its name
-    # as a string to the Estimator constructor. e.g., obj_function= SSE
-    # Sum of squared error function
-    # def SSE(model):
-    #     expr = (model.experiment_outputs[model.y] - model.y) ** 2
-    #     return expr
-
     # Create an experiment list
     exp_list = []
     for i in range(data.shape[0]):
         exp_list.append(RooneyBieglerExperiment(data.loc[i, :]))
 
-    # View one model
-    # exp0_model = exp_list[0].get_labeled_model()
-    # exp0_model.pprint()
-
     # Create an instance of the parmest estimator
     pest = parmest.Estimator(exp_list, obj_function="SSE")
 

From edf0c64aae1a000f09e5d33ccaad8c41488af521 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Fri, 12 Dec 2025 16:06:42 -0500
Subject: [PATCH 16/17] Delete the label_model() in `rooney_biegler.py` and
 remove the import in test_parmest.

---
 .../parmest/examples/rooney_biegler/rooney_biegler.py        | 5 -----
 pyomo/contrib/parmest/tests/test_parmest.py                  | 1 -
 2 files changed, 6 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index 6ce70a56dc0..d71c6113d8a 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -90,12 +90,7 @@ def label_model(self):
         m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.experiment_inputs.update([(m.hour, self.data['hour'])])
 
-    def finalize_model(self):
-        pass
-
     def get_labeled_model(self):
         self.create_model()
         self.label_model()
-        self.finalize_model()
-
         return self.model
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index 7f71a56e1d2..ab5c49cfe16 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -9,7 +9,6 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
-from pyexpat import model
 import sys
 import os
 import subprocess

From 86d6b68fe5ebfc07e8c8e5e69836476368286319 Mon Sep 17 00:00:00 2001
From: Shuvashish Mondal <smondal@nd.edu>
Date: Fri, 12 Dec 2025 16:13:36 -0500
Subject: [PATCH 17/17] Correct typo

Co-authored-by: Bethany Nicholson <blnicho@users.noreply.github.com>
---
 pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
index d71c6113d8a..50d8cd72861 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py
@@ -68,7 +68,7 @@ def label_model(self):
         m = self.model
 
         # Add experiment outputs as a suffix
-        # Experiment outputs suffix is required for parest
+        # Experiment outputs suffix is required for parmest
         m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
         m.experiment_outputs.update([(m.y, self.data['y'])])