From f46b9d48cc62b39bba5d15ff8cea4ce9e0f18150 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Wed, 30 Jun 2021 12:52:11 -0700
Subject: [PATCH 01/16] Allow for indexed variables as theta in parmest model

---
 pyomo/contrib/parmest/parmest.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py
index 3fab4406a68..4f00d4cd0ab 100644
--- a/pyomo/contrib/parmest/parmest.py
+++ b/pyomo/contrib/parmest/parmest.py
@@ -45,6 +45,7 @@
 
 from pyomo.opt import SolverFactory
 from pyomo.environ import Block, ComponentUID
+from pyomo.environ import Block, ComponentUID
 
 import pyomo.contrib.parmest.mpi_utils as mpiu
 import pyomo.contrib.parmest.ipopt_solver_wrapper as ipopt_solver_wrapper
@@ -337,7 +338,11 @@ def _create_parmest_model(self, data):
                     # If the component that was found is not a variable,
                     # this will generate an exception (and the warning
                     # in the 'except')
-                    var_validate.fixed = False
+                    if var_validate.is_indexed():
+                        for v in var_validate.values():
+                            v.fixed = False
+                    else:
+                        var_validate.fixed = False
                     # We want to standardize on the CUID string
                     # representation
                     self.theta_names[i] = repr(var_cuid)

From f2ef5ff97dca3d4aca6e82008fd7e766aab7f0c7 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Wed, 30 Jun 2021 12:52:33 -0700
Subject: [PATCH 02/16] Updated RB example to use romodel

---
 .../reactor_design/parmest_RO_example.py      | 80 +++++++++++++++++++
 .../examples/reactor_design/reactor_design.py | 56 +++++++++++++
 2 files changed, 136 insertions(+)
 create mode 100644 pyomo/contrib/parmest/examples/reactor_design/parmest_RO_example.py

diff --git a/pyomo/contrib/parmest/examples/reactor_design/parmest_RO_example.py b/pyomo/contrib/parmest/examples/reactor_design/parmest_RO_example.py
new file mode 100644
index 00000000000..3a65e6af72a
--- /dev/null
+++ b/pyomo/contrib/parmest/examples/reactor_design/parmest_RO_example.py
@@ -0,0 +1,80 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright 2017 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 pandas as pd
+import pyomo.environ as pyo
+import pyomo.contrib.parmest.parmest as parmest
+from pyomo.contrib.parmest.examples.reactor_design.reactor_design import reactor_design_model, reactor_design_model_RO
+
+if __name__ == "__main__":
+    # Note, if this is not called from main, you get multiprocessing warnings (windows issue)
+    
+    # Data
+    data = pd.read_excel('reactor_data.xlsx') 
+    
+    ### Reactor design model with data from a single experiment
+    m = reactor_design_model(data.iloc[-1])
+    solver = pyo.SolverFactory('ipopt')
+    solver.solve(m)
+    
+    m.pprint()
+    
+    print('k1', m.k1.value)
+    print('k2', m.k2.value)
+    print('k3', m.k3.value)
+    print('ca', m.ca.value)
+    print('cb', m.cb.value)
+    print('cc', m.cc.value)
+    print('cd', m.cd.value)
+    
+    ### Estimate k1, k2, k3 using data from all experiments
+    # Vars to estimate
+    theta_names = ['k1', 'k2', 'k3']
+
+    # Sum of squared error function
+    def SSE(model, data): 
+        expr = (float(data['ca']) - model.ca)**2 + \
+               (float(data['cb']) - model.cb)**2 + \
+               (float(data['cc']) - model.cc)**2 + \
+               (float(data['cd']) - model.cd)**2
+        return expr
+    
+    pest = parmest.Estimator(reactor_design_model, data.loc[0:17,:], theta_names, SSE)
+    obj, theta, cov = pest.theta_est(calc_cov=True)
+    print(obj)
+    print(theta)
+    print(cov)
+    
+    parmest.graphics.pairwise_plot((theta, cov, 1000), theta_star=theta, alpha=0.8, 
+                               distributions=['MVN'])
+    
+    
+    ### Reactor design model with data from a single experiment using RO with theta and cov
+    m = reactor_design_model_RO(data.iloc[-1], theta, cov)
+    
+    solver = pyo.SolverFactory('romodel.cuts')
+    solver.solve(m)
+    
+    m.pprint()
+    
+    theta_ro = {'k1': m.k['k1'].value,
+                'k2': m.k['k2'].value,
+                'k3': m.k['k3'].value}
+    print(theta_ro)
+    print('ca', m.ca.value)
+    print('cb', m.cb.value)
+    print('cc', m.cc.value)
+    print('cd', m.cd.value)
+    print(data.iloc[-1])
+
+    parmest.graphics.pairwise_plot((theta_ro, cov, 1000), theta_star=theta_ro, alpha=0.8, 
+                               distributions=['MVN'])
+
+        
\ No newline at end of file
diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
index 4b243fc9c69..ba4609a138c 100644
--- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
+++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
@@ -11,8 +11,10 @@
 Continuously stirred tank reactor model, based on
 pyomo\examples\doc\pyomobook\nonlinear-ch\react_design\ReactorDesign.py
 """
+import numpy as np
 import pandas as pd
 from pyomo.environ import ConcreteModel, Var, PositiveReals, Objective, Constraint, maximize, SolverFactory
+from romodel import UncSet, UncParam
 
 def reactor_design_model(data):
     
@@ -60,6 +62,60 @@ def reactor_design_model(data):
     
     return model
 
+
+
+def reactor_design_model_RO(data, theta, cov):
+    
+    # Create the concrete model
+    model = ConcreteModel()
+        
+    # Rate constants defined as an uncertain parameter
+    theta_names = theta.keys()
+    model.E = UncSet()
+    model.k = UncParam(theta_names, uncset=model.E, nominal=theta) # {'k1': 5.0/6.0, 'k2': 5.0/3.0, 'k3': 1.0/6000.0}) #
+    # Ellipsoidal set defined using covariance informed by data/parmest
+    mu = list(theta.values())
+    cov = cov.values
+    A = np.linalg.inv(cov).tolist()
+    lhs = 0
+    for i, ki in enumerate(theta_names):
+        for j, kj in enumerate(theta_names):
+            lhs += (model.k[ki] - mu[i])*A[i][j]*(model.k[kj] - mu[j])
+    model.E.cons = Constraint(expr=lhs <= 1)
+        
+    # Inlet concentration of A, gmol/m^3
+    model.caf = Var(initialize = float(data['caf']), within=PositiveReals)
+    model.caf.fixed = True
+        
+    # Space velocity (flowrate/volume)
+    model.sv = Var(initialize = float(data['sv']), within=PositiveReals)
+    model.sv.fixed = True
+        
+    # Outlet concentration of each component
+    model.ca = Var(initialize = float(data['ca']), within=PositiveReals) 
+    model.cb = Var(initialize = float(data['cb']), within=PositiveReals) 
+    model.cc = Var(initialize = float(data['cc']), within=PositiveReals) 
+    model.cd = Var(initialize = float(data['cd']), within=PositiveReals)
+
+    # Objective
+    model.obj = Objective(expr = model.cb, sense=maximize)
+        
+    # Constraints
+    model.ca_bal = Constraint(expr = (0 == model.sv * model.caf \
+                     - model.sv * model.ca - model.k['k1'] * model.ca \
+                     -  2.0 * model.k['k3'] * model.ca ** 2.0))
+        
+    model.cb_bal = Constraint(expr=(0 == -model.sv * model.cb \
+                     + model.k['k1'] * model.ca - model.k['k2'] * model.cb))
+    
+    model.cc_bal = Constraint(expr=(0 == -model.sv * model.cc \
+                     + model.k['k2'] * model.cb))
+        
+    model.cd_bal = Constraint(expr=(0 == -model.sv * model.cd \
+                     + model.k['k3'] * model.ca ** 2.0))
+        
+    return model
+    
 if __name__ == "__main__":
     
     # For a range of sv values, return ca, cb, cc, and cd

From 0f3865ac5edcf77d0b42b05791182a105ae29e76 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Mon, 19 Jul 2021 15:49:38 -0700
Subject: [PATCH 03/16] Added indexed reactor design example

---
 .../examples/reactor_design/reactor_design.py | 50 ++++++++++++++++++-
 1 file changed, 48 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
index ba4609a138c..9f50f381308 100644
--- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
+++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
@@ -13,7 +13,7 @@
 """
 import numpy as np
 import pandas as pd
-from pyomo.environ import ConcreteModel, Var, PositiveReals, Objective, Constraint, maximize, SolverFactory
+from pyomo.environ import ConcreteModel, RangeSet, Var, PositiveReals, Objective, Constraint, maximize, SolverFactory
 from romodel import UncSet, UncParam
 
 def reactor_design_model(data):
@@ -63,6 +63,52 @@ def reactor_design_model(data):
     return model
 
 
+def reactor_design_model_indexed(data):
+    
+    # Create the concrete model
+    model = ConcreteModel()
+    
+    # Rate constants
+    model.rxn = RangeSet(1,3)
+    
+    model.k = Var(model.rxn, initialize = 5.0/6.0, within=PositiveReals) 
+    model.k[1].fixed = True
+    model.k[2].fixed = True
+    model.k[3].fixed = True
+    
+    # Inlet concentration of A, gmol/m^3
+    model.caf = Var(initialize = float(data['caf']), within=PositiveReals)
+    model.caf.fixed = True
+    
+	# Space velocity (flowrate/volume)
+    model.sv = Var(initialize = float(data['sv']), within=PositiveReals)
+    model.sv.fixed = True
+    
+    # Outlet concentration of each component
+    model.ca = Var(initialize = 5000.0, within=PositiveReals) 
+    model.cb = Var(initialize = 2000.0, within=PositiveReals) 
+    model.cc = Var(initialize = 2000.0, within=PositiveReals) 
+    model.cd = Var(initialize = 1000.0, within=PositiveReals)
+    
+    # Objective
+    model.obj = Objective(expr = model.cb, sense=maximize)
+    
+    # Constraints
+    model.ca_bal = Constraint(expr = (0 == model.sv * model.caf \
+                     - model.sv * model.ca - model.k[1] * model.ca \
+                     -  2.0 * model.k[3] * model.ca ** 2.0))
+    
+    model.cb_bal = Constraint(expr=(0 == -model.sv * model.cb \
+                     + model.k[1] * model.ca - model.k[2] * model.cb))
+    
+    model.cc_bal = Constraint(expr=(0 == -model.sv * model.cc \
+                     + model.k[2] * model.cb))
+    
+    model.cd_bal = Constraint(expr=(0 == -model.sv * model.cd \
+                     + model.k[3] * model.ca ** 2.0))
+    
+    return model
+
 
 def reactor_design_model_RO(data, theta, cov):
     
@@ -82,7 +128,7 @@ def reactor_design_model_RO(data, theta, cov):
         for j, kj in enumerate(theta_names):
             lhs += (model.k[ki] - mu[i])*A[i][j]*(model.k[kj] - mu[j])
     model.E.cons = Constraint(expr=lhs <= 1)
-        
+
     # Inlet concentration of A, gmol/m^3
     model.caf = Var(initialize = float(data['caf']), within=PositiveReals)
     model.caf.fixed = True

From da22286bda2bedc18134afbb3fd7e867282188f3 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Tue, 20 Jul 2021 09:01:17 -0700
Subject: [PATCH 04/16] removed romodel import

---
 .../examples/reactor_design/reactor_design.py | 53 -------------------
 1 file changed, 53 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
index 9f50f381308..b22cac944be 100644
--- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
+++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
@@ -14,7 +14,6 @@
 import numpy as np
 import pandas as pd
 from pyomo.environ import ConcreteModel, RangeSet, Var, PositiveReals, Objective, Constraint, maximize, SolverFactory
-from romodel import UncSet, UncParam
 
 def reactor_design_model(data):
     
@@ -109,58 +108,6 @@ def reactor_design_model_indexed(data):
     
     return model
 
-
-def reactor_design_model_RO(data, theta, cov):
-    
-    # Create the concrete model
-    model = ConcreteModel()
-        
-    # Rate constants defined as an uncertain parameter
-    theta_names = theta.keys()
-    model.E = UncSet()
-    model.k = UncParam(theta_names, uncset=model.E, nominal=theta) # {'k1': 5.0/6.0, 'k2': 5.0/3.0, 'k3': 1.0/6000.0}) #
-    # Ellipsoidal set defined using covariance informed by data/parmest
-    mu = list(theta.values())
-    cov = cov.values
-    A = np.linalg.inv(cov).tolist()
-    lhs = 0
-    for i, ki in enumerate(theta_names):
-        for j, kj in enumerate(theta_names):
-            lhs += (model.k[ki] - mu[i])*A[i][j]*(model.k[kj] - mu[j])
-    model.E.cons = Constraint(expr=lhs <= 1)
-
-    # Inlet concentration of A, gmol/m^3
-    model.caf = Var(initialize = float(data['caf']), within=PositiveReals)
-    model.caf.fixed = True
-        
-    # Space velocity (flowrate/volume)
-    model.sv = Var(initialize = float(data['sv']), within=PositiveReals)
-    model.sv.fixed = True
-        
-    # Outlet concentration of each component
-    model.ca = Var(initialize = float(data['ca']), within=PositiveReals) 
-    model.cb = Var(initialize = float(data['cb']), within=PositiveReals) 
-    model.cc = Var(initialize = float(data['cc']), within=PositiveReals) 
-    model.cd = Var(initialize = float(data['cd']), within=PositiveReals)
-
-    # Objective
-    model.obj = Objective(expr = model.cb, sense=maximize)
-        
-    # Constraints
-    model.ca_bal = Constraint(expr = (0 == model.sv * model.caf \
-                     - model.sv * model.ca - model.k['k1'] * model.ca \
-                     -  2.0 * model.k['k3'] * model.ca ** 2.0))
-        
-    model.cb_bal = Constraint(expr=(0 == -model.sv * model.cb \
-                     + model.k['k1'] * model.ca - model.k['k2'] * model.cb))
-    
-    model.cc_bal = Constraint(expr=(0 == -model.sv * model.cc \
-                     + model.k['k2'] * model.cb))
-        
-    model.cd_bal = Constraint(expr=(0 == -model.sv * model.cd \
-                     + model.k['k3'] * model.ca ** 2.0))
-        
-    return model
     
 if __name__ == "__main__":
     

From 3c0107d1cf1a369a0404d03de6f5cddb3ca19d93 Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Tue, 20 Jul 2021 12:17:19 -0400
Subject: [PATCH 05/16] Added example from Paul.

---
 .../simple_reaction_parmest_example.py        | 90 +++++++++++++++++++
 1 file changed, 90 insertions(+)
 create mode 100644 pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py

diff --git a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
new file mode 100644
index 00000000000..a07c48c4c9c
--- /dev/null
+++ b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
@@ -0,0 +1,90 @@
+import pandas as pd
+from pandas import DataFrame
+from os import path
+
+from pyomo.environ import (ConcreteModel, Param, Var, PositiveReals, Objective,
+                           Constraint, RangeSet, Expression, minimize, exp, value)
+
+#from idaes.core.util import get_default_solver
+import pyomo.contrib.parmest.parmest as parmest
+
+# =======================================================================
+data = [{'experiment': 1, 'x1': 0.1, 'x2': 100, 'y': 0.98},
+        {'experiment': 2, 'x1': 0.2, 'x2': 100, 'y': 0.983},
+        {'experiment': 3, 'x1': 0.3, 'x2': 100, 'y': 0.955},
+        {'experiment': 4, 'x1': 0.4, 'x2': 100, 'y': 0.979},
+        {'experiment': 5, 'x1': 0.5, 'x2': 100, 'y': 0.993},
+        {'experiment': 6, 'x1': 0.05, 'x2': 200, 'y': 0.626},
+        {'experiment': 7, 'x1': 0.1, 'x2': 200, 'y': 0.544},
+        {'experiment': 8, 'x1': 0.15, 'x2': 200, 'y': 0.455},
+        {'experiment': 9, 'x1': 0.2, 'x2': 200, 'y': 0.225},
+        {'experiment': 10, 'x1': 0.25, 'x2': 200, 'y': 0.167},
+        {'experiment': 11, 'x1': 0.02, 'x2': 300, 'y': 0.566},
+        {'experiment': 12, 'x1': 0.04, 'x2': 300, 'y': 0.317},
+        {'experiment': 13, 'x1': 0.06, 'x2': 300, 'y': 0.034},
+        {'experiment': 14, 'x1': 0.08, 'x2': 300, 'y': 0.016},
+        {'experiment': 15, 'x1': 0.1, 'x2': 300, 'y': 0.006}]
+
+# =======================================================================
+
+def simple_reaction_model(data):
+
+    # Create the concrete model
+    model = ConcreteModel()
+
+    model.x1 = Param(initialize=float(data['x1']))
+    model.x2 = Param(initialize=float(data['x2']))
+
+    # Rate constants
+    model.rxn = RangeSet(2)
+    initial_guess = {1: 750, 2: 1200}
+    model.k = Var(model.rxn, initialize=initial_guess, within=PositiveReals)
+
+    # reaction product
+    model.y = Expression(expr=exp(-model.k[1] *
+                                  model.x1 * exp(-model.k[2] / model.x2)))
+                                  
+    # fix the rate constants
+    model.k.fix()
+
+    # linked variables and constraints
+    #model.k1 = Var(initialize=750, within=PositiveReals)
+    #model.k2 = Var(initialize=1200, within=PositiveReals)
+    #model.eq_L1 = Constraint(expr = model.k1 == model.k[1])
+    #model.eq_L2 = Constraint(expr = model.k2 == model.k[2])
+    #===================================================================
+    # 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
+
+if __name__ == "__main__":
+
+    # =======================================================================
+    # Parameter estimation without covariance estimate
+    #solver = get_default_solver
+    theta_names = ['k[1]']
+    pest = parmest.Estimator(simple_reaction_model, data, theta_names)
+    obj, theta = pest.theta_est()
+    print(obj)
+    print(theta)
+    #=======================================================================
+    # Parameter estimation covariance estimate
+    
+    obj, theta, cov = pest.theta_est(calc_cov=True)
+    print(obj)
+    print(theta)
+    print(cov)
+    
+    #=======================================================================
\ No newline at end of file

From e8701a6a5e21b862893158587918bea55ba68147 Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Tue, 20 Jul 2021 12:23:25 -0400
Subject: [PATCH 06/16] Added comments to the example

---
 .../simple_reaction_parmest_example.py        | 26 ++++++++++++++-----
 1 file changed, 19 insertions(+), 7 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
index a07c48c4c9c..2be571cc685 100644
--- a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
+++ b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
@@ -8,6 +8,17 @@
 #from idaes.core.util import get_default_solver
 import pyomo.contrib.parmest.parmest as parmest
 
+''' Example from Y. Bard's "Nonlinear Parameter Estimation"
+
+This example shows:
+1. How to define the unknown (to be regressed parameters) with an index
+2. How to call parmest to only estimate some of the parameters (and fix the rest)
+
+Code provided by Paul Akula.
+
+'''
+
+
 # =======================================================================
 data = [{'experiment': 1, 'x1': 0.1, 'x2': 100, 'y': 0.98},
         {'experiment': 2, 'x1': 0.2, 'x2': 100, 'y': 0.983},
@@ -44,14 +55,10 @@ def simple_reaction_model(data):
     model.y = Expression(expr=exp(-model.k[1] *
                                   model.x1 * exp(-model.k[2] / model.x2)))
                                   
-    # fix the rate constants
+    # fix all of the regressed parameters
     model.k.fix()
 
-    # linked variables and constraints
-    #model.k1 = Var(initialize=750, within=PositiveReals)
-    #model.k2 = Var(initialize=1200, within=PositiveReals)
-    #model.eq_L1 = Constraint(expr = model.k1 == model.k[1])
-    #model.eq_L2 = Constraint(expr = model.k2 == model.k[2])
+
     #===================================================================
     # Stage-specific cost computations
     def ComputeFirstStageCost_rule(model):
@@ -74,13 +81,18 @@ def total_cost_rule(m):
     # =======================================================================
     # Parameter estimation without covariance estimate
     #solver = get_default_solver
+    
+    # Only estimate the parameter k[1]. The parameter k[2] will remain fixed
+    # at its initial value
     theta_names = ['k[1]']
+    
     pest = parmest.Estimator(simple_reaction_model, data, theta_names)
     obj, theta = pest.theta_est()
     print(obj)
     print(theta)
+    
     #=======================================================================
-    # Parameter estimation covariance estimate
+    # Estimate the covariance matrix
     
     obj, theta, cov = pest.theta_est(calc_cov=True)
     print(obj)

From 92207529748c37bf1cbd79d84403574befeb7c3e Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Tue, 20 Jul 2021 09:25:24 -0700
Subject: [PATCH 07/16] removed RO example

---
 .../reactor_design/parmest_RO_example.py      | 80 -------------------
 1 file changed, 80 deletions(-)
 delete mode 100644 pyomo/contrib/parmest/examples/reactor_design/parmest_RO_example.py

diff --git a/pyomo/contrib/parmest/examples/reactor_design/parmest_RO_example.py b/pyomo/contrib/parmest/examples/reactor_design/parmest_RO_example.py
deleted file mode 100644
index 3a65e6af72a..00000000000
--- a/pyomo/contrib/parmest/examples/reactor_design/parmest_RO_example.py
+++ /dev/null
@@ -1,80 +0,0 @@
-#  ___________________________________________________________________________
-#
-#  Pyomo: Python Optimization Modeling Objects
-#  Copyright 2017 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 pandas as pd
-import pyomo.environ as pyo
-import pyomo.contrib.parmest.parmest as parmest
-from pyomo.contrib.parmest.examples.reactor_design.reactor_design import reactor_design_model, reactor_design_model_RO
-
-if __name__ == "__main__":
-    # Note, if this is not called from main, you get multiprocessing warnings (windows issue)
-    
-    # Data
-    data = pd.read_excel('reactor_data.xlsx') 
-    
-    ### Reactor design model with data from a single experiment
-    m = reactor_design_model(data.iloc[-1])
-    solver = pyo.SolverFactory('ipopt')
-    solver.solve(m)
-    
-    m.pprint()
-    
-    print('k1', m.k1.value)
-    print('k2', m.k2.value)
-    print('k3', m.k3.value)
-    print('ca', m.ca.value)
-    print('cb', m.cb.value)
-    print('cc', m.cc.value)
-    print('cd', m.cd.value)
-    
-    ### Estimate k1, k2, k3 using data from all experiments
-    # Vars to estimate
-    theta_names = ['k1', 'k2', 'k3']
-
-    # Sum of squared error function
-    def SSE(model, data): 
-        expr = (float(data['ca']) - model.ca)**2 + \
-               (float(data['cb']) - model.cb)**2 + \
-               (float(data['cc']) - model.cc)**2 + \
-               (float(data['cd']) - model.cd)**2
-        return expr
-    
-    pest = parmest.Estimator(reactor_design_model, data.loc[0:17,:], theta_names, SSE)
-    obj, theta, cov = pest.theta_est(calc_cov=True)
-    print(obj)
-    print(theta)
-    print(cov)
-    
-    parmest.graphics.pairwise_plot((theta, cov, 1000), theta_star=theta, alpha=0.8, 
-                               distributions=['MVN'])
-    
-    
-    ### Reactor design model with data from a single experiment using RO with theta and cov
-    m = reactor_design_model_RO(data.iloc[-1], theta, cov)
-    
-    solver = pyo.SolverFactory('romodel.cuts')
-    solver.solve(m)
-    
-    m.pprint()
-    
-    theta_ro = {'k1': m.k['k1'].value,
-                'k2': m.k['k2'].value,
-                'k3': m.k['k3'].value}
-    print(theta_ro)
-    print('ca', m.ca.value)
-    print('cb', m.cb.value)
-    print('cc', m.cc.value)
-    print('cd', m.cd.value)
-    print(data.iloc[-1])
-
-    parmest.graphics.pairwise_plot((theta_ro, cov, 1000), theta_star=theta_ro, alpha=0.8, 
-                               distributions=['MVN'])
-
-        
\ No newline at end of file

From f4f1c2b79568290f662f05afbb560cabbd844a06 Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Tue, 20 Jul 2021 12:28:11 -0400
Subject: [PATCH 08/16] Added reference for example.

---
 .../reaction_kinetics/simple_reaction_parmest_example.py     | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
index 2be571cc685..831d0ac9e6f 100644
--- a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
+++ b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
@@ -8,7 +8,7 @@
 #from idaes.core.util import get_default_solver
 import pyomo.contrib.parmest.parmest as parmest
 
-''' Example from Y. Bard's "Nonlinear Parameter Estimation"
+''' Example from Y. Bard, "Nonlinear Parameter Estimation", (pg. 124)
 
 This example shows:
 1. How to define the unknown (to be regressed parameters) with an index
@@ -20,6 +20,9 @@
 
 
 # =======================================================================
+''' Data from Table 5.2 in  Y. Bard, "Nonlinear Parameter Estimation", (pg. 124)
+'''
+
 data = [{'experiment': 1, 'x1': 0.1, 'x2': 100, 'y': 0.98},
         {'experiment': 2, 'x1': 0.2, 'x2': 100, 'y': 0.983},
         {'experiment': 3, 'x1': 0.3, 'x2': 100, 'y': 0.955},

From e90e1eac0e1a33de7957e8f89708de85221cb79a Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Tue, 20 Jul 2021 09:36:55 -0700
Subject: [PATCH 09/16] doc update

---
 pyomo/contrib/parmest/parmest.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py
index 413e9208390..c4091ed0de9 100644
--- a/pyomo/contrib/parmest/parmest.py
+++ b/pyomo/contrib/parmest/parmest.py
@@ -651,7 +651,7 @@ def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None, calc_cov
         Parameters
         ----------
         solver: string, optional
-            "ef_ipopt" or "k_aug". Default is "ef_ipopt".
+            Currently only "ef_ipopt" is supported. Default is "ef_ipopt".
         return_values: list, optional
             List of Variable names used to return values from the model
         bootlist: list, optional
@@ -667,8 +667,6 @@ def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None, calc_cov
             A dictionary of all values for theta
         variable values: pd.DataFrame
             Variable values for each variable name in return_values (only for solver='ef_ipopt')
-        Hessian: dict
-            A dictionary of dictionaries for the Hessian (only for solver='k_aug')
         cov: pd.DataFrame
             Covariance matrix of the fitted parameters (only for solver='ef_ipopt')
         """

From 45ff401de4a94ae3fc64047eb38a04e981445ba9 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Fri, 27 Aug 2021 08:11:37 -0700
Subject: [PATCH 10/16] removed indexed model and added to the test

---
 .../examples/reactor_design/reactor_design.py | 48 -------------------
 1 file changed, 48 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
index b22cac944be..7d990cf35f7 100644
--- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
+++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
@@ -60,54 +60,6 @@ def reactor_design_model(data):
                      + model.k3 * model.ca ** 2.0))
     
     return model
-
-
-def reactor_design_model_indexed(data):
-    
-    # Create the concrete model
-    model = ConcreteModel()
-    
-    # Rate constants
-    model.rxn = RangeSet(1,3)
-    
-    model.k = Var(model.rxn, initialize = 5.0/6.0, within=PositiveReals) 
-    model.k[1].fixed = True
-    model.k[2].fixed = True
-    model.k[3].fixed = True
-    
-    # Inlet concentration of A, gmol/m^3
-    model.caf = Var(initialize = float(data['caf']), within=PositiveReals)
-    model.caf.fixed = True
-    
-	# Space velocity (flowrate/volume)
-    model.sv = Var(initialize = float(data['sv']), within=PositiveReals)
-    model.sv.fixed = True
-    
-    # Outlet concentration of each component
-    model.ca = Var(initialize = 5000.0, within=PositiveReals) 
-    model.cb = Var(initialize = 2000.0, within=PositiveReals) 
-    model.cc = Var(initialize = 2000.0, within=PositiveReals) 
-    model.cd = Var(initialize = 1000.0, within=PositiveReals)
-    
-    # Objective
-    model.obj = Objective(expr = model.cb, sense=maximize)
-    
-    # Constraints
-    model.ca_bal = Constraint(expr = (0 == model.sv * model.caf \
-                     - model.sv * model.ca - model.k[1] * model.ca \
-                     -  2.0 * model.k[3] * model.ca ** 2.0))
-    
-    model.cb_bal = Constraint(expr=(0 == -model.sv * model.cb \
-                     + model.k[1] * model.ca - model.k[2] * model.cb))
-    
-    model.cc_bal = Constraint(expr=(0 == -model.sv * model.cc \
-                     + model.k[2] * model.cb))
-    
-    model.cd_bal = Constraint(expr=(0 == -model.sv * model.cd \
-                     + model.k[3] * model.ca ** 2.0))
-    
-    return model
-
     
 if __name__ == "__main__":
     

From c56cd244935ee2be3a8a7580eb279d3513c64b22 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Fri, 27 Aug 2021 08:12:04 -0700
Subject: [PATCH 11/16] cleaned up tests

---
 pyomo/contrib/parmest/tests/test_parmest.py   | 79 ++++++++-----------
 .../parmest/tests/test_scenariocreator.py     |  4 +-
 2 files changed, 33 insertions(+), 50 deletions(-)

diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index bac451ee419..cfc31874986 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -41,13 +41,14 @@
 @unittest.skipIf(not parmest.parmest_available,
                  "Cannot test parmest: required dependencies are missing")
 @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
-class parmest_object_Tester_RB(unittest.TestCase):
+class TestRooneyBiegler(unittest.TestCase):
 
     def setUp(self):
         from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import rooney_biegler_model
 
+        # Note, the data used in this test has been corrected to use data.loc[5,'hour'] = 7 (instead of 6)
         data = pd.DataFrame(data=[[1,8.3],[2,10.3],[3,19.0],
-                                  [4,16.0],[5,15.6],[6,19.8]], columns=['hour', 'y'])
+                                  [4,16.0],[5,15.6],[7,19.8]], columns=['hour', 'y'])
 
         theta_names = ['asymptote', 'rate_constant']
 
@@ -65,9 +66,9 @@ def SSE(model, data):
     def test_theta_est(self):
         objval, thetavals = self.pest.theta_est()
 
-        self.assertAlmostEqual(objval, 4.4675, places=2)
-        self.assertAlmostEqual(thetavals['asymptote'], 19.2189, places=2) # 19.1426 from the paper
-        self.assertAlmostEqual(thetavals['rate_constant'], 0.5312, places=2) # 0.5311 from the paper
+        self.assertAlmostEqual(objval, 4.3317112, places=2)
+        self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2) 
+        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2) 
 
     @unittest.skipIf(not graphics.imports_available,
                      "parmest.graphics imports are unavailable")
@@ -109,8 +110,8 @@ def test_likelihood_ratio(self):
         LR = self.pest.likelihood_ratio_test(obj_at_theta, objval, [0.8, 0.9, 1.0])
 
         self.assertTrue(set(LR.columns) >= set([0.8, 0.9, 1.0]))
-        self.assertTrue(LR[0.8].sum() == 7)
-        self.assertTrue(LR[0.9].sum() == 11)
+        self.assertTrue(LR[0.8].sum() == 6)
+        self.assertTrue(LR[0.9].sum() == 10)
         self.assertTrue(LR[1.0].sum() == 60) # all true
 
         graphics.pairwise_plot(LR, thetavals, 0.8)
@@ -185,37 +186,6 @@ def test_theta_k_aug_for_Hessian(self):
         objval, thetavals, Hessian = self.pest.theta_est(solver="k_aug")
         self.assertAlmostEqual(objval, 4.4675, places=2)
 
-
-'''
-The test cases above were developed with a transcription mistake in the dataset.
-This test works with the correct dataset.
-'''
-@unittest.skipIf(not parmest.parmest_available,
-                 "Cannot test parmest: required dependencies are missing")
-@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
-class parmest_object_Tester_RB_match_paper(unittest.TestCase):
-
-    def setUp(self):
-        from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import rooney_biegler_model
-
-        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'])
-
-        theta_names = ['asymptote', 'rate_constant']
-
-        def SSE(model, data):
-            expr = sum((data.y[i] - model.response_function[data.hour[i]])**2 for i in data.index)
-            return expr
-
-        self.pest = parmest.Estimator(rooney_biegler_model, data, theta_names, SSE)
-
-    def test_theta_est(self):
-        objval, thetavals = self.pest.theta_est(calc_cov=False)
-
-        self.assertAlmostEqual(objval, 4.3317112, places=2)
-        self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2) # 19.1426 from the paper
-        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2) # 0.5311 from the paper
-
     @unittest.skipIf(not pynumero_ASL_available, "pynumero ASL is not available")
     @unittest.skipIf(not parmest.inverse_reduced_hessian_available,
                      "Cannot test covariance matrix: required ASL dependency is missing")
@@ -223,8 +193,8 @@ def test_theta_est_cov(self):
         objval, thetavals, cov = self.pest.theta_est(calc_cov=True)
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
-        self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2) # 19.1426 from the paper
-        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2) # 0.5311 from the paper
+        self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2) 
+        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2)
 
         # Covariance matrix
         self.assertAlmostEqual(cov.iloc[0,0], 6.30579403, places=2) # 6.22864 from paper
@@ -244,7 +214,7 @@ def test_theta_est_cov(self):
 @unittest.skipIf(not parmest.parmest_available,
                  "Cannot test parmest: required dependencies are missing")
 @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
-class Test_parmest_indexed_variables(unittest.TestCase):
+class TestIndexedVariables(unittest.TestCase):
 
     def make_model(self, theta_names):
 
@@ -262,7 +232,9 @@ def rooney_biegler_model_alternate(data):
             model.var_names = pyo.Set(initialize=['asymptote','rate_constant'])
 
             model.theta = pyo.Var(model.var_names, initialize={'asymptote':15, 'rate_constant':0.5})
-
+            model.theta['asymptote'].fixed = True # parmest will unfix theta variables, even when they are indexed
+            model.theta['rate_constant'].fixed = True
+            
             def response_rule(m, h):
                 expr = m.theta['asymptote'] * (1 - pyo.exp(-m.theta['rate_constant'] * h))
                 return expr
@@ -279,7 +251,18 @@ def SSE(model, data):
             return expr
 
         return parmest.Estimator(rooney_biegler_model_alternate, data, theta_names, SSE)
+    
+    def test_theta_est(self):
 
+        theta_names = ["theta"]
+
+        pest = self.make_model(theta_names)
+        objval, thetavals = pest.theta_est(calc_cov=False)
+
+        self.assertAlmostEqual(objval, 4.3317112, places=2)
+        self.assertAlmostEqual(thetavals["theta[asymptote]"], 19.1426, places=2) 
+        self.assertAlmostEqual(thetavals["theta[rate_constant]"], 0.5311, places=2) 
+        
     def test_theta_est_quotedIndex(self):
 
         theta_names = ["theta['asymptote']", "theta['rate_constant']"]
@@ -288,8 +271,8 @@ def test_theta_est_quotedIndex(self):
         objval, thetavals = pest.theta_est(calc_cov=False)
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
-        self.assertAlmostEqual(thetavals["theta[asymptote]"], 19.1426, places=2) # 19.1426 from the paper
-        self.assertAlmostEqual(thetavals['theta[rate_constant]'], 0.5311, places=2) # 0.5311 from the paper
+        self.assertAlmostEqual(thetavals["theta[asymptote]"], 19.1426, places=2) 
+        self.assertAlmostEqual(thetavals["theta[rate_constant]"], 0.5311, places=2) 
 
     def test_theta_est_impliedStrIndex(self):
 
@@ -300,7 +283,7 @@ def test_theta_est_impliedStrIndex(self):
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
         self.assertAlmostEqual(thetavals["theta[asymptote]"], 19.1426, places=2) # 19.1426 from the paper
-        self.assertAlmostEqual(thetavals['theta[rate_constant]'], 0.5311, places=2) # 0.5311 from the paper
+        self.assertAlmostEqual(thetavals["theta[rate_constant]"], 0.5311, places=2) # 0.5311 from the paper
 
 
     @unittest.skipIf(not pynumero_ASL_available, "pynumero ASL is not available")
@@ -308,7 +291,7 @@ def test_theta_est_impliedStrIndex(self):
         not parmest.inverse_reduced_hessian_available,
         "Cannot test covariance matrix: required ASL dependency is missing")
     def test_theta_est_cov(self):
-        theta_names = ["theta[asymptote]", "theta[rate_constant]"]
+        theta_names = ["theta"]
 
         pest = self.make_model(theta_names)
         objval, thetavals, cov = pest.theta_est(calc_cov=True)
@@ -329,7 +312,7 @@ def test_theta_est_cov(self):
                  "Cannot test parmest: required dependencies are missing")
 @unittest.skipIf(not ipopt_available,
                  "The 'ipopt' solver is not available")
-class parmest_object_Tester_reactor_design(unittest.TestCase):
+class TestReactorDesign(unittest.TestCase):
 
     def setUp(self):
         from pyomo.contrib.parmest.examples.reactor_design.reactor_design import reactor_design_model
@@ -389,7 +372,7 @@ def test_return_values(self):
 @unittest.skipIf(not graphics.imports_available,
                  "parmest.graphics imports are unavailable")
 @unittest.skipIf(is_osx, "Disabling graphics tests on OSX due to issue in Matplotlib, see Pyomo PR #1337")
-class parmest_graphics(unittest.TestCase):
+class TestGraphics(unittest.TestCase):
 
     def setUp(self):
         self.A = pd.DataFrame(np.random.randint(0,100,size=(100,4)), columns=list('ABCD'))
diff --git a/pyomo/contrib/parmest/tests/test_scenariocreator.py b/pyomo/contrib/parmest/tests/test_scenariocreator.py
index ee58ebb8f4b..92794d334f4 100644
--- a/pyomo/contrib/parmest/tests/test_scenariocreator.py
+++ b/pyomo/contrib/parmest/tests/test_scenariocreator.py
@@ -31,7 +31,7 @@
 
 @unittest.skipIf(not parmest.parmest_available, "Cannot test parmest: required dependencies are missing")
 @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
-class pamest_Scenario_creator_reactor_design(unittest.TestCase):
+class TestScenarioReactorDesign(unittest.TestCase):
     
     def setUp(self):
         from pyomo.contrib.parmest.examples.reactor_design.reactor_design import reactor_design_model
@@ -99,7 +99,7 @@ def test_no_csv_if_empty(self):
 
 @unittest.skipIf(not parmest.parmest_available, "Cannot test parmest: required dependencies are missing")
 @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
-class  pamest_Scenario_creator_semibatch(unittest.TestCase):
+class TestScenarioSemibatch(unittest.TestCase):
     
     def setUp(self):
         import pyomo.contrib.parmest.examples.semibatch.semibatch as sb

From 0144b45efb198ac20a13fbf35a0f5aa24bd8a978 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Fri, 27 Aug 2021 09:11:33 -0700
Subject: [PATCH 12/16] cleanup

---
 .../simple_reaction_parmest_example.py        | 41 ++++++++++---------
 .../examples/reactor_design/reactor_design.py |  5 +--
 pyomo/contrib/parmest/parmest.py              |  1 -
 pyomo/contrib/parmest/tests/test_parmest.py   |  8 ++--
 4 files changed, 28 insertions(+), 27 deletions(-)

diff --git a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
index 831d0ac9e6f..c17cd7fc22e 100644
--- a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
+++ b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py
@@ -1,3 +1,21 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright 2017 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.
+#  ___________________________________________________________________________
+''' 
+Example from Y. Bard, "Nonlinear Parameter Estimation", (pg. 124)
+
+This example shows:
+1. How to define the unknown (to be regressed parameters) with an index
+2. How to call parmest to only estimate some of the parameters (and fix the rest)
+
+Code provided by Paul Akula.
+'''
 import pandas as pd
 from pandas import DataFrame
 from os import path
@@ -8,17 +26,6 @@
 #from idaes.core.util import get_default_solver
 import pyomo.contrib.parmest.parmest as parmest
 
-''' Example from Y. Bard, "Nonlinear Parameter Estimation", (pg. 124)
-
-This example shows:
-1. How to define the unknown (to be regressed parameters) with an index
-2. How to call parmest to only estimate some of the parameters (and fix the rest)
-
-Code provided by Paul Akula.
-
-'''
-
-
 # =======================================================================
 ''' Data from Table 5.2 in  Y. Bard, "Nonlinear Parameter Estimation", (pg. 124)
 '''
@@ -83,23 +90,19 @@ def total_cost_rule(m):
 
     # =======================================================================
     # Parameter estimation without covariance estimate
-    #solver = get_default_solver
-    
     # Only estimate the parameter k[1]. The parameter k[2] will remain fixed
     # at its initial value
     theta_names = ['k[1]']
-    
     pest = parmest.Estimator(simple_reaction_model, data, theta_names)
     obj, theta = pest.theta_est()
     print(obj)
     print(theta)
-    
+    print()
     #=======================================================================
-    # Estimate the covariance matrix
-    
+    # Estimate both k1 and k2 and compute the covariance matrix
+    theta_names = ['k']
+    pest = parmest.Estimator(simple_reaction_model, data, theta_names)
     obj, theta, cov = pest.theta_est(calc_cov=True)
     print(obj)
     print(theta)
     print(cov)
-    
-    #=======================================================================
\ No newline at end of file
diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
index 7d990cf35f7..4b243fc9c69 100644
--- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
+++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py
@@ -11,9 +11,8 @@
 Continuously stirred tank reactor model, based on
 pyomo\examples\doc\pyomobook\nonlinear-ch\react_design\ReactorDesign.py
 """
-import numpy as np
 import pandas as pd
-from pyomo.environ import ConcreteModel, RangeSet, Var, PositiveReals, Objective, Constraint, maximize, SolverFactory
+from pyomo.environ import ConcreteModel, Var, PositiveReals, Objective, Constraint, maximize, SolverFactory
 
 def reactor_design_model(data):
     
@@ -60,7 +59,7 @@ def reactor_design_model(data):
                      + model.k3 * model.ca ** 2.0))
     
     return model
-    
+
 if __name__ == "__main__":
     
     # For a range of sv values, return ca, cb, cc, and cd
diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py
index c4091ed0de9..f958c02cc43 100644
--- a/pyomo/contrib/parmest/parmest.py
+++ b/pyomo/contrib/parmest/parmest.py
@@ -45,7 +45,6 @@
 
 from pyomo.opt import SolverFactory
 from pyomo.environ import Block, ComponentUID
-from pyomo.environ import Block, ComponentUID
 
 import pyomo.contrib.parmest.mpi_utils as mpiu
 import pyomo.contrib.parmest.ipopt_solver_wrapper as ipopt_solver_wrapper
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index cfc31874986..dfbd48e7e75 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -67,8 +67,8 @@ def test_theta_est(self):
         objval, thetavals = self.pest.theta_est()
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
-        self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2) 
-        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2) 
+        self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2)  # 19.1426 from the paper
+        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2)  # 19.1426 from the paper
 
     @unittest.skipIf(not graphics.imports_available,
                      "parmest.graphics imports are unavailable")
@@ -193,8 +193,8 @@ def test_theta_est_cov(self):
         objval, thetavals, cov = self.pest.theta_est(calc_cov=True)
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
-        self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2) 
-        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2)
+        self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2) # 19.1426 from the paper
+        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2) # 19.1426 from the paper
 
         # Covariance matrix
         self.assertAlmostEqual(cov.iloc[0,0], 6.30579403, places=2) # 6.22864 from paper

From c56b5c4fcb606b9c806c7dc99c9d71352027e5a2 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Fri, 27 Aug 2021 09:31:23 -0700
Subject: [PATCH 13/16] added init to example, cleanup

---
 .../parmest/examples/reaction_kinetics/__init__.py    |  9 +++++++++
 .../parmest/examples/reactor_design/__init__.py       | 11 ++++++++++-
 .../parmest/examples/rooney_biegler/__init__.py       |  2 +-
 pyomo/contrib/parmest/examples/semibatch/__init__.py  |  2 +-
 pyomo/contrib/parmest/tests/test_parmest.py           |  4 ++--
 5 files changed, 23 insertions(+), 5 deletions(-)
 create mode 100644 pyomo/contrib/parmest/examples/reaction_kinetics/__init__.py

diff --git a/pyomo/contrib/parmest/examples/reaction_kinetics/__init__.py b/pyomo/contrib/parmest/examples/reaction_kinetics/__init__.py
new file mode 100644
index 00000000000..cd6b0b75748
--- /dev/null
+++ b/pyomo/contrib/parmest/examples/reaction_kinetics/__init__.py
@@ -0,0 +1,9 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright 2017 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.
+#  ___________________________________________________________________________
diff --git a/pyomo/contrib/parmest/examples/reactor_design/__init__.py b/pyomo/contrib/parmest/examples/reactor_design/__init__.py
index 8d1c8b69c3f..6b39dd18d6a 100644
--- a/pyomo/contrib/parmest/examples/reactor_design/__init__.py
+++ b/pyomo/contrib/parmest/examples/reactor_design/__init__.py
@@ -1 +1,10 @@
- 
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright 2017 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.
+#  ___________________________________________________________________________
+
diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/__init__.py b/pyomo/contrib/parmest/examples/rooney_biegler/__init__.py
index d9f70706c29..cd6b0b75748 100644
--- a/pyomo/contrib/parmest/examples/rooney_biegler/__init__.py
+++ b/pyomo/contrib/parmest/examples/rooney_biegler/__init__.py
@@ -1,4 +1,4 @@
- #  ___________________________________________________________________________
+#  ___________________________________________________________________________
 #
 #  Pyomo: Python Optimization Modeling Objects
 #  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
diff --git a/pyomo/contrib/parmest/examples/semibatch/__init__.py b/pyomo/contrib/parmest/examples/semibatch/__init__.py
index 5296dafcc78..6b39dd18d6a 100644
--- a/pyomo/contrib/parmest/examples/semibatch/__init__.py
+++ b/pyomo/contrib/parmest/examples/semibatch/__init__.py
@@ -1,4 +1,4 @@
- #  ___________________________________________________________________________
+#  ___________________________________________________________________________
 #
 #  Pyomo: Python Optimization Modeling Objects
 #  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index dfbd48e7e75..c424c9303d2 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -68,7 +68,7 @@ def test_theta_est(self):
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
         self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2)  # 19.1426 from the paper
-        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2)  # 19.1426 from the paper
+        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2)  # 0.5311 from the paper
 
     @unittest.skipIf(not graphics.imports_available,
                      "parmest.graphics imports are unavailable")
@@ -194,7 +194,7 @@ def test_theta_est_cov(self):
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
         self.assertAlmostEqual(thetavals['asymptote'], 19.1426, places=2) # 19.1426 from the paper
-        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2) # 19.1426 from the paper
+        self.assertAlmostEqual(thetavals['rate_constant'], 0.5311, places=2) # 0.5311 from the paper
 
         # Covariance matrix
         self.assertAlmostEqual(cov.iloc[0,0], 6.30579403, places=2) # 6.22864 from paper

From 9858238ed5a879ad187abba03dd8d6447cee0650 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Fri, 27 Aug 2021 09:33:08 -0700
Subject: [PATCH 14/16] Changed use of dictionary for IO to pd.Series

---
 pyomo/contrib/parmest/parmest.py | 35 +++++++++++++++++---------------
 1 file changed, 19 insertions(+), 16 deletions(-)

diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py
index f958c02cc43..320ee130aea 100644
--- a/pyomo/contrib/parmest/parmest.py
+++ b/pyomo/contrib/parmest/parmest.py
@@ -273,7 +273,7 @@ class Estimator(object):
     model_function: function
         Function that generates an instance of the Pyomo model using 'data' 
         as the input argument
-    data: pandas DataFrame, list of dictionaries, or list of json file names
+    data: pd.DataFrame, list of dictionaries, or list of json file names
         Data that is used to build an instance of the Pyomo model and build 
         the objective function
     theta_names: list of strings
@@ -506,6 +506,8 @@ def _Q_opt(self, ThetaVals=None, solver="ef_ipopt",
                 cov = 2 * sse / (n - l) * inv_red_hes
                 cov = pd.DataFrame(cov, index=thetavals.keys(), columns=thetavals.keys())
             
+            thetavals = pd.Series(thetavals)
+            
             if len(return_values) > 0:
                 var_values = []
                 for exp_i in self.ef_instance.component_objects(Block, descend_into=False):
@@ -662,8 +664,8 @@ def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None, calc_cov
         -------
         objectiveval: float
             The objective function value
-        thetavals: dict
-            A dictionary of all values for theta
+        thetavals: pd.Series
+            Estimated values for theta
         variable values: pd.DataFrame
             Variable values for each variable name in return_values (only for solver='ef_ipopt')
         cov: pd.DataFrame
@@ -698,7 +700,7 @@ def theta_est_bootstrap(self, bootstrap_samples, samplesize=None,
         
         Returns
         -------
-        bootstrap_theta: DataFrame 
+        bootstrap_theta: pd.DataFrame 
             Theta values for each sample and (if return_samples = True) 
             the sample numbers used in each estimation
         """
@@ -760,7 +762,7 @@ def theta_est_leaveNout(self, lNo, lNo_samples=None, seed=None,
         
         Returns
         -------
-        lNo_theta: DataFrame 
+        lNo_theta: pd.DataFrame 
             Theta values for each sample and (if return_samples = True) 
             the sample numbers left out of each estimation
         """
@@ -889,12 +891,12 @@ def objective_at_theta(self, theta_values):
 
         Parameters
         ----------
-        theta_values: DataFrame, columns=theta_names
+        theta_values: pd.DataFrame, columns=theta_names
             Values of theta used to compute the objective
             
         Returns
         -------
-        obj_at_theta: DataFrame
+        obj_at_theta: pd.DataFrame
             Objective value for each theta (infeasible solutions are 
             omitted).
         """
@@ -929,7 +931,7 @@ def likelihood_ratio_test(self, obj_at_theta, obj_value, alphas,
         
         Parameters
         ----------
-        obj_at_theta: DataFrame, columns = theta_names + 'obj'
+        obj_at_theta: pd.DataFrame, columns = theta_names + 'obj'
             Objective values for each theta value (returned by 
             objective_at_theta)
         obj_value: int or float
@@ -941,10 +943,10 @@ def likelihood_ratio_test(self, obj_at_theta, obj_value, alphas,
             
         Returns
         -------
-        LR: DataFrame 
+        LR: pd.DataFrame 
             Objective values for each theta value along with True or False for 
             each alpha
-        thresholds: dictionary
+        thresholds: pd.Series
             If return_threshold = True, the thresholds are also returned.
         """
         assert isinstance(obj_at_theta, pd.DataFrame)
@@ -960,6 +962,8 @@ def likelihood_ratio_test(self, obj_at_theta, obj_value, alphas,
             thresholds[a] = obj_value * ((chi2_val / (S - 2)) + 1)
             LR[a] = LR['obj'] < thresholds[a]
         
+        thresholds = pd.Series(thresholds)
+        
         if return_thresholds:
             return LR, thresholds
         else:
@@ -974,7 +978,7 @@ def confidence_region_test(self, theta_values, distribution, alphas,
         
         Parameters
         ----------
-        theta_values: DataFrame, columns = theta_names
+        theta_values: pd.DataFrame, columns = theta_names
             Theta values used to generate a confidence region 
             (generally returned by theta_est_bootstrap)
         distribution: string
@@ -984,23 +988,22 @@ def confidence_region_test(self, theta_values, distribution, alphas,
         alphas: list
             List of alpha values used to determine if theta values are inside 
             or outside the region.
-        test_theta_values: dictionary or DataFrame, keys/columns = theta_names, optional
+        test_theta_values: pd.Series or pd.DataFrame, keys/columns = theta_names, optional
             Additional theta values that are compared to the confidence region
             to determine if they are inside or outside.
         
         Returns
-        -------
-        training_results: DataFrame 
+        training_results: pd.DataFrame 
             Theta value used to generate the confidence region along with True 
             (inside) or False (outside) for each alpha
-        test_results: DataFrame 
+        test_results: pd.DataFrame 
             If test_theta_values is not None, returns test theta value along 
             with True (inside) or False (outside) for each alpha
         """
         assert isinstance(theta_values, pd.DataFrame)
         assert distribution in ['Rect', 'MVN', 'KDE']
         assert isinstance(alphas, list)
-        assert isinstance(test_theta_values, (type(None), dict, pd.DataFrame))
+        assert isinstance(test_theta_values, (type(None), dict, pd.Series, pd.DataFrame))
         
         if isinstance(test_theta_values, dict):
             test_theta_values = pd.Series(test_theta_values).to_frame().transpose()

From c59c9ceb5c4e84bc12c152fa85ad9fc06958736e Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Fri, 27 Aug 2021 09:55:00 -0700
Subject: [PATCH 15/16] test_theta_values must be a dataframe

---
 pyomo/contrib/parmest/parmest.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py
index 320ee130aea..ae694af141d 100644
--- a/pyomo/contrib/parmest/parmest.py
+++ b/pyomo/contrib/parmest/parmest.py
@@ -1005,7 +1005,7 @@ def confidence_region_test(self, theta_values, distribution, alphas,
         assert isinstance(alphas, list)
         assert isinstance(test_theta_values, (type(None), dict, pd.Series, pd.DataFrame))
         
-        if isinstance(test_theta_values, dict):
+        if isinstance(test_theta_values, (dict, pd.Series)):
             test_theta_values = pd.Series(test_theta_values).to_frame().transpose()
             
         training_results = theta_values.copy()

From 173ebf080ed1724f8378f7429cf578f54aa9b6d6 Mon Sep 17 00:00:00 2001
From: Katherine Klise <kaklise@sandia.gov>
Date: Tue, 31 Aug 2021 14:34:11 -0700
Subject: [PATCH 16/16] Updated to use unfix()

---
 pyomo/contrib/parmest/parmest.py | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py
index ae694af141d..127b7faa39d 100644
--- a/pyomo/contrib/parmest/parmest.py
+++ b/pyomo/contrib/parmest/parmest.py
@@ -337,11 +337,7 @@ def _create_parmest_model(self, data):
                     # If the component that was found is not a variable,
                     # this will generate an exception (and the warning
                     # in the 'except')
-                    if var_validate.is_indexed():
-                        for v in var_validate.values():
-                            v.fixed = False
-                    else:
-                        var_validate.fixed = False
+                    var_validate.unfix()
                     # We want to standardize on the CUID string
                     # representation
                     self.theta_names[i] = repr(var_cuid)