implemented window and fixed short window, also debug pyomo opt mode

src/dispatch/Dispatcher.py

@@ -112,14 +112,15 @@ def validate(self, components, activity, times, meta):
 
   # ---------------------------------------------
   # UTILITY METHODS
-  def _compute_cashflows(self, components, activity, times, meta, state_args=None):
+  def _compute_cashflows(self, components, activity, times, meta, state_args=None, time_offset=0):
     """
       Method to compute CashFlow evaluations given components and their activity.
       @ In, components, list, HERON components whose cashflows should be evaluated
       @ In, activity, DispatchState instance, activity by component/resources/time
       @ In, times, np.array(float), time values to evaluate; may be length 1 or longer
       @ In, meta, dict, additional info to be passed through to functional evaluations
       @ In, state_args, dict, optional, additional arguments to pass while getting activity state
+      @ In, time_offset, int, optional, increase time index tracker by this value if provided
       @ Out, total, float, total cashflows for given components
     """
     if state_args is None:
@@ -138,7 +139,7 @@ def _compute_cashflows(self, components, activity, times, meta, state_args=None)
         specific_activity = {}
         for resource in resource_indexer[comp]:
           specific_activity[resource] = activity.get_activity(comp, resource, time, **state_args)
-        specific_meta['HERON']['time_index'] = t
+        specific_meta['HERON']['time_index'] = t + time_offset
         specific_meta['HERON']['time_value'] = time
         cfs = comp.get_state_cost(specific_activity, specific_meta)
         time_subtotal = sum(cfs.values())

src/dispatch/pyomo_dispatch.py

@@ -57,6 +57,16 @@ def get_input_specs(cls):
       @ Out, specs, InputData, specs
     """
     specs = InputData.parameterInputFactory('pyomo', ordered=False, baseNode=None)
+    specs.addSub(InputData.parameterInputFactory('rolling_window_length', contentType=InputTypes.IntegerType,
+        descr=r"""Sets the length of the rolling window that the Pyomo optimization algorithm
+        uses to break down histories. Longer window lengths will minimize boundary effects, such as
+        nonoptimal storage dispatch, at the cost of slower optimization solves.
+        Note that if the rolling window results in a window of length 1 (such as at the end of a history),
+        this can cause problems for pyomo.
+        \default{24}"""))
+    specs.addSub(InputData.parameterInputFactory('debug_mode', contentType=InputTypes.BoolType,
+        descr=r"""Enables additional printing in the pyomo dispatcher. Highly discouraged for production runs.
+              \default{False}."""))
     # TODO specific for pyomo dispatcher
     return specs
 
@@ -67,6 +77,7 @@ def __init__(self):
       @ Out, None
     """
     self.name = 'PyomoDispatcher' # identifying name
+    self.debug_mode = False       # whether to print additional information
     self._window_len = 24         # time window length to dispatch at a time # FIXME user input
 
   def read_input(self, specs):
@@ -75,7 +86,15 @@ def read_input(self, specs):
       @ In, specs, RAVEN InputData, specifications
       @ Out, None
     """
-    print('DEBUGG specs:', specs)
+    super().read_input(specs)
+
+    window_len_node = specs.findFirst('rolling_window_length')
+    if window_len_node is not None:
+      self._window_len = window_len_node.value
+
+    debug_node = specs.findFirst('debug_mode')
+    if debug_node is not None:
+      self.debug_mode = debug_node.value
 
   ### API
   def dispatch(self, case, components, sources, meta):
@@ -101,13 +120,26 @@ def dispatch(self, case, components, sources, meta):
     while start_index < final_index:
       end_index = start_index + self._window_len
       if end_index > final_index:
-        end_index = final_index # + 1?
+        end_index = final_index
+      if end_index - start_index == 1:
+        # TODO custom error raise for catching in DispatchManager?
+        raise IOError("A rolling window of length 1 was requested, but this causes crashes in pyomo. " +
+                      "Change the length of the rolling window to avoid length 1 histories.")
       specific_time = time[start_index:end_index]
       print('DEBUGG starting window {} to {}'.format(start_index, end_index))
       start = time_mod.time()
-      subdisp = self.dispatch_window(specific_time,
+      # set initial storage levels
+      initial_levels = {}
+      for comp in components:
+        if comp.get_interaction().is_type('Storage'):
+          if start_index == 0:
+            initial_levels[comp] = comp.get_interaction().get_initial_level(meta)
+          else:
+            initial_levels[comp] = subdisp[comp.name][comp.get_interaction().get_resource()][-1]
+      # dispatch
+      subdisp = self.dispatch_window(specific_time, start_index,
                                      case, components, sources, resources,
-                                     meta)
+                                     initial_levels, meta)
       end = time_mod.time()
       print('DEBUGG solve time: {} s'.format(end-start))
       # store result in corresponding part of dispatch
@@ -118,23 +150,26 @@ def dispatch(self, case, components, sources, meta):
     return dispatch
 
   ### INTERNAL
-  def dispatch_window(self, time,
+  def dispatch_window(self, time, time_offset,
                       case, components, sources, resources,
-                      meta):
+                      initial_storage, meta):
     """
       Dispatches one part of a rolling window.
       @ In, time, np.array, value of time to evaluate
+      @ In, time_offset, int, offset of the time index in the greater history
       @ In, case, HERON Case, Case that this dispatch is part of
       @ In, components, list, HERON components available to the dispatch
       @ In, sources, list, HERON source (placeholders) for signals
       @ In, resources, list, sorted list of all resources in problem
+      @ In, initial_storage, dict, initial storage levels if any
       @ In, meta, dict, additional variables passed through
       @ Out, result, dict, results of window dispatch
     """
     # build the Pyomo model
     # TODO abstract this model as much as possible BEFORE, then concrete initialization per window
     m = pyo.ConcreteModel()
     # indices
+    # XXX FIXME the time offset index is needed for the stuff coming from ARMAs!!!!
     C = np.arange(0, len(components), dtype=int) # indexes component
     R = np.arange(0, len(resources), dtype=int) # indexes resources
     # T = np.arange(start_index, end_index, dtype=int) # indexes resources
@@ -143,6 +178,7 @@ def dispatch_window(self, time,
     m.R = pyo.Set(initialize=R)
     m.T = pyo.Set(initialize=T)
     m.Times = time
+    m.time_offset = time_offset
     m.resource_index_map = meta['HERON']['resource_indexer'] # maps the resource to its index WITHIN APPLICABLE components (sparse matrix)
                                                              #   e.g. component: {resource: local index}, ... etc}
     # properties
@@ -160,13 +196,14 @@ def dispatch_window(self, time,
       self._create_capacity(m, comp, prod_name, meta)    # capacity constraints
       self._create_transfer(m, comp, prod_name)    # transfer functions (constraints)
       # ramp rates TODO ## INCLUDING previous-time boundary condition TODO
-    self._create_conservation(m, resources, meta) # conservation of resources (e.g. production == consumption)
+    self._create_conservation(m, resources, initial_storage, meta) # conservation of resources (e.g. production == consumption)
     self._create_objective(meta, m) # objective
     # start a solution search
     done_and_checked = False
     attempts = 0
     # DEBUGG show variables, bounds
-    # self._debug_pyomo_print(m)
+    if self.debug_mode:
+      self._debug_pyomo_print(m)
     while not done_and_checked:
       attempts += 1
       print(f'DEBUGG solve attempt {attempts} ...:')
@@ -202,8 +239,9 @@ def dispatch_window(self, time,
         done_and_checked = True
       if attempts > 100:
         raise RuntimeError('Exceeded validation attempt limit!')
-    # soln.write() # DEBUGG
-    # self._debug_print_soln(m) # DEBUGG
+    if self.debug_mode:
+      soln.write()
+      self._debug_print_soln(m)
     # return dict of numpy arrays
     result = self._retrieve_solution(m)
     return result
@@ -281,7 +319,7 @@ def _create_capacity(self, m, comp, prod_name, meta):
     caps = []
     mins = []
     for t, time in enumerate(m.Times):
-      meta['HERON']['time_index'] = t
+      meta['HERON']['time_index'] = t + m.time_offset
       cap = comp.get_capacity(meta)[0][cap_res] # value of capacity limit (units of governing resource)
       caps.append(cap)
       minimum = comp.get_minimum(meta)[0][cap_res]
@@ -328,16 +366,17 @@ def _create_transfer(self, m, comp, prod_name):
       constr = pyo.Constraint(m.T, rule=rule)
       setattr(m, rule_name, constr)
 
-  def _create_conservation(self, m, resources, meta):
+  def _create_conservation(self, m, resources, initial_storage, meta):
     """
       Creates pyomo conservation constraints
       @ In, m, pyo.ConcreteModel, associated model
       @ In, resources, list, list of resources in problem
+      @ In, initial_storage, dict, initial storage levels
       @ In, meta, dict, dictionary of state variables
       @ Out, None
     """
     for res, resource in enumerate(resources):
-      rule = partial(self._conservation_rule, meta, resource)
+      rule = partial(self._conservation_rule, initial_storage, meta, resource)
       constr = pyo.Constraint(m.T, rule=rule)
       setattr(m, '{r}_conservation'.format(r=resource), constr)
 
@@ -458,12 +497,14 @@ def _cashflow_rule(self, meta, m):
     """
     activity = m.Activity # dict((comp, getattr(m, f"{comp.name}_production")) for comp in m.Components)
     state_args = {'valued': False}
-    total = self._compute_cashflows(m.Components, activity, m.Times, meta, state_args=state_args)
+    total = self._compute_cashflows(m.Components, activity, m.Times, meta,
+                                    state_args=state_args, time_offset=m.time_offset)
     return total
 
-  def _conservation_rule(self, meta, res, m, t):
+  def _conservation_rule(self, initial_storage, meta, res, m, t):
     """
       Constructs conservation constraints.
+      @ In, initial_storage, dict, initial storage levels at t==0 (not t+offset==0)
       @ In, meta, dict, dictionary of state variables
       @ In, res, str, name of resource
       @ In, m, pyo.ConcreteModel, associated model
@@ -484,7 +525,7 @@ def _conservation_rule(self, meta, res, m, t):
             dt = m.Times[t] - m.Times[t-1]
           else:
             # FIXME check this with a variety of ValuedParams
-            previous = comp.get_interaction().get_initial_level(meta)
+            previous = initial_storage[comp] # comp.get_interaction().get_initial_level(meta)
             dt = m.Times[1] - m.Times[0]
           new = var[r, t]
           production = -1 * (new - previous) / dt # swap sign b/c negative is absorbing, positive is emitting
@@ -560,9 +601,9 @@ def _debug_print_soln(self, m):
       print('  component:', c, name)
       for res, r in m.resource_index_map[comp].items():
         print('    resource:', r, res)
-        for t, time_index in enumerate(m.T):
+        for t, time in enumerate(m.Times):
           prod = getattr(m, '{n}_production'.format(n=name))
-          print('      time:', t, time_index, prod[r, time_index].value)
+          print('      time:', t + m.time_offset, time, prod[r, t].value)
     print('*'*80)
 
 

tests/integration_tests/mechanics/pyomo_options/gold/Sweep_Runs_o/sweep.csv

@@ -0,0 +1,3 @@
+steamer_capacity,steam_storage_capacity,generator_capacity,electr_market_capacity,electr_flex_capacity,mean_NPV,std_NPV,med_NPV,max_NPV,min_NPV,perc_5_NPV,perc_95_NPV,samp_NPV,var_NPV,ProbabilityWeight,PointProbability,prefix,ProbabilityWeight-steamer_capacity
+1.0,100.0,-99.0,-2.0,-1e+200,28.0280044015,1.63813697062e-08,28.0280043932,28.0280044203,28.0280043909,28.0280043909,28.0280044203,3.0,2.68349273452e-16,0.5,0.010101010101,1,0.5
+100.0,100.0,-99.0,-2.0,-1e+200,463.491305175,2.89609407576e-07,463.491305156,463.491305474,463.491304895,463.491304895,463.491305474,3.0,8.38736089568e-14,0.5,0.010101010101,2,0.5

tests/integration_tests/mechanics/pyomo_options/heron_input.xml

@@ -0,0 +1,129 @@
+<HERON>
+  <TestInfo>
+    <name>storage</name>
+    <author>talbpaul</author>
+    <created>2020-09-28</created>
+    <description>
+      Tests including a storage unit as part of a flexible analysis case. Extends from
+      workflows/storage case. Note with the shorter rolling window, the optimizer cannot
+      make as optimal decisions as with the full window, resulting in differences especially
+      in the 100 steam production case.
+    </description>
+    <classesTested>HERON</classesTested>
+  </TestInfo>
+
+  <Case name="Sweep_Runs">
+    <mode>sweep</mode>
+    <num_arma_samples>3</num_arma_samples>
+    <time_discretization>
+      <time_variable>Time</time_variable>
+      <end_time>2</end_time>
+      <num_steps>21</num_steps>
+    </time_discretization>
+    <economics>
+      <ProjectTime>3</ProjectTime>
+      <DiscountRate>0.08</DiscountRate>
+      <tax>0.0</tax>
+      <inflation>0.0</inflation>
+      <verbosity>50</verbosity>
+    </economics>
+    <dispatcher>
+      <pyomo>
+        <rolling_window_length>8</rolling_window_length>
+        <debug_mode>True</debug_mode>
+      </pyomo>
+    </dispatcher>
+  </Case>
+
+  <Components>
+    <Component name="steamer">
+      <produces resource="steam" dispatch="fixed">
+        <capacity resource="steam">
+          <sweep_values>1, 100</sweep_values>
+        </capacity>
+      </produces>
+      <economics>
+        <lifetime>27</lifetime>
+      </economics>
+    </Component>
+
+    <Component name="steam_storage">
+      <stores resource="steam" dispatch="independent">
+        <capacity resource="steam">
+          <fixed_value>100</fixed_value>
+        </capacity>
+        <initial_stored>
+          <fixed_value>1</fixed_value>
+        </initial_stored>
+      </stores>
+      <economics>
+        <lifetime>10</lifetime>
+      </economics>
+    </Component>
+
+    <Component name="generator">
+      <produces resource="electricity" dispatch="independent">
+        <consumes>steam</consumes>
+        <capacity resource="steam">
+          <fixed_value>-99</fixed_value>
+        </capacity>
+        <transfer>
+          <linear>
+            <rate resource="steam">-1</rate>
+            <rate resource="electricity">0.5</rate>
+          </linear>
+        </transfer>
+      </produces>
+      <economics>
+        <lifetime>27</lifetime>
+      </economics>
+    </Component>
+
+    <Component name="electr_market">
+      <demands resource="electricity" dispatch="dependent"><!---->
+        <capacity>
+          <fixed_value>-2</fixed_value>
+        </capacity>
+      </demands>
+      <economics>
+        <lifetime>30</lifetime>
+        <CashFlow name="e_sales" type="repeating" taxable='True' inflation='none' mult_target='False'>
+          <driver>
+            <Function method="electric_consume">transfers</Function>
+            <!-- <resource multiplier="-1" absolute="True">electricity</resource> -->
+          </driver>
+          <reference_price>
+            <fixed_value>0.5</fixed_value>
+          </reference_price>
+        </CashFlow>
+      </economics>
+    </Component>
+
+
+    <Component name="electr_flex">
+      <demands resource="electricity" dispatch="dependent">
+        <capacity>
+          <fixed_value>-1e200</fixed_value>
+        </capacity>
+      </demands>
+      <economics>
+        <lifetime>30</lifetime>
+        <CashFlow name="e_sales" type="repeating" taxable='True' inflation='none' mult_target='False'>
+          <driver>
+            <Function method="electric_consume">transfers</Function>
+          </driver>
+          <reference_price>
+            <Function method="flex_price">transfers</Function>
+          </reference_price>
+        </CashFlow>
+      </economics>
+    </Component>
+
+  </Components>
+
+  <DataGenerators>
+    <ARMA name='Price' variable="Signal">%HERON%/tests/integration_tests/ARMA/Sine/arma.pk</ARMA>
+    <Function name="transfers">transfers.py</Function>
+  </DataGenerators>
+
+</HERON>

tests/integration_tests/mechanics/pyomo_options/tests

@@ -0,0 +1,14 @@
+[Tests]
+  [./PyomoOptions]
+    type = HeronIntegration
+    input = heron_input.xml
+    # prereq = SineArma
+    [./csv]
+      type = OrderedCSV
+      output = 'Sweep_Runs_o/sweep.csv'
+      zero_threshold = 1e-6
+      rel_err = 1e-6
+    [../]
+  [../]
+
+[]