From 8495d93d6b167c65d523e5147c4c1ebe18935006 Mon Sep 17 00:00:00 2001 From: Dylan McDowell Date: Wed, 24 Jan 2024 09:51:43 -0700 Subject: [PATCH] Refactor Pyomo Dispatch (#326) * Refactor Pyomo Dispatch * Refactor dispatch method to be more clear * More refactoring of pyomo_dispatch * Testing out PyomoModelHandler * Fix PyomoModelHandler * Working PMH for now... * Remove Commented Code and Refactor Dispatcher * Add Docstrings * Get Tests Passing * Trim trailing whitespace * Address Comments + Other Minor Changes * Fix more docstrings * Trim trailing whitespace * Fix Docstrings --- src/_utils.py | 11 - src/dispatch/DispatchState.py | 43 + src/dispatch/Dispatcher.py | 114 +-- src/dispatch/PyomoModelHandler.py | 558 +++++++++++++ src/dispatch/PyomoRuleLibrary.py | 302 +++++++ src/dispatch/putils.py | 183 ++++ src/dispatch/pyomo_dispatch.py | 1296 +++++------------------------ 7 files changed, 1290 insertions(+), 1217 deletions(-) create mode 100644 src/dispatch/PyomoModelHandler.py create mode 100644 src/dispatch/PyomoRuleLibrary.py create mode 100644 src/dispatch/putils.py diff --git a/src/_utils.py b/src/_utils.py index 58e2cf13..0d33f573 100644 --- a/src/_utils.py +++ b/src/_utils.py @@ -92,17 +92,6 @@ def get_farm_loc(raven_path=None): # Added by Haoyu Wang, May 25, 2022 farm_loc = plugin_handler.getPluginLocation('FARM') return farm_loc -def get_all_resources(components): - """ - Provides a set of all resources used among all components - @ In, components, list, HERON component objects - @ Out, resources, list, resources used in case - """ - res = set() - for comp in components: - res.update(comp.get_resources()) - return res - def get_project_lifetime(case, components): """ obtains the project lifetime diff --git a/src/dispatch/DispatchState.py b/src/dispatch/DispatchState.py index 18b133b3..11cacb15 100644 --- a/src/dispatch/DispatchState.py +++ b/src/dispatch/DispatchState.py @@ -199,3 +199,46 @@ def set_activity_vector(self, comp, res, values, tracker='production', start_idx r = self._resources[comp][res] self._data[f'{comp.name}_{tracker}'][r, start_idx:end_idx] = values + +# DispatchState for Pyomo dispatcher +class PyomoState(DispatchState): + def __init__(self): + """ + Constructor. + @ In, None + @ Out, None + """ + DispatchState.__init__(self) + self._model = None # Pyomo model object + + def initialize(self, components, resources_map, times, model): + """ + Connect information about this State to other objects + @ In, components, list, HERON components + @ In, resources_map, dict, map of component names to resources used + @ In, times, np.array, values of "time" this state represents + @ In, model, pyomo.Model, associated model for this state + @ Out, None + """ + DispatchState.initialize(self, components, resources_map, times) + self._model = model + + def get_activity_indexed(self, comp, activity, r, t, valued=True, **kwargs): + """ + Getter for activity level. + @ In, comp, HERON Component, component whose information should be retrieved + @ In, activity, str, tracking variable name for activity subset + @ In, r, int, index of resource to retrieve (as given by meta[HERON][resource_indexer]) + @ In, t, int, index of time at which activity should be provided + @ In, valued, bool, optional, if True then get float value instead of pyomo expression + @ In, kwargs, dict, additional pass-through keyword arguments + @ Out, activity, float, amount of resource "res" produced/consumed by "comp" at time "time"; + note positive is producting, negative is consuming + """ + prod = getattr(self._model, f'{comp.name}_{activity}')[r, t] + if valued: + return prod() + return prod + + def set_activity_indexed(self, comp, r, t, value, valued=False): + raise NotImplementedError diff --git a/src/dispatch/Dispatcher.py b/src/dispatch/Dispatcher.py index dd386453..5d06128d 100644 --- a/src/dispatch/Dispatcher.py +++ b/src/dispatch/Dispatcher.py @@ -4,9 +4,16 @@ """ Base class for dispatchers. """ -from ravenframework.utils import InputData, InputTypes +from ravenframework.utils import InputData from ravenframework.BaseClasses import MessageUser, InputDataUser +class DispatchError(Exception): + """ + Custom exception for dispatch errors. + """ + pass + + class Dispatcher(MessageUser, InputDataUser): """ Base class for strategies for consecutive dispatching in a continuous period. @@ -122,108 +129,3 @@ def validate(self, components, activity, times, meta): # no validator, nothing needs to be changed return {} - # --------------------------------------------- - # UTILITY METHODS - 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: - state_args = {} - - if meta['HERON']['Case'].use_levelized_inner: - total = self._compute_levelized_cashflows(components, activity, times, meta, state_args, time_offset) - return total - - total = 0 - specific_meta = dict(meta) # TODO what level of copying do we need here? - resource_indexer = meta['HERON']['resource_indexer'] - - #print('DEBUGG computing cashflows!') - for comp in components: - #print(f'DEBUGG ... comp {comp.name}') - specific_meta['HERON']['component'] = comp - comp_subtotal = 0 - for t, time in enumerate(times): - #print(f'DEBUGG ... ... time {t}') - # NOTE care here to assure that pyomo-indexed variables work here too - specific_activity = {} - for tracker in comp.get_tracking_vars(): - specific_activity[tracker] = {} - for resource in resource_indexer[comp]: - specific_activity[tracker][resource] = activity.get_activity(comp, tracker, resource, time, **state_args) - specific_meta['HERON']['time_index'] = t + time_offset - specific_meta['HERON']['time_value'] = time - cfs = comp.get_state_cost(specific_activity, specific_meta, marginal=True) - time_subtotal = sum(cfs.values()) - comp_subtotal += time_subtotal - total += comp_subtotal - return total - - def _compute_levelized_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 - """ - total = 0 - specific_meta = dict(meta) # TODO what level of copying do we need here? - resource_indexer = meta['HERON']['resource_indexer'] - - # How does this work? - # The general equation looks like: - # - # SUM(Non-Multiplied Terms) + x * SUM(Multiplied Terms) = Target - # - # and we are solving for `x`. Target is 0 by default. Terms here are marginal cashflows. - # Summations here occur over: components, time steps, tracking variables, and resources. - # Typically, there is only 1 multiplied term/cash flow. - - multiplied = 0 - non_multiplied = 0 - - for comp in components: - specific_meta['HERON']['component'] = comp - multiplied_comp = 0 - non_multiplied_comp = 0 - for t, time in enumerate(times): - # NOTE care here to assure that pyomo-indexed variables work here too - specific_activity = {} - for tracker in comp.get_tracking_vars(): - specific_activity[tracker] = {} - for resource in resource_indexer[comp]: - specific_activity[tracker][resource] = activity.get_activity(comp, tracker, resource, time, **state_args) - specific_meta['HERON']['time_index'] = t + time_offset - specific_meta['HERON']['time_value'] = time - cfs = comp.get_state_cost(specific_activity, specific_meta, marginal=True) - - # there is an assumption here that if a component has a levelized cost, marginal cashflow - # then it is the only marginal cashflow - if comp.levelized_meta: - for cf in comp.levelized_meta.keys(): - lcf = cfs.pop(cf) # this should be ok as long as HERON init checks are successful - multiplied_comp += lcf - else: - time_subtotal = sum(cfs.values()) - non_multiplied_comp += time_subtotal - - multiplied += multiplied_comp - non_multiplied += non_multiplied_comp - - # at this point, there should be a not None NPV Target - multiplied += self._eps - total = (meta['HERON']['Case'].npv_target - non_multiplied) / multiplied - total *= -1 - return total diff --git a/src/dispatch/PyomoModelHandler.py b/src/dispatch/PyomoModelHandler.py new file mode 100644 index 00000000..a23b77d4 --- /dev/null +++ b/src/dispatch/PyomoModelHandler.py @@ -0,0 +1,558 @@ +# Copyright 2020, Battelle Energy Alliance, LLC +# ALL RIGHTS RESERVED +""" + This module constructs the dispatch optimization model used by HERON. +""" +import numpy as np +import pyomo.environ as pyo + +from . import PyomoRuleLibrary as prl +from . import putils +from .DispatchState import PyomoState + +class PyomoModelHandler: + """ + Class for constructing the pyomo model, populate with objective/constraints, and evaluate. + """ + + _eps = 1e-9 + + def __init__(self, time, time_offset, case, components, resources, initial_storage, meta) -> None: + """ + Initializes a PyomoModelHandler instance. + @ In, time, np.array(float), time values to evaluate; may be length 1 or longer + @ In, time_offset, int, optional, increase time index tracker by this value if provided + @ In, case, HERON Case, case to evaluate + @ In, components, list, HERON components to evaluate + @ In, resources, list, HERON resources to evaluate + @ In, initial_storage, dict, initial storage levels + @ In, meta, dict, additional state information + @ Out, None + """ + self.time = time + self.time_offset = time_offset + self.case = case + self.components = components + self.resources = resources + self.initial_storage = initial_storage + self.meta = meta + self.model = self.build_model() + + + def build_model(self): + """ + Construct the skeleton of the pyomo model. + @ In, None + @ Out, model, pyo.ConcreteModel, model + """ + model = pyo.ConcreteModel() + C = np.arange(0, len(self.components), dtype=int) # indexes component + R = np.arange(0, len(self.resources), dtype=int) # indexes resources + T = np.arange(0, len(self.time), dtype=int) # indexes resources + model.C = pyo.Set(initialize=C) + model.R = pyo.Set(initialize=R) + model.T = pyo.Set(initialize=T) + model.Times = self.time + model.time_offset = self.time_offset + # maps the resource to its index WITHIN APPLICABLE components (sparse matrix) + # e.g. component: {resource: local index}, ... etc} + model.resource_index_map = self.meta['HERON']['resource_indexer'] + # properties + model.Case = self.case + model.Components = self.components + model.Activity = PyomoState() + model.Activity.initialize(model.Components, model.resource_index_map, model.Times, model) + return model + + + def populate_model(self): + """ + Populate the pyomo model with generated objectives/contraints. + @ In, None + @ Out, None + """ + for comp in self.components: + self._process_component(comp) + self._create_conservation() # conservation of resources (e.g. production == consumption) + self._create_objective() # objective function + + + def _process_component(self, component): + """ + Determine what kind of component this is and process it accordingly. + @ In, component, HERON Component, component to process + @ Out, None + """ + interaction = component.get_interaction() + if interaction.is_governed(): + self._process_governed_component(component, interaction) + elif interaction.is_type("Storage"): + self._create_storage(component) + else: + self._create_production(component) + + + def _process_governed_component(self, component, interaction): + """ + Process a component that is governed since it requires special attention. + @ In, component, HERON Component, component to process + @ In, interaction, HERON Interaction, interaction to process + @ Out, None + """ + self.meta["request"] = {"component": component, "time": self.time} + if interaction.is_type("Storage"): + self._process_storage_component(component, interaction) + else: + activity = interaction.get_strategy().evaluate(self.meta)[0]['level'] + self._create_production_param(component, activity) + + + def _process_storage_component(self, component, interaction): + """ + Process a storage component. + @ In, component, HERON Component, component to process + @ In, interaction, HERON Interaction, interaction to process + """ + activity = interaction.get_strategy().evaluate(self.meta)[0]["level"] + self._create_production_param(component, activity, tag="level") + dt = self.model.Times[1] - self.model.Times[0] + rte2 = component.get_sqrt_RTE() + deltas = np.zeros(len(activity)) + deltas[1:] = activity[1:] - activity[:-1] + deltas[0] = activity[0] - interaction.get_initial_level(self.meta) + charge = np.where(deltas > 0, -deltas / dt / rte2, 0) + discharge = np.where(deltas < 0, -deltas / dt * rte2, 0) + self._create_production_param(component, charge, tag="charge") + self._create_production_param(component, discharge, tag="discharge") + + + def _create_production_limit(self, validation): + """ + Creates pyomo production constraint given validation errors + @ In, validation, dict, information from Validator about limit violation + @ Out, None + """ + # TODO could validator write a symbolic expression on request? That'd be sweet. + comp = validation['component'] + resource = validation['resource'] + r = self.model.resource_index_map[comp][resource] + t = validation['time_index'] + limit = validation['limit'] + limits = np.zeros(len(self.model.Times)) + limits[t] = limit + limit_type = validation['limit_type'] + prod_name = f'{comp.name}_production' + rule = lambda mod: prl.prod_limit_rule(prod_name, r, limits, limit_type, t, mod) + constr = pyo.Constraint(rule=rule) + counter = 1 + name_template = f'{comp.name}_{resource}_{t}_vld_limit_constr_{{i}}' + # make sure we get a unique name for this constraint + name = name_template.format(i=counter) + while getattr(self.model, name, None) is not None: + counter += 1 + name = name_template.format(i=counter) + setattr(self.model, name, constr) + print(f'DEBUGG added validation constraint "{name}"') + + + def _create_production_param(self, comp, values, tag=None): + """ + Creates production pyomo fixed parameter object for a component + @ In, comp, HERON Component, component to make production variables for + @ In, values, np.array(float), values to set for param + @ In, tag, str, optional, if not None then name will be component_[tag] + @ Out, prod_name, str, name of production variable + """ + name = comp.name + if tag is None: + tag = 'production' + # create pyomo indexer for this component's resources + res_indexer = pyo.Set(initialize=range(len(self.model.resource_index_map[comp]))) + setattr(self.model, f'{name}_res_index_map', res_indexer) + prod_name = f'{name}_{tag}' + init = (((0, t), values[t]) for t in self.model.T) + prod = pyo.Param(res_indexer, self.model.T, initialize=dict(init)) + setattr(self.model, prod_name, prod) + return prod_name + + + def _create_production(self, comp): + """ + Creates all pyomo variable objects for a non-storage component + @ In, comp, HERON Component, component to make production variables for + @ Out, prod_name, str, name of the production variable + """ + prod_name = self._create_production_variable(comp) + ## if you cannot set limits directly in the production variable, set separate contraint: + ## Method 1: set variable bounds directly --> TODO more work needed, but would be nice + # lower, upper = self._get_prod_bounds(m, comp) + # limits should be None unless specified, so use "getters" from dictionaries + # bounds = lambda m, r, t: (lower.get(r, None), upper.get(r, None)) + ## Method 2: set variable bounds directly --> TODO more work needed, but would be nice + # self._create_capacity(m, comp, prod_name, meta) # capacity constraints + # transfer function governs input -> output relationship + self._create_transfer(comp, prod_name) + # ramp rates + if comp.ramp_limit is not None: + self._create_ramp_limit(comp, prod_name) + return prod_name + + + def _create_production_variable(self, comp, tag=None, add_bounds=True, **kwargs): + """ + Creates production pyomo variable object for a component + @ In, comp, HERON Component, component to make production variables for + @ In, tag, str, optional, if not None then name will be component_[tag]; otherwise "production" + @ In, add_bounds, bool, optional, if True then determine and set bounds for variable + @ In, kwargs, dict, optional, passalong kwargs to pyomo variable + @ Out, prod_name, str, name of production variable + """ + if tag is None: + tag = 'production' + name = comp.name + cap_res = comp.get_capacity_var() # name of resource that defines capacity + limit_r = self.model.resource_index_map[comp][cap_res] # production index of the governing resource + # create pyomo indexer for this component's resources + indexer_name = f'{name}_res_index_map' + indexer = getattr(self.model, indexer_name, None) + if indexer is None: + indexer = pyo.Set(initialize=range(len(self.model.resource_index_map[comp]))) + setattr(self.model, indexer_name, indexer) + prod_name = f'{name}_{tag}' + caps, mins = self._find_production_limits(comp) + if min(caps) < 0: + # quick check that capacities signs are consistent #FIXME: revisit, this is an assumption + assert max(caps) <= 0, \ + 'Capacities are inconsistent: mix of positive and negative values not currently supported.' + # we have a unit that's consuming, so we need to flip the variables to be sensible + mins, caps = caps, mins + inits = caps + else: + inits = mins + if add_bounds: + # create bounds based in min, max operation + bounds = lambda m, r, t: (mins[t] if r == limit_r else None, caps[t] if r == limit_r else None) + initial = lambda m, r, t: inits[t] if r == limit_r else 0 + else: + bounds = (None, None) + initial = 0 + # production variable depends on resources, time + #FIXME initials! Should be lambda with mins for tracking var! + prod = pyo.Var(indexer, self.model.T, initialize=initial, bounds=bounds, **kwargs) + # TODO it may be that we need to set variable values to avoid problems in some solvers. + # if comp.is_dispatchable() == 'fixed': + # for t, _ in enumerate(m.Times): + # prod[limit_r, t].fix(caps[t]) + setattr(self.model, prod_name, prod) + return prod_name + + + def _create_ramp_limit(self, comp, prod_name): + """ + Creates ramping limitations for a producing component + @ In, comp, HERON Component, component to make ramping limits for + @ In, prod_name, str, name of production variable + @ Out, None + """ + # ramping is defined in terms of the capacity variable + cap_res = comp.get_capacity_var() # name of resource that defines capacity + cap = comp.get_capacity(self.meta)[0][cap_res] + r = self.model.resource_index_map[comp][cap_res] # production index of the governing resource + # NOTE: this includes the built capacity * capacity factor, if any, which assumes + # the ramp rate depends on the available capacity, not the built capacity. + limit_delta = comp.ramp_limit * cap # NOTE: if cap is negative, then this is negative. + if limit_delta < 0: + neg_cap = True + else: + neg_cap = False + # if we're limiting ramp frequency, make vars and rules for that + if comp.ramp_freq: + # create binaries for tracking ramping + up = pyo.Var(self.model.T, initialize=0, domain=pyo.Binary) + down = pyo.Var(self.model.T, initialize=0, domain=pyo.Binary) + steady = pyo.Var(self.model.T, initialize=1, domain=pyo.Binary) + setattr(self.model, f'{comp.name}_up_ramp_tracker', up) + setattr(self.model, f'{comp.name}_down_ramp_tracker', down) + setattr(self.model, f'{comp.name}_steady_ramp_tracker', steady) + ramp_trackers = (down, up, steady) + else: + ramp_trackers = None + # limit production changes when ramping down + ramp_rule_down = lambda mod, t: prl.ramp_rule_down(prod_name, r, limit_delta, neg_cap, t, mod, bins=ramp_trackers) + constr = pyo.Constraint(self.model.T, rule=ramp_rule_down) + setattr(self.model, f'{comp.name}_ramp_down_constr', constr) + # limit production changes when ramping up + ramp_rule_up = lambda mod, t: prl.ramp_rule_up(prod_name, r, limit_delta, neg_cap, t, mod, bins=ramp_trackers) + constr = pyo.Constraint(self.model.T, rule=ramp_rule_up) + setattr(self.model, f'{comp.name}_ramp_up_constr', constr) + # if ramping frequency limit, impose binary constraints + if comp.ramp_freq: + # binaries rule, for exclusive choice up/down/steady + binaries_rule = lambda mod, t: prl.ramp_freq_bins_rule(down, up, steady, t, mod) + constr = pyo.Constraint(self.model.T, rule=binaries_rule) + setattr(self.model, f'{comp.name}_ramp_freq_binaries', constr) + # limit frequency of ramping + # TODO calculate "tao" window using ramp freq and dt + # -> for now, just use the integer for number of windows + freq_rule = lambda mod, t: prl.ramp_freq_rule(down, up, comp.ramp_freq, t, mod) + constr = pyo.Constraint(self.model.T, rule=freq_rule) + setattr(self.model, f'{comp.name}_ramp_freq_constr', constr) + + + def _create_capacity_constraints(self, comp, prod_name): + """ + Creates pyomo capacity constraints + @ In, comp, HERON Component, component to make variables for + @ In, prod_name, str, name of production variable + @ Out, None + """ + cap_res = comp.get_capacity_var() # name of resource that defines capacity + r = self.model.resource_index_map[comp][cap_res] # production index of the governing resource + caps, mins = self._find_production_limits(comp) + # capacity + max_rule = lambda mod, t: prl.capacity_rule(prod_name, r, caps, mod, t) + constr = pyo.Constraint(self.model.T, rule=max_rule) + setattr(self.model, f'{comp.name}_{cap_res}_capacity_constr', constr) + # minimum + min_rule = lambda mod, t: prl.min_prod_rule(prod_name, r, caps, mins, mod, t) + constr = pyo.Constraint(self.model.T, rule=min_rule) + # set initial conditions + for t, time in enumerate(self.model.Times): + cap = caps[t] + if cap == mins[t]: + # initialize values so there's no boundary errors + var = getattr(self.model, prod_name) + values = var.get_values() + for k in values: + values[k] = cap + var.set_values(values) + setattr(self.model, f'{comp.name}_{cap_res}_minprod_constr', constr) + + + def _find_production_limits(self, comp): + """ + Determines the capacity limits of a unit's operation, in time. + @ In, comp, HERON Component, component to make variables for + @ Out, caps, array, max production values by time + @ Out, mins, array, min production values by time + """ + cap_res = comp.get_capacity_var() # name of resource that defines capacity + r = self.model.resource_index_map[comp][cap_res] # production index of the governing resource + # production is always lower than capacity + ## NOTE get_capacity returns (data, meta) and data is dict + ## TODO does this work with, e.g., ARMA-based capacities? + ### -> "time" is stored on "m" and could be used to correctly evaluate the capacity + caps = [] + mins = [] + for t, time in enumerate(self.model.Times): + self.meta['HERON']['time_index'] = t + self.model.time_offset + cap = comp.get_capacity(self.meta)[0][cap_res] # value of capacity limit (units of governing resource) + caps.append(cap) + if (comp.is_dispatchable() == 'fixed'): + minimum = cap + else: + minimum = comp.get_minimum(self.meta)[0][cap_res] + mins.append(minimum) + return caps, mins + + + def _create_transfer(self, comp, prod_name): + """ + Creates pyomo transfer function constraints + @ In, comp, HERON Component, component to make variables for + @ In, prod_name, str, name of production variable + @ Out, None + """ + name = comp.name + # transfer functions + # e.g. 2A + 3B -> 1C + 2E + # get linear coefficients + # TODO this could also take a transfer function from an external Python function assuming + # we're careful about how the expression-vs-float gets used + # and figure out how to handle multiple ins, multiple outs + ratios = putils.get_transfer_coeffs(self.model, comp) + ref_r, ref_name, _ = ratios.pop('__reference', (None, None, None)) + for resource, ratio in ratios.items(): + r = self.model.resource_index_map[comp][resource] + rule_name = f'{name}_{resource}_{ref_name}_transfer' + rule = lambda mod, t: prl.transfer_rule(ratio, r, ref_r, prod_name, mod, t) + constr = pyo.Constraint(self.model.T, rule=rule) + setattr(self.model, rule_name, constr) + + + def _create_storage(self, comp): + """ + Creates storage pyomo variable objects for a storage component + Similar to create_production, but for storages + @ In, comp, HERON Component, component to make production variables for + @ Out, None + """ + prefix = comp.name + # what resource index? Isn't it always 0? # assumption + r = 0 # NOTE this is only true if each storage ONLY uses 1 resource + # storages require a few variables: + # (1) a level tracker, + level_name = self._create_production_variable(comp, tag='level') + # -> set operational limits + # self._create_capacity(m, comp, level_name, meta) + # (2, 3) separate charge/discharge trackers, so we can implement round-trip efficiency and ramp rates + charge_name = self._create_production_variable(comp, tag='charge', add_bounds=False, within=pyo.NonPositiveReals) + discharge_name = self._create_production_variable(comp, tag='discharge', add_bounds=False, within=pyo.NonNegativeReals) + # balance level, charge/discharge + level_rule_name = prefix + '_level_constr' + rule = lambda mod, t: prl.level_rule(comp, level_name, charge_name, discharge_name, self.initial_storage, r, mod, t) + setattr(self.model, level_rule_name, pyo.Constraint(self.model.T, rule=rule)) + # periodic boundary condition for storage level + if comp.get_interaction().apply_periodic_level: + periodic_rule_name = prefix + '_level_periodic_constr' + rule = lambda mod, t: prl.periodic_level_rule(comp, level_name, self.initial_storage, r, mod, t) + setattr(self.model, periodic_rule_name, pyo.Constraint(self.model.T, rule=rule)) + + # (4) a binary variable to track whether we're charging or discharging, to prevent BOTH happening + # -> 0 is charging, 1 is discharging + # -> TODO make this a user-based option to disable, if they want to allow dual operation + # -> -> but they should really think about if that's what they want! + # FIXME currently introducing the bigM strategy also makes solves numerically unstable, + # and frequently results in spurious errors. For now, disable it. + allow_both = True # allow simultaneous charging and discharging + if not allow_both: + bin_name = self._create_production_variable(comp, tag='dcforcer', add_bounds=False, within=pyo.Binary) + # we need a large epsilon, but not so large that addition stops making sense + # -> we don't know what any values for this component will be! How do we choose? + # -> NOTE that choosing this value has VAST impact on solve stability!! + large_eps = 1e8 #0.01 * sys.float_info.max + # charging constraint: don't charge while discharging (note the sign matters) + charge_rule_name = prefix + '_charge_constr' + rule = lambda mod, t: prl.charge_rule(charge_name, bin_name, large_eps, r, mod, t) + setattr(self.model, charge_rule_name, pyo.Constraint(self.model.T, rule=rule)) + discharge_rule_name = prefix + '_discharge_constr' + rule = lambda mod, t: prl.discharge_rule(discharge_name, bin_name, large_eps, r, mod, t) + setattr(self.model, discharge_rule_name, pyo.Constraint(self.model.T, rule=rule)) + + + def _create_conservation(self): + """ + Creates pyomo conservation constraints + @ In, None + @ Out, None + """ + for resource in self.resources: + rule = lambda mod, t: prl.conservation_rule(resource, mod, t) + constr = pyo.Constraint(self.model.T, rule=rule) + setattr(self.model, f'{resource}_conservation', constr) + + + def _create_objective(self): + """ + Creates pyomo objective function + @ In, None + @ Out, None + """ + # cashflow eval + rule = lambda mod: prl.cashflow_rule(self._compute_cashflows, self.meta, mod) + self.model.obj = pyo.Objective(rule=rule, sense=pyo.maximize) + + + 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: + state_args = {} + + if meta['HERON']['Case'].use_levelized_inner: + total = self._compute_levelized_cashflows(components, activity, times, meta, state_args, time_offset) + return total + + total = 0 + specific_meta = dict(meta) # TODO what level of copying do we need here? + resource_indexer = meta['HERON']['resource_indexer'] + + #print('DEBUGG computing cashflows!') + for comp in components: + #print(f'DEBUGG ... comp {comp.name}') + specific_meta['HERON']['component'] = comp + comp_subtotal = 0 + for t, time in enumerate(times): + #print(f'DEBUGG ... ... time {t}') + # NOTE care here to assure that pyomo-indexed variables work here too + specific_activity = {} + for tracker in comp.get_tracking_vars(): + specific_activity[tracker] = {} + for resource in resource_indexer[comp]: + specific_activity[tracker][resource] = activity.get_activity(comp, tracker, resource, time, **state_args) + specific_meta['HERON']['time_index'] = t + time_offset + specific_meta['HERON']['time_value'] = time + cfs = comp.get_state_cost(specific_activity, specific_meta, marginal=True) + time_subtotal = sum(cfs.values()) + comp_subtotal += time_subtotal + total += comp_subtotal + return total + + def _compute_levelized_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 + """ + total = 0 + specific_meta = dict(meta) # TODO what level of copying do we need here? + resource_indexer = meta['HERON']['resource_indexer'] + + # How does this work? + # The general equation looks like: + # + # SUM(Non-Multiplied Terms) + x * SUM(Multiplied Terms) = Target + # + # and we are solving for `x`. Target is 0 by default. Terms here are marginal cashflows. + # Summations here occur over: components, time steps, tracking variables, and resources. + # Typically, there is only 1 multiplied term/cash flow. + + multiplied = 0 + non_multiplied = 0 + + for comp in components: + specific_meta['HERON']['component'] = comp + multiplied_comp = 0 + non_multiplied_comp = 0 + for t, time in enumerate(times): + # NOTE care here to assure that pyomo-indexed variables work here too + specific_activity = {} + for tracker in comp.get_tracking_vars(): + specific_activity[tracker] = {} + for resource in resource_indexer[comp]: + specific_activity[tracker][resource] = activity.get_activity(comp, tracker, resource, time, **state_args) + specific_meta['HERON']['time_index'] = t + time_offset + specific_meta['HERON']['time_value'] = time + cfs = comp.get_state_cost(specific_activity, specific_meta, marginal=True) + + # there is an assumption here that if a component has a levelized cost, marginal cashflow + # then it is the only marginal cashflow + if comp.levelized_meta: + for cf in comp.levelized_meta.keys(): + lcf = cfs.pop(cf) # this should be ok as long as HERON init checks are successful + multiplied_comp += lcf + else: + time_subtotal = sum(cfs.values()) + non_multiplied_comp += time_subtotal + + multiplied += multiplied_comp + non_multiplied += non_multiplied_comp + + # at this point, there should be a not None NPV Target + multiplied += self._eps + total = (meta['HERON']['Case'].npv_target - non_multiplied) / multiplied + total *= -1 + return total diff --git a/src/dispatch/PyomoRuleLibrary.py b/src/dispatch/PyomoRuleLibrary.py new file mode 100644 index 00000000..3bd26d53 --- /dev/null +++ b/src/dispatch/PyomoRuleLibrary.py @@ -0,0 +1,302 @@ +# Copyright 2020, Battelle Energy Alliance, LLC +# ALL RIGHTS RESERVED +""" +Library of Pyomo rules for HERON dispatch. +""" +import pyomo.environ as pyo + +def charge_rule(charge_name, bin_name, large_eps, r, m, t) -> bool: + """ + Constructs pyomo don't-charge-while-discharging constraints. + For storage units specificially. + @ In, charge_name, str, name of charging variable + @ In, bin_name, str, name of forcing binary variable + @ In, large_eps, float, a large-ish number w.r.t. storage activity + @ In, r, int, index of stored resource (is this always 0?) + @ In, m, pyo.ConcreteModel, associated model + @ In, t, int, time index for capacity rule + @ Out, rule, bool, inequality used to limit charge behavior + """ + charge_var = getattr(m, charge_name) + bin_var = getattr(m, bin_name) + return -charge_var[r, t] <= (1 - bin_var[r, t]) * large_eps + +def discharge_rule(discharge_name, bin_name, large_eps, r, m, t) -> bool: + """ + Constructs pyomo don't-discharge-while-charging constraints. + For storage units specificially. + @ In, discharge_name, str, name of discharging variable + @ In, bin_name, str, name of forcing binary variable + @ In, large_eps, float, a large-ish number w.r.t. storage activity + @ In, r, int, index of stored resource (is this always 0?) + @ In, m, pyo.ConcreteModel, associated model + @ In, t, int, time index for capacity rule + @ Out, rule, bool, inequality used to limit discharge behavior + """ + discharge_var = getattr(m, discharge_name) + bin_var = getattr(m, bin_name) + return discharge_var[r, t] <= bin_var[r, t] * large_eps + +def level_rule(comp, level_name, charge_name, discharge_name, initial_storage, r, m, t) -> bool: + """ + Constructs pyomo charge-discharge-level balance constraints. + For storage units specificially. + @ In, comp, Component, storage component of interest + @ In, level_name, str, name of level-tracking variable + @ In, charge_name, str, name of charging variable + @ In, discharge_name, str, name of discharging variable + @ In, initial_storage, dict, initial storage levels by component + @ In, r, int, index of stored resource (is this always 0?) + @ In, m, pyo.ConcreteModel, associated model + @ In, t, int, time index for capacity rule + @ Out, rule, bool, inequality used to limit level behavior + """ + level_var = getattr(m, level_name) + charge_var = getattr(m, charge_name) + discharge_var = getattr(m, discharge_name) + if t > 0: + previous = level_var[r, t-1] + dt = m.Times[t] - m.Times[t-1] + else: + previous = initial_storage[comp] + dt = m.Times[1] - m.Times[0] + rte2 = comp.get_sqrt_RTE() # square root of the round-trip efficiency + production = - rte2 * charge_var[r, t] - discharge_var[r, t] / rte2 + return level_var[r, t] == previous + production * dt + +def periodic_level_rule(comp, level_name, initial_storage, r, m, t) -> bool: + """ + Mandates storage units end with the same level they start with, which prevents + "free energy" or "free sink" due to initial starting levels. + For storage units specificially. + @ In, comp, Component, storage component of interest + @ In, level_name, str, name of level-tracking variable + @ In, initial_storage, dict, initial storage levels by component + @ In, r, int, index of stored resource (is this always 0?) + @ In, m, pyo.ConcreteModel, associated model + @ In, t, int, time index for capacity rule + @ Out, rule, bool, inequality used to limit level behavior + """ + return getattr(m, level_name)[r, m.T[-1]] == initial_storage[comp] + +def capacity_rule(prod_name, r, caps, m, t) -> bool: + """ + Constructs pyomo capacity constraints. + @ In, prod_name, str, name of production variable + @ In, r, int, index of resource for capacity constraining + @ In, caps, list(float), value to constrain resource at in time + @ In, m, pyo.ConcreteModel, associated model + @ In, t, int, time index for capacity rule + """ + kind = 'lower' if min(caps) < 0 else 'upper' + return prod_limit_rule(prod_name, r, caps, kind, t, m) + +def prod_limit_rule(prod_name, r, limits, kind, t, m) -> bool: + """ + Constructs pyomo production constraints. + @ In, prod_name, str, name of production variable + @ In, r, int, index of resource for capacity constraining + @ In, limits, list(float), values in time at which to constrain resource production + @ In, kind, str, either 'upper' or 'lower' for limiting production + @ In, t, int, time index for production rule (NOTE not pyomo index, rather fixed index) + @ In, m, pyo.ConcreteModel, associated model + @ Out, rule, bool, pyomo expression contraint for production limits + """ + prod = getattr(m, prod_name) + if kind == 'lower': + # production must exceed value + return prod[r, t] >= limits[t] + elif kind == 'upper': + return prod[r, t] <= limits[t] + else: + raise TypeError('Unrecognized production limit "kind":', kind) + +def ramp_rule_down(prod_name, r, limit, neg_cap, t, m, bins=None) -> bool: + """ + Constructs pyomo production ramping constraints for reducing production level. + Note that this is number-getting-less-positive for positive-defined capacity, while + it is number-getting-less-negative for negative-defined capacity. + This means that dQ is negative for positive-defined capacity, but positive for vice versa + @ In, prod_name, str, name of production variable + @ In, r, int, index of resource for capacity constraining + @ In, limit, float, limiting change in production level across time steps. NOTE: negative for negative-defined capacity. + @ In, neg_cap, bool, True if capacity is expressed as negative (consumer) + @ In, t, int, time index for ramp limit rule (NOTE not pyomo index, rather fixed index) + @ In, m, pyo.ConcreteModel, associated model + @ In, bins, tuple, optional, (lower, upper, steady) binaries if limiting ramp frequency + @ Out, rule, expression, evaluation for Pyomo constraint + """ + prod = getattr(m, prod_name) + if t == 0: + return pyo.Constraint.Skip + delta = prod[r, t] - prod[r, t-1] + # special treatment if we have frequency-limiting binaries available + if bins is None: + if neg_cap: + # NOTE change in production should be "less positive" than the max + # "negative decrease" in production (decrease is positive when defined by consuming) + return delta <= - limit + else: + # dq is negative, - limit is negative + return delta >= - limit + else: + eps = 1.0 # aux parameter to force binaries to behave, TODO needed? + down = bins[0][t] + up = bins[1][t] + # NOTE we're following the convention that "less negative" is ramping "down" + # for capacity defined by consumption + # e.g. consuming 100 ramps down to consuming 70 is (-100 -> -70), dq = 30 + if neg_cap: + # dq <= limit * dt * Bu + eps * Bd, if limit <= 0 + # dq is positive, - limit is positive + return delta <= - limit * down - eps * up + else: + # dq <= limit * dt * Bu - eps * Bd, if limit >= 0 + # dq is negative, - limit is negative + return delta >= - limit * down + eps * up + +def ramp_rule_up(prod_name, r, limit, neg_cap, t, m, bins=None) -> bool: + """ + Constructs pyomo production ramping constraints. + @ In, prod_name, str, name of production variable + @ In, r, int, index of resource for capacity constraining + @ In, limit, float, limiting change in production level across time steps + @ In, neg_cap, bool, True if capacity is expressed as negative (consumer) + @ In, t, int, time index for ramp limit rule (NOTE not pyomo index, rather fixed index) + @ In, m, pyo.ConcreteModel, associated model + @ In, bins, tuple, optional, (lower, steady, upper) binaries if limiting ramp frequency + @ Out, rule, expression, evaluation for Pyomo constraint + """ + prod = getattr(m, prod_name) + if t == 0: + return pyo.Constraint.Skip + delta = prod[r, t] - prod[r, t-1] + if bins is None: + if neg_cap: + # NOTE change in production should be "more positive" than the max + # "negative increase" in production (increase is negative when defined by consuming) + return delta >= limit + else: + # change in production should be less than the max production increase + return delta <= limit + else: + # special treatment if we have frequency-limiting binaries available + eps = 1.0 # aux parameter to force binaries to behave, TODO needed? + down = bins[0][t] + up = bins[1][t] + # NOTE we're following the convention that "more negative" is ramping "up" + # for capacity defined by consumption + # e.g. consuming 100 ramps up to consuming 130 is (-100 -> -130), dq = -30 + if neg_cap: + # dq >= limit * dt * Bu + eps * Bd, if limit <= 0 + return delta >= limit * up + eps * down + else: + # dq <= limit * dt * Bu - eps * Bd, if limit >= 0 + return delta <= limit * up - eps * down + +def ramp_freq_rule(Bd, Bu, tao, t, m) -> bool: + """ + Constructs pyomo frequency-of-ramp constraints. + @ In, Bd, bool var, binary tracking down-ramp events + @ In, Bu, bool var, binary tracking up-ramp events + @ In, tao, int, number of time steps to look back + @ In, t, int, time step indexer + @ In, m, pyo.ConcreteModel, pyomo model + @ Out, rule, expression, evaluation for Pyomo constraint + """ + if t == 0: + return pyo.Constraint.Skip + # looking-back-window shouldn't be longer than existing time + tao = min(t, tao) + # how many ramp-down events in backward window? + tally = sum(1 - Bd[tm] for tm in range(t - tao, t)) + # only allow ramping up if no rampdowns in back window + ## but we can't use if statements, so use binary math + return Bu[t] <= 1 / tao * tally + +def ramp_freq_bins_rule(Bd, Bu, Bn, t, m) -> bool: + """ + Constructs pyomo constraint for ramping event tracking variables. + This forces choosing between ramping up, ramping down, and steady state operation. + @ In, Bd, bool var, binary tracking down-ramp events + @ In, Bu, bool var, binary tracking up-ramp events + @ In, Bn, bool var, binary tracking no-ramp events + @ In, t, int, time step indexer + @ In, m, pyo.ConcreteModel, pyomo model + @ Out, rule, expression, evaluation for Pyomo constraint + """ + return Bd[t] + Bu[t] + Bn[t] == 1 + +def cashflow_rule(compute_cashflows, meta, m) -> float: + """ + Objective function rule. + @ In, compute_cashflows, function, function to compute cashflows + @ In, meta, dict, additional variable passthrough + @ In, m, pyo.ConcreteModel, associated model + @ Out, total, float, evaluation of cost + """ + activity = m.Activity + state_args = {'valued': False} + total = compute_cashflows(m.Components, activity, m.Times, meta, state_args=state_args, time_offset=m.time_offset) + return total + +def conservation_rule(res, m, t) -> bool: + """ + Constructs conservation constraints. + @ In, res, str, name of resource + @ In, m, pyo.ConcreteModel, associated model + @ In, t, int, index of time variable + @ Out, conservation, bool, balance check + """ + balance = 0 # sum of production rates, which needs to be zero + for comp, res_dict in m.resource_index_map.items(): + if res in res_dict: + # activity information depends on if storage or component + r = res_dict[res] + intr = comp.get_interaction() + if intr.is_type('Storage'): + # Storages have 3 variables: level, charge, and discharge + # -> so calculate activity + charge = getattr(m, f'{comp.name}_charge') + discharge = getattr(m, f'{comp.name}_discharge') + # note that "charge" is negative (as it's consuming) and discharge is positive + # -> so the intuitive |discharge| - |charge| becomes discharge + charge + production = discharge[r, t] + charge[r, t] + else: + var = getattr(m, f'{comp.name}_production') + # TODO move to this? balance += m._activity.get_activity(comp, res, t) + production = var[r, t] + balance += production + return balance == 0 # TODO tol? + +def min_prod_rule(prod_name, r, caps, minimums, m, t) -> bool: + """ + Constructs minimum production constraint + @ In, prod_name, str, name of production variable + @ In, r, int, index of resource for capacity constraining + @ In, caps, list(float), capacity(t) value for component + @ In, minimums, list(float), minimum allowable production in time + @ In, m, pyo.ConcreteModel, associated model + @ In, t, int, index of time variable + @ Out, minimum, bool, min check + """ + prod = getattr(m, prod_name) + # negative capacity means consuming instead of producing + if max(caps) > 0: + return prod[r, t] >= minimums[t] + else: + return prod[r, t] <= minimums[t] + +def transfer_rule(ratio, r, ref_r, prod_name, m, t) -> bool: + """ + Constructs transfer function constraints + @ In, ratio, float, ratio for resource to nominal first resource + @ In, r, int, index of transfer resource + @ In, ref_r, int, index of reference resource + @ In, prod_name, str, name of production variable + @ In, m, pyo.ConcreteModel, associated model + @ In, t, int, index of time variable + @ Out, transfer, bool, transfer ratio check + """ + prod = getattr(m, prod_name) + return prod[r, t] == prod[ref_r, t] * ratio # TODO tolerance?? diff --git a/src/dispatch/putils.py b/src/dispatch/putils.py new file mode 100644 index 00000000..940656d7 --- /dev/null +++ b/src/dispatch/putils.py @@ -0,0 +1,183 @@ +# Copyright 2020, Battelle Energy Alliance, LLC +# ALL RIGHTS RESERVED +""" +Pyomo Utilities for Model Dispatch. +""" +import platform +import numpy as np +import pyomo.environ as pyo +from pyomo.common.errors import ApplicationError + +def check_solver_availability(requested_solver: str) -> str: + """ + Check if any of the requested solvers are available. If not, display available options. + @ In, requested_solver, str, requested solver (e.g. 'cbc', 'glpk', 'ipopt') + @ Out, solver, str, name of solver that is available to use. + """ + # Choose solver; CBC is a great choice unless we're on Windows + if platform.system() == 'Windows': + platform_solvers = ['glpk', 'cbc', 'ipopt'] + else: + platform_solvers = ['cbc', 'glpk', 'ipopt'] + + solvers_to_check = platform_solvers if requested_solver is None else [requested_solver] + for solver in solvers_to_check: + if is_solver_available(solver): + # Early return if everything is a-ok + return solver + + # Otherwise raise an error + all_options = pyo.SolverFactory._cls.keys() + available_solvers = [op for op in all_options if not op.startswith('_') and is_solver_available(op)] + raise RuntimeError( + f'Requested solver "{requested_solver}" not found. Available options may include: {available_solvers}.' + ) + +def is_solver_available(solver: str) -> bool: + """ + Check if specified soler is available on the system. + @ In, solver, str, name of solver to check. + @ Out, is_available, bool, True if solver is available. + """ + try: + return pyo.SolverFactory(solver).available() + except (ApplicationError, NameError, ImportError): + return False + +def get_all_resources(components): + """ + Provides a set of all resources used among all components + @ In, components, list, HERON component objects + @ Out, resources, list, resources used in case + """ + res = set() + for comp in components: + res.update(comp.get_resources()) + return res + + +def get_initial_storage_levels(components: list, meta: dict, start_index: int) -> dict: + """ + Return initial storage levels for 'Storage' component types. + @ In, components, list, HERON components available to the dispatch. + @ In, meta, dict, additional variables passed through. + @ In, start_index, int, index of the start of the window. + @ Out, initial_levels, dict, initial storage levels for 'Storage' component types. + """ + 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) + # NOTE: There used to be an else conditional here that depended on the + # variable `subdisp` which was not defined yet. Leaving an unreachable + # branch of code, thus, I removed it. So currently, this function assumes + # start_index will always be zero, otherwise it will return an empty dict. + # Here was the line in case we need it in the future: + # else: initial_levels[comp] = subdisp[comp.name]['level'][comp.get_interaction().get_resource()][-1] + return initial_levels + + +def get_transfer_coeffs(m, comp) -> dict: + """ + Obtains transfer function ratios (assuming Linear ValuedParams) + Form: 1A + 3B -> 2C + 4D + Ratios are calculated with respect to first resource listed, so e.g. B = 3/1 * A + TODO I think we can handle general external functions, maybe? + @ In, m, pyo.ConcreteModel, associated model + @ In, comp, HERON component, component to get coefficients of + @ Out, ratios, dict, ratios of transfer function variables + """ + transfer = comp.get_interaction().get_transfer() + if transfer is None: + return {} + + # linear transfer coefficients, dict as {resource: coeff}, SIGNS MATTER + # it's all about ratios -> store as ratio of resource / first resource (arbitrary) + coeffs = transfer.get_coefficients() + coeffs_iter = iter(coeffs.items()) + first_name, first_coef = next(coeffs_iter) + first_r = m.resource_index_map[comp][first_name] + ratios = {'__reference': (first_r, first_name, first_coef)} + + for resource, coef in coeffs_iter: + ratios[resource] = coef / first_coef + + return ratios + +def retrieve_solution(m) -> dict: + """ + Extracts solution from Pyomo optimization + @ In, m, pyo.ConcreteModel, associated (solved) model + @ Out, result, dict, {comp: {activity: {resource: [production]}} e.g. generator[production][steam] + """ + return { + component.name: { + tag: retrieve_value_from_model(m, component, tag) + for tag in component.get_tracking_vars() + } + for component in m.Components + } + +def retrieve_value_from_model(m, comp, tag) -> dict: + """ + Retrieve values of a series from the pyomo model. + @ In, m, pyo.ConcreteModel, associated (solved) model + @ In, comp, Component, relevant component + @ In, tag, str, particular history type to retrieve + @ Out, result, dict, {resource: [array], etc} + """ + result = {} + prod = getattr(m, f'{comp.name}_{tag}') + kind = 'Var' if isinstance(prod, pyo.Var) else 'Param' + for res, comp_r in m.resource_index_map[comp].items(): + if kind == 'Var': + result[res] = np.array([prod[comp_r, t].value for t in m.T], dtype=float) + elif kind == 'Param': + result[res] = np.array([prod[comp_r, t] for t in m.T], dtype=float) + return result + +### DEBUG +def debug_pyomo_print(m) -> None: + """ + Prints the setup pieces of the Pyomo model + @ In, m, pyo.ConcreteModel, model to interrogate + @ Out, None + """ + print('/' + '='*80) + print('DEBUGG model pieces:') + print(' -> objective:') + print(' ', m.obj.pprint()) + print(' -> variables:') + for var in m.component_objects(pyo.Var): + print(' ', var.pprint()) + print(' -> constraints:') + for constr in m.component_objects(pyo.Constraint): + print(' ', constr.pprint()) + print('\\' + '='*80) + print('') + +def debug_print_soln(m) -> None: + """ + Prints the solution from the Pyomo model + @ In, m, pyo.ConcreteModel, model to interrogate + @ Out, None + """ + output = ['*' * 80, "DEBUGG solution:", f' objective value: {m.obj()}'] + for c, comp in enumerate(m.Components): + name = comp.name + output.append(f' component: {c} {name}') + for tracker in comp.get_tracking_vars(): + prod = getattr(m, f'{name}_{tracker}') + kind = 'Var' if isinstance(prod, pyo.Var) else 'Param' + for res, r in m.resource_index_map[comp].items(): + output.append(f' tracker: {tracker} resource {r}: {res}') + for t, time in enumerate(m.Times): + if kind == 'Var': + value = prod[r, t].value + elif kind == 'Param': + value = prod[r, t] + output.append(f' time: {t + m.time_offset} {time} {value}') + + output.append('*' * 80) + print('\n'.join(output)) \ No newline at end of file diff --git a/src/dispatch/pyomo_dispatch.py b/src/dispatch/pyomo_dispatch.py index 88477867..7841c5f4 100644 --- a/src/dispatch/pyomo_dispatch.py +++ b/src/dispatch/pyomo_dispatch.py @@ -1,46 +1,29 @@ - # Copyright 2020, Battelle Energy Alliance, LLC # ALL RIGHTS RESERVED """ pyomo-based dispatch strategy """ - -import os -import sys import time as time_mod -import platform -from itertools import compress import pprint - import numpy as np +import pyutilib.subprocess.GlobalData + import pyomo.environ as pyo from pyomo.opt import SolverStatus, TerminationCondition - from ravenframework.utils import InputData, InputTypes +from . import putils +from .PyomoModelHandler import PyomoModelHandler +from .Dispatcher import Dispatcher, DispatchError +from .DispatchState import NumpyState + # allows pyomo to solve on threaded processes -import pyutilib.subprocess.GlobalData pyutilib.subprocess.GlobalData.DEFINE_SIGNAL_HANDLERS_DEFAULT = False -from pyomo.common.errors import ApplicationError - -from .Dispatcher import Dispatcher -from .DispatchState import DispatchState, NumpyState -try: - import _utils as hutils -except (ModuleNotFoundError, ImportError): - sys.path.append(os.path.join(os.path.dirname(__file__), '..')) - import _utils as hutils - -# Choose solver; CBC is a great choice unless we're on Windows -if platform.system() == 'Windows': - SOLVERS = ['glpk', 'cbc', 'ipopt'] -else: - SOLVERS = ['cbc', 'glpk', 'ipopt'] # different solvers express "tolerance" for converging solution in different # ways. Further, they mean different things for different solvers. This map # just tracks the "nominal" argument that we should pass through pyomo. -solver_tol_map = { +SOLVER_TOL_MAP = { 'ipopt': 'tol', 'cbc': 'primalTolerance', 'glpk': 'mipgap', @@ -50,13 +33,14 @@ class Pyomo(Dispatcher): """ Dispatches using rolling windows in Pyomo """ - naming_template = {'comp prod': '{comp}|{res}|prod', - 'comp transfer': '{comp}|{res}|trans', - 'comp max': '{comp}|{res}|max', - 'comp ramp up': '{comp}|{res}|rampup', - 'comp ramp down': '{comp}|{res}|rampdown', - 'conservation': '{res}|consv', - } + # naming_template = { + # 'comp prod': '{comp}|{res}|prod', + # 'comp transfer': '{comp}|{res}|trans', + # 'comp max': '{comp}|{res}|max', + # 'comp ramp up': '{comp}|{res}|rampup', + # 'comp ramp down': '{comp}|{res}|rampdown', + # 'conservation': '{res}|consv', + # } ### INITIALIZATION @classmethod @@ -66,30 +50,54 @@ def get_input_specs(cls): @ In, None @ Out, specs, InputData, specs """ - specs = InputData.parameterInputFactory('pyomo', ordered=False, baseNode=None, - descr=r"""The \texttt{pyomo} dispatcher uses analytic modeling and rolling windows to - solve dispatch optimization with perfect information via the pyomo optimization library.""") - 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}.""")) - specs.addSub(InputData.parameterInputFactory('solver', contentType=InputTypes.StringType, - descr=r"""Indicates which solver should be used by pyomo. Options depend on individual installation. - \default{'glpk' for Windows, 'cbc' otherwise}.""")) - specs.addSub(InputData.parameterInputFactory('tol', contentType=InputTypes.FloatType, - descr=r"""Relative tolerance for converging final optimal dispatch solutions. Specific implementation - depends on the solver selected. Changing this value could have significant impacts on the - dispatch optimization time and quality. - \default{solver dependent, often 1e-6}.""")) + specs = InputData.parameterInputFactory( + 'pyomo', ordered=False, baseNode=None, + descr=r"""The \texttt{pyomo} dispatcher uses analytic modeling and rolling + windows to solve dispatch optimization with perfect information via the + pyomo optimization library.""" + ) + + 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}.""" + ) + ) + + specs.addSub( + InputData.parameterInputFactory( + 'solver', contentType=InputTypes.StringType, + descr=r"""Indicates which solver should be used by pyomo. Options depend + on individual installation. \default{'glpk' for Windows, 'cbc' otherwise}.""" + ) + ) + + specs.addSub( + InputData.parameterInputFactory( + 'tol', contentType=InputTypes.FloatType, + descr=r"""Relative tolerance for converging final optimal dispatch solutions. + Specific implementation depends on the solver selected. Changing this value + could have significant impacts on the dispatch optimization time and quality. + \default{solver dependent, often 1e-6}.""" + ) + ) # TODO specific for pyomo dispatcher return specs + def __init__(self): """ Constructor. @@ -104,7 +112,8 @@ def __init__(self): self._solver = None # overwrite option for solver self._picard_limit = 10 # iterative solve limit - def read_input(self, specs): + + def read_input(self, specs) -> None: """ Read in input specifications. @ In, specs, RAVEN InputData, specifications @@ -124,56 +133,31 @@ def read_input(self, specs): if solver_node is not None: self._solver = solver_node.value - # the tolerance value needs to be saved for after the solver is set - # since we don't know what key to assign it to otherwise tol_node = specs.findFirst('tol') if tol_node is not None: solver_tol = tol_node.value else: solver_tol = None - if self._solver is None: - solvers_to_check = SOLVERS - else: - solvers_to_check = [self._solver] - # check solver exists - for solver in solvers_to_check: - self._solver = solver - found_solver = True - try: - if not pyo.SolverFactory(self._solver).available(): - found_solver = False - else: - break - except ApplicationError: - found_solver = False - # NOTE: we probably need a consistent way to test and check viable solvers, - # maybe through a unit test that mimics the model setup here. For now, I assume - # that anything that shows as not available or starts with an underscore is not - # viable, and it will crash if the solver can't solve our kinds of models. - # This should only come up if the user is specifically requesting a solver, though, - # the default glpk and cbc are tested. - if not found_solver: - all_options = pyo.SolverFactory._cls.keys() # TODO shorten to list of tested options? - solver_filter = [] - for op in all_options: - if op.startswith('_'): # These don't seem like legitimate options, based on testing - solver_filter.append(False) - continue - try: - solver_filter.append(pyo.SolverFactory(op).available()) - except (ApplicationError, NameError, ImportError): - solver_filter.append(False) - available = list(compress(all_options, solver_filter)) - msg = f'Requested solver "{self._solver}" was not found for pyomo dispatcher!' - msg += f' Options MAY include: {available}' - raise RuntimeError(msg) + self._solver = putils.check_solver_availability(self._solver) if solver_tol is not None: - key = solver_tol_map[self._solver] - self.solve_options[key] = solver_tol + key = SOLVER_TOL_MAP.get(self._solver, None) + if key is not None: + self.solve_options[key] = solver_tol + else: + raise ValueError(f"Tolerance setting not available for solver '{self._solver}'.") + + + def get_solver(self): + """ + Retrieves the solver information (if applicable) + @ In, None + @ Out, solver, str, name of solver used + """ + return self._solver + - ### API def dispatch(self, case, components, sources, meta): """ Performs dispatch. @@ -184,77 +168,74 @@ def dispatch(self, case, components, sources, meta): @ Out, disp, DispatchScenario, resulting dispatch """ t_start, t_end, t_num = self.get_time_discr() - time = np.linspace(t_start, t_end, t_num) # Note we don't care about segment/cluster here - resources = sorted(list(hutils.get_all_resources(components))) # list of all active resources - # pre-build results structure - ## we can use NumpyState here so we don't need to worry about a Pyomo model object - dispatch = NumpyState()# dict((comp.name, dict((res, np.zeros(len(time))) for res in comp.get_resources())) for comp in components) + time = np.linspace(t_start, t_end, t_num) + resources = sorted(putils.get_all_resources(components)) + dispatch = NumpyState() dispatch.initialize(components, meta['HERON']['resource_indexer'], time) - # rolling window + start_index = 0 final_index = len(time) - # TODO window overlap! ( )[ ] -> ( [ ) ] + initial_levels = putils.get_initial_storage_levels(components, meta, start_index) + while start_index < final_index: - end_index = start_index + self._window_len - if end_index > final_index: - end_index = final_index + end_index = min(start_index + self._window_len, 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.") + raise DispatchError("Window length of 1 detected, which is not supported.") + specific_time = time[start_index:end_index] - print(f'DEBUGG starting window {start_index} to {end_index}') - start = time_mod.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]['level'][comp.get_interaction().get_resource()][-1] - # allow for converging solution iteratively - converged = False - conv_counter = 0 - previous = None - while not converged: - conv_counter += 1 - if conv_counter > self._picard_limit: - break - # dispatch - subdisp = self.dispatch_window(specific_time, start_index, - case, components, sources, resources, - initial_levels, meta) - # do we need a convergence criteria? Check now. - if self.needs_convergence(components): - print(f'DEBUGG iteratively solving window, iteration {conv_counter}/{self._picard_limit} ...') - converged = self.check_converged(subdisp, previous, components) - previous = subdisp - else: - converged = True + print(f"Start: {start_index} End: {end_index}") + subdisp, solve_time = self._handle_dispatch_window_solve( + specific_time, start_index, case, components, sources, resources, initial_levels, meta + ) + print(f'DEBUGG solve time: {solve_time} s') - end = time_mod.time() - print(f'DEBUGG solve time: {end-start} s') - # store result in corresponding part of dispatch + # Store Results of optimization into dispatch container for comp in components: for tag in comp.get_tracking_vars(): for res, values in subdisp[comp.name][tag].items(): dispatch.set_activity_vector(comp, res, values, tracker=tag, start_idx=start_index, end_idx=end_index) start_index = end_index + return dispatch - def get_solver(self): + + def _handle_dispatch_window_solve(self, specific_time, start_index, case, components, sources, resources, initial_levels, meta): """ - Retrieves the solver information (if applicable) - @ In, None - @ Out, solver, str, name of solver used + Set up convergence criteria and collect results from a dispatch window solve. + @ In, specific_time, np.array, value of time to evaluate. + @ In, start_index, int, index of the start of the window. + @ 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_levels, dict, initial storage levels if any. + @ In, meta, dict, additional variables passed through. + @ Out, subdisp, dict, results of window dispatch. """ - return self._solver + start = time_mod.time() + subdisp = self._dispatch_window(specific_time, start_index, case, components, resources, initial_levels, meta) + + if self._needs_convergence(components): + conv_counter = 0 + converged = False + previous = None + + while not converged and conv_counter < self._picard_limit: + conv_counter += 1 + print(f'DEBUGG iteratively solving window, iteration {conv_counter}/{self._picard_limit} ...') + subdisp = self._dispatch_window(specific_time, start_index, case, components, resources, initial_levels, meta) + converged = self._check_if_converged(subdisp, previous, components) + previous = subdisp - ### INTERNAL - def dispatch_window(self, time, time_offset, - case, components, sources, resources, - initial_storage, meta): + if conv_counter >= self._picard_limit and not converged: + raise DispatchError(f"Convergence not reached after {self._picard_limit} iterations.") + + end = time_mod.time() + solve_time = end - start + return subdisp, solve_time + + + def _dispatch_window(self, time, time_offset, case, components, resources, initial_storage, meta): """ Dispatches one part of a rolling window. @ In, time, np.array, value of time to evaluate @@ -267,991 +248,106 @@ def dispatch_window(self, time, time_offset, @ 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 - 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 - T = np.arange(0, len(time), dtype=int) # indexes resources - m.C = pyo.Set(initialize=C) - m.R = pyo.Set(initialize=R) - m.T = pyo.Set(initialize=T) - m.Times = time - dt = m.Times[1] - m.Times[0] # TODO assumes consistent step sizing - 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 - m.Case = case - m.Components = components - m.Activity = PyomoState() - m.Activity.initialize(m.Components, m.resource_index_map, m.Times, m) - # constraints and variables - for comp in components: - # components using a governing strategy (not opt) are Parameters, not Variables - # TODO should this come BEFORE or AFTER each dispatch opt solve? - # -> responsive or proactive? - intr = comp.get_interaction() - if intr.is_governed(): - meta['request'] = {'component': comp, 'time': time} - activity = intr.get_strategy().evaluate(meta)[0]['level'] - if intr.is_type('Storage'): - self._create_production_param(m, comp, activity, tag='level') - # set up "activity" rates (change in level over time, plus efficiency) - rte2 = comp.get_sqrt_RTE() # square root of the round-trip efficiency - L = len(activity) - deltas = np.zeros(L) - deltas[1:] = activity[1:] - activity[:-1] - deltas[0] = activity[0] - intr.get_initial_level(meta) - # rate of charge - # change sign, since increasing level means absorbing energy from system - # also scale by RTE, since to get level increase you have to over-absorb - charge = np.zeros(L) - charge_mask = np.where(deltas > 0) - charge[charge_mask] = - deltas[charge_mask] / dt / rte2 - self._create_production_param(m, comp, charge, tag='charge') - # rate of discharge - # change sign, since decreasing level means emitting energy into system - # also scale by RTE, since level decrease yields less to system - discharge = np.zeros(L) - discharge_mask = np.where(deltas < 0) - discharge[discharge_mask] = - deltas[discharge_mask] / dt * rte2 - self._create_production_param(m, comp, discharge, tag='discharge') - else: - self._create_production_param(m, comp, activity) - continue - # NOTE: "fixed" components could hypothetically be treated differently - ## however, in order for the "production" variable for components to be treatable as the - ## same as other production variables, we create components with limitation - ## lowerbound == upperbound == capacity (just for "fixed" dispatch components) - if intr.is_type('Storage'): - self._create_storage(m, comp, initial_storage, meta) - else: - self._create_production(m, comp, meta) # variables - 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 - if self.debug_mode: - self._debug_pyomo_print(m) - while not done_and_checked: - attempts += 1 - print(f'DEBUGG solve attempt {attempts} ...:') - # solve - soln = pyo.SolverFactory(self._solver).solve(m, options=self.solve_options) - # check solve status - if soln.solver.status == SolverStatus.ok and soln.solver.termination_condition == TerminationCondition.optimal: - print('DEBUGG ... solve was successful!') - else: - print('DEBUGG ... solve was unsuccessful!') - print('DEBUGG ... status:', soln.solver.status) - print('DEBUGG ... termination:', soln.solver.termination_condition) - self._debug_pyomo_print(m) - print('Resource Map:') - pprint.pprint(m.resource_index_map) - raise RuntimeError(f"Solve was unsuccessful! Status: {soln.solver.status} Termination: {soln.solver.termination_condition}") - # try validating - print('DEBUGG ... validating ...') - validation_errs = self.validate(m.Components, m.Activity, m.Times, meta) - if validation_errs: - done_and_checked = False - print('DEBUGG ... validation concerns raised:') - for e in validation_errs: - print(f"DEBUGG ... ... Time {e['time_index']} ({e['time']}) \n" + - f"Component \"{e['component'].name}\" Resource \"{e['resource']}\": {e['msg']}") - self._create_production_limit(m, e) - # go back and solve again - # raise NotImplementedError('Validation failed, but idk how to handle that yet') - else: - print('DEBUGG Solve successful and no validation concerns raised.') - done_and_checked = True - if attempts > 100: - raise RuntimeError('Exceeded validation attempt limit!') - if self.debug_mode: - soln.write() - self._debug_print_soln(m) - # return dict of numpy arrays - result = self._retrieve_solution(m) + model = PyomoModelHandler(time, time_offset, case, components, resources, initial_storage, meta) + model.populate_model() + result = self._solve_dispatch(model, meta) return result - def check_converged(self, new, old, components): + + def _check_if_converged(self, new, old, components, tol=1e-4): """ Checks convergence of consecutive dispatch solves @ In, new, dict, results of dispatch # TODO should this be the model rather than dict? @ In, old, dict, results of previous dispatch @ In, components, list, HERON component list + @ In, tol, float, optional, tolerance for convergence @ Out, converged, bool, True if convergence is met """ - tol = 1e-4 # TODO user option if old is None: return False - converged = True + for comp in components: intr = comp.get_interaction() - name = comp.name - tracker = comp.get_tracking_vars()[0] if intr.is_governed(): # by "is_governed" we mean "isn't optimized in pyomo" # check activity L2 norm as a differ # TODO this may be specific to storage right now + name = comp.name + tracker = comp.get_tracking_vars()[0] res = intr.get_resource() scale = np.max(old[name][tracker][res]) - diff = np.linalg.norm(new[name][tracker][res] - old[name][tracker][res]) / (scale if scale != 0 else 1) + # Avoid division by zero + if scale == 0: + diff = np.linalg.norm(new[name][tracker][res] - old[name][tracker][res]) + else: + diff = np.linalg.norm(new[name][tracker][res] - old[name][tracker][res]) / scale if diff > tol: - converged = False - return converged + return False + return True + - def needs_convergence(self, components): + def _needs_convergence(self, components): """ Determines whether the current setup needs convergence to solve. @ In, components, list, HERON component list @ Out, needs_convergence, bool, True if iteration is needed """ - for comp in components: - intr = comp.get_interaction() - # storages with a prescribed strategy MAY need iteration - if intr.is_governed(): - return True - # if we get here, no iteration is needed - return False - - - ### PYOMO Element Constructors - def _create_production_limit(self, m, validation): - """ - Creates pyomo production constraint given validation errors - @ In, m, pyo.ConcreteModel, associated model - @ In, validation, dict, information from Validator about limit violation - @ Out, None - """ - # TODO could validator write a symbolic expression on request? That'd be sweet. - comp = validation['component'] - resource = validation['resource'] - r = m.resource_index_map[comp][resource] - t = validation['time_index'] - limit = validation['limit'] - limits = np.zeros(len(m.Times)) - limits[t] = limit - limit_type = validation['limit_type'] - prod_name = f'{comp.name}_production' - rule = lambda mod: self._prod_limit_rule(prod_name, r, limits, limit_type, t, mod) - constr = pyo.Constraint(rule=rule) - counter = 1 - name_template = f'{comp.name}_{resource}_{t}_vld_limit_constr_{{i}}' - # make sure we get a unique name for this constraint - name = name_template.format(i=counter) - while getattr(m, name, None) is not None: - counter += 1 - name = name_template.format(i=counter) - setattr(m, name, constr) - print(f'DEBUGG added validation constraint "{name}"') - - def _create_production_param(self, m, comp, values, tag=None): - """ - Creates production pyomo fixed parameter object for a component - @ In, m, pyo.ConcreteModel, associated model - @ In, comp, HERON Component, component to make production variables for - @ In, values, np.array(float), values to set for param - @ In, tag, str, optional, if not None then name will be component_[tag] - @ Out, prod_name, str, name of production variable - """ - name = comp.name - if tag is None: - tag = 'production' - # create pyomo indexer for this component's resources - res_indexer = pyo.Set(initialize=range(len(m.resource_index_map[comp]))) - setattr(m, f'{name}_res_index_map', res_indexer) - prod_name = f'{name}_{tag}' - init = (((0, t), values[t]) for t in m.T) - prod = pyo.Param(res_indexer, m.T, initialize=dict(init)) - setattr(m, prod_name, prod) - return prod_name - - def _create_production(self, m, comp, meta): - """ - Creates all pyomo variable objects for a non-storage component - @ In, m, pyo.ConcreteModel, associated model - @ In, comp, HERON Component, component to make production variables for - @ In, meta, dict, dictionary of state variables - @ Out, None - """ - prod_name = self._create_production_variable(m, comp, meta) - ## if you cannot set limits directly in the production variable, set separate contraint: - ## Method 1: set variable bounds directly --> TODO more work needed, but would be nice - # lower, upper = self._get_prod_bounds(m, comp) - # limits should be None unless specified, so use "getters" from dictionaries - # bounds = lambda m, r, t: (lower.get(r, None), upper.get(r, None)) - ## Method 2: set variable bounds directly --> TODO more work needed, but would be nice - # self._create_capacity(m, comp, prod_name, meta) # capacity constraints - # transfer function governs input -> output relationship - self._create_transfer(m, comp, prod_name) - # ramp rates - if comp.ramp_limit is not None: - self._create_ramp_limit(m, comp, prod_name, meta) - return prod_name - - def _create_production_variable(self, m, comp, meta, tag=None, add_bounds=True, **kwargs): - """ - Creates production pyomo variable object for a component - @ In, m, pyo.ConcreteModel, associated model - @ In, comp, HERON Component, component to make production variables for - @ In, tag, str, optional, if not None then name will be component_[tag]; otherwise "production" - @ In, add_bounds, bool, optional, if True then determine and set bounds for variable - @ In, kwargs, dict, optional, passalong kwargs to pyomo variable - @ Out, prod_name, str, name of production variable - """ - if tag is None: - tag = 'production' - name = comp.name - cap_res = comp.get_capacity_var() # name of resource that defines capacity - limit_r = m.resource_index_map[comp][cap_res] # production index of the governing resource - # create pyomo indexer for this component's resources - indexer_name = f'{name}_res_index_map' - indexer = getattr(m, indexer_name, None) - if indexer is None: - indexer = pyo.Set(initialize=range(len(m.resource_index_map[comp]))) - setattr(m, indexer_name, indexer) - prod_name = f'{name}_{tag}' - caps, mins = self._find_production_limits(m, comp, meta) - if min(caps) < 0: - # quick check that capacities signs are consistent #FIXME: revisit, this is an assumption - assert max(caps) <= 0, \ - 'Capacities are inconsistent: mix of positive and negative values not currently supported.' - # we have a unit that's consuming, so we need to flip the variables to be sensible - mins, caps = caps, mins - inits = caps - else: - inits = mins - if add_bounds: - # create bounds based in min, max operation - bounds = lambda m, r, t: (mins[t] if r == limit_r else None, caps[t] if r == limit_r else None) - initial = lambda m, r, t: inits[t] if r == limit_r else 0 - else: - bounds = (None, None) - initial = 0 - # production variable depends on resources, time - #FIXME initials! Should be lambda with mins for tracking var! - prod = pyo.Var(indexer, m.T, initialize=initial, bounds=bounds, **kwargs) - # TODO it may be that we need to set variable values to avoid problems in some solvers. - # if comp.is_dispatchable() == 'fixed': - # for t, _ in enumerate(m.Times): - # prod[limit_r, t].fix(caps[t]) - setattr(m, prod_name, prod) - return prod_name - - def _create_ramp_limit(self, m, comp, prod_name, meta): - """ - Creates ramping limitations for a producing component - @ In, m, pyo.ConcreteModel, associated model - @ In, comp, HERON Component, component to make ramping limits for - @ In, prod_name, str, name of production variable - @ In, meta, dict, dictionary of state variables - @ Out, None - """ - # ramping is defined in terms of the capacity variable - cap_res = comp.get_capacity_var() # name of resource that defines capacity - cap = comp.get_capacity(meta)[0][cap_res] - r = m.resource_index_map[comp][cap_res] # production index of the governing resource - # NOTE: this includes the built capacity * capacity factor, if any, which assumes - # the ramp rate depends on the available capacity, not the built capacity. - limit_delta = comp.ramp_limit * cap # NOTE: if cap is negative, then this is negative. - if limit_delta < 0: - neg_cap = True - else: - neg_cap = False - # if we're limiting ramp frequency, make vars and rules for that - if comp.ramp_freq: - # create binaries for tracking ramping - up = pyo.Var(m.T, initialize=0, domain=pyo.Binary) - down = pyo.Var(m.T, initialize=0, domain=pyo.Binary) - steady = pyo.Var(m.T, initialize=1, domain=pyo.Binary) - setattr(m, f'{comp.name}_up_ramp_tracker', up) - setattr(m, f'{comp.name}_down_ramp_tracker', down) - setattr(m, f'{comp.name}_steady_ramp_tracker', steady) - ramp_trackers = (down, up, steady) - else: - ramp_trackers = None - # limit production changes when ramping down - ramp_rule_down = lambda mod, t: \ - self._ramp_rule_down(prod_name, r, limit_delta, neg_cap, t, mod, - bins=ramp_trackers) - constr = pyo.Constraint(m.T, rule=ramp_rule_down) - setattr(m, f'{comp.name}_ramp_down_constr', constr) - # limit production changes when ramping up - ramp_rule_up = lambda mod, t: \ - self._ramp_rule_up(prod_name, r, limit_delta, neg_cap, t, mod, - bins=ramp_trackers) - constr = pyo.Constraint(m.T, rule=ramp_rule_up) - setattr(m, f'{comp.name}_ramp_up_constr', constr) - # if ramping frequency limit, impose binary constraints - if comp.ramp_freq: - # binaries rule, for exclusive choice up/down/steady - binaries_rule = lambda mod, t: \ - self._ramp_freq_bins_rule(down, up, steady, t, mod) - constr = pyo.Constraint(m.T, rule=binaries_rule) - setattr(m, f'{comp.name}_ramp_freq_binaries', constr) - # limit frequency of ramping - # TODO calculate "tao" window using ramp freq and dt - # -> for now, just use the integer for number of windows - freq_rule = lambda mod, t: self._ramp_freq_rule(down, up, comp.ramp_freq, t, m) - constr = pyo.Constraint(m.T, rule=freq_rule) - setattr(m, f'{comp.name}_ramp_freq_constr', constr) - - def _create_capacity_constraints(self, m, comp, prod_name, meta): - """ - Creates pyomo capacity constraints - @ In, m, pyo.ConcreteModel, associated model - @ In, comp, HERON Component, component to make variables for - @ In, prod_name, str, name of production variable - @ In, meta, dict, additional state information - @ Out, None - """ - cap_res = comp.get_capacity_var() # name of resource that defines capacity - r = m.resource_index_map[comp][cap_res] # production index of the governing resource - caps, mins = self._find_production_limits(m, comp, meta) - # capacity - max_rule = lambda mod, t: self._capacity_rule(prod_name, r, caps, mod, t) - constr = pyo.Constraint(m.T, rule=max_rule) - setattr(m, f'{comp.name}_{cap_res}_capacity_constr', constr) - # minimum - min_rule = lambda mod, t: self._min_prod_rule(prod_name, r, caps, mins, mod, t) - constr = pyo.Constraint(m.T, rule=min_rule) - # set initial conditions - for t, time in enumerate(m.Times): - cap = caps[t] - if cap == mins[t]: - # initialize values so there's no boundary errors - var = getattr(m, prod_name) - values = var.get_values() - for k in values: - values[k] = cap - var.set_values(values) - setattr(m, f'{comp.name}_{cap_res}_minprod_constr', constr) - - def _find_production_limits(self, m, comp, meta): - """ - Determines the capacity limits of a unit's operation, in time. - @ In, m, pyo.ConcreteModel, associated model - @ In, comp, HERON Component, component to make variables for - @ In, meta, dict, additional state information - @ Out, caps, array, max production values by time - @ Out, mins, array, min production values by time - """ - cap_res = comp.get_capacity_var() # name of resource that defines capacity - r = m.resource_index_map[comp][cap_res] # production index of the governing resource - # production is always lower than capacity - ## NOTE get_capacity returns (data, meta) and data is dict - ## TODO does this work with, e.g., ARMA-based capacities? - ### -> "time" is stored on "m" and could be used to correctly evaluate the capacity - caps = [] - mins = [] - for t, time in enumerate(m.Times): - 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) - if (comp.is_dispatchable() == 'fixed'): - minimum = cap - else: - minimum = comp.get_minimum(meta)[0][cap_res] - mins.append(minimum) - return caps, mins - - def _create_transfer(self, m, comp, prod_name): - """ - Creates pyomo transfer function constraints - @ In, m, pyo.ConcreteModel, associated model - @ In, comp, HERON Component, component to make variables for - @ In, prod_name, str, name of production variable - @ Out, None - """ - name = comp.name - # transfer functions - # e.g. 2A + 3B -> 1C + 2E - # get linear coefficients - # TODO this could also take a transfer function from an external Python function assuming - # we're careful about how the expression-vs-float gets used - # and figure out how to handle multiple ins, multiple outs - ratios = self._get_transfer_coeffs(m, comp) - ref_r, ref_name, _ = ratios.pop('__reference', (None, None, None)) - for resource, ratio in ratios.items(): - r = m.resource_index_map[comp][resource] - rule_name = f'{name}_{resource}_{ref_name}_transfer' - rule = lambda mod, t: self._transfer_rule(ratio, r, ref_r, prod_name, mod, t) - constr = pyo.Constraint(m.T, rule=rule) - setattr(m, rule_name, constr) + return any(comp.get_interaction().is_governed() for comp in components) - def _create_storage(self, m, comp, initial_storage, meta): - """ - Creates storage pyomo variable objects for a storage component - Similar to create_production, but for storages - @ In, m, pyo.ConcreteModel, associated model - @ In, comp, HERON Component, component to make production variables for - @ In, initial_storage, dict, initial storage levels - @ In, meta, dict, additional state information - @ Out, level_name, str, name of storage level variable - """ - prefix = comp.name - # what resource index? Isn't it always 0? # assumption - r = 0 # NOTE this is only true if each storage ONLY uses 1 resource - # storages require a few variables: - # (1) a level tracker, - level_name = self._create_production_variable(m, comp, meta, tag='level') - # -> set operational limits - # self._create_capacity(m, comp, level_name, meta) - # (2, 3) separate charge/discharge trackers, so we can implement round-trip efficiency and ramp rates - charge_name = self._create_production_variable(m, comp, meta, tag='charge', add_bounds=False, within=pyo.NonPositiveReals) - discharge_name = self._create_production_variable(m, comp, meta, tag='discharge', add_bounds=False, within=pyo.NonNegativeReals) - # balance level, charge/discharge - level_rule_name = prefix + '_level_constr' - rule = lambda mod, t: self._level_rule(comp, level_name, charge_name, discharge_name, - initial_storage, r, mod, t) - setattr(m, level_rule_name, pyo.Constraint(m.T, rule=rule)) - # periodic boundary condition for storage level - if comp.get_interaction().apply_periodic_level: - periodic_rule_name = prefix + '_level_periodic_constr' - rule = lambda mod, t: self._periodic_level_rule(comp, level_name, initial_storage, r, mod, t) - setattr(m, periodic_rule_name, pyo.Constraint(m.T, rule=rule)) - - # (4) a binary variable to track whether we're charging or discharging, to prevent BOTH happening - # -> 0 is charging, 1 is discharging - # -> TODO make this a user-based option to disable, if they want to allow dual operation - # -> -> but they should really think about if that's what they want! - # FIXME currently introducing the bigM strategy also makes solves numerically unstable, - # and frequently results in spurious errors. For now, disable it. - allow_both = True # allow simultaneous charging and discharging - if not allow_both: - bin_name = self._create_production_variable(m, comp, meta, tag='dcforcer', add_bounds=False, within=pyo.Binary) - # we need a large epsilon, but not so large that addition stops making sense - # -> we don't know what any values for this component will be! How do we choose? - # -> NOTE that choosing this value has VAST impact on solve stability!! - large_eps = 1e8 #0.01 * sys.float_info.max - # charging constraint: don't charge while discharging (note the sign matters) - charge_rule_name = prefix + '_charge_constr' - rule = lambda mod, t: self._charge_rule(charge_name, bin_name, large_eps, r, mod, t) - setattr(m, charge_rule_name, pyo.Constraint(m.T, rule=rule)) - discharge_rule_name = prefix + '_discharge_constr' - rule = lambda mod, t: self._discharge_rule(discharge_name, bin_name, large_eps, r, mod, t) - setattr(m, discharge_rule_name, pyo.Constraint(m.T, rule=rule)) - - 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 resource in resources: - rule = lambda mod, t: self._conservation_rule(initial_storage, meta, resource, mod, t) - constr = pyo.Constraint(m.T, rule=rule) - setattr(m, f'{resource}_conservation', constr) - - def _create_objective(self, meta, m): - """ - Creates pyomo objective function - @ In, meta, dict, additional variables to pass through - @ In, m, pyo.ConcreteModel, associated model - @ Out, None - """ - # cashflow eval - rule = lambda mod: self._cashflow_rule(meta, mod) - m.obj = pyo.Objective(rule=rule, sense=pyo.maximize) - - ### UTILITIES for general use - def _get_prod_bounds(self, m, comp, meta): - """ - Determines the production limits of the given component - @ In, comp, HERON component, component to get bounds of - @ In, meta, dict, additional state information - @ Out, lower, dict, lower production limit (might be negative for consumption) - @ Out, upper, dict, upper production limit (might be 0 for consumption) - """ - raise NotImplementedError - cap_res = comp.get_capacity_var() # name of resource that defines capacity - r = m.resource_index_map[comp][cap_res] - maxs = [] - mins = [] - for t, time in enumerate(m.Times): - meta['HERON']['time_index'] = t + m.time_offset - cap = comp.get_capacity(meta)[0][cap_res] - low = comp.get_minimum(meta)[0][cap_res] - maxs.append(cap) - if (comp.is_dispatchable() == 'fixed') or (low == cap): - low = cap - # initialize values to avoid boundary errors - var = getattr(m, prod_name) - values = var.get_values() - for k in values: - values[k] = cap - var.set_values(values) - mins.append(low) - maximum = comp.get_capacity(None, None, None, None)[0][cap_res] - # TODO minimum! - # producing or consuming the defining resource? - # -> put these in dictionaries so we can "get" limits or return None - if maximum > 0: - lower = {r: 0} - upper = {r: maximum} - else: - lower = {r: maximum} - upper = {r: 0} - return lower, upper - - def _get_transfer_coeffs(self, m, comp): - """ - Obtains transfer function ratios (assuming Linear ValuedParams) - Form: 1A + 3B -> 2C + 4D - Ratios are calculated with respect to first resource listed, so e.g. B = 3/1 * A - TODO I think we can handle general external functions, maybe? - @ In, m, pyo.ConcreteModel, associated model - @ In, comp, HERON component, component to get coefficients of - @ Out, ratios, dict, ratios of transfer function variables - """ - name = comp.name - transfer = comp.get_interaction().get_transfer() # get the transfer ValuedParam, if any - if transfer is None: - return {} - coeffs = transfer.get_coefficients() # linear transfer coefficients, dict as {resource: coeff}, SIGNS MATTER - # it's all about ratios -> store as ratio of resource / first resource (arbitrary) - first_r = None - first_name = None - first_coef = None - ratios = {} - for resource, coef in coeffs.items(): - r = m.resource_index_map[comp][resource] - if first_r is None: - # set up nominal resource to compare to - first_r = r - first_name = resource - first_coef = coef - ratios['__reference'] = (first_r, first_name, first_coef) - continue - ratio = coef / first_coef - ratios[resource] = ratio - return ratios - - def _retrieve_solution(self, m): - """ - Extracts solution from Pyomo optimization - @ In, m, pyo.ConcreteModel, associated (solved) model - @ Out, result, dict, {comp: {activity: {resource: [production]}} e.g. generator[production][steam] - """ - result = {} # {component: {activity_tag: {resource: production}}} - for comp in m.Components: - result[comp.name] = {} - for tag in comp.get_tracking_vars(): - tag_values = self._retrieve_value_from_model(m, comp, tag) - result[comp.name][tag] = tag_values - return result - - def _retrieve_value_from_model(self, m, comp, tag): - """ - Retrieve values of a series from the pyomo model. - @ In, m, pyo.ConcreteModel, associated (solved) model - @ In, comp, Component, relevant component - @ In, tag, str, particular history type to retrieve - @ Out, result, dict, {resource: [array], etc} - """ - result = {} - prod = getattr(m, f'{comp.name}_{tag}') - if isinstance(prod, pyo.Var): - kind = 'Var' - elif isinstance(prod, pyo.Param): - kind = 'Param' - for res, comp_r in m.resource_index_map[comp].items(): - if kind == 'Var': - result[res] = np.fromiter((prod[comp_r, t].value for t in m.T), dtype=float, count=len(m.T)) - elif kind == 'Param': - result[res] = np.fromiter((prod[comp_r, t] for t in m.T), dtype=float, count=len(m.T)) - return result - - ### RULES for lambda function calls - # these get called using lambda functions to make Pyomo constraints, vars, objectives, etc - - def _charge_rule(self, charge_name, bin_name, large_eps, r, m, t): - """ - Constructs pyomo don't-charge-while-discharging constraints. - For storage units specificially. - @ In, charge_name, str, name of charging variable - @ In, bin_name, str, name of forcing binary variable - @ In, large_eps, float, a large-ish number w.r.t. storage activity - @ In, r, int, index of stored resource (is this always 0?) - @ In, m, pyo.ConcreteModel, associated model - @ In, t, int, time index for capacity rule - @ Out, rule, bool, inequality used to limit charge behavior - """ - charge_var = getattr(m, charge_name) - bin_var = getattr(m, bin_name) - return -charge_var[r, t] <= (1 - bin_var[r, t]) * large_eps - - def _discharge_rule(self, discharge_name, bin_name, large_eps, r, m, t): - """ - Constructs pyomo don't-discharge-while-charging constraints. - For storage units specificially. - @ In, discharge_name, str, name of discharging variable - @ In, bin_name, str, name of forcing binary variable - @ In, large_eps, float, a large-ish number w.r.t. storage activity - @ In, r, int, index of stored resource (is this always 0?) - @ In, m, pyo.ConcreteModel, associated model - @ In, t, int, time index for capacity rule - @ Out, rule, bool, inequality used to limit discharge behavior - """ - discharge_var = getattr(m, discharge_name) - bin_var = getattr(m, bin_name) - return discharge_var[r, t] <= bin_var[r, t] * large_eps - - def _level_rule(self, comp, level_name, charge_name, discharge_name, initial_storage, r, m, t): - """ - Constructs pyomo charge-discharge-level balance constraints. - For storage units specificially. - @ In, comp, Component, storage component of interest - @ In, level_name, str, name of level-tracking variable - @ In, charge_name, str, name of charging variable - @ In, discharge_name, str, name of discharging variable - @ In, initial_storage, dict, initial storage levels by component - @ In, r, int, index of stored resource (is this always 0?) - @ In, m, pyo.ConcreteModel, associated model - @ In, t, int, time index for capacity rule - @ Out, rule, bool, inequality used to limit level behavior - """ - level_var = getattr(m, level_name) - charge_var = getattr(m, charge_name) - discharge_var = getattr(m, discharge_name) - if t > 0: - previous = level_var[r, t-1] - dt = m.Times[t] - m.Times[t-1] - else: - previous = initial_storage[comp] - dt = m.Times[1] - m.Times[0] - rte2 = comp.get_sqrt_RTE() # square root of the round-trip efficiency - production = - rte2 * charge_var[r, t] - discharge_var[r, t] / rte2 - return level_var[r, t] == previous + production * dt - - def _periodic_level_rule(self, comp, level_name, initial_storage, r, m, t): - """ - Mandates storage units end with the same level they start with, which prevents - "free energy" or "free sink" due to initial starting levels. - For storage units specificially. - @ In, comp, Component, storage component of interest - @ In, level_name, str, name of level-tracking variable - @ In, initial_storage, dict, initial storage levels by component - @ In, r, int, index of stored resource (is this always 0?) - @ In, m, pyo.ConcreteModel, associated model - @ In, t, int, time index for capacity rule - @ Out, rule, bool, inequality used to limit level behavior - """ - return getattr(m, level_name)[r, m.T[-1]] == initial_storage[comp] - def _capacity_rule(self, prod_name, r, caps, m, t): + def _solve_dispatch(self, m, meta): """ - Constructs pyomo capacity constraints. - @ In, prod_name, str, name of production variable - @ In, r, int, index of resource for capacity constraining - @ In, caps, list(float), value to constrain resource at in time - @ In, m, pyo.ConcreteModel, associated model - @ In, t, int, time index for capacity rule + Solves the dispatch problem. + @ In, m, PyomoModelHandler, model object to solve + @ In, meta, dict, additional variables passed through + @ Out, result, dict, results of solved dispatch """ - kind = 'lower' if min(caps) < 0 else 'upper' - return self._prod_limit_rule(prod_name, r, caps, kind, t, m) + # start a solution search + done_and_checked = False + attempts = 0 + # DEBUGG show variables, bounds + if self.debug_mode: + putils.debug_pyomo_print(m.model) - def _prod_limit_rule(self, prod_name, r, limits, kind, t, m): - """ - Constructs pyomo production constraints. - @ In, prod_name, str, name of production variable - @ In, r, int, index of resource for capacity constraining - @ In, limits, list(float), values in time at which to constrain resource production - @ In, kind, str, either 'upper' or 'lower' for limiting production - @ In, t, int, time index for production rule (NOTE not pyomo index, rather fixed index) - @ In, m, pyo.ConcreteModel, associated model - """ - prod = getattr(m, prod_name) - if kind == 'lower': - # production must exceed value - return prod[r, t] >= limits[t] - elif kind == 'upper': - return prod[r, t] <= limits[t] - else: - raise TypeError('Unrecognized production limit "kind":', kind) + while not done_and_checked: + attempts += 1 + print(f'DEBUGG solve attempt {attempts} ...:') + soln = pyo.SolverFactory(self._solver).solve(m.model, options=self.solve_options) - def _ramp_rule_down(self, prod_name, r, limit, neg_cap, t, m, bins=None): - """ - Constructs pyomo production ramping constraints for reducing production level. - Note that this is number-getting-less-positive for positive-defined capacity, while - it is number-getting-less-negative for negative-defined capacity. - This means that dQ is negative for positive-defined capacity, but positive for vice versa - @ In, prod_name, str, name of production variable - @ In, r, int, index of resource for capacity constraining - @ In, limit, float, limiting change in production level across time steps. NOTE: negative for negative-defined capacity. - @ In, neg_cap, bool, True if capacity is expressed as negative (consumer) - @ In, t, int, time index for ramp limit rule (NOTE not pyomo index, rather fixed index) - @ In, m, pyo.ConcreteModel, associated model - @ In, bins, tuple, optional, (lower, upper, steady) binaries if limiting ramp frequency - @ Out, rule, expression, evaluation for Pyomo constraint - """ - prod = getattr(m, prod_name) - if t == 0: - return pyo.Constraint.Skip - delta = prod[r, t] - prod[r, t-1] - # special treatment if we have frequency-limiting binaries available - if bins is None: - if neg_cap: - # NOTE change in production should be "less positive" than the max - # "negative decrease" in production (decrease is positive when defined by consuming) - return delta <= - limit - else: - # dq is negative, - limit is negative - return delta >= - limit - else: - eps = 1.0 # aux parameter to force binaries to behave, TODO needed? - down = bins[0][t] - up = bins[1][t] - # NOTE we're following the convention that "less negative" is ramping "down" - # for capacity defined by consumption - # e.g. consuming 100 ramps down to consuming 70 is (-100 -> -70), dq = 30 - if neg_cap: - # dq <= limit * dt * Bu + eps * Bd, if limit <= 0 - # dq is positive, - limit is positive - return delta <= - limit * down - eps * up + # check solve status + if soln.solver.status == SolverStatus.ok and soln.solver.termination_condition == TerminationCondition.optimal: + print('DEBUGG ... solve was successful!') else: - # dq <= limit * dt * Bu - eps * Bd, if limit >= 0 - # dq is negative, - limit is negative - return delta >= - limit * down + eps * up + putils.debug_pyomo_print(m.model) + print('Resource Map:') + pprint.pprint(m.model.resource_index_map) + raise DispatchError( + f"Solve was unsuccessful! Status: {soln.solver.status} Termination: {soln.solver.termination_condition}" + ) - def _ramp_rule_up(self, prod_name, r, limit, neg_cap, t, m, bins=None): - """ - Constructs pyomo production ramping constraints. - @ In, prod_name, str, name of production variable - @ In, r, int, index of resource for capacity constraining - @ In, limit, float, limiting change in production level across time steps - @ In, neg_cap, bool, True if capacity is expressed as negative (consumer) - @ In, t, int, time index for ramp limit rule (NOTE not pyomo index, rather fixed index) - @ In, m, pyo.ConcreteModel, associated model - @ In, bins, tuple, optional, (lower, steady, upper) binaries if limiting ramp frequency - """ - prod = getattr(m, prod_name) - if t == 0: - return pyo.Constraint.Skip - delta = prod[r, t] - prod[r, t-1] - if bins is None: - if neg_cap: - # NOTE change in production should be "more positive" than the max - # "negative increase" in production (increase is negative when defined by consuming) - return delta >= limit - else: - # change in production should be less than the max production increase - return delta <= limit - else: - # special treatment if we have frequency-limiting binaries available - eps = 1.0 # aux parameter to force binaries to behave, TODO needed? - down = bins[0][t] - up = bins[1][t] - # NOTE we're following the convention that "more negative" is ramping "up" - # for capacity defined by consumption - # e.g. consuming 100 ramps up to consuming 130 is (-100 -> -130), dq = -30 - if neg_cap: - # dq >= limit * dt * Bu + eps * Bd, if limit <= 0 - return delta >= limit * up + eps * down + # try validating + print('DEBUGG ... validating ...') + validation_errs = self.validate(m.model.Components, m.model.Activity, m.model.Times, meta) + if validation_errs: + done_and_checked = False + print('DEBUGG ... validation concerns raised:') + for e in validation_errs: + print(f"DEBUGG ... ... Time {e['time_index']} ({e['time']}) \n" + + f"Component \"{e['component'].name}\" Resource \"{e['resource']}\": {e['msg']}") + m._create_production_limit(e) + # TODO go back and solve again + # raise DispatchError('Validation failed, but idk how to handle that yet') else: - # dq <= limit * dt * Bu - eps * Bd, if limit >= 0 - return delta <= limit * up - eps * down - - def _ramp_freq_rule(self, Bd, Bu, tao, t, m): - """ - Constructs pyomo frequency-of-ramp constraints. - @ In, Bd, bool var, binary tracking down-ramp events - @ In, Bu, bool var, binary tracking up-ramp events - @ In, tao, int, number of time steps to look back - @ In, t, int, time step indexer - @ In, m, pyo.ConcreteModel, pyomo model - """ - if t == 0: - return pyo.Constraint.Skip - # looking-back-window shouldn't be longer than existing time - tao = min(t, tao) - # how many ramp-down events in backward window? - tally = sum(1 - Bd[tm] for tm in range(t - tao, t)) - # only allow ramping up if no rampdowns in back window - ## but we can't use if statements, so use binary math - return Bu[t] <= 1 / tao * tally - - def _ramp_freq_bins_rule(self, Bd, Bu, Bn, t, m): - """ - Constructs pyomo constraint for ramping event tracking variables. - This forces choosing between ramping up, ramping down, and steady state operation. - @ In, Bd, bool var, binary tracking down-ramp events - @ In, Bu, bool var, binary tracking up-ramp events - @ In, Bn, bool var, binary tracking no-ramp events - @ In, t, int, time step indexer - @ In, m, pyo.ConcreteModel, pyomo model - """ - return Bd[t] + Bu[t] + Bn[t] == 1 - - def _cashflow_rule(self, meta, m): - """ - Objective function rule. - @ In, meta, dict, additional variable passthrough - @ In, m, pyo.ConcreteModel, associated model - @ Out, total, float, evaluation of cost - """ - 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, time_offset=m.time_offset) - return total - - 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 - @ In, t, int, index of time variable - @ Out, conservation, bool, balance check - """ - balance = 0 # sum of production rates, which needs to be zero - for comp, res_dict in m.resource_index_map.items(): - if res in res_dict: - # activity information depends on if storage or component - r = res_dict[res] - intr = comp.get_interaction() - if intr.is_type('Storage'): - # Storages have 3 variables: level, charge, and discharge - # -> so calculate activity - charge = getattr(m, f'{comp.name}_charge') - discharge = getattr(m, f'{comp.name}_discharge') - # note that "charge" is negative (as it's consuming) and discharge is positive - # -> so the intuitive |discharge| - |charge| becomes discharge + charge - production = discharge[r, t] + charge[r, t] - else: - var = getattr(m, f'{comp.name}_production') - # TODO move to this? balance += m._activity.get_activity(comp, res, t) - production = var[r, t] - balance += production - return balance == 0 # TODO tol? - - def _min_prod_rule(self, prod_name, r, caps, minimums, m, t): - """ - Constructs minimum production constraint - @ In, prod_name, str, name of production variable - @ In, r, int, index of resource for capacity constraining - @ In, caps, list(float), capacity(t) value for component - @ In, minimums, list(float), minimum allowable production in time - @ In, m, pyo.ConcreteModel, associated model - @ In, t, int, index of time variable - @ Out, minimum, bool, min check - """ - prod = getattr(m, prod_name) - # negative capacity means consuming instead of producing - if max(caps) > 0: - return prod[r, t] >= minimums[t] - else: - return prod[r, t] <= minimums[t] - - def _transfer_rule(self, ratio, r, ref_r, prod_name, m, t): - """ - Constructs transfer function constraints - @ In, ratio, float, ratio for resource to nominal first resource - @ In, r, int, index of transfer resource - @ In, ref_r, int, index of reference resource - @ In, prod_name, str, name of production variable - @ In, m, pyo.ConcreteModel, associated model - @ In, t, int, index of time variable - @ Out, transfer, bool, transfer ratio check - """ - prod = getattr(m, prod_name) - return prod[r, t] == prod[ref_r, t] * ratio # TODO tolerance?? - - ### DEBUG - def _debug_pyomo_print(self, m): - """ - Prints the setup pieces of the Pyomo model - @ In, m, pyo.ConcreteModel, model to interrogate - @ Out, None - """ - print('/' + '='*80) - print('DEBUGG model pieces:') - print(' -> objective:') - print(' ', m.obj.pprint()) - print(' -> variables:') - for var in m.component_objects(pyo.Var): - print(' ', var.pprint()) - print(' -> constraints:') - for constr in m.component_objects(pyo.Constraint): - print(' ', constr.pprint()) - print('\\' + '='*80) - print('') - - def _debug_print_soln(self, m): - """ - Prints the solution from the Pyomo model - @ In, m, pyo.ConcreteModel, model to interrogate - @ Out, None - """ - print('*'*80) - print('DEBUGG solution:') - print(' objective value:', m.obj()) - for c, comp in enumerate(m.Components): - name = comp.name - print(' component:', c, name) - for tracker in comp.get_tracking_vars(): - for res, r in m.resource_index_map[comp].items(): - print(f' tracker: {tracker} resource {r}: {res}') - for t, time in enumerate(m.Times): - prod = getattr(m, f'{name}_{tracker}') - if isinstance(prod, pyo.Var): - print(' time:', t + m.time_offset, time, prod[r, t].value) - elif isinstance(prod, pyo.Param): - print(' time:', t + m.time_offset, time, prod[r, t]) - print('*'*80) - - - - - -# DispatchState for Pyomo dispatcher -class PyomoState(DispatchState): - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - DispatchState.__init__(self) - self._model = None # Pyomo model object + print('DEBUGG Solve successful and no validation concerns raised.') + done_and_checked = True - def initialize(self, components, resources_map, times, model): - """ - Connect information about this State to other objects - @ In, components, list, HERON components - @ In, resources_map, dict, map of component names to resources used - @ In, times, np.array, values of "time" this state represents - @ In, model, pyomo.Model, associated model for this state - @ Out, None - """ - DispatchState.initialize(self, components, resources_map, times) - self._model = model + # In the words of Charles Bukowski, "Don't Try [too many times]" + if attempts > 100: + raise DispatchError('Exceeded validation attempt limit!') - def get_activity_indexed(self, comp, activity, r, t, valued=True, **kwargs): - """ - Getter for activity level. - @ In, comp, HERON Component, component whose information should be retrieved - @ In, activity, str, tracking variable name for activity subset - @ In, r, int, index of resource to retrieve (as given by meta[HERON][resource_indexer]) - @ In, t, int, index of time at which activity should be provided - @ In, valued, bool, optional, if True then get float value instead of pyomo expression - @ In, kwargs, dict, additional pass-through keyword arguments - @ Out, activity, float, amount of resource "res" produced/consumed by "comp" at time "time"; - note positive is producting, negative is consuming - """ - prod = getattr(self._model, f'{comp.name}_{activity}')[r, t] - if valued: - return prod() - return prod + if self.debug_mode: + soln.write() + putils.debug_print_soln(m.model) - def set_activity_indexed(self, comp, r, t, value, valued=False): - raise NotImplementedError + # return dict of numpy arrays + result = putils.retrieve_solution(m.model) + return result