Add mapping keys between nominal and differential values for the Diff…

…erentialParameterSweep object (watertap-org#1079)

* added mapping information for differntial samples.

* all tests pass in serial.

* mapping indices are now stored in the case of differential parameter sweep and tests pass. Cleanup pending.

* run black.

* cleanup and run black.

* fixed input for regular parameter sweep. recursive and differential next.

* fixed recursive parameter sweep outputs.

* updated tests to reflect how the new implementation works.

* minor cleanup.

* commit midway of debugging PR.

* backtracked storing inputs within the output key.

* fix differential ps tests.

* fix parameter_sweep_optimize.

* added documentation on how to use differential parameter sweep.

---------

Co-authored-by: Kinshuk <[email protected]>

docs/how_to_guides/how_to_run_differential_parameter_sweep.rst

@@ -0,0 +1,217 @@
+How to Run Differential Parameter Sweep
+============================================
+
+Overview
+--------
+
+This guide explains how to run a differential parameter sweep for conducting, e.g., value of innovation (VOI), analysis.
+
+How To
+------
+
+As before, we begin by importing or explicitly programming any functions relating to flowsheet building/specification, simulation, and optimization setup steps.  We will use the same RO with energy recovery flowsheet for this example.
+
+.. testsetup::
+
+    # quiet idaes logs
+    import idaes.logger as idaeslogger
+    idaeslogger.getLogger('ideas.core').setLevel('CRITICAL')
+    idaeslogger.getLogger('ideas.core.util.scaling').setLevel('CRITICAL')
+    idaeslogger.getLogger('idaes.init').setLevel('CRITICAL')
+
+.. testcode::
+
+    # replace this with your own flowsheet module, e.g.
+    # import my_flowsheet_module as mfm
+    import watertap.examples.flowsheets.RO_with_energy_recovery.RO_with_energy_recovery as RO_flowsheet
+
+Once this is done, we import the differential parameter sweep tool and sampling classes.
+
+.. testcode::
+
+    from watertap.tools.parameter_sweep import differential_parameter_sweep, UniformSample, NormalSample
+
+We will use the :ref:`same setup steps as before<how_to_use_parameter_sweep>` to set up the generating functions for our model, sweep params, and outputs:
+
+.. testcode::
+
+    def build_model(**kwargs):
+
+        # replace these function calls with
+        # those in your own flowsheet module
+
+        # set up system
+        m = RO_flowsheet.build()
+        RO_flowsheet.set_operating_conditions(m)
+        RO_flowsheet.initialize_system(m)
+
+        # simulate
+        RO_flowsheet.solve(m)
+
+        # set up the model for optimization
+        RO_flowsheet.optimize_set_up(m)
+
+        return m
+
+Once the model has been setup, we specify the variables to sample using a dictionary
+
+.. testcode::
+
+    def build_sweep_params(model, num_samples=5):
+        sweep_params = dict()
+        sweep_params["A_comp"] = UniformSample(model.fs.RO.A_comp, 4.2e-12, 2.1e-11, num_samples)
+        sweep_params["membrane_cost"] = UniformSample(
+            model.fs.costing.reverse_osmosis.membrane_cost, 30, 10, num_samples
+        )
+        return sweep_params
+
+where the ``A_comp`` attribute will be randomly selected from a uniform distribution of values in the range :math:`[4.2e-12, 2.1e-11]` and ``membrane_cost`` will be drawn from a uniform distribution between :math:`[30, 10]`.  For this example, we will extract flowsheet outputs associated with cost, the levelized cost of water (LCOW) and energy consumption (EC), defined via another dictionary.
+Next we will construct a specification dictionary to run the differential parameter sweep.
+
+.. testcode::
+
+    def build_diff_sweep_param_specs(model):
+        differential_sweep_specs = dict()
+
+        differential_sweep_specs["A_comp"] = {
+            "diff_sample_type": NormalSample,
+            "std_dev": 0.3e-12,
+            "pyomo_object": model.fs.RO.A_comp,
+        }
+
+        differential_sweep_specs["membrane_cost"] = {
+            "diff_sample_type": UniformSample,
+            "diff_mode": "percentile",
+            "nominal_lb" : sweep_params["membrane_cost"].lower_limit,
+            "nominal_ub" : sweep_params["membrane_cost"].upper_limit,
+            "relative_lb" : -0.05,
+            "relative_ub" : -0.05,
+            "pyomo_object": model.fs.costing.reverse_osmosis.membrane_cost,
+        }
+
+        return differential_sweep_specs
+
+``differential_sweep_specs`` is a specification dictionary that contains details for how to construct the parameter sweep dictionary for differential sweep. This is a nested dictionary where the first level denotes the variable names for which the differential sweep needs to be carried out. The second level denotes various options to be used for each variable. The number of samples for each differential sweep is specified while initializing the ``DifferentialParameterSweep`` object using the keyword ``num_diff_samples``. There are 4 modes of setting up a variable to undergo differential sweep:
+
+#. ``NormalSample`` : Uses the nominal value as the mean and expects ``std_dev`` key for the differential sweep sampling. It looks like the following:
+
+    .. code-block:: python
+
+        differential_sweep_specs["A_comp"] = {
+                "diff_sample_type": NormalSample,
+                "std_dev": 0.3e-12,
+                "pyomo_object": model.fs.RO.A_comp,
+            }
+
+    This differential mode is unique to variables that expect normal sampling. *All other sampling types expect one of the other 3 differential modes below.*
+
+#. ``sum`` : Perturbs the nominal value by a certain absolute percentage to create an upper and lower bound for the differential solve. The logic in the code looks as follows:
+
+    .. code-block:: python
+
+        lower_bound = nominal_val * (1 - relative_lb)
+        upper_bound = nominal_val * (1 + relative_ub)
+
+#. ``product``: Perturbs the nominal value by a scaling factor to create upper and lower bounds for the differential sweep. It uses the following logic
+
+    .. code-block:: python
+
+        lower_bound = nominal_val * relative_lb
+        upper_bound = nominal_val * relative_ub  
+
+#. ``percentile``: Perturbs the nominal value by a percentage of the difference between the nominal upper and lower bound values. The logic is 
+
+    .. code-block:: python
+
+        delta_nominal = abs(upper_nominal - lower_nominal)
+        lower_bound = nominal_val + delta_nominal * relative_lb
+        upper_bound = nominal_val + delta_nominal * relative_ub
+
+An example differential sweep spec dictionary may look like the following:
+
+.. code-block:: python
+
+    differential_sweep_specs = dict()
+    differential_sweep_specs["membrane_cost"] = {
+            "diff_sample_type": UniformSample,
+            "diff_mode": "percentile",
+            "nominal_lb" : sweep_params["membrane_cost"].lower_limit,
+            "nominal_ub" : sweep_params["membrane_cost"].upper_limit,
+            "relative_lb" : -0.05,
+            "relative_ub" : -0.05,
+            "pyomo_object": model.fs.costing.reverse_osmosis.membrane_cost,
+        }
+    differential_sweep_specs["px_cost"] = {
+        "diff_sample_type": LinearSample,
+        "diff_mode": "sum",
+        "relative_lb" : -0.05,
+        "relative_ub" : -0.05,
+        "pyomo_object": m.fs.costing.pressure_exchanger.cost,
+    }
+    differential_sweep_specs["px_efficiency"] = {
+        "diff_sample_type": UniformSample,
+        "diff_mode": "product",
+        "relative_lb" : 0.001,
+        "relative_ub" : 0.001,
+        "pyomo_object": m.fs.PXR.efficiency_pressure_exchanger,
+    }
+
+.. important:: The user can only conduct differential sweeps for variables specified with ``sweep_params``.
+
+Continuing with the example test code from above, we will use the following function for building the outputs.
+
+.. testcode::
+
+    def build_outputs(model, sweep_params):
+        outputs = dict()
+        outputs['EC'] = model.fs.costing.specific_energy_consumption
+        outputs['LCOW'] = model.fs.costing.LCOW
+        return outputs
+
+With the flowsheet defined and suitably initialized, along with the definitions for ``sweep_params``, ``differential_sweep_specs``, and ``outputs`` on hand, we can call the ``differential_parameter_sweep`` function as before.
+
+.. note:: This documentation currently uses the older API for calling the differential parameter sweep. This API will be deprecated in the near future. The documentation will be changed to reflect this accordingly. We recommend running the differential parameter sweep in serial or with MPI only.
+
+.. testcode::
+
+    # Define the local results directory, num_samples, and seed (if desired)
+    num_samples = 5
+    seed = None
+
+    model = build_model()
+    sweep_params = build_sweep_params(model, num_samples=num_samples)
+    differential_sweep_specs = build_diff_sweep_param_specs(model)
+    outputs = build_outputs(model, sweep_params)
+
+    # Run the parameter sweep
+    global_results = differential_parameter_sweep(
+            model, 
+            sweep_params, 
+            differential_sweep_specs,
+            outputs, 
+            h5_results_file_name='monte_carlo_results.h5',
+            optimize_function=RO_flowsheet.optimize,
+            debugging_data_dir=None,
+            num_samples=num_samples,
+            num_diff_samples=2,
+            seed=seed,
+        )
+
+.. testoutput::
+
+    ...
+
+.. testcleanup::
+
+    import os
+    import shutil
+    try:
+        os.remove('monte_carlo_results.h5')
+        os.remove('monte_carlo_results.h5.txt')
+    except:
+        print("monte_carlo_results.h5 does not exist, nothing to delete.")
+
+Module Documentation
+--------------------
+
+* :mod:`watertap.tools.parameter_sweep`

docs/how_to_guides/index.rst

@@ -18,6 +18,7 @@ How To Guides
    how_to_scale_chemical_process_energy_balance
    how_to_use_parameter_sweep
    how_to_use_parameter_sweep_monte_carlo
+   how_to_run_differential_parameter_sweep
    how_to_install_electrolyte_database
    how_to_use_EDB
    how_to_use_ui_api

watertap/tools/parameter_sweep/parameter_sweep.py

@@ -361,8 +361,8 @@ def _create_local_output_skeleton(self, model, sweep_params, outputs, num_sample
         sweep_param_objs = ComponentSet()
 
         # Store the inputs
-        for param_name, sweep_param in sweep_params.items():
-            var = sweep_param.pyomo_object
+        for param_name, sampling_obj in sweep_params.items():
+            var = sampling_obj.pyomo_object
             sweep_param_objs.add(var)
             output_dict["sweep_params"][
                 param_name
@@ -373,9 +373,16 @@ def _create_local_output_skeleton(self, model, sweep_params, outputs, num_sample
             for pyo_obj in model.component_data_objects(
                 (pyo.Var, pyo.Expression, pyo.Objective, pyo.Param), active=True
             ):
-                output_dict["outputs"][
-                    pyo_obj.name
-                ] = self._create_component_output_skeleton(pyo_obj, num_samples)
+                # We do however need to make sure that the short name for the inputs is used here
+                for param_name, sampling_obj in sweep_params.items():
+                    if pyo_obj.name == sampling_obj.pyomo_object.name:
+                        output_dict["outputs"][
+                            param_name
+                        ] = self._create_component_output_skeleton(pyo_obj, num_samples)
+                    else:
+                        output_dict["outputs"][
+                            pyo_obj.name
+                        ] = self._create_component_output_skeleton(pyo_obj, num_samples)
 
         else:
             # Save only the outputs specified in the outputs dictionary
@@ -452,7 +459,7 @@ def _create_global_output(self, local_output_dict, req_num_samples=None):
             global_output_dict = copy.deepcopy(local_output_dict)
             # Create a global value array of inputs in the dictionary
             for key, item in global_output_dict.items():
-                if key != "solve_successful":
+                if key in ["sweep_params", "outputs"]:
                     for subkey, subitem in item.items():
                         subitem["value"] = np.zeros(num_total_samples, dtype=float)
 
@@ -461,7 +468,7 @@ def _create_global_output(self, local_output_dict, req_num_samples=None):
 
         # Finally collect the values
         for key, item in local_output_dict.items():
-            if key != "solve_successful":
+            if key in ["sweep_params", "outputs"]:
                 for subkey, subitem in item.items():
                     self.parallel_manager.gather_arrays_to_root(
                         sendbuf=subitem["value"],
@@ -711,7 +718,7 @@ def _combine_gather_results(self, all_results):
     one process's run.
     """
 
-    def _combine_outputs(self, gathered_results):
+    def _combine_output_array(self, gathered_results):
         outputs = gathered_results["outputs"]
         if len(outputs) == 0:
             return []
@@ -850,8 +857,8 @@ def parameter_sweep(
             all_parameter_combinations,
         )
 
-        global_sweep_results = self._combine_gather_results(all_results)
-        combined_outputs = self._combine_outputs(global_sweep_results)
+        global_sweep_results_dict = self._combine_gather_results(all_results)
+        combined_output_arr = self._combine_output_array(global_sweep_results_dict)
 
         # save the results for all simulations run by this process and its children
         for results in self.parallel_manager.results_from_local_tree(all_results):
@@ -860,14 +867,14 @@ def parameter_sweep(
                 results.parameters,
                 all_parameter_combinations,
                 results.results,
-                global_sweep_results,
-                combined_outputs,
+                global_sweep_results_dict,
+                combined_output_arr,
                 process_number=results.process_number,
             )
 
-        global_save_data = np.hstack((all_parameter_combinations, combined_outputs))
+        global_save_data = np.hstack((all_parameter_combinations, combined_output_arr))
 
-        return global_save_data, global_sweep_results
+        return global_save_data, global_sweep_results_dict
 
 
 class RecursiveParameterSweep(_ParameterSweepBase):