[Depends on #3575] ParmEst: Add profile likelihood analysis for model identifiability

pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py

@@ -52,23 +52,6 @@ def simple_reaction_model(data):
     # fix all of the regressed parameters
     model.k.fix()
 
-    # ===================================================================
-    # Stage-specific cost computations
-    def ComputeFirstStageCost_rule(model):
-        return 0
-
-    model.FirstStageCost = Expression(rule=ComputeFirstStageCost_rule)
-
-    def AllMeasurements(m):
-        return (float(data['y']) - m.y) ** 2
-
-    model.SecondStageCost = Expression(rule=AllMeasurements)
-
-    def total_cost_rule(m):
-        return m.FirstStageCost + m.SecondStageCost
-
-    model.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize)
-
     return model
 
 
@@ -90,6 +73,8 @@ def label_model(self):
         m.experiment_outputs.update(
             [(m.x1, self.data['x1']), (m.x2, self.data['x2']), (m.y, self.data['y'])]
         )
+        m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.measurement_error.update([(m.y, None), (m.x1, None), (m.x2, None)])
 
         return m
 
@@ -161,7 +146,7 @@ def main():
     # Only estimate the parameter k[1]. The parameter k[2] will remain fixed
     # at its initial value
 
-    pest = parmest.Estimator(exp_list)
+    pest = parmest.Estimator(exp_list, obj_function="SSE")
     obj, theta = pest.theta_est()
     print(obj)
     print(theta)
@@ -174,9 +159,11 @@ def main():
 
     # =======================================================================
     # Estimate both k1 and k2 and compute the covariance matrix
-    pest = parmest.Estimator(exp_list)
-    n = 15  # total number of data points used in the objective (y in 15 scenarios)
-    obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=n)
+    pest = parmest.Estimator(exp_list, obj_function="SSE")
+    # Calculate the objective value and estimated parameters
+    obj, theta = pest.theta_est()
+    # Compute the covariance matrix using the reduced Hessian method
+    cov = pest.cov_est(method="reduced_hessian")
     print(obj)
     print(theta)
     print(cov)

pyomo/contrib/parmest/examples/reactor_design/multistart_example.py

@@ -0,0 +1,56 @@
+# ____________________________________________________________________________________
+#
+# 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.
+# ____________________________________________________________________________________
+
+from pyomo.common.dependencies import numpy as np, pandas as pd
+from itertools import product
+from os.path import join, abspath, dirname
+import pyomo.contrib.parmest.parmest as parmest
+from pyomo.contrib.parmest.examples.reactor_design.reactor_design import (
+    ReactorDesignExperiment,
+)
+
+
+def main():
+
+    # Read in data
+    file_dirname = dirname(abspath(str(__file__)))
+    file_name = abspath(join(file_dirname, "reactor_data.csv"))
+    data = pd.read_csv(file_name)
+
+    # Create an experiment list
+    exp_list = [ReactorDesignExperiment(data, i) for i in range(data.shape[0])]
+
+    # Solver options belong here (Ipopt options shown as example)
+    solver_options = {"max_iter": 1000, "tol": 1e-6}
+
+    pest = parmest.Estimator(
+        exp_list, obj_function="SSE", solver_options=solver_options
+    )
+
+    # Single-start estimation
+    obj, theta = pest.theta_est()
+    print("Single-start objective:", obj)
+    print("Single-start theta:\n", theta)
+
+    # Multistart estimation
+    results_df, best_theta, best_obj = pest.theta_est_multistart(
+        n_restarts=10,
+        multistart_sampling_method="uniform_random",
+        seed=42,
+        save_results=False,  # True if you want CSV via file_name=
+    )
+
+    print("\nMultistart best objective:", best_obj)
+    print("Multistart best theta:", best_theta)
+    print("\nAll multistart results:")
+    print(results_df)
+
+
+if __name__ == "__main__":
+    main()

pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py

@@ -33,11 +33,17 @@ def main():
 
     pest = parmest.Estimator(exp_list, obj_function='SSE')
 
-    # Parameter estimation with covariance
-    obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=17)
-    print(obj)
+    # Parameter estimation
+    obj, theta = pest.theta_est()
+    print("Least squares objective value:", obj)
+    print("Estimated parameters (theta):\n")
     print(theta)
 
+    # Compute the covariance matrix at the estimated parameter
+    cov = pest.cov_est()
+    print("Covariance matrix:\n")
+    print(cov)
+
 
 if __name__ == "__main__":
     main()

pyomo/contrib/parmest/examples/reactor_design/reactor_design.py

@@ -22,15 +22,15 @@ def reactor_design_model():
     # Create the concrete model
     model = pyo.ConcreteModel()
 
-    # Rate constants
-    model.k1 = pyo.Param(
-        initialize=5.0 / 6.0, within=pyo.PositiveReals, mutable=True
+    # Rate constants, make unknown parameters variables
+    model.k1 = pyo.Var(
+        initialize=5.0 / 6.0, within=pyo.PositiveReals, bounds=(0.1, 10.0)
     )  # min^-1
-    model.k2 = pyo.Param(
-        initialize=5.0 / 3.0, within=pyo.PositiveReals, mutable=True
+    model.k2 = pyo.Var(
+        initialize=5.0 / 3.0, within=pyo.PositiveReals, bounds=(0.1, 10.0)
     )  # min^-1
-    model.k3 = pyo.Param(
-        initialize=1.0 / 6000.0, within=pyo.PositiveReals, mutable=True
+    model.k3 = pyo.Var(
+        initialize=1.0 / 6000.0, within=pyo.PositiveReals, bounds=(1e-5, 1e-3)
     )  # m^3/(gmol min)
 
     # Inlet concentration of A, gmol/m^3
@@ -119,6 +119,11 @@ def label_model(self):
             (k, pyo.ComponentUID(k)) for k in [m.k1, m.k2, m.k3]
         )
 
+        m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.measurement_error.update(
+            [(m.ca, None), (m.cb, None), (m.cc, None), (m.cd, None)]
+        )
+
         return m
 
     def get_labeled_model(self):

pyomo/contrib/parmest/examples/rooney_biegler/multistart_example.py

@@ -0,0 +1,62 @@
+# ____________________________________________________________________________________
+#
+# 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.
+# ____________________________________________________________________________________
+
+from pyomo.common.dependencies import numpy as np, pandas as pd
+from itertools import product
+import pyomo.contrib.parmest.parmest as parmest
+from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
+    RooneyBieglerExperiment,
+)
+
+
+def main():
+
+    # Data
+    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'],
+    )
+
+    # 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()
+
+    # Solver options belong here (Ipopt options shown as example)
+    solver_options = {"max_iter": 1000, "tol": 1e-6}
+
+    pest = parmest.Estimator(
+        exp_list, obj_function="SSE", solver_options=solver_options
+    )
+
+    # Single-start estimation
+    obj, theta = pest.theta_est()
+    print("Single-start objective:", obj)
+    print("Single-start theta:\n", theta)
+
+    # Multistart estimation
+    results_df, best_theta, best_obj = pest.theta_est_multistart(
+        n_restarts=10,
+        multistart_sampling_method="uniform_random",
+        seed=42,
+        save_results=False,  # True if you want CSV via file_name=
+    )
+
+    print("\nMultistart best objective:", best_obj)
+    print("Multistart best theta:", best_theta)
+    print("\nAll multistart results:")
+    print(results_df)
+
+
+if __name__ == "__main__":
+    main()

pyomo/contrib/parmest/examples/rooney_biegler/profile_likelihood_example.py

@@ -0,0 +1,74 @@
+# ____________________________________________________________________________________
+#
+# 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.
+# ____________________________________________________________________________________
+
+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 (
+    RooneyBieglerExperiment,
+)
+
+
+def main():
+    # Data
+    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"],
+    )
+
+    # Build experiment list
+    exp_list = [RooneyBieglerExperiment(data.loc[i, :]) for i in range(data.shape[0])]
+
+    # Create estimator
+    pest = parmest.Estimator(exp_list, obj_function="SSE")
+
+    # Compute profile likelihood for both unknown parameters.
+    # Use a small grid for quick terminal runs.
+    profile_results = pest.profile_likelihood(
+        profiled_theta=["asymptote", "rate_constant"],
+        n_grid=9,
+        solver="ef_ipopt",
+        warmstart="neighbor",
+        # Demonstrate baseline from multistart integration:
+        use_multistart_for_baseline=True,
+        baseline_multistart_kwargs={
+            "n_restarts": 5,
+            "multistart_sampling_method": "uniform_random",
+            "seed": 7,
+        },
+    )
+
+    # Display a compact summary table
+    profiles = profile_results["profiles"]
+    print("\nBaseline:")
+    print(profile_results["baseline"])
+    print("\nProfile results (first 12 rows):")
+    print(
+        profiles[
+            [
+                "profiled_theta",
+                "theta_value",
+                "obj",
+                "delta_obj",
+                "lr_stat",
+                "status",
+                "success",
+            ]
+        ].head(12)
+    )
+
+    # Plot profile curves to file for terminal/non-GUI usage
+    out_file = "rooney_biegler_profile_likelihood.png"
+    parmest.graphics.profile_likelihood_plot(
+        profile_results, alpha=0.95, filename=out_file
+    )
+    print(f"\nSaved plot to: {out_file}")
+
+
+if __name__ == "__main__":
+    main()

pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py

@@ -43,8 +43,8 @@ def rooney_biegler_model(data, theta=None):
     if theta is None:
         theta = {'asymptote': 15, 'rate_constant': 0.5}
 
-    model.asymptote = pyo.Var(initialize=theta['asymptote'])
-    model.rate_constant = pyo.Var(initialize=theta['rate_constant'])
+    model.asymptote = pyo.Var(initialize=theta['asymptote'], bounds=(0.1, 100))
+    model.rate_constant = pyo.Var(initialize=theta['rate_constant'], bounds=(0, 10))
 
     # Fix the unknown parameters
     model.asymptote.fix()