Parmest utils submodule, added support for Params and Vars

doc/OnlineDocs/bibliography.rst

@@ -15,6 +15,11 @@ Bibliography
    engineering. AIChE J. 2021; 67:e17175. DOI `10.1002/aic.17175
    <https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.17175>`_
 
+.. [mpisppy] Bernard Knueven, David Mildebrath, Christopher Muir,
+   John D Siirola, Jean-Paul Watson, and David L Woodruff, A Parallel 
+   Hub-and-Spoke System for Large-Scale Scenario-Based Optimization 
+   Under Uncertainty, pre-print, 2020
+
 .. [ParmestPaper] Katherine A. Klise, Bethany L. Nicholson, Andrea
    Staid, David L.Woodruff.  Parmest: Parameter Estimation Via Pyomo.
    Computer Aided Chemical Engineering, 47 (2019): 41-46.

doc/OnlineDocs/contributed_packages/parmest/driver.rst

@@ -5,10 +5,11 @@ Parameter Estimation
 
 Parameter Estimation using parmest requires a Pyomo model, experimental
 data which defines multiple scenarios, and a list of parameter names
-(thetas) to estimate.  parmest uses PySP [PyomoBookII]_ to solve a
+(thetas) to estimate.  parmest uses Pyomo [PyomoBookII]_ and (optionally) 
+mpi-sppy [mpisppy]_ to solve a
 two-stage stochastic programming problem, where the experimental data is
 used to create a scenario tree.  The objective function needs to be
-written in PySP form with the Pyomo Expression for first stage cost
+written with the Pyomo Expression for first stage cost
 (named "FirstStageCost") set to zero and the Pyomo Expression for second
 stage cost (named "SecondStageCost") defined as the deviation between
 the model and the observations (typically defined as the sum of squared
@@ -17,8 +18,8 @@ deviation between model values and observed values).
 If the Pyomo model is not formatted as a two-stage stochastic
 programming problem in this format, the user can supply a custom
 function to use as the second stage cost and the Pyomo model will be
-modified within parmest to match the specifications required by PySP.
-The PySP callback function is also defined within parmest.  The callback
+modified within parmest to match the required specifications.
+The stochastic programming callback function is also defined within parmest.  The callback
 function returns a populated and initialized model for each scenario.
 
 To use parmest, the user creates a :class:`~pyomo.contrib.parmest.parmest.Estimator` object 
@@ -88,11 +89,16 @@ Model function
 
 The first argument is a function which uses data for a single scenario
 to return a populated and initialized Pyomo model for that scenario.
-Parameters that the user would like to estimate must be defined as
-variables (Pyomo `Var`).  The variables can be fixed (parmest unfixes
-variables that will be estimated).  The model does not have to be
-specifically written for parmest. That is, parmest can modify the
-objective for PySP, see :ref:`ObjFunction` below.
+
+Parameters that the user would like to estimate can be defined as
+**mutable parameters (Pyomo `Param`) or variables (Pyomo `Var`)**.  
+Within parmest, any parameters that are to be estimated are converted to unfixed variables. 
+Variables that are to be estimated are also unfixed.
+
+The model does not have to be specifically written as a 
+two-stage stochastic programming problem for parmest. 
+That is, parmest can modify the
+objective, see :ref:`ObjFunction` below.
 
 Data
 ----
@@ -104,6 +110,8 @@ model.  Supported data formats include:
   names refer to observed quantities.  Pandas DataFrames are easily
   stored and read in from csv, excel, or databases, or created directly
   in Python.
+* **List of Pandas Dataframe** where each entry in the list is a separate scenario. 
+  Dataframes store observed quantities, referenced by index and column.
 * **List of dictionaries** where each entry in the list is a separate
   scenario and the keys (or nested keys) refer to observed quantities.
   Dictionaries are often preferred over DataFrames when using static and
@@ -122,8 +130,8 @@ functions, see :ref:`ObjFunction` below.
 Theta names
 -----------
 
-The third argument is a list of variable names that the user wants to
-estimate.  The list contains strings with `Var` names from the Pyomo
+The third argument is a list of parameters or variable names that the user wants to
+estimate.  The list contains strings with `Param` and/or `Var` names from the Pyomo
 model.
 
 .. _ObjFunction:
@@ -132,11 +140,14 @@ Objective function
 ------------------
 
 The fourth argument is an optional argument which defines the
-optimization objective function to use in parameter estimation.  If no
-objective function is specified, the Pyomo model is used "as is" and
+optimization objective function to use in parameter estimation.  
+
+If no objective function is specified, the Pyomo model is used "as is" and
 should be defined with "FirstStageCost" and "SecondStageCost"
-expressions that are used to build an objective for PySP.  If the Pyomo
-model is not written as a two stage stochastic programming problem in
+expressions that are used to build an objective for the two-stage 
+stochastic programming problem.  
+
+If the Pyomo model is not written as a two-stage stochastic programming problem in
 this format, and/or if the user wants to use an objective that is
 different than the original model, a custom objective function can be
 defined for parameter estimation.  The objective function arguments

doc/OnlineDocs/contributed_packages/parmest/index.rst

@@ -3,7 +3,7 @@ Parameter Estimation with ``parmest``
 
 ``parmest`` is a Python package built on the Pyomo optimization modeling
 language ([PyomoJournal]_, [PyomoBookII]_) to support parameter estimation using experimental data along with
-confidence regions and subsequent creation of scenarios for PySP.
+confidence regions and subsequent creation of scenarios for stochastic programming.
 
 Citation for parmest
 ^^^^^^^^^^^^^^^^^^^^

doc/OnlineDocs/contributed_packages/parmest/installation.rst

@@ -12,6 +12,7 @@ Python package dependencies
 #. numpy
 #. pandas
 #. pyomo
+#. mpisppy (optional)
 #. matplotlib (optional)
 #. scipy.stats (optional)
 #. seaborn (optional)

pyomo/contrib/parmest/__init__.py

@@ -8,3 +8,28 @@
 #  rights in this software.
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
+
+from pyomo.common.deprecation import relocated_module_attribute
+
+relocated_module_attribute(
+    'create_ef', 
+    'pyomo.contrib.parmest.utils.create_ef', 
+    '6.4.2', 
+    '6.5.0')
+relocated_module_attribute(
+    'ipopt_solver_wrapper ', 
+    'pyomo.contrib.parmest.utils.ipopt_solver_wrapper ', 
+    '6.4.2', 
+    '6.5.0')
+relocated_module_attribute(
+    'mpi_utils ', 
+    'pyomo.contrib.parmest.utils.mpi_utils ', 
+    '6.4.2', 
+    '6.5.0')
+relocated_module_attribute(
+    'scenario_tree ', 
+    'pyomo.contrib.parmest.utils.scenario_tree ', 
+    '6.4.2', 
+    '6.5.0')
+
+

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

@@ -14,30 +14,26 @@
 """
 import pandas as pd
 from pyomo.environ import (
-    ConcreteModel, Var, PositiveReals, Objective, Constraint, maximize,
+    ConcreteModel, Param, Var, PositiveReals, Objective, Constraint, maximize,
     SolverFactory
 )
 
+
 def reactor_design_model(data):
 
     # Create the concrete model
     model = ConcreteModel()
 
     # Rate constants
-    model.k1 = Var(initialize = 5.0/6.0, within=PositiveReals) # min^-1
-    model.k2 = Var(initialize = 5.0/3.0, within=PositiveReals) # min^-1
-    model.k3 = Var(initialize = 1.0/6000.0, within=PositiveReals) # m^3/(gmol min)
-    model.k1.fixed = True
-    model.k2.fixed = True
-    model.k3.fixed = True
-
+    model.k1 = Param(initialize = 5.0/6.0, within=PositiveReals, mutable=True) # min^-1
+    model.k2 = Param(initialize = 5.0/3.0, within=PositiveReals, mutable=True) # min^-1
+    model.k3 = Param(initialize = 1.0/6000.0, within=PositiveReals, mutable=True) # m^3/(gmol min)
+
     # Inlet concentration of A, gmol/m^3
-    model.caf = Var(initialize = float(data['caf']), within=PositiveReals)
-    model.caf.fixed = True
+    model.caf = Param(initialize = float(data['caf']), within=PositiveReals)
 
 	# Space velocity (flowrate/volume)
-    model.sv = Var(initialize = float(data['sv']), within=PositiveReals)
-    model.sv.fixed = True
+    model.sv = Param(initialize = float(data['sv']), within=PositiveReals)
 
     # Outlet concentration of each component
     model.ca = Var(initialize = 5000.0, within=PositiveReals) 

pyomo/contrib/parmest/parmest.py

@@ -25,8 +25,8 @@
     import mpisppy.opt.ef as st
     import mpisppy.scenario_tree as scenario_tree
 else:
-    import pyomo.contrib.parmest.create_ef as local_ef
-    import pyomo.contrib.parmest.scenario_tree as scenario_tree
+    import pyomo.contrib.parmest.utils.create_ef as local_ef
+    import pyomo.contrib.parmest.utils.scenario_tree as scenario_tree
 
 import re
 import importlib as im
@@ -47,8 +47,7 @@
 from pyomo.opt import SolverFactory
 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
+import pyomo.contrib.parmest.utils as utils
 import pyomo.contrib.parmest.graphics as graphics
 
 parmest_available = numpy_available & pandas_available & scipy_available
@@ -170,19 +169,19 @@ def _experiment_instance_creation_callback(scenario_name, node_names=None, cb_da
                                                 nonant_list=nonant_list,
                                                 scen_model=instance)]
 
-
     if "ThetaVals" in outer_cb_data:
         thetavals = outer_cb_data["ThetaVals"]
 
         # dlw august 2018: see mea code for more general theta
         for vstr in thetavals:
-            object = instance.find_component(vstr)
+            theta_cuid = ComponentUID(vstr)
+            theta_object = theta_cuid.find_component_on(instance)
             if thetavals[vstr] is not None:
                 #print("Fixing",vstr,"at",str(thetavals[vstr]))
-                object.fix(thetavals[vstr])
+                theta_object.fix(thetavals[vstr])
             else:
                 #print("Freeing",vstr)
-                object.fixed = False
+                theta_object.unfix()
 
     return instance
 
@@ -327,17 +326,33 @@ def _create_parmest_model(self, data):
         """
         Modify the Pyomo model for parameter estimation
         """
-        from pyomo.core import Objective
-
         model = self.model_function(data)
 
         if (len(self.theta_names) == 1) and (self.theta_names[0] == 'parmest_dummy_var'):
             model.parmest_dummy_var = pyo.Var(initialize = 1.0)
+
+        # Add objective function (optional)
+        if self.obj_function:
+            for obj in model.component_objects(pyo.Objective):
+                if obj.name in ["Total_Cost_Objective"]:
+                    raise RuntimeError("Parmest will not override the existing model Objective named "+ obj.name)
+                obj.deactivate()
+
+            for expr in model.component_data_objects(pyo.Expression):
+                if expr.name in ["FirstStageCost", "SecondStageCost"]:
+                    raise RuntimeError("Parmest will not override the existing model Expression named "+ expr.name)
+            model.FirstStageCost = pyo.Expression(expr=0)
+            model.SecondStageCost = pyo.Expression(rule=_SecondStageCostExpr(self.obj_function, data))
 
+            def TotalCost_rule(model):
+                return model.FirstStageCost + model.SecondStageCost
+            model.Total_Cost_Objective = pyo.Objective(rule=TotalCost_rule, sense=pyo.minimize)
+
+        # Convert theta Params to Vars, and unfix theta Vars
+        model = utils.convert_params_to_vars(model, self.theta_names)
+
+        # Update theta names list to use CUID string representation
         for i, theta in enumerate(self.theta_names):
-            # First, leverage the parser in ComponentUID to locate the
-            # component.  If that fails, fall back on the original
-            # (insecure) use of 'eval'
             var_cuid = ComponentUID(theta)
             var_validate = var_cuid.find_component_on(model)
             if var_validate is None:
@@ -346,29 +361,14 @@ def _create_parmest_model(self, data):
                     (i, theta))
             else:
                 try:
-                    # If the component that was found is not a variable,
+                    # If the component is not a variable,
                     # this will generate an exception (and the warning
                     # in the 'except')
                     var_validate.unfix()
-                    # We want to standardize on the CUID string
-                    # representation
                     self.theta_names[i] = repr(var_cuid)
                 except:
                     logger.warning(theta + ' is not a variable')
-
-        if self.obj_function:
-            for obj in model.component_objects(Objective):
-                obj.deactivate()
-
-            def FirstStageCost_rule(model):
-                return 0
-            model.FirstStageCost = pyo.Expression(rule=FirstStageCost_rule)
-            model.SecondStageCost = pyo.Expression(rule=_SecondStageCostExpr(self.obj_function, data))
-
-            def TotalCost_rule(model):
-                return model.FirstStageCost + model.SecondStageCost
-            model.Total_Cost_Objective = pyo.Objective(rule=TotalCost_rule, sense=pyo.minimize)
-
+
         self.parmest_model = model
 
         return model
@@ -405,8 +405,8 @@ def _Q_opt(self, ThetaVals=None, solver="ef_ipopt",
 
         # (Bootstrap scenarios will use indirection through the bootlist)
         if bootlist is None:
-            senario_numbers = list(range(len(self.callback_data)))
-            scen_names = ["Scenario{}".format(i) for i in senario_numbers]
+            scenario_numbers = list(range(len(self.callback_data)))
+            scen_names = ["Scenario{}".format(i) for i in scenario_numbers]
         else:
             scen_names = ["Scenario{}".format(i) for i in range(len(bootlist))]
 
@@ -574,13 +574,12 @@ def _Q_at_theta(self, thetavals):
 
         # start block of code to deal with models with no constraints
         # (ipopt will crash or complain on such problems without special care)
-        instance = _experiment_instance_creation_callback("FOO1", None, dummy_cb)
+        instance = _experiment_instance_creation_callback("FOO0", None, dummy_cb)
         try: # deal with special problems so Ipopt will not crash
             first = next(instance.component_objects(pyo.Constraint, active=True))
+            active_constraints = True
         except:
-            sillylittle = True 
-        else:
-            sillylittle = False
+            active_constraints = False 
         # end block of code to deal with models with no constraints
 
         WorstStatus = pyo.TerminationCondition.optimal
@@ -589,12 +588,12 @@ def _Q_at_theta(self, thetavals):
         for snum in senario_numbers:
             sname = "scenario_NODE"+str(snum)
             instance = _experiment_instance_creation_callback(sname, None, dummy_cb)
-            if not sillylittle:
+            if active_constraints:
                 if self.diagnostic_mode:
                     print('      Experiment = ',snum)
                     print('     First solve with with special diagnostics wrapper')
                     status_obj, solved, iters, time, regu \
-                        = ipopt_solver_wrapper.ipopt_solve_with_stats(instance, optimizer, max_iter=500, max_cpu_time=120)
+                        = utils.ipopt_solve_with_stats(instance, optimizer, max_iter=500, max_cpu_time=120)
                     print("   status_obj, solved, iters, time, regularization_stat = ",
                            str(status_obj), str(solved), str(iters), str(time), str(regu))
 
@@ -608,10 +607,11 @@ def _Q_at_theta(self, thetavals):
                     # DLW: Aug2018: not distinguishing "middlish" conditions
                     if WorstStatus != pyo.TerminationCondition.infeasible:
                         WorstStatus = results.solver.termination_condition
-                    
+
             objobject = getattr(instance, self._second_stage_cost_exp)
             objval = pyo.value(objobject)
             totobj += objval
+
         retval = totobj / len(senario_numbers) # -1??
 
         return retval, thetavals, WorstStatus
@@ -728,7 +728,7 @@ def theta_est_bootstrap(self, bootstrap_samples, samplesize=None,
         global_list = self._get_sample_list(samplesize, bootstrap_samples, 
                                             replacement)
 
-        task_mgr = mpiu.ParallelTaskManager(bootstrap_samples)
+        task_mgr = utils.ParallelTaskManager(bootstrap_samples)
         local_list = task_mgr.global_to_local_data(global_list)
 
         bootstrap_theta = list()
@@ -781,7 +781,7 @@ def theta_est_leaveNout(self, lNo, lNo_samples=None, seed=None,
 
         global_list = self._get_sample_list(samplesize, lNo_samples, replacement=False)
 
-        task_mgr = mpiu.ParallelTaskManager(len(global_list))
+        task_mgr = utils.ParallelTaskManager(len(global_list))
         local_list = task_mgr.global_to_local_data(global_list)
 
         lNo_theta = list()
@@ -901,7 +901,7 @@ def objective_at_theta(self, theta_values):
         # for parallel code we need to use lists and dicts in the loop
         theta_names = theta_values.columns
         all_thetas = theta_values.to_dict('records')
-        task_mgr = mpiu.ParallelTaskManager(len(all_thetas))
+        task_mgr = utils.ParallelTaskManager(len(all_thetas))
         local_thetas = task_mgr.global_to_local_data(all_thetas)
 
         # walk over the mesh, return objective function