diff --git a/__init__.py b/__init__.py index fc233cf6..dd78e9ac 100644 --- a/__init__.py +++ b/__init__.py @@ -1,5 +1,4 @@ # Copyright 2020, Battelle Energy Alliance, LLC # ALL RIGHTS RESERVED -from HERON.src import DispatchPlot from HERON.src import DispatchManager diff --git a/src/Cases.py b/src/Cases.py index dbe38e43..ce142d07 100644 --- a/src/Cases.py +++ b/src/Cases.py @@ -14,15 +14,16 @@ from HERON.src.base import Base -from HERON.src.dispatch.Factory import known as known_dispatchers -from HERON.src.dispatch.Factory import get_class as get_dispatcher - from HERON.src.ValuedParams import factory as vp_factory from HERON.src.ValuedParamHandler import ValuedParamHandler from HERON.src.validators.Factory import known as known_validators from HERON.src.validators.Factory import get_class as get_validator +from DOVE.src.Dispatch.Factory import known as known_dispatchers +from DOVE.src.Dispatch.Factory import get_class as get_dispatcher + + from collections import OrderedDict import HERON.src._utils as hutils diff --git a/src/Components.py b/src/Components.py index e4d350f4..2daf45ec 100644 --- a/src/Components.py +++ b/src/Components.py @@ -12,36 +12,38 @@ import xml.etree.ElementTree as ET from HERON.src.Economics import CashFlowUser from HERON.src.ValuedParams import factory as vp_factory -from HERON.src.TransferFuncs import factory as tf_factory + from HERON.src.ValuedParamHandler import ValuedParamHandler from HERON.src import _utils as hutils +from DOVE.src.Components import Component as DoveComponent +from DOVE.src.TransferFuncs import factory as tf_factory +from DOVE.src.Interactions import Interaction as DoveInteraction +from DOVE.src.Interactions import Producer as DoveProducer +from DOVE.src.Interactions import Demand as DoveDemand +from DOVE.src.Interactions import Storage as DoveStorage + try: import ravenframework except ModuleNotFoundError: framework_path = hutils.get_raven_loc() sys.path.append(framework_path) -from ravenframework.utils import InputData, xmlUtils,InputTypes +from ravenframework.utils import InputData, xmlUtils, InputTypes -# TODO can we use EntityFactory from RAVEN? -def factory(xml, method='sweep'): - """ - Tool for constructing compnents without the input_loader - TODO can this be set up so the input_loader calls it instead of the methods directly? - @ In, xml, ET.Element, node from which to read component settings - @ In, method, string, optional, operational mode for case - @ Out, comp, Component instance, component constructed - """ - comp = Component() - comp.read_input(xml, method) - comp.finalize_init() - return comp -class Component(Base, CashFlowUser): +class HeronComponent(DoveComponent): """ Represents a unit in the grid analysis. Each component has a single "interaction" that describes what it can do (produce, store, demand) """ + def __repr__(self): + """ + String representation. + @ In, None + @ Out, __repr__, string representation + """ + return f'' + @classmethod def get_input_specs(cls): """ @@ -49,48 +51,111 @@ def get_input_specs(cls): @ In, None @ Out, input_specs, InputData, specs """ - input_specs = InputData.parameterInputFactory('Component', ordered=False, baseNode=None, - descr=r"""defines a component as an element of the grid system. Components are defined by the action they - perform such as \xmlNode{produces} or \xmlNode{consumes}; see details below.""") - input_specs.addParam('name', param_type=InputTypes.StringType, required=True, - descr=r"""identifier for the component. This identifier will be used to generate variables - and relate signals to this component throughout the HERON analysis.""") - # production - ## this unit may be able to make stuff, possibly from other stuff - input_specs.addSub(Producer.get_input_specs()) - # storage - ## this unit may be able to store stuff - input_specs.addSub(Storage.get_input_specs()) - # demands - ## this unit may have a certain demand that must be met - input_specs.addSub(Demand.get_input_specs()) - # this unit probably has some economics - input_specs = CashFlowUser.get_input_specs(input_specs) + ## DEVELOPER NOTE: + ## You should NOT add new subspecs to this method unless they have nothing to + ## do with DOVE (In which case maybe rethink its applicability to HERON). + ## All InputSpecs for Interactions are defined within the DOVE.src.Interactions. + ## If a new VP node needs to be added, find someway to add a fixed-value version + ## to DOVE first! Then modify it in here. This will keep feature parity with DOVE. + ## YOU HAVE BEEN WARNED! + + # Grab all the DOVE input specs -- these input specs have no ValuedParams in + # them, so we need to modify the input spec to allow for those VPs + input_specs = super().get_input_specs() + + # Define the subs to modify along with their configurations -- if we need to + # modify more subs later, just add them to this dict and they'll be added. + interact_subs_to_modify = { + "capacity": { + "add_params": [("resource", "resource")], + "allowed": None + }, + "capacity_factor": { + "add_params": [], + "allowed": ['ARMA', 'CSV'] + }, + "minimum": { + "add_params": [("resource", "resource")], + "allowed": None, + }, + "initial_stored": { + "add_params": [], + "allowed": None, + }, + "strategy": { + "add_params": [], + "allowed": ['Function'], + }, + } + + econ_subs_to_modify = { + "driver":{ + "add_params": [], + "allowed": None, + }, + "reference_price": { + "add_params": [], + "allowed": None, + }, + "reference_driver": { + "add_params": [], + "allowed": None, + }, + "scaling_factor_x": { + "add_params": [], + "allowed": None, + }, + } + + # Iterate over the subs to modify + for sub in input_specs.subs: + for sub_name, config in interact_subs_to_modify.items(): + current_sub = sub.getSub(sub_name) + if current_sub is not None: + print(f"INSIDE HERE -- {sub.getName()} - {sub_name}") + new_sub = vp_factory.make_input_specs(sub_name, descr=sub.description, allowed=config["allowed"]) + # Add parameters if any + for param_name, param_key in config["add_params"]: + new_sub.addParam(param_name, descr=current_sub.parameters[param_key]['description']) + # Replace the old sub with the new one + _ = sub.popSub(sub_name) + sub.addSub(new_sub) + + for sub in input_specs.subs: + if sub.getName() == "economics": + for econ_sub in sub.subs: + if econ_sub.getName() == "CashFlow": + for sub_name, config in econ_subs_to_modify.items(): + current_sub = econ_sub.getSub(sub_name) + if current_sub is not None: + print(f"INSIDE HERE -- {econ_sub.getName()} - {sub_name}") + new_sub = vp_factory.make_input_specs(sub_name, descr=sub.description, allowed=config["allowed"]) + # Add parameters if any + for param_name, param_key in config["add_params"]: + new_sub.addParam(param_name, descr=current_sub.parameters[param_key]['description']) + # Replace the old sub with the new one + econ_sub.popSub(sub_name) + econ_sub.addSub(new_sub) + return input_specs - def __init__(self, **kwargs): - """ - Constructor - @ In, kwargs, dict, optional, arguments to pass to other constructors - @ Out, None - """ - Base.__init__(self, **kwargs) - CashFlowUser.__init__(self) - self.name = None - self._produces = [] - self._stores = [] - self._demands = [] - self.levelized_meta = {} + # def __init__(self, **kwargs): + # """ + # Constructor + # @ In, kwargs, dict, optional, arguments to pass to other constructors + # @ Out, None + # """ + # super().__init__(**kwargs) + # # Base.__init__(self, **kwargs) + # # CashFlowUser.__init__(self) + # self.name = None + # self._produces = [] + # self._stores = [] + # self._demands = [] + # self.levelized_meta = {} - def __repr__(self): - """ - String representation. - @ In, None - @ Out, __repr__, string representation - """ - return f'' - def read_input(self, xml, mode): + def read_input(self, xml, mode="opt"): """ Sets settings from input file @ In, xml, xml.etree.ElementTree.Element, input from user @@ -98,7 +163,10 @@ def read_input(self, xml, mode): @ Out, None """ # get specs for allowable inputs + print("STARRTING READ INPUT") specs = self.get_input_specs()() + print(specs.generateLatex(recDepth=2)) + print("PARSENODE") specs.parseNode(xml) self.name = specs.parameterValues['name'] self.raiseADebug(f'Loading component "{self.name}"') @@ -132,208 +200,7 @@ def read_input(self, xml, mode): self.raiseAnError(IOError, f' node missing from component "{self.name}"!') CashFlowUser.read_input(self, econ_node) - def finalize_init(self): - """ - Finalizes the initialization of the component, and checks input settings. - @ In, None - @ Out, None - """ - self.get_interaction().finalize_init() - - def get_crossrefs(self): - """ - Collect the required value entities needed for this component to function. - @ In, None - @ Out, crossrefs, dict, mapping of dictionaries with information about the entities required. - """ - inter = self.get_interaction() - crossrefs = {inter: inter.get_crossrefs()} - crossrefs.update(self._economics.get_crossrefs()) - return crossrefs - - def set_crossrefs(self, refs): - """ - Connect cross-reference material from other entities to the ValuedParams in this component. - @ In, refs, dict, dictionary of entity information - @ Out, None - """ - try_match = self.get_interaction() - for interaction in list(refs.keys()): - # find associated interaction - if try_match == interaction: - try_match.set_crossrefs(refs.pop(interaction)) - break - # send what's left to the economics - self._economics.set_crossrefs(refs) - # if anything left, there's an issue - assert not refs - - def get_interaction(self): - """ - Return the interactions this component uses. - TODO could this just return the only non-empty one, since there can only be one? - @ In, None - @ Out, interactions, list, list of Interaction instances - """ - try: - return (self._produces + self._stores + self._demands)[0] - except IndexError: # there are no interactions! - return None - - def get_sqrt_RTE(self): - """ - Provide the square root of the round-trip efficiency for this component. - Note we use the square root due to splitting loss across the input and output. - @ In, None - @ Out, RTE, float, round-trip efficiency as a multiplier - """ - return self.get_interaction().get_sqrt_RTE() - - def print_me(self, tabs=0, tab=' '): - """ - Prints info about self - @ In, tabs, int, optional, number of tabs to insert before prints - @ In, tab, str, optional, characters to use to denote hierarchy - @ Out, None - """ - pre = tab*tabs - self.raiseADebug(pre+'Component:') - self.raiseADebug(pre+' name:', self.name) - self.get_interaction().print_me(tabs=tabs+1, tab=tab) - - def get_inputs(self): - """ - returns list of all resources consumed here - @ In, None - @ Out, inputs, set, set of input resources as strings (resources that are taken/consumed/stored) - """ - inputs = set() - # simply combine the inputs for the interaction - inputs.update(self.get_interaction().get_inputs()) - return inputs - - def get_outputs(self): - """ - returns list of all resources producable here - @ In, None - @ Out, outputs, set, set of output resources as strings (resources that are produced/provided) - """ - outputs = set() - outputs.update(self.get_interaction().get_outputs()) - return outputs - - def get_resources(self): - """ - Provides the full set of resources used by this component. - @ In, None - @ Out, res, set, set(str) of resource names - """ - res = set() - res.update(self.get_inputs()) - res.update(self.get_outputs()) - return res - - def get_capacity(self, meta, raw=False): - """ - returns the capacity of the interaction of this component - @ In, meta, dict, arbitrary metadata from EGRET - @ In, raw, bool, optional, if True then return the ValuedParam instance for capacity, instead of the evaluation - @ Out, capacity, float (or ValuedParam), the capacity of this component's interaction - """ - return self.get_interaction().get_capacity(meta, raw=raw) - - def get_minimum(self, meta, raw=False): - """ - returns the minimum of the interaction of this component - @ In, meta, dict, arbitrary metadata from EGRET - @ In, raw, bool, optional, if True then return the ValuedParam instance for capacity, instead of the evaluation - @ Out, capacity, float (or ValuedParam), the capacity of this component's interaction - """ - return self.get_interaction().get_minimum(meta, raw=raw) - - def get_capacity_var(self): - """ - Returns the variable that is used to define this component's capacity. - @ In, None - @ Out, var, str, name of capacity resource - """ - return self.get_interaction().get_capacity_var() - - def get_tracking_vars(self): - """ - Provides the variables used by this component to track dispatch - @ In, None - @ Out, get_tracking_vars, list, variable name list - """ - return self.get_interaction().get_tracking_vars() - - def is_dispatchable(self): - """ - Returns the dispatchability indicator of this component. - TODO Note that despite the name, this is NOT boolean, but a string indicator. - @ In, None - @ Out, dispatchable, str, dispatchability (e.g. independent, dependent, fixed) - """ - return self.get_interaction().is_dispatchable() - - def is_governed(self): - """ - Determines if this component is optimizable or governed by some function. - @ In, None - @ Out, is_governed, bool, whether this component is governed. - """ - return self.get_interaction().is_governed() - - def set_capacity(self, cap): - """ - Set the float value of the capacity of this component's interaction - @ In, cap, float, value - @ Out, None - """ - return self.get_interaction().set_capacity(cap) - - @property - def ramp_limit(self): - """ - Accessor for ramp limits on interactions. - @ In, None - @ Out, limit, float, limit - """ - return self.get_interaction().ramp_limit - - @property - def ramp_freq(self): - """ - Accessor for ramp frequency limits on interactions. - @ In, None - @ Out, limit, float, limit - """ - return self.get_interaction().ramp_freq - - def get_capacity_param(self): - """ - Provides direct access to the ValuedParam for the capacity of this component. - @ In, None - @ Out, cap, ValuedParam, capacity valued param - """ - intr = self.get_interaction() - return intr.get_capacity(None, None, None, None, raw=True) - - def set_levelized_cost_meta(self, cashflows): - """ - Create a dictionary for determining correct resource to use per cashflow if using levelized - inner objective (only an option when selecting LC as an econ metric) - @ In, cashflows, list, list of Interaction instances - @ Out, None - """ - for cf in cashflows: - tracker = cf.get_driver()._vp.get_tracking_var() - resource = cf.get_driver()._vp.get_resource() - self.levelized_meta[cf.name] = {tracker:resource} - - - -class Interaction(Base): +class HeronInteraction(DoveInteraction): """ Base class for component interactions (e.g. Producer, Storage, Demand) """ @@ -346,51 +213,50 @@ def get_input_specs(cls): @ In, None @ Out, input_specs, InputData, specs """ - if cls.tag == 'produces': - desc = r"""indicates that this component produces one or more resources by consuming other resources.""" - resource_desc = r"""the resource produced by this component's activity.""" - elif cls.tag == 'stores': - desc = r"""indicates that this component stores one resource, potentially absorbing or providing that resource.""" - resource_desc = r"""the resource stored by this component.""" - elif cls.tag == "demands": - desc = r"""indicates that this component exclusively consumes a resource.""" - resource_desc = r"""the resource consumed by this component.""" - specs = InputData.parameterInputFactory(cls.tag, ordered=False, descr=desc) - specs.addParam('resource', param_type=InputTypes.StringListType, required=True, - descr=resource_desc) - dispatch_opts = InputTypes.makeEnumType('dispatch_opts', 'dispatch_opts', ['fixed', 'independent', 'dependent']) - specs.addParam('dispatch', param_type=dispatch_opts, - descr=r"""describes the way this component should be dispatched, or its flexibility. - \texttt{fixed} indicates the component always fully dispatched at its maximum level. - \texttt{independent} indicates the component is fully dispatchable by the dispatch optimization algorithm. - \texttt{dependent} indicates that while this component is not directly controllable by the dispatch - algorithm, it can however be flexibly dispatched in response to other units changing dispatch level. - For example, when attempting to increase profitability, the \texttt{fixed} components are not adjustable, - but the \texttt{independent} components can be adjusted to attempt to improve the economic metric. - In response to the \texttt{independent} component adjustment, the \texttt{dependent} components - may respond to balance the resource usage from the changing behavior of other components.""") - - descr = r"""the maximum value at which this component can act, in units corresponding to the indicated resource. """ - cap = vp_factory.make_input_specs('capacity', descr=descr) - cap.addParam('resource', param_type=InputTypes.StringType, - descr=r"""indicates the resource that defines the capacity of this component's operation. For example, - if a component consumes steam and electricity to produce hydrogen, the capacity of the component - can be defined by the maximum steam consumable, maximum electricity consumable, or maximum - hydrogen producable. Any choice should be nominally equivalent, but determines the units - of the value of this node.""") - specs.addSub(cap) - - descr = r"""the actual value at which this component can act, as a unitless fraction of total rated capacity. - Note that these factors are applied within the dispatch optimization; we assume that the capacity factor - is not a variable in the outer optimization.""" - capfactor = vp_factory.make_input_specs('capacity_factor', descr=descr, allowed=['ARMA', 'CSV']) - specs.addSub(capfactor) + ## DEVELOPER NOTE: + ## You should NOT add new subspecs to this method unless they have nothing to + ## do with DOVE (In which case maybe rethink its applicability to HERON). + ## All InputSpecs for Interactions are defined within the DOVE.src.Interactions. + ## If a new VP node needs to be added, find someway to add a fixed-value version + ## to DOVE first! Then modify it in here. This will keep feature parity with DOVE. + ## YOU HAVE BEEN WARNED! + + # Grab all the DOVE input specs -- these input specs have no ValuedParams in + # them, so we need to modify the input spec to allow for those VPs + input_specs = super().get_input_specs() + + # Define the subs to modify along with their configurations -- if we need to + # modify more subs later, just add them to this dict and they'll be added. + subs_to_modify = { + "capacity": { + "add_params": [("resource", "resource")], + "allowed": None + }, + "capacity_factor": { + "add_params": [], + "allowed": ['ARMA', 'CSV'] + }, + "minimum": { + "add_params": [("resource", "resource")], + "allowed": None + }, + } + + # Iterate over the subs to modify + for sub_name, config in subs_to_modify.items(): + # Get the sub from input_specs + sub = input_specs.getSub(sub_name) + if sub is not None: + # Create a new input spec using vp_factory + new_sub = vp_factory.make_input_specs(sub_name, descr=sub.description, allowed=config["allowed"]) + # Add parameters if any + for param_name, param_key in config["add_params"]: + new_sub.addParam(param_name, descr=sub.parameters[param_key]['description']) + # Replace the old sub with the new one + input_specs.popSub(sub_name) + input_specs.addSub(new_sub) - descr = r"""provides the minimum value at which this component can act, as a percentage of the installed capacity. """ - minn = vp_factory.make_input_specs('minimum', descr=descr) - specs.addSub(minn) - - return specs + return input_specs def __init__(self, **kwargs): """ @@ -669,7 +535,7 @@ def get_transfer(self): -class Producer(Interaction): +class HeronProducer(HeronInteraction, DoveProducer): """ Explains a particular interaction, where resources are consumed to produce other resources """ @@ -723,383 +589,383 @@ def get_input_specs(cls): ) return specs - - def __init__(self, **kwargs): - """ - Constructor - @ In, None - @ Out, None - """ - Interaction.__init__(self, **kwargs) - self._produces = [] # the resource(s) produced by this interaction - self._consumes = [] # the resource(s) consumed by this interaction - self._tracking_vars = ['production'] - - def read_input(self, specs, mode, comp_name): - """ - Sets settings from input file - @ In, specs, InputData, specs - @ In, mode, string, case mode to operate in (e.g. 'sweep' or 'opt') - @ In, comp_name, string, name of component this Interaction belongs to - @ Out, None - """ - # specs were already checked in Component - Interaction.read_input(self, specs, mode, comp_name) - self._produces = specs.parameterValues['resource'] - for item in specs.subparts: - if item.getName() == 'consumes': - self._consumes = item.value - elif item.getName() == 'transfer': - self._set_transfer_func('_transfer', comp_name, item) - elif item.getName() == 'ramp_limit': - self.ramp_limit = item.value - elif item.getName() == 'ramp_freq': - self.ramp_freq = item.value - - # input checking - ## if a transfer function not given, can't be consuming a resource - if self._transfer is None: - if self._consumes: - self.raiseAnError(IOError, 'Any component that consumes a resource must have a transfer function describing the production process!') - ## transfer elements are all in IO list - if self._transfer is not None: - self._transfer.check_io(self.get_inputs(), self.get_outputs(), comp_name) - self._transfer.set_io_signs(self.get_inputs(), self.get_outputs()) - ## ramp limit is (0, 1] - if self.ramp_limit is not None and not 0 < self.ramp_limit <= 1: - self.raiseAnError(IOError, f'Ramp limit must be (0, 1] but got "{self.ramp_limit}"') - - def _set_transfer_func(self, name, comp, spec): - """ - Sets up a Transfer Function - @ In, name, str, name of member of this class - @ In, comp, str, name of associated component - @ In, spec, inputparam, input specifications - @ Out, None - """ - known = tf_factory.knownTypes() - found = False - for sub in spec.subparts: - if sub.getName() in known: - if found: - self.raiseAnError(IOError, f'Received multiple Transfer Functions for component "{name}"!') - self._transfer = tf_factory.returnInstance(sub.getName()) - self._transfer.read(comp, spec) - found = True - - def get_inputs(self): - """ - Returns the set of resources that are inputs to this interaction. - @ In, None - @ Out, inputs, set, set of inputs - """ - inputs = Interaction.get_inputs(self) - inputs.update(np.atleast_1d(self._consumes)) - return inputs - - def get_outputs(self): - """ - Returns the set of resources that are outputs to this interaction. - @ In, None - @ Out, outputs, set, set of outputs - """ - outputs = set(np.atleast_1d(self._produces)) - return outputs - - def print_me(self, tabs: int=0, tab: str=' ') -> None: - """ - Prints info about self - @ In, tabs, int, optional, number of tabs to insert before prints - @ In, tab, str, optional, characters to use to denote hierarchy - @ Out, None - """ - pre = tab*tabs - self.raiseADebug(pre+'Producer:') - self.raiseADebug(pre+' produces:', self._produces) - self.raiseADebug(pre+' consumes:', self._consumes) - self.raiseADebug(pre+' transfer:', self._transfer) - self.raiseADebug(pre+' capacity:', self._capacity) - - - -class Storage(Interaction): - """ - Explains a particular interaction, where a resource is stored and released later - """ - tag = 'stores' # node name in input file - - @classmethod - def get_input_specs(cls): - """ - Collects input specifications for this class. - @ In, None - @ Out, input_specs, InputData, specs - """ - specs = super().get_input_specs() - # TODO unused, please implement ... : - # descr = r"""the limiting charge/discharge rate of this storage. """ - # specs.addSub(ValuedParam.get_input_specs('rate')) - # initial stored - descr=r"""indicates what percent of the storage unit is full at the start of each optimization sequence, - from 0 to 1. Overwritten if using periodic level conditions, in which case the initial level is - solved as part of the optimization, but the initial and final levels must match. \default{0.0}. """ - sub = vp_factory.make_input_specs('initial_stored', descr=descr) - specs.addSub(sub) - # periodic level boundary condition - descr=r"""indicates whether the level of the storage should be required to return to its initial level - at the end of each modeling window. If True, replaces the \xmlNode{initial_stored} with an optimization - variable. If False, this increases the flexibility of the storage at the cost of potentially - violating conservation of resources. \default{True}. """ - sub = InputData.parameterInputFactory('periodic_level', contentType=InputTypes.BoolType, descr=descr) - specs.addSub(sub) - # control strategy - descr=r"""control strategy for operating the storage. If not specified, uses a perfect foresight strategy. """ - specs.addSub(vp_factory.make_input_specs('strategy', allowed=['Function'], descr=descr)) - # round trip efficiency - descr = r"""round-trip efficiency for this component as a scalar multiplier. \default{1.0}""" - specs.addSub(InputData.parameterInputFactory('RTE', contentType=InputTypes.FloatType, descr=descr)) - return specs - - def __init__(self, **kwargs): - """ - Constructor - @ In, kwargs, dict, passthrough args - @ Out, None - """ - Interaction.__init__(self, **kwargs) - self.apply_periodic_level = True # whether to apply periodic boundary conditions for the level of the storage - self._stores = None # the resource stored by this interaction - self._rate = None # the rate at which this component can store up or discharge - self._initial_stored = None # how much resource does this component start with stored? - self._strategy = None # how to operate storage unit - self._tracking_vars = ['level', 'charge', 'discharge'] # stored quantity, charge activity, discharge activity - - def read_input(self, specs, mode, comp_name): - """ - Sets settings from input file - @ In, specs, InputData, specs - @ In, mode, string, case mode to operate in (e.g. 'sweep' or 'opt') - @ In, comp_name, string, name of component this Interaction belongs to - @ Out, None - """ - # specs were already checked in Component - Interaction.read_input(self, specs, mode, comp_name) - self._stores = specs.parameterValues['resource'] - for item in specs.subparts: - if item.getName() == 'rate': - self._set_valued_param('_rate', comp_name, item, mode) - elif item.getName() == 'initial_stored': - self._set_valued_param('_initial_stored', comp_name, item, mode) - elif item.getName() == 'periodic_level': - self.apply_periodic_level = item.value - elif item.getName() == 'strategy': - self._set_valued_param('_strategy', comp_name, item, mode) - elif item.getName() == 'RTE': - self._sqrt_rte = np.sqrt(item.value) - assert len(self._stores) == 1, f'Multiple storage resources given for component "{comp_name}"' - self._stores = self._stores[0] - # checks and defaults - if self._initial_stored is None: - self.raiseAWarning(f'Initial storage level for "{comp_name}" was not provided! Defaulting to 0%.') - # make a fake reader node for a 0 value - vp = ValuedParamHandler('initial_stored') - vp.set_const_VP(0.0) - self._initial_stored = vp - # the capacity is limited by the stored resource. - self._capacity_var = self._stores - - def get_inputs(self): - """ - Returns the set of resources that are inputs to this interaction. - @ In, None - @ Out, inputs, set, set of inputs - """ - inputs = Interaction.get_inputs(self) - inputs.update(np.atleast_1d(self._stores)) - return inputs - - def get_outputs(self): - """ - Returns the set of resources that are outputs to this interaction. - @ In, None - @ Out, outputs, set, set of outputs - """ - outputs = Interaction.get_outputs(self) - outputs.update(np.atleast_1d(self._stores)) - return outputs - - def get_resource(self): - """ - Returns the resource this unit stores. - @ In, None - @ Out, stores, str, resource stored - """ - return self._stores - - def get_strategy(self): - """ - Returns the resource this unit stores. - @ In, None - @ Out, stores, str, resource stored - """ - return self._strategy - - def is_governed(self): - """ - Determines if this interaction is optimizable or governed by some function. - @ In, None - @ Out, is_governed, bool, whether this component is governed. - """ - return self._strategy is not None - - def print_me(self, tabs=0, tab=' '): - """ - Prints info about self - @ In, tabs, int, optional, number of tabs to insert before prints - @ In, tab, str, optional, characters to use to denote hierarchy - @ Out, None - """ - pre = tab*tabs - self.raiseADebug(pre+'Storage:') - self.raiseADebug(pre+' stores:', self._stores) - self.raiseADebug(pre+' rate:', self._rate) - self.raiseADebug(pre+' capacity:', self._capacity) - - def _check_capacity_limit(self, res, amt, balance, meta, raven_vars, dispatch, t, level): - """ - Check to see if capacity limits of this component have been violated. - overloads Interaction method, since units for storage are "res" not "res per second" - @ In, res, str, name of capacity-limiting resource - @ In, amt, float, requested amount of resource used in interaction - @ In, balance, dict, results of requested interaction - @ In, meta, dict, additional variable passthrough - @ In, raven_vars, dict, TODO part of meta! consolidate! - @ In, dispatch, dict, TODO part of meta! consolidate! - @ In, t, int, TODO part of meta! consolidate! - @ In, level, float, current level of storage - @ Out, balance, dict, new results of requested action, possibly modified if capacity hit - @ Out, meta, dict, additional variable passthrough - """ - # note "amt" has units of AMOUNT not RATE (resource, not resource per second) - sign = np.sign(amt) - # are we storing or providing? - #print('DEBUGG supposed current level:', level) - if sign < 0: - # we are being asked to consume some - cap, meta = self.get_capacity(meta, raven_vars, dispatch, t) - available_amount = cap[res] - level - #print('Supposed Capacity, Only calculated ins sign<0 (being asked to consumer)',cap) - else: - # we are being asked to produce some - available_amount = level - # the amount we can consume is the minimum of the requested or what's available - delta = sign * min(available_amount, abs(amt)) - return {res: delta}, meta - - def _check_rate_limit(self, res, amt, balance, meta, raven_vars, dispatch, t): - """ - Determines the limiting rate of in/out production for storage - @ In, res, str, name of capacity-limiting resource - @ In, amt, float, requested amount of resource used in interaction - @ In, balance, dict, results of requested interaction - @ In, meta, dict, additional variable passthrough - @ In, raven_vars, dict, TODO part of meta! consolidate! - @ In, dispatch, dict, TODO part of meta! consolidate! - @ In, t, int, TODO part of meta! consolidate! - @ Out, balance, dict, new results of requested action, possibly modified if capacity hit - @ Out, meta, dict, additional variable passthrough - """ - # TODO distinct up/down rates - # check limiting rate for resource flow in/out, if any - if self._rate: - request = {res: None} - inputs = {'request': request, - 'meta': meta, - 'raven_vars': raven_vars, - 'dispatch': dispatch, - 't': t} - max_rate = self._rate.evaluate(inputs, target_var=res)[0][res] - delta = np.sign(amt) * min(max_rate, abs(amt)) - print('max_rate in _check_rate_limit',max_rate, 'delta (min of maxrate and abs(amt)',delta) - return {res: delta}, meta - return {res: amt}, meta - - def get_initial_level(self, meta): - """ - Find initial level of the storage - @ In, meta, dict, additional variable passthrough - @ Out, initial, float, initial level - """ - res = self.get_resource() - request = {res: None} - meta['request'] = request - pct = self._initial_stored.evaluate(meta, target_var=res)[0][res] - if not (0 <= pct <= 1): - self.raiseAnError(ValueError, f'While calculating initial storage level for storage "{self.tag}", ' + - f'an invalid percent was provided/calculated ({pct}). Initial levels should be between 0 and 1, inclusive.') - amt = pct * self.get_capacity(meta)[0][res] - return amt - - - - -class Demand(Interaction): - """ - Explains a particular interaction, where a resource is demanded - """ - tag = 'demands' # node name in input file - - @classmethod - def get_input_specs(cls): - """ - Collects input specifications for this class. - @ In, None - @ Out, input_specs, InputData, specs - """ - specs = super().get_input_specs() - return specs - - def __init__(self, **kwargs): - """ - Constructor - @ In, kwargs, dict, arguments - @ Out, None - """ - Interaction.__init__(self, **kwargs) - self._demands = None # the resource demanded by this interaction - self._tracking_vars = ['production'] - - def read_input(self, specs, mode, comp_name): - """ - Sets settings from input file - @ In, specs, InputData, specs - @ In, mode, string, case mode to operate in (e.g. 'sweep' or 'opt') - @ In, comp_name, string, name of component this Interaction belongs to - @ Out, None - """ - # specs were already checked in Component - # must set demands first, so that "capacity" can access it - self._demands = specs.parameterValues['resource'] - Interaction.read_input(self, specs, mode, comp_name) - - def get_inputs(self): - """ - Returns the set of resources that are inputs to this interaction. - @ In, None - @ Out, inputs, set, set of inputs - """ - inputs = Interaction.get_inputs(self) - inputs.update(np.atleast_1d(self._demands)) - return inputs - - def print_me(self, tabs: int=0, tab: str=' ') -> None: - """ - Prints info about self - @ In, tabs, int, optional, number of tabs to insert before prints - @ In, tab, str, optional, characters to use to denote hierarchy - @ Out, None - """ - pre = tab*tabs - self.raiseADebug(pre+'Demand/Load:') - self.raiseADebug(pre+' demands:', self._demands) - self.raiseADebug(pre+' capacity:', self._capacity) +# +# def __init__(self, **kwargs): +# """ +# Constructor +# @ In, None +# @ Out, None +# """ +# Interaction.__init__(self, **kwargs) +# self._produces = [] # the resource(s) produced by this interaction +# self._consumes = [] # the resource(s) consumed by this interaction +# self._tracking_vars = ['production'] +# +# def read_input(self, specs, mode, comp_name): +# """ +# Sets settings from input file +# @ In, specs, InputData, specs +# @ In, mode, string, case mode to operate in (e.g. 'sweep' or 'opt') +# @ In, comp_name, string, name of component this Interaction belongs to +# @ Out, None +# """ +# # specs were already checked in Component +# Interaction.read_input(self, specs, mode, comp_name) +# self._produces = specs.parameterValues['resource'] +# for item in specs.subparts: +# if item.getName() == 'consumes': +# self._consumes = item.value +# elif item.getName() == 'transfer': +# self._set_transfer_func('_transfer', comp_name, item) +# elif item.getName() == 'ramp_limit': +# self.ramp_limit = item.value +# elif item.getName() == 'ramp_freq': +# self.ramp_freq = item.value +# +# # input checking +# ## if a transfer function not given, can't be consuming a resource +# if self._transfer is None: +# if self._consumes: +# self.raiseAnError(IOError, 'Any component that consumes a resource must have a transfer function describing the production process!') +# ## transfer elements are all in IO list +# if self._transfer is not None: +# self._transfer.check_io(self.get_inputs(), self.get_outputs(), comp_name) +# self._transfer.set_io_signs(self.get_inputs(), self.get_outputs()) +# ## ramp limit is (0, 1] +# if self.ramp_limit is not None and not 0 < self.ramp_limit <= 1: +# self.raiseAnError(IOError, f'Ramp limit must be (0, 1] but got "{self.ramp_limit}"') +# +# def _set_transfer_func(self, name, comp, spec): +# """ +# Sets up a Transfer Function +# @ In, name, str, name of member of this class +# @ In, comp, str, name of associated component +# @ In, spec, inputparam, input specifications +# @ Out, None +# """ +# known = tf_factory.knownTypes() +# found = False +# for sub in spec.subparts: +# if sub.getName() in known: +# if found: +# self.raiseAnError(IOError, f'Received multiple Transfer Functions for component "{name}"!') +# self._transfer = tf_factory.returnInstance(sub.getName()) +# self._transfer.read(comp, spec) +# found = True +# +# def get_inputs(self): +# """ +# Returns the set of resources that are inputs to this interaction. +# @ In, None +# @ Out, inputs, set, set of inputs +# """ +# inputs = Interaction.get_inputs(self) +# inputs.update(np.atleast_1d(self._consumes)) +# return inputs +# +# def get_outputs(self): +# """ +# Returns the set of resources that are outputs to this interaction. +# @ In, None +# @ Out, outputs, set, set of outputs +# """ +# outputs = set(np.atleast_1d(self._produces)) +# return outputs +# +# def print_me(self, tabs: int=0, tab: str=' ') -> None: +# """ +# Prints info about self +# @ In, tabs, int, optional, number of tabs to insert before prints +# @ In, tab, str, optional, characters to use to denote hierarchy +# @ Out, None +# """ +# pre = tab*tabs +# self.raiseADebug(pre+'Producer:') +# self.raiseADebug(pre+' produces:', self._produces) +# self.raiseADebug(pre+' consumes:', self._consumes) +# self.raiseADebug(pre+' transfer:', self._transfer) +# self.raiseADebug(pre+' capacity:', self._capacity) +# +# +# +# class Storage(Interaction): +# """ +# Explains a particular interaction, where a resource is stored and released later +# """ +# tag = 'stores' # node name in input file +# +# @classmethod +# def get_input_specs(cls): +# """ +# Collects input specifications for this class. +# @ In, None +# @ Out, input_specs, InputData, specs +# """ +# specs = super().get_input_specs() +# # TODO unused, please implement ... : +# # descr = r"""the limiting charge/discharge rate of this storage. """ +# # specs.addSub(ValuedParam.get_input_specs('rate')) +# # initial stored +# descr=r"""indicates what percent of the storage unit is full at the start of each optimization sequence, +# from 0 to 1. Overwritten if using periodic level conditions, in which case the initial level is +# solved as part of the optimization, but the initial and final levels must match. \default{0.0}. """ +# sub = vp_factory.make_input_specs('initial_stored', descr=descr) +# specs.addSub(sub) +# # periodic level boundary condition +# descr=r"""indicates whether the level of the storage should be required to return to its initial level +# at the end of each modeling window. If True, replaces the \xmlNode{initial_stored} with an optimization +# variable. If False, this increases the flexibility of the storage at the cost of potentially +# violating conservation of resources. \default{True}. """ +# sub = InputData.parameterInputFactory('periodic_level', contentType=InputTypes.BoolType, descr=descr) +# specs.addSub(sub) +# # control strategy +# descr=r"""control strategy for operating the storage. If not specified, uses a perfect foresight strategy. """ +# specs.addSub(vp_factory.make_input_specs('strategy', allowed=['Function'], descr=descr)) +# # round trip efficiency +# descr = r"""round-trip efficiency for this component as a scalar multiplier. \default{1.0}""" +# specs.addSub(InputData.parameterInputFactory('RTE', contentType=InputTypes.FloatType, descr=descr)) +# return specs +# +# def __init__(self, **kwargs): +# """ +# Constructor +# @ In, kwargs, dict, passthrough args +# @ Out, None +# """ +# Interaction.__init__(self, **kwargs) +# self.apply_periodic_level = True # whether to apply periodic boundary conditions for the level of the storage +# self._stores = None # the resource stored by this interaction +# self._rate = None # the rate at which this component can store up or discharge +# self._initial_stored = None # how much resource does this component start with stored? +# self._strategy = None # how to operate storage unit +# self._tracking_vars = ['level', 'charge', 'discharge'] # stored quantity, charge activity, discharge activity +# +# def read_input(self, specs, mode, comp_name): +# """ +# Sets settings from input file +# @ In, specs, InputData, specs +# @ In, mode, string, case mode to operate in (e.g. 'sweep' or 'opt') +# @ In, comp_name, string, name of component this Interaction belongs to +# @ Out, None +# """ +# # specs were already checked in Component +# Interaction.read_input(self, specs, mode, comp_name) +# self._stores = specs.parameterValues['resource'] +# for item in specs.subparts: +# if item.getName() == 'rate': +# self._set_valued_param('_rate', comp_name, item, mode) +# elif item.getName() == 'initial_stored': +# self._set_valued_param('_initial_stored', comp_name, item, mode) +# elif item.getName() == 'periodic_level': +# self.apply_periodic_level = item.value +# elif item.getName() == 'strategy': +# self._set_valued_param('_strategy', comp_name, item, mode) +# elif item.getName() == 'RTE': +# self._sqrt_rte = np.sqrt(item.value) +# assert len(self._stores) == 1, f'Multiple storage resources given for component "{comp_name}"' +# self._stores = self._stores[0] +# # checks and defaults +# if self._initial_stored is None: +# self.raiseAWarning(f'Initial storage level for "{comp_name}" was not provided! Defaulting to 0%.') +# # make a fake reader node for a 0 value +# vp = ValuedParamHandler('initial_stored') +# vp.set_const_VP(0.0) +# self._initial_stored = vp +# # the capacity is limited by the stored resource. +# self._capacity_var = self._stores +# +# def get_inputs(self): +# """ +# Returns the set of resources that are inputs to this interaction. +# @ In, None +# @ Out, inputs, set, set of inputs +# """ +# inputs = Interaction.get_inputs(self) +# inputs.update(np.atleast_1d(self._stores)) +# return inputs +# +# def get_outputs(self): +# """ +# Returns the set of resources that are outputs to this interaction. +# @ In, None +# @ Out, outputs, set, set of outputs +# """ +# outputs = Interaction.get_outputs(self) +# outputs.update(np.atleast_1d(self._stores)) +# return outputs +# +# def get_resource(self): +# """ +# Returns the resource this unit stores. +# @ In, None +# @ Out, stores, str, resource stored +# """ +# return self._stores +# +# def get_strategy(self): +# """ +# Returns the resource this unit stores. +# @ In, None +# @ Out, stores, str, resource stored +# """ +# return self._strategy +# +# def is_governed(self): +# """ +# Determines if this interaction is optimizable or governed by some function. +# @ In, None +# @ Out, is_governed, bool, whether this component is governed. +# """ +# return self._strategy is not None +# +# def print_me(self, tabs=0, tab=' '): +# """ +# Prints info about self +# @ In, tabs, int, optional, number of tabs to insert before prints +# @ In, tab, str, optional, characters to use to denote hierarchy +# @ Out, None +# """ +# pre = tab*tabs +# self.raiseADebug(pre+'Storage:') +# self.raiseADebug(pre+' stores:', self._stores) +# self.raiseADebug(pre+' rate:', self._rate) +# self.raiseADebug(pre+' capacity:', self._capacity) +# +# def _check_capacity_limit(self, res, amt, balance, meta, raven_vars, dispatch, t, level): +# """ +# Check to see if capacity limits of this component have been violated. +# overloads Interaction method, since units for storage are "res" not "res per second" +# @ In, res, str, name of capacity-limiting resource +# @ In, amt, float, requested amount of resource used in interaction +# @ In, balance, dict, results of requested interaction +# @ In, meta, dict, additional variable passthrough +# @ In, raven_vars, dict, TODO part of meta! consolidate! +# @ In, dispatch, dict, TODO part of meta! consolidate! +# @ In, t, int, TODO part of meta! consolidate! +# @ In, level, float, current level of storage +# @ Out, balance, dict, new results of requested action, possibly modified if capacity hit +# @ Out, meta, dict, additional variable passthrough +# """ +# # note "amt" has units of AMOUNT not RATE (resource, not resource per second) +# sign = np.sign(amt) +# # are we storing or providing? +# #print('DEBUGG supposed current level:', level) +# if sign < 0: +# # we are being asked to consume some +# cap, meta = self.get_capacity(meta, raven_vars, dispatch, t) +# available_amount = cap[res] - level +# #print('Supposed Capacity, Only calculated ins sign<0 (being asked to consumer)',cap) +# else: +# # we are being asked to produce some +# available_amount = level +# # the amount we can consume is the minimum of the requested or what's available +# delta = sign * min(available_amount, abs(amt)) +# return {res: delta}, meta +# +# def _check_rate_limit(self, res, amt, balance, meta, raven_vars, dispatch, t): +# """ +# Determines the limiting rate of in/out production for storage +# @ In, res, str, name of capacity-limiting resource +# @ In, amt, float, requested amount of resource used in interaction +# @ In, balance, dict, results of requested interaction +# @ In, meta, dict, additional variable passthrough +# @ In, raven_vars, dict, TODO part of meta! consolidate! +# @ In, dispatch, dict, TODO part of meta! consolidate! +# @ In, t, int, TODO part of meta! consolidate! +# @ Out, balance, dict, new results of requested action, possibly modified if capacity hit +# @ Out, meta, dict, additional variable passthrough +# """ +# # TODO distinct up/down rates +# # check limiting rate for resource flow in/out, if any +# if self._rate: +# request = {res: None} +# inputs = {'request': request, +# 'meta': meta, +# 'raven_vars': raven_vars, +# 'dispatch': dispatch, +# 't': t} +# max_rate = self._rate.evaluate(inputs, target_var=res)[0][res] +# delta = np.sign(amt) * min(max_rate, abs(amt)) +# print('max_rate in _check_rate_limit',max_rate, 'delta (min of maxrate and abs(amt)',delta) +# return {res: delta}, meta +# return {res: amt}, meta +# +# def get_initial_level(self, meta): +# """ +# Find initial level of the storage +# @ In, meta, dict, additional variable passthrough +# @ Out, initial, float, initial level +# """ +# res = self.get_resource() +# request = {res: None} +# meta['request'] = request +# pct = self._initial_stored.evaluate(meta, target_var=res)[0][res] +# if not (0 <= pct <= 1): +# self.raiseAnError(ValueError, f'While calculating initial storage level for storage "{self.tag}", ' + +# f'an invalid percent was provided/calculated ({pct}). Initial levels should be between 0 and 1, inclusive.') +# amt = pct * self.get_capacity(meta)[0][res] +# return amt +# +# +# +# +# class Demand(Interaction): +# """ +# Explains a particular interaction, where a resource is demanded +# """ +# tag = 'demands' # node name in input file +# +# @classmethod +# def get_input_specs(cls): +# """ +# Collects input specifications for this class. +# @ In, None +# @ Out, input_specs, InputData, specs +# """ +# specs = super().get_input_specs() +# return specs +# +# def __init__(self, **kwargs): +# """ +# Constructor +# @ In, kwargs, dict, arguments +# @ Out, None +# """ +# Interaction.__init__(self, **kwargs) +# self._demands = None # the resource demanded by this interaction +# self._tracking_vars = ['production'] +# +# def read_input(self, specs, mode, comp_name): +# """ +# Sets settings from input file +# @ In, specs, InputData, specs +# @ In, mode, string, case mode to operate in (e.g. 'sweep' or 'opt') +# @ In, comp_name, string, name of component this Interaction belongs to +# @ Out, None +# """ +# # specs were already checked in Component +# # must set demands first, so that "capacity" can access it +# self._demands = specs.parameterValues['resource'] +# Interaction.read_input(self, specs, mode, comp_name) +# +# def get_inputs(self): +# """ +# Returns the set of resources that are inputs to this interaction. +# @ In, None +# @ Out, inputs, set, set of inputs +# """ +# inputs = Interaction.get_inputs(self) +# inputs.update(np.atleast_1d(self._demands)) +# return inputs +# +# def print_me(self, tabs: int=0, tab: str=' ') -> None: +# """ +# Prints info about self +# @ In, tabs, int, optional, number of tabs to insert before prints +# @ In, tab, str, optional, characters to use to denote hierarchy +# @ Out, None +# """ +# pre = tab*tabs +# self.raiseADebug(pre+'Demand/Load:') +# self.raiseADebug(pre+' demands:', self._demands) +# self.raiseADebug(pre+' capacity:', self._capacity) diff --git a/src/DispatchPlot.py b/src/DispatchPlot.py deleted file mode 100644 index 16559ad2..00000000 --- a/src/DispatchPlot.py +++ /dev/null @@ -1,350 +0,0 @@ -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" - Author: dylanjm - Date: 2021-05-18 -""" -import itertools as it -from collections import defaultdict - -import matplotlib as mpl -mpl.use('Agg') # Prevents the script from blocking while plotting -import matplotlib.pyplot as plt - -from typing import List, Dict -import numpy as np - -try: - from ravenframework.PluginBaseClasses.OutStreamPlotPlugin import PlotPlugin, InputTypes, InputData -except ModuleNotFoundError: - import sys - from . import _utils - sys.path.append(_utils.get_raven_loc()) - from ravenframework.PluginBaseClasses.OutStreamPlotPlugin import PlotPlugin, InputTypes, InputData - - -# default color cycler, hatches -colormap = plt.get_cmap('tab10').colors -hatchmap = [None, '..', '\\\\', 'xx', '--', 'oo', '++', '**', 'OO'] - -# Matplotlib Global Settings -plt.rc("figure", figsize=(12, 8), titleweight='bold') # type: ignore -plt.rc("axes", titleweight="bold", labelsize=12, axisbelow=True, grid=True) # type: ignore -plt.rc("savefig", bbox="tight") # type: ignore -plt.rc("legend", fontsize=12) # type:ignore -plt.rc(["xtick", "ytick"], labelsize=10) # type: ignore - - -class DispatchPlot(PlotPlugin): - - @classmethod - def getInputSpecification(cls): - """ - Define the acceptable user inputs for this class. - @ In, None - @ Out, specs, InputData.ParameterInput, - """ - specs = super().getInputSpecification() - specs.addSub(InputData.parameterInputFactory('source', contentType=InputTypes.StringType)) - specs.addSub(InputData.parameterInputFactory('macro_variable', contentType=InputTypes.StringType)) - specs.addSub(InputData.parameterInputFactory('micro_variable', contentType=InputTypes.StringType)) - specs.addSub(InputData.parameterInputFactory('signals', contentType=InputTypes.StringListType)) - return specs - - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - super().__init__() - self.printTag = 'HERON.DispatchPlot' - self._sourceName = None - self._source = None - self._macroName = None - self._microName = None - self._addSignals = [] - - def handleInput(self, spec): - """ - Reads in data from the input file - @ In, spec, InputData.ParameterInput, input information - @ Out, None - """ - super().handleInput(spec) - for node in spec.subparts: - if node.getName() == 'source': - self._sourceName = node.value - elif node.getName() == 'macro_variable': - self._macroName = node.value - elif node.getName() == 'micro_variable': - self._microName = node.value - elif node.getName() == 'signals': - self._addSignals = node.value - - def initialize(self, stepEntities): - """ - Set up plotter for each run - @ In, stepEntities, dict, entities from the Step - @ Out, None - """ - super().initialize(stepEntities) - src = self.findSource(self._sourceName, stepEntities) - if src is None: - self.raiseAnError(IOError, f'Source DataObject "{self._sourceName}" was not found in the Step!') - self._source = src - - @staticmethod - def _group_by(iterable: List[str], idx: int) -> Dict[str, List[str]]: - """ - Return a dictionary containing grouped dispatch variables. - @ In, iterable, List[str], a list of variable names to group-by. - @ In, idx, int, the index of the variable to group-by. - @ Out, gr, Dict[str, List[str]], a dictionary mapping of grouped variable names. - """ - gr = {} - for var in iterable: - # var is expected to have the form: 'Dispatch__component__tracker__resource' - key = var.split('__')[idx] - if key in gr.keys(): - gr[key].append(var) - else: - gr[key] = [var] - return gr - - def plot_dispatch(self, figs, axes, df, grp_vars, sid, mstep, cid, cdict) -> None: - """ - Plot and output the optimized dispatch for a specific sample, year, and cluster. - @ In, figs, matplotlib.figure.Figure, current figures used for plotting. - @ In, axes, List[List[matplotlib.Axes]], a list of axes across figures to plot each variable. - @ In, df, pandas.DataFrame, a dataframe containing data to plot. - @ In, grp_vars, Dict[str, List[str]], a dictionary mapping components to variables. - @ In, sid, int, the sample ID. - @ In, mstep, int, the macro step. - @ In, cid, int, the cluster ID. - @ In, cdict, Dict[Dict[str, Tuple[Float]]], a dictionary contains color code to variables - @ Out, None - """ - alpha = 0.7 - time = df[self._microName].to_numpy() - for (key, group), ax in zip(grp_vars.items(), axes.flat): - # Define list for data, label, and color. Seperate 'level'(line plot) with other variables (stack plot) - positive_dat = [] - positive_label = [] - positive_color = [] - positive_hatch = [] - - negative_dat = [] - negative_label = [] - negative_color = [] - negative_hatch = [] - - level_dat = [] - level_label = [] - level_color = [] - - # Secondary y axis for storage levels - ax2 = ax.twinx() - - # Fill the lists - for var in group: - _, comp_name, tracker, _ = var.split('__') - comp_label = comp_name.replace('_', ' ')#.title() - var_label = f'{comp_label}, {tracker.title()}' - ls = '-' - # Fill the positive, negative, and level lists - info = cdict[comp_name]#[res][tracker] - color = info['color'] - hatch = info['hatch'] - if (df[var] != 0).any(): # no plotting variables that have all zeros values - if tracker == 'level': - # this is the level of a storage component - level_dat.append(var) - level_label.append(var_label) - level_color.append(color) - else: - if (df[var] > 0).any(): - # these are production tracking variables - positive_dat.append(var) - positive_label.append(var_label) - positive_color.append(tuple([*color, alpha])) - positive_hatch.append(hatch) - else: - # these are consumption tracking variables - negative_dat.append(var) - negative_label.append(var_label) - negative_color.append(tuple([*color, alpha])) - negative_hatch.append(hatch) - - # Plot the micro-step variable on the x-axis (i.e Time) - # center (0) line - ax.plot([time[0],time[-1]], [0,0], 'k-') - # Production - if len(positive_dat) > 0: - pos_stacks = ax.stackplot(time, - *[df[key] for key in positive_dat], - labels=positive_label, - baseline='zero', - colors=positive_color) - for stack, hatch in zip(pos_stacks, positive_hatch): - stack.set_hatch(hatch) - - # Consumption - if len(negative_dat) > 0: - neg_stacks = ax.stackplot(time, - *[df[key] for key in negative_dat], - labels= negative_label, - baseline='zero', - colors=negative_color) - for stack, hatch in zip(neg_stacks, negative_hatch): - stack.set_hatch(hatch) - - # Levels - if(len(level_dat) > 0): - for key, c, llabel in zip(level_dat, level_color, level_label): - ax2.plot(df[self._microName], df[key], linestyle=ls, label=llabel, color=c) - - # Set figure title, legend, and grid - ax.set_title(key.title().split('_')[-1]) - ax.set_xlabel(self._microName) - if(len(positive_label) > 0 or len(negative_label) > 0): - ax.legend(loc='upper left', bbox_to_anchor=(1.1, 0.6), fontsize = 10) - if(len(level_label) > 0): - ax2.legend(loc='lower left', bbox_to_anchor=(1.1, 0.6), fontsize = 10) - # Add the label and adjust location - ax.set_ylabel('Activity', fontsize=10, rotation=0) - ax2.set_ylabel('Level', fontsize=10, rotation=0) - ax.yaxis.set_label_coords(-0.01,1.02) - ax2.yaxis.set_label_coords(1,1.07) - ax.grid(None) - ax2.grid(None) - # Output and save the image - for f, fig in enumerate(figs): - file_name = f"dispatch_id{sid}_y{mstep}_c{cid}_f{f+1}.png" - fig.suptitle(f'Dispatch ID {sid} Year {mstep} Cluster {cid},\nFigure {f+1}/{len(figs)}') - fig.tight_layout() - fig.savefig(file_name) - self.raiseAMessage(f'Saved figure {f+1}/{len(figs)} to "{file_name}"') - plt.clf() - - def plot_signal(self, fig, axes, df, sid, mstep, cid) -> None: - """ - Plot and output the synthetic history for a specific sample, year, and cluster. - @ In, fig, matplotlib.figure.Figure, a current figure used for plotting. - @ In, axes, List[List[matplotlib.Axes]], a list of axes to plot each variable. - @ In, df, pandas.DataFrame, a dataframe containing data to plot. - @ In, sid, int, the sample ID. - @ In, mstep, int, the macro step. - @ In, cid, int, the cluster ID. - """ - for name, ax in zip(self._addSignals, axes.flat): - var = df.get(name, None) - if var is None: - self.raiseAnError(RuntimeError, f'Requested signal variable "{name}" but variable not in data!') - ax.plot(df[self._microName], var, marker='.', linestyle='-', label=name) - ax.set_title(name.title()) - ax.set_xlabel(self._microName) - ax.legend(loc='center left', bbox_to_anchor=(1.03, 0.5)) - - signal_file_name = f"dispatch_id{sid}_y{mstep}_c{cid}_SIGNAL.png" - fig.tight_layout() - fig.savefig(signal_file_name) - self.raiseAMessage(f'Saved figure to "{signal_file_name}"') - plt.clf() - - def color_style(self, grp_vars): - """ - Set the coloring scheme for each of the variables that will be plotted. - @ In, grp_vars, Dict[str, List[str]], a dictionary mapping components to variables. - @ Out, comp_colors, Dict[Dict[str, Tuple[Float]]], contains rgb color code and hatching for components - """ - # DESIGN: - # -> components should be clearly different in color, and consistent across resource plots - # get the components, resources they use, and trackers per resource - # TODO this is just a rearrangement of the data, is it really useful, or is there another way? - comps = defaultdict(dict) - for res, group in grp_vars.items(): - for variable in group: - _, comp, tracker, res = variable.split('__') - if res not in comps[comp]: - comps[comp][res] = [tracker] - else: - comps[comp][res].append(tracker) - n_comps = len(comps) - n_colors = len(colormap) - n_hatches = len(hatchmap) - n_uniques = n_colors * n_hatches - if n_comps > n_uniques: - self.raiseAWarning(f'A total of {n_comps} exist to plot, but only {n_uniques} unique identifying ' + - 'colors and patterns are available! This may lead to dispatch plot confusion.') - print('Assigning colors ...') - comp_colors = {} - # kept for easy debugging - #print('DEBUGG | c | ci | hi | comp | hatch | color_r | color_g | color_b') - for c, (comp, ress) in enumerate(comps.items()): - hatch_index, color_index = divmod(c, n_colors) - color = colormap[color_index] - hatch = hatchmap[hatch_index] - # kept for easy debugging - #print(f'DEBUGG | {c:2d} | {color_index:1d} | {hatch_index:1d} | {comp:20s} | '+ - # f'{hatch if hatch is not None else "None":4s} | '+ - # f'{color[0]*255:1.8f} | {color[1]*255:1.8f} | {color[2]*255:1.8f}') - comp_colors[comp] = {'color': color, 'hatch': hatch} - return comp_colors - - def run(self): - """ - Generate the plot - @ In, None - @ Out, None - """ - assert self._source is not None - ds = self._source.asDataset() - if ds is None: - self.raiseAWarning(f'No data in "{self._source.name}" data object; nothing to plot!') - return - df = ds.to_dataframe().reset_index() - dispatch_vars = list(filter(lambda x: "Dispatch__" in x, df.columns)) - grouped_vars = self._group_by(dispatch_vars, -1) - - # Dimension variables to plot - sample_ids = df[self._source.sampleTag].unique() - cluster_ids = df['_ROM_Cluster'].unique() # TODO: find way to not hardcode name - macro_steps = df[self._macroName].unique() - - # Assign colors - cdict = self.color_style(grouped_vars) - - resources = set([x for x in grouped_vars]) - for sample_id, macro_step, cluster_id in it.product(sample_ids, macro_steps, cluster_ids): - # Filter data to plot correct values for current dimension - dat = df[ - (df[self._source.sampleTag] == sample_id) & - (df[self._macroName] == macro_step) & - (df['_ROM_Cluster'] == cluster_id) - ] - - # TODO: find a way to combine both plots into one output. - # Currently, this is difficult to do because of the nested - # nature of the subplots, as well as the dynamic number of - # components and signals to plot (i.e. dynamically nested subplots) - - # If only 3 resources, make one figure; otherwise, 2 resources per figure - if len(resources) <= 3: - fig, res_axs = plt.subplots(len(resources), 1, sharex=True, squeeze=False) - res_figs = [fig] - else: - res_figs = [] - res_axs = [] - for _ in range(int(np.ceil(len(resources) / 2))): - fig, axs = plt.subplots(2, 1, sharex=True, squeeze=False) - res_figs.append(fig) - res_axs.extend(axs) - - # Output optimized component dispatch for current dimension. - res_axs = np.asarray(res_axs) - self.plot_dispatch(res_figs, res_axs, dat, grouped_vars, sample_id, macro_step, cluster_id, cdict) - - # Output synthetic time series signal for current dimension. - sig_figs, sig_axs = plt.subplots(len(self._addSignals), 1, sharex=True, squeeze=False) - self.plot_signal(sig_figs, sig_axs, dat, sample_id, macro_step, cluster_id) diff --git a/src/TransferFuncs/Factory.py b/src/TransferFuncs/Factory.py deleted file mode 100644 index 8479a714..00000000 --- a/src/TransferFuncs/Factory.py +++ /dev/null @@ -1,40 +0,0 @@ - -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED - -from ravenframework.utils import InputData, InputTypes -from ravenframework.EntityFactoryBase import EntityFactory - -from .Ratio import Ratio -from .Polynomial import Polynomial - -class TransferFuncFactory(EntityFactory): - """ - Factory for Transfer Functions - """ - def make_input_specs(self, name, descr=None): - """ - Fill input specs for the provided name and description. - @ In, name, str, name of new spec - @ In, descr, str, optional, description of spec - @ Out, spec, InputData, specification - """ - add_descr = '' - if descr is None: - description = add_descr - else: - description = descr + r"""\\ \\""" + add_descr - - spec = InputData.parameterInputFactory(name, descr=description) - for name, klass in self._registeredTypes.items(): - sub_spec = klass.get_input_specs() - sub_spec.name = name - spec.addSub(sub_spec) - return spec - -factory = TransferFuncFactory('TransferFunc') - -# fixed in inner -factory.registerType('ratio', Ratio) -factory.registerType('linear', Ratio) -factory.registerType('poly', Polynomial) diff --git a/src/TransferFuncs/Polynomial.py b/src/TransferFuncs/Polynomial.py deleted file mode 100644 index 2d0dd86f..00000000 --- a/src/TransferFuncs/Polynomial.py +++ /dev/null @@ -1,86 +0,0 @@ - -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" - Transfer fucntions that are expressed as a polynomial relationship. - For example, ax^2 + bxy + cy^2 + dx + ey = fm^2 + gmn + hn^2 + im + jn + k -""" -from collections import defaultdict - -from .TransferFunc import TransferFunc, InputData, InputTypes - -class Polynomial(TransferFunc): - """ - Represents a ValuedParam that is a polynomial relationship. - """ - - @classmethod - def get_input_specs(cls): - """ - Define parameters for a polynomial transfer function. - @ In, None - @ Out, spec, InputData, value-based spec - """ - spec = InputData.parameterInputFactory('poly', contentType=InputTypes.StringType, - descr=r"""indicates this transfer function is expressed by a polynomial relationship of arbitrary order. - Note the polynomial must be specified in weak form, with all terms on one side of the equation - set equal to zero. For instance, the equation $ax^2 + bx + c = dy^2 + fy + g$ should be reformulated - as $ax^2 + bx + (c-g) - dy^2 - fy = 0$.""") - coeff = InputData.parameterInputFactory('coeff', contentType=InputTypes.FloatType, - descr=r"""one coefficient for one poloynomial term of the specified \xmlAttr{resources}. - Care should be taken to assure the sign of the coefficient is working as expected, - as consumed resources have a negative sign while produced resources have a positive - sign, and the full equation should have the form 0 = ... .""") - coeff.addParam('resource', param_type=InputTypes.StringListType, - descr=r"""indicates the resource(s) for which the polynomial coefficient is being provided in this node. - Note that the order of the resources matters for specifying the polynomial \xmlAttr{order}.""") - coeff.addParam('order', param_type=InputTypes.IntegerListType, - descr=r"""indicates the orders of the polynomial for each resource specified, in order. - For example, if \xmlAttr{resources} is ``x, y'', then order ``2,3'' would mean - the specified coefficient is for $x^{2}y^{3}$.""") - spec.addSub(coeff) - return spec - - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - super().__init__() - # coeffs, stored as { (resources): { (order): coeff, }, } - # e.g., {(water, flour): {(1,1): 42, (1,2): 3.14}, - # (water): {(1): 1.61} } - self._coefficients = defaultdict(dict) - - def read(self, comp_name, spec): - """ - Used to read valued param from XML input - @ In, comp_name, str, name of component that this valued param will be attached to; only used for print messages - @ In, spec, InputData params, input specifications - @ Out, None - """ - super().read(comp_name, spec) - for coeff_node in spec.findFirst('poly').findAll('coeff'): - resource = coeff_node.parameterValues['resource'] - order = coeff_node.parameterValues['order'] - self._coefficients[tuple(resource)][tuple(order)] = coeff_node.value - - def get_resources(self): - """ - Provides the resources used in this transfer function. - @ In, None - @ Out, resources, set, set of resources used - """ - res_set = set() - for res_tup, ord_dict in self._coefficients.items(): - res_set = res_set.union(set(res_tup)) - return res_set - - def get_coefficients(self): - """ - Returns linear coefficients. - @ In, None - @ Out, coeffs, dict, coefficient mapping - """ - return self._coefficients diff --git a/src/TransferFuncs/Ratio.py b/src/TransferFuncs/Ratio.py deleted file mode 100644 index 741baf51..00000000 --- a/src/TransferFuncs/Ratio.py +++ /dev/null @@ -1,99 +0,0 @@ - -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" - Values that are expressed as linear ratios of one another. - Primarily intended for transfer functions. -""" -import numpy as np - -from .TransferFunc import TransferFunc, InputData, InputTypes - -# class for custom dynamically-evaluated quantities -class Ratio(TransferFunc): - """ - Represents a TransferFunc that uses a linear balance of resources, such as 3a + 7b -> 2c. - This means the ratios of the resources must be maintained, NOT 3a + 7b = 2c! - """ - - @classmethod - def get_input_specs(cls): - """ - Input specification for this class. - @ In, None - @ Out, spec, InputData, value-based spec - """ - spec = InputData.parameterInputFactory('ratio', contentType=InputTypes.StringType, - descr=r"""indicates this transfer function is a constant linear combination of resources. For example, - a balance equation might be written as 3a + 7b -> 2c, implying that to make 2c, it always takes - 3 parts a and 7 parts b, or the balance ratio (3a, 7b, 2c). This means that the ratio of (3, 7, 2) must be - maintained between (a, b, c) for all production levels. Note that the coefficient signs are automatically fixed - internally to be negative for consumed quantities and positive for produced resources, regardless of signs used - by the user. For an equation-based transfer function instead of balance ratio, see Polynomial.""") - rate = InputData.parameterInputFactory('rate', contentType=InputTypes.FloatType, - descr=r"""linear coefficient for the indicated \xmlAttr{resource}.""") - rate.addParam('resource', param_type=InputTypes.StringType, - descr=r"""indicates the resource for which the linear transfer ratio is being provided in this node.""") - spec.addSub(rate) - return spec - - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - super().__init__() - self._coefficients = None # ratios, stored in a dict as key: value - - def read(self, comp_name, spec): - """ - Used to read transfer func from XML input - @ In, comp_name, str, name of component that this valued param will be attached to; only used for print messages - @ In, spec, InputData params, input specifications - @ Out, None - """ - super().read(comp_name, spec) - self._coefficients = {} - node = spec.findFirst('ratio') - # ALIAS SUPPORT - if node is None: - node = spec.findFirst('linear') - if node is None: - self.raiseAnError(IOError, f'Unrecognized transfer function for component "{comp_name}": "{spec.name}"') - self.raiseAWarning('"linear" has been deprecated and will be removed in the future; see "ratio" transfer function!') - for rate_node in node.findAll('rate'): - resource = rate_node.parameterValues['resource'] - self._coefficients[resource] = rate_node.value - - def get_resources(self): - """ - Provides the resources used in this transfer function. - @ In, None - @ Out, resources, set, set of resources used - """ - return set(self._coefficients.keys()) - - def get_coefficients(self): - """ - Returns linear coefficients. - @ In, None - @ Out, coeffs, dict, coefficient mapping - """ - return self._coefficients - - def set_io_signs(self, consumed, produced): - """ - Fix up input/output signs, if interpretable - @ In, consumed, list, list of resources consumed in the transfer - @ In, produced, list, list of resources produced in the transfer - @ Out, None - """ - for res, coef in self.get_coefficients().items(): - if res in consumed: - self._coefficients[res] = - np.abs(coef) - elif res in produced: - self._coefficients[res] = np.abs(coef) - else: - # should not be able to get here after IO gets checked! - raise RuntimeError('While checking transfer coefficient, resource "{res}" was neither consumed nor produced!') diff --git a/src/TransferFuncs/TransferFunc.py b/src/TransferFuncs/TransferFunc.py deleted file mode 100644 index 3cd92309..00000000 --- a/src/TransferFuncs/TransferFunc.py +++ /dev/null @@ -1,93 +0,0 @@ - -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" - Defines the TransferFunc entity. - These define the transfer functions for generating Components. -""" -import sys -from HERON.src import _utils as hutils -try: - import ravenframework -except ModuleNotFoundError: - framework_path = hutils.get_raven_loc() - sys.path.append(framework_path) -from ravenframework.utils import InputData, InputTypes -from ravenframework.BaseClasses import MessageUser - -# class for potentially dynamically-evaluated quantities -class TransferFunc(MessageUser): - """ - These define the transfer functions for generating Components. - """ - - @classmethod - def get_input_specs(cls, name): - """ - Define inputs for this VP. - @ In, name, string, name for spec (tag) - @ Out, spec, InputData, value-based spec - """ - spec = InputData.parameterInputFactory(name) - return spec - - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - super().__init__() - self.type = self.__class__.__name__ # class type, for easy checking - - def __repr__(self) -> str: - """ - Return Object Representation String - @ In, None - @ Out, None - """ - return "" - - def read(self, comp_name, spec): - """ - Used to read transfer function from XML input - @ In, comp_name, str, name of component that this transfer function is describing - @ In, spec, InputData params, input specifications - @ Out, needs, list, signals needed to evaluate this ValuedParam at runtime - """ - - def check_io(self, inputs, outputs, comp_name): - """ - Checks that the transfer function uses all and only the resources used by the component. - @ In, inputs, list, list of input resources to check against. - @ In, outputs, list, list of output resources to check against. - @ In, comp_name, str, name of component that this transfer function is describing - @ Out, None - """ - used = self.get_resources() - inps = set(inputs) - outs = set(outputs) - excess_inputs = inps - used - excess_outputs = outs - used - unrecog = used - inps.union(outs) - if excess_inputs or excess_outputs or unrecog: - msg = f'Transfer function for Component "{comp_name}" has a mismatch with consumed and produced!' - msg += f'\n... All Consumed: {inps}' - msg += f'\n... All Produced: {outs}' - msg += f'\n... All in Transfer Function: {used}' - if excess_inputs: - msg += f'\n... Consumed but not used in transfer: {excess_inputs}' - if excess_outputs: - msg += f'\n... Produced but not used in transfer: {excess_outputs}' - if unrecog: - msg += f'\n... In transfer but not consumed or produced: {unrecog}' - self.raiseAnError(IOError, msg) - - def set_io_signs(self, consumed, produced): - """ - Fix up input/output signs, if interpretable - @ In, consumed, list, list of resources consumed in the transfer - @ In, produced, list, list of resources produced in the transfer - @ Out, None - """ - # nothing to do by default diff --git a/src/TransferFuncs/__init__.py b/src/TransferFuncs/__init__.py deleted file mode 100644 index b8422e5d..00000000 --- a/src/TransferFuncs/__init__.py +++ /dev/null @@ -1,15 +0,0 @@ -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" - Transfer functions describe the balance between consumed and produced - resources for generating components. This module defines the templates - that can be used to describe transfer functions. -""" -# only type references here, as needed -from .TransferFunc import TransferFunc -from .Ratio import Ratio -from .Polynomial import Polynomial - -# provide easy name access to module -from .Factory import factory - diff --git a/src/dispatch/CustomDispatcher.py b/src/dispatch/CustomDispatcher.py deleted file mode 100644 index 7a0b9bbe..00000000 --- a/src/dispatch/CustomDispatcher.py +++ /dev/null @@ -1,107 +0,0 @@ - -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" - Interface for user-provided dispatching strategies. -""" -import os -import inspect -import numpy as np - -from ravenframework.utils import utils, InputData, InputTypes -from .Dispatcher import Dispatcher -from .DispatchState import NumpyState - -class Custom(Dispatcher): - """ - Base class for strategies for consecutive dispatching in a continuous period. - """ - # --------------------------------------------- - # INITIALIZATION - @classmethod - def get_input_specs(cls): - """ - Set acceptable input specifications. - @ In, None - @ Out, specs, InputData, specs - """ - specs = InputData.parameterInputFactory('custom', ordered=False, baseNode=None) - specs.addSub(InputData.parameterInputFactory('location', contentType=InputTypes.StringType, - descr=r"""The hard drive location of the custom dispatcher. Relative paths are taken with - respect to the HERON run location. Custom dispatchers must implement - a \texttt{dispatch} method that accepts the HERON case, components, and sources; this - method must return the activity for each resource of each component.""")) - return specs - - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - Dispatcher.__init__(self) - self.name = 'CustomDispatcher' - self._usr_loc = None # user-provided path to custom dispatcher module - self._file = None # resolved pathlib.Path to the custom dispatcher module - - def read_input(self, inputs): - """ - Loads settings based on provided inputs - @ In, inputs, InputData.InputSpecs, input specifications - @ Out, None - """ - usr_loc = inputs.findFirst('location') - if usr_loc is None: - raise RuntimeError('No provided for dispatch strategy in !') - # assure python extension, omitting it is a user convenience - if not usr_loc.value.endswith('.py'): - usr_loc.value += '.py' - self._usr_loc = os.path.abspath(os.path.expanduser(usr_loc.value)) - - def initialize(self, case, components, sources, **kwargs): - """ - Initialize dispatcher properties. - @ In, case, Case, HERON case instance - @ In, components, list, HERON components - @ In, sources, list, HERON sources - @ In, kwargs, dict, keyword arguments - @ Out, None - """ - start_loc = case.run_dir - file_loc = os.path.abspath(os.path.join(start_loc, self._usr_loc)) - # check that it exists - if not os.path.isfile(file_loc): - raise IOError(f'Custom dispatcher not found at "{file_loc}"! (input dir "{start_loc}", provided path "{self._usr_loc}"') - self._file = file_loc - print(f'Loading custom dispatch at "{self._file}"') - # load user module - load_string, _ = utils.identifyIfExternalModelExists(self, self._file, '') - module = utils.importFromPath(load_string, True) - # check it works as intended - if not 'dispatch' in dir(module): - raise IOError(f'Custom Dispatch at "{self._file}" does not have a "dispatch" method!') - # TODO other checks? - - def dispatch(self, case, components, sources, meta): - """ - Performs technoeconomic dispatch. - @ In, case, Case, HERON case - @ In, components, list, HERON components - @ In, sources, list, HERON sources - @ Out, results, dict, economic and production metrics - """ - # load up time indices - 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 - # load user module - load_string, _ = utils.identifyIfExternalModelExists(self, self._file, '') - module = utils.importFromPath(load_string, True) - state = NumpyState() - indexer = meta['HERON']['resource_indexer'] - state.initialize(components, indexer, time) - # run dispatch - results = module.dispatch(meta, state) - # TODO: Check to make sure user has uploaded all activity data. - return state - - diff --git a/src/dispatch/DispatchState.py b/src/dispatch/DispatchState.py deleted file mode 100644 index 11cacb15..00000000 --- a/src/dispatch/DispatchState.py +++ /dev/null @@ -1,244 +0,0 @@ - -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -import numpy as np -from io import StringIO - -class DispatchState: - """ utility that expresses the activity (i.e. production level) of all the components in the system """ - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - self._components = None # list of HERON Component objects - self._resources = None # Map of resources to indices for components, as {comp.name: {res, r}} - self._times = None # numpy array of time values, monotonically increasing - - def initialize(self, components, resources_map, times): - """ - Set up dispatch state to hold data - @ In, components, list, HERON components to be stored - @ In, resources_map, dict, map of resources to indices for each component - @ In, times, list, float times to store activity - @ Out, None - """ - self._components = components - self._resources = resources_map - self._times = times - - def __repr__(self): - """ - Compiles string representation of object. - @ In, None - @ Out, repr, str, string representation - """ - return '' - - def get_activity(self, comp, activity, res, time, **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, res, string, name of resource to retrieve - @ In, time, float, time at which activity should be provided - @ Out, activity, float, amount of resource "res" produced/consumed by "comp" at time "time"; - note positive is producting, negative is consuming - """ - r = self._resources[comp][res] - t = np.searchsorted(self._times, time) # TODO protect against value not present - return self.get_activity_indexed(comp, activity, r, t, **kwargs) - - def set_activity(self, comp, activity, res, time, value, **kwargs): - """ - Setter for activity level. - @ In, comp, HERON Component, component whose information should be set - @ In, activity, str, tracking variable name for activity subset - @ In, res, string, name of resource to retrieve - @ In, time, float, time at which activity should be provided - @ In, value, float, activity level; note positive is producting, negative is consuming - @ In, kwargs, dict, additional pass-through keyword arguments - @ Out, None - """ - r = self._resources[comp][res] - t = np.searchsorted(self._times, time) # TODO protect against value not present - self.set_activity_indexed(comp, activity, r, t, value, **kwargs) - - def get_activity_indexed(self, *args, **kwargs): - """ - Getter for activity level, using indexes instead of values for r and t - @ In, comp, HERON Component, component whose information should be retrieved - @ 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, 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 - """ - # to be overwritten by implementing classes - raise NotImplementedError - - def set_activity_indexed(self, *args, **kwargs): - """ - Getter for activity level, using indexes instead of values for r and t - @ In, args, list, additional pass-through arguments - @ 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 - """ - # to be overwritten by implementing classes - raise NotImplementedError - - def create_raven_vars(self, template): - """ - Writes out RAVEN variables as expected - @ In, template, str, formating string for variable names (using {comp}, {res}) - @ Out, data, dict, map of raven var names to numpy array data - """ - #template = 'Dispatch__{c}__{r}' # standardized via input - data = {} - for comp in self._components: - for tracker in comp.get_tracking_vars(): - for res, r in self._resources[comp].items(): - result = np.empty(len(self._times)) - for t in range(len(self._times)): - result[t] = self.get_activity_indexed(comp, tracker, r, t) - data[template.format(comp=comp.name, tracker=tracker, res=res)] = result - return data - -# NumpyState is the nominal DispatchState implementation -class NumpyState(DispatchState): - """ implemenatation of DispatchState using Numpy. A good nominal choice if additional functionality isn't needed. """ - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - DispatchState.__init__(self) - self._data = None # dict of numpy ND array for data - - def initialize(self, components, resources_map, times): - """ - Set up dispatch state to hold data - @ In, components, list, HERON components to be stored - @ In, resources_map, dict, map of resources to indices for each component - @ In, time, list, float times to store - @ Out, None - """ - DispatchState.initialize(self, components, resources_map, times) - self._data = {} - for comp in components: - for tag in comp.get_tracking_vars(): - self._data[f'{comp.name}_{tag}'] = np.zeros((len(self._resources[comp]), len(times))) - - def __repr__(self): - """ - Compiles string representation of object. - @ In, None - @ Out, repr, str, string representation - """ - msg = StringIO() - msg.write('') - return msg.getvalue() - - def get_activity_indexed(self, comp, activity, r, t, **kwargs): - """ - Getter for activity level. - Note, if any of the arguments are "None" it is assumed that means "all" - @ 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, 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 - """ - return self._data[f'{comp.name}_{activity}'][r, t] - - def set_activity_indexed(self, comp, activity, r, t, value, **kwargs): - """ - Setter for activity level. - @ In, comp, HERON Component, component whose information should be set - @ In, activity, str, tracking variable name for activity subset - @ In, res, string, name of resource to retrieve - @ In, time, float, time at which activity should be provided - @ In, value, float, activity level; note positive is producting, negative is consuming - @ In, kwargs, dict, additional pass-through keyword arguments - @ Out, None - """ - self._data[f'{comp.name}_{activity}'][r, t] = value - - # def set_activity_vector(self, comp, tracker, res, start_time, end_time, values): - def set_activity_vector(self, comp, res, values, tracker='production', start_idx=0, end_idx=None): - """ - Shortcut utility for setting values all-at-once in a vector. - @ In, comp, HERON Component, component whose information should be set - @ In, res, string, name of resource to retrieve - @ In, values, np.array, activity level; note positive is producting, negative is consuming - @ In, tracker, str, optional, tracking variable name for activity subset, Default: 'production' - @ In, start_idx, int, optional, first time index at which activity is provided, Default: 0 - @ In, end_idx, int, optional, last time index at which activity is provided, Default: None - @ Out, None - """ - if comp.get_interaction().is_type('Storage') and tracker == 'production': - raise RuntimeError(f'Component "{comp.name}" is Storage and the provided tracker is "production"!') - if end_idx is None: - end_idx = len(self._times) - - 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 deleted file mode 100644 index 5d06128d..00000000 --- a/src/dispatch/Dispatcher.py +++ /dev/null @@ -1,131 +0,0 @@ - -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" - Base class for dispatchers. -""" -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. - """ - # --------------------------------------------- - # INITIALIZATION - @classmethod - def get_input_specs(cls): - """ - Set acceptable input specifications. - @ In, None - @ Out, specs, InputData, specs - """ - specs = InputData.parameterInputFactory('Dispatcher', ordered=False, baseNode=None) - # only place specifications that are true for ALL dispatchers here - return specs - - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - super().__init__() - self.name = 'BaseDispatcher' - self._time_discretization = None # (start, end, num_steps) to build time discretization - self._validator = None # can be used to validate activity - self._solver = None - self._eps = 1e-9 # small constant to add to denominator - - def read_input(self, inputs): - """ - Loads settings based on provided inputs - @ In, inputs, InputData.InputSpecs, input specifications - @ Out, None - """ - pass # add here if true for ALL dispatchers - - def initialize(self, case, components, sources, **kwargs): - """ - Initialize dispatcher properties. - @ In, case, Case, HERON case instance - @ In, components, list, HERON components - @ In, sources, list, HERON sources - @ In, kwargs, dict, keyword arguments - @ Out, None - """ - pass - - # --------------------------------------------- - # GETTER AND SETTERS - def get_time_discr(self): - """ - Retrieves the time discretization information. - @ In, None - @ Out, info, tuple, (start, end, number of steps) for time discretization - """ - return self._time_discretization - - def set_time_discr(self, info): - """ - Stores the time discretization information. - @ In, info, tuple, (start, end, number of steps) for time discretization - @ Out, None - """ - assert info is not None - # TODO is this the right idea? - # don't expand into linspace right now, just store the pieces - self._time_discretization = info - - def set_validator(self, validator): - """ - Sets the dispatch validation instance to use in dispatching. - @ In, validator, HERON Validator, instance of validator - @ Out, None - """ - self._validator = validator - - def get_solver(self): - """ - Retrieves the solver information (if applicable) - @ In, None - @ Out, solver, str, name of solver used - """ - return self._solver - - # --------------------------------------------- - # API - # TODO make this a virtual method? - def dispatch(self, case, components, sources): - """ - Performs technoeconomic dispatch. - @ In, case, Case, HERON case - @ In, components, list, HERON components - @ In, sources, list, HERON sources - @ Out, results, dict, economic and production metrics - """ - raise NotImplementedError # must be implemented by inheriting classes - - def validate(self, components, activity, times, meta): - """ - Method to validate a dispatch 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, extra information needed for validation - @ Out, validation, dict, information about validation - """ - # default implementation - if self._validator is not None: - return self._validator.validate(components, activity, times, meta) - else: - # no validator, nothing needs to be changed - return {} - diff --git a/src/dispatch/Factory.py b/src/dispatch/Factory.py deleted file mode 100644 index 75baabd1..00000000 --- a/src/dispatch/Factory.py +++ /dev/null @@ -1,18 +0,0 @@ - -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -from .pyomo_dispatch import Pyomo -from .CustomDispatcher import Custom - -known = { - 'pyomo': Pyomo, - 'custom': Custom, -} - -def get_class(typ): - """ - Returns the requested dispatcher type. - @ In, typ, str, name of one of the dispatchers - @ Out, class, object, class object - """ - return known.get(typ, None) diff --git a/src/dispatch/PyomoModelHandler.py b/src/dispatch/PyomoModelHandler.py deleted file mode 100644 index 4d37397f..00000000 --- a/src/dispatch/PyomoModelHandler.py +++ /dev/null @@ -1,591 +0,0 @@ -# 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 - """ - transfer = comp.get_interaction().get_transfer() - if transfer is None: - return - if transfer.type == 'Ratio': - self._create_transfer_ratio(transfer, comp, prod_name) - elif transfer.type == 'Polynomial': - self._create_transfer_poly(transfer, comp, prod_name) - else: - raise NotImplementedError(f'Transfer function type "{transfer.type}" not implemented for PyomoModelHandler!') - - def _create_transfer_ratio(self, transfer, comp, prod_name): - """ - Create a balance ratio-based transfer function. - This comes in the form of a balance expression (not an equality). - @ In, transfer, TransferFunc, Ratio transfer function - @ In, comp, Component, component object for this transfer - @ In, prod_name, str, name of production element - @ Out, None - """ - name = comp.name - coeffs = transfer.get_coefficients() - coeffs_iter = iter(coeffs.items()) - first_name, first_coef = next(coeffs_iter) - first_r = self.model.resource_index_map[comp][first_name] - for resource, coef in coeffs_iter: - ratio = coef / first_coef - r = self.model.resource_index_map[comp][resource] - rule_name = f'{name}_{resource}_{first_name}_transfer' - rule = lambda mod, t: prl.ratio_transfer_rule(ratio, r, first_r, prod_name, mod, t) - constr = pyo.Constraint(self.model.T, rule=rule) - setattr(self.model, rule_name, constr) - - def _create_transfer_poly(self, transfer, comp, prod_name): - """ - Create a polynomial transfer function. This comes in the form of an equality expression. - @ In, transfer, TransferFunc, Ratio transfer function - @ In, comp, Component, component object for this transfer - @ In, prod_name, str, name of production element - @ Out, None - """ - rule_name = f'{comp.name}_transfer_func' - coeffs = transfer.get_coefficients() - # dict of form {(r1, r2): {(o1, o2): n}} - # where: - # r1, r2 are resource names - # o1, o2 are polynomial orders (may not be integers?) - # n is the float polynomial coefficient for the term - rule = lambda mod, t: prl.poly_transfer_rule(coeffs, self.model.resource_index_map[comp], 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' - if comp.get_interaction().apply_periodic_level: - level_var = getattr(self.model, level_name) - initial = level_var[(r, self.model.T[-1])] - else: - initial = self.initial_storage[comp] - rule = lambda mod, t: prl.level_rule(comp, level_name, charge_name, discharge_name, initial, r, mod, t) - setattr(self.model, level_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 deleted file mode 100644 index 1c7c687a..00000000 --- a/src/dispatch/PyomoRuleLibrary.py +++ /dev/null @@ -1,330 +0,0 @@ -# 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, float or var, initial storage level - @ 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 - 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: - # special treatment if we have frequency-limiting binaries available - # eps is an aux parameter (small number) to ensure that binary vars for ramp-up and ramp- - # down events behave correctly (e.g., are zero when activity remains constant) - eps = 1e-6*limit - 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 is an aux parameter (small number) to ensure that binary vars for ramp-up and ramp- - # down events behave correctly (e.g., are zero when activity remains constant) - eps = 1e-6*limit - 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 ratio_transfer_rule(ratio: float, r: int, ref_r: int,prod_name: str, m, t) -> bool: - """ - Constructs transfer function constraints - @ In, ratio, float, balanced ratio of this resource to the reference resource - @ In, r, int, index for this resource in the activity map - @ In, ref_r, int, index of the reference resource in the activity map - @ 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 - """ - activity = getattr(m, prod_name) - return activity[r, t] == activity[ref_r, t] * ratio - -def poly_transfer_rule(coeffs, r_map, prod_name, m, t) -> bool: - """ - Constructs transfer function constraints - @ In, coeffs, dict, nested mapping of resources and polynomial orders to coefficients - as {(r1, r2): {(o1, o2): n}} - @ In, r_map, dict, mapping of resources to activity indices for this component - @ 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 - """ - activity = getattr(m, prod_name) - eqn = 0 - for resources, ord_dict in coeffs.items(): - for orders, coeff in ord_dict.items(): - term = coeff - for r, res in enumerate(resources): - map_index = r_map[res] - prod = activity[map_index, t] - term *= prod ** orders[r] - eqn += term - return eqn == 0 diff --git a/src/dispatch/__init__.py b/src/dispatch/__init__.py deleted file mode 100644 index 5e3ed630..00000000 --- a/src/dispatch/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ - -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED diff --git a/src/dispatch/putils.py b/src/dispatch/putils.py deleted file mode 100644 index 876487ab..00000000 --- a/src/dispatch/putils.py +++ /dev/null @@ -1,183 +0,0 @@ -# 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)) diff --git a/src/dispatch/pyomo_dispatch.py b/src/dispatch/pyomo_dispatch.py deleted file mode 100644 index 9d6d9695..00000000 --- a/src/dispatch/pyomo_dispatch.py +++ /dev/null @@ -1,357 +0,0 @@ -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" - pyomo-based dispatch strategy -""" -import os -import time as time_mod -import pprint -import numpy as np -import pyutilib.subprocess.GlobalData -import logging - -import pyomo.environ as pyo -from pyomo.opt import SolverStatus, TerminationCondition -from pyomo.util.infeasible import log_infeasible_constraints -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 -pyutilib.subprocess.GlobalData.DEFINE_SIGNAL_HANDLERS_DEFAULT = False - -# 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 = { - 'ipopt': 'tol', - 'cbc': 'primalTolerance', - 'glpk': 'mipgap', -} - -class DispatchError(Exception): - """ - Custom exception for dispatch errors. - """ - pass - -class Pyomo(Dispatcher): - """ - Dispatches using rolling windows in Pyomo - """ - ### INITIALIZATION - @classmethod - def get_input_specs(cls): - """ - Set acceptable input specifications. - @ 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}.""" - ) - ) - # TODO specific for pyomo dispatcher - return specs - - - def __init__(self): - """ - Constructor. - @ In, None - @ Out, None - """ - super().__init__() - self.name = 'PyomoDispatcher' # identifying name - self.debug_mode = False # whether to print additional information - self.solve_options = {} # options passed from Pyomo to the solver - self._window_len = 24 # time window length to dispatch at a time # FIXME user input - self._solver = None # overwrite option for solver - self._picard_limit = 10 # iterative solve limit - - - def read_input(self, specs) -> None: - """ - Read in input specifications. - @ In, specs, RAVEN InputData, specifications - @ Out, None - """ - super().read_input(specs) - - window_len_node = specs.findFirst('rolling_window_length') - if window_len_node is not None: - self._window_len = window_len_node.value - - debug_node = specs.findFirst('debug_mode') - if debug_node is not None: - self.debug_mode = debug_node.value - - solver_node = specs.findFirst('solver') - if solver_node is not None: - self._solver = solver_node.value - - tol_node = specs.findFirst('tol') - if tol_node is not None: - solver_tol = tol_node.value - else: - solver_tol = None - - self._solver = putils.check_solver_availability(self._solver) - - if solver_tol is not None: - 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 - - - def dispatch(self, case, components, sources, meta): - """ - Performs dispatch. - @ 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, meta, dict, additional variables passed through - @ Out, disp, DispatchScenario, resulting dispatch - """ - t_start, t_end, t_num = self.get_time_discr() - 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) - - start_index = 0 - final_index = len(time) - initial_levels = putils.get_initial_storage_levels(components, meta, start_index) - - while start_index < final_index: - end_index = min(start_index + self._window_len, final_index) - if end_index - start_index == 1: - raise DispatchError("Window length of 1 detected, which is not supported.") - - specific_time = time[start_index:end_index] - 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') - - # 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 _handle_dispatch_window_solve(self, specific_time, start_index, case, components, sources, resources, initial_levels, meta): - """ - 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. - """ - 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 - - 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 - @ In, time_offset, int, offset of the time index in the greater history - @ In, case, HERON Case, Case that this dispatch is part of - @ In, components, list, HERON components available to the dispatch - @ In, sources, list, HERON source (placeholders) for signals - @ In, resources, list, sorted list of all resources in problem - @ In, initial_storage, dict, initial storage levels if any - @ In, meta, dict, additional variables passed through - @ Out, result, dict, results of window dispatch - """ - model = PyomoModelHandler(time, time_offset, case, components, resources, initial_storage, meta) - model.populate_model() - result = self._solve_dispatch(model, meta) - return result - - - 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 - """ - if old is None: - return False - - for comp in components: - intr = comp.get_interaction() - 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]) - # 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: - return False - return True - - - 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 - """ - return any(comp.get_interaction().is_governed() for comp in components) - - - def _solve_dispatch(self, m, meta): - """ - 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 - """ - # start a solution search - done_and_checked = False - attempts = 0 - # DEBUGG show variables, bounds - if self.debug_mode: - putils.debug_pyomo_print(m.model) - - while not done_and_checked: - attempts += 1 - print(f'DEBUGG using solver: {self._solver}') - print(f'DEBUGG solve attempt {attempts} ...:') - soln = pyo.SolverFactory(self._solver).solve(m.model, 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: - putils.debug_pyomo_print(m.model) - print('Resource Map:') - pprint.pprint(m.model.resource_index_map) - log_infeasible_constraints(m.model, log_expression=True, log_variables=True) - log_name = 'constraint_violations.log' - logging.basicConfig(filename=log_name, encoding='utf-8', level=logging.INFO) - raise DispatchError( - f'Solve was unsuccessful, see log file located at: "{os.getcwd()}/{log_name}" for more details! Status: {soln.solver.status} Termination: {soln.solver.termination_condition}' - ) - - # 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: - print('DEBUGG Solve successful and no validation concerns raised.') - done_and_checked = True - - # In the words of Charles Bukowski, "Don't Try [too many times]" - if attempts > 100: - raise DispatchError('Exceeded validation attempt limit!') - - if self.debug_mode: - soln.write() - putils.debug_print_soln(m.model) - - # return dict of numpy arrays - result = putils.retrieve_solution(m.model) - return result diff --git a/src/dispatch/twin_pyomo_limited_ramp.py b/src/dispatch/twin_pyomo_limited_ramp.py deleted file mode 100644 index d2dcfecf..00000000 --- a/src/dispatch/twin_pyomo_limited_ramp.py +++ /dev/null @@ -1,354 +0,0 @@ -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" -Enables quick testing of Pyomo problem solves, intending to set up a system -similar to how it is set up in the Pyomo dispatcher. This allows rapid testing -of different configurations for the rolling window optimization. -""" -import platform # only for choosing solver based on HERON installation - -import numpy as np -import matplotlib.pyplot as plt - -import pyomo.environ as pyo - -# This isn't strictly required, just how HERON installs default options -if platform.system() == 'Windows': - SOLVER = 'glpk' -else: - SOLVER = 'cbc' - -# -------------------------------------------------------------------------------- -# set up physical system, resources, time discretization -# -components = ['NPP', 'BOP', 'Grid'] -resources = ['steam', 'electricity'] -time = np.linspace(0, 24, 24) # from @1 to @2 in @3 steps -dt = time[1] - time[0] -resource_map = {'NPP': {'steam': 0}, - 'BOP': {'steam': 0, 'electricity': 1}, - 'Grid': {'electricity': 0}, - } - -# sizing specifications -steam_produced = 100 # MWt/h of steam at NPP -gen_consume_limit = 110 # consumes at most 110 MWt steam at BOP -sink_limit = 1000 # MWh = MW of electricity at Grid - -ramp_up_limit = 100 # MWt/h steam NPP ramp up limit -ramp_down_limit = 100 # MWt/h steam NPP ramp down limit -delta = 1 # aux var TODO good value? -time_between_ramp_up = 10 # hours before unit can ramp up again - -# -------------------------------------------------------------------------------- -# general economics -# -# marginal cost of NPP, 5 $/MWe = 1.67 $ / MWt -margin_NPP = 1.67 -# sales price of electricity: use a sine function for demo -max_price = 10 -freq = 2 * np.pi / 12 -elec_price = max_price * np.sin(freq * time) ** 2 -# debug plot of price curve -# import matplotlib.pyplot as plt -# plt.plot(time, elec_price, 'o-') -# plt.show() - -# -------------------------------------------------------------------------------- -# set up tracker for how components are dispatched ("activity") -# Map is activity: component: resource -# -> e.g. activity[NPP][steam] -# -activity = {} -for comp in components: - activity[comp] = np.zeros((len(resources), len(time)), dtype=float) - -# -------------------------------------------------------------------------------- -# Set up method to construct the Pyomo concrete model. -# -def make_concrete_model(): - """ - Test writing a simple concrete model with terms typical to the pyomo dispatcher. - @ In, None - @ Out, m, pyo.ConcreteModel, instance of the model to solve - """ - 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(0, len(time), dtype=int) # indexes time - # move onto model - m.C = pyo.Set(initialize=C) - m.R = pyo.Set(initialize=R) - m.T = pyo.Set(initialize=T) - #******************* - # FOCUS: create equations for limiting ramp frequency - m.ramp_up = pyo.Var(m.T, initialize=0, within=pyo.Binary) - m.ramp_down = pyo.Var(m.T, initialize=0, within=pyo.Binary) - m.ramp_none = pyo.Var(m.T, initialize=0, within=pyo.Binary) - # store some stuff for reference -> NOT NOTICED by Pyomo, we hope - #******************* - m.Times = time - m.Components = components - m.resource_index_map = resource_map - m.Activity = activity - #******************* - # set up optimization variables - # -> for now we just do this manually - # NPP - m.NPP_index_map = pyo.Set(initialize=range(len(m.resource_index_map['NPP']))) - m.NPP_production = pyo.Var(m.NPP_index_map, m.T, initialize=0) - # BOP - m.BOP_index_map = pyo.Set(initialize=range(len(m.resource_index_map['BOP']))) - m.BOP_production = pyo.Var(m.BOP_index_map, m.T, initialize=0) - # Grid - m.Grid_index_map = pyo.Set(initialize=range(len(m.resource_index_map['Grid']))) - m.Grid_production = pyo.Var(m.Grid_index_map, m.T, initialize=0) - #******************* - # set up lower, upper bounds - # -> for testing we just do this manually - # -> consuming is negative sign by convention! - # -> producing is positive sign by convention! - # steam source produces between 0 and # steam - m.NPP_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.NPP_production[0, t] >= 0) - m.NPP_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.NPP_production[0, t] <= steam_produced) - # elec generator can consume steam to produce electricity; 0 < consumed steam < 1000 - # -> this effectively limits electricity production, but we're defining capacity in steam terms for fun - # -> therefore signs are negative, -1000 < consumed steam < 0! - m.BOP_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.BOP_production[0, t] >= -gen_consume_limit) - m.BOP_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.BOP_production[0, t] <= 0) - # elec sink can take any amount of electricity - # -> consuming, so -10000 < consumed elec < 0 - m.Grid_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.Grid_production[0, t] >= -sink_limit) - m.Grid_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.Grid_production[0, t] <= 0) - #******************* - # create transfer function - # 2 steam make 1 electricity (sure, why not) - m.BOP_transfer = pyo.Constraint(m.T, rule=_generator_transfer) - #******************* - # create conservation rules - # steam - m.steam_conservation = pyo.Constraint(m.T, rule=_conserve_steam) - # electricity - m.elec_conservation = pyo.Constraint(m.T, rule=_conserve_electricity) - #******************* - # FOCUS: ramping rules - m.npp_ramp_up_limit = pyo.Constraint(m.T, rule=_ramp_up_NPP) - m.npp_ramp_down_limit = pyo.Constraint(m.T, rule=_ramp_down_NPP) - m.npp_ramp_binaries = pyo.Constraint(m.T, rule=_ramp_binaries_NPP) - m.npp_ramp_freq_limit = pyo.Constraint(m.T, rule=_ramp_time_limit_NPP) - #******************* - # create objective function - m.OBJ = pyo.Objective(sense=pyo.maximize, rule=_economics) - ####### - # return - return m - -####### -# -# Callback Functions -# -def _generator_transfer(m, t): - """ - Constraint rule for electricity generation in generator - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - # 1 MWt steam -> 0.33 MWe power - return - m.BOP_production[0, t] == 3.0 * m.BOP_production[1, t] - -def _conserve_steam(m, t): - """ - Constraint rule for conserving steam - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - # signs are tricky here, consumption is negative and production is positive - sources = m.NPP_production[0, t] - sinks = m.BOP_production[0, t] - return sources + sinks == 0 - -def _conserve_electricity(m, t): - """ - Constraint rule for conserving electricity - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - sources = m.BOP_production[1, t] - sinks = m.Grid_production[0, t] - return sources + sinks == 0 - -def _ramp_up_NPP(m, t): - """ - Constraining rule for ramping up the NPP. - @ In, m, pyomo.ConcreteModel, model - @ In, t, int, relevant time index - @ Out, rule, expression, evaluatable for Pyomo constraint - """ - if t > 0: - ineq = m.NPP_production[0, t] - m.NPP_production[0, t-1] <= ramp_up_limit * m.ramp_up[t] - delta * m.ramp_down[t] - else: - ineq = pyo.Constraint.Skip - return ineq - -def _ramp_down_NPP(m, t): - """ - Constraining rule for ramping down the NPP. - @ In, m, pyomo.ConcreteModel, model - @ In, t, int, relevant time index - @ Out, rule, expression, evaluatable for Pyomo constraint - """ - if t > 0: - ineq = m.NPP_production[0, t-1] - m.NPP_production[0, t] <= ramp_down_limit * m.ramp_down[t] - delta * m.ramp_up[t] - else: - ineq = pyo.Constraint.Skip - return ineq - -def _ramp_binaries_NPP(m, t): - """ - Binaries to assure only one of (down, steady, up) can occur per time step - @ In, m, pyomo.ConcreteModel, model - @ In, t, int, relevant time index - @ Out, rule, expression, evaluatable for Pyomo constraint - """ - return m.ramp_up[t] + m.ramp_down[t] + m.ramp_none[t] == 1 - -def _ramp_time_limit_NPP(m, t): - """ - Limits time between sequential ramp events - @ In, m, pyomo.ConcreteModel, model - @ In, t, int, relevant time index - @ Out, rule, expression, evaluatable for Pyomo constraint - """ - if t == 0: - return pyo.Constraint.Skip - tao = min(t, time_between_ramp_up) - limit = 0 - for tm in range(t-tao, t): - # equivalent and slightly easier to understand? - limit += 1 - m.ramp_down[tm] - return m.ramp_up[t] <= 1/tao * limit - -def _economics(m): - """ - Constraint rule for optimization target - @ In, m, pyo.ConcreteModel, model containing problem - @ Out, objective, float, constraining evaluation - """ - # marginal cost of operating NPP - opex = sum(m.BOP_production[0, t] for t in m.T) * margin_NPP # will be negative b/c consumed - # grid sales, calculated above; negative sign because we're consuming - sales = - sum((m.Grid_production[0, t] * elec_price[t]) for t in m.T) # net positive because consumed - return opex + sales - -####### -# -# Debug printing functions -# -def print_setup(m): - """ - Debug printing for pre-solve model setup - @ In, m, pyo.ConcreteModel, model containing problem - @ 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 extract_soln(m): - """ - Extracts final solution from model evaluation - @ In, m, pyo.ConcreteModel, model - @ Out, res, dict, results dictionary for dispatch - """ - res = {} - T = len(m.T) - res['prices'] = np.zeros(T) - res['NPP_steam'] = np.zeros(T) - res['BOP_steam'] = np.zeros(T) - res['BOP_elec'] = np.zeros(T) - res['Grid_elec'] = np.zeros(T) - for t in m.T: - res['prices'][t] = elec_price[t] - res['NPP_steam'][t] = m.NPP_production[0, t].value - res['BOP_steam'][t] = m.BOP_production[0, t].value - res['BOP_elec'][t] = m.BOP_production[1, t].value - res['Grid_elec'][t] = m.Grid_production[0, t].value - return res - -def plot_solution(m): - """ - Plots solution from optimized model - @ In, m, pyo.ConcreteModel, model - @ Out, None - """ - res = extract_soln(m) - fig, axs = plt.subplots(3, 1, sharex=True) - axs[0].set_ylabel(r'Steam MW$_t$') - axs[1].set_ylabel(r'Elec MW$_e$') - axs[2].set_ylabel(r'Prices \$/MW$_e$') - axs[2].set_xlabel('Time (h)') - axs[0].plot(time, res['NPP_steam'], 'o-', label='NPP_steam') - axs[0].plot(time, res['BOP_steam'], 'o-', label='BOP_steam') - axs[0].legend() - axs[1].plot(time, res['BOP_elec'], 'o-', label='BOP_elec') - axs[1].plot(time, res['Grid_elec'], 'o-', label='Grid_elec') - axs[1].legend() - axs[2].plot(time, res['prices'], 'o-', label='Prices') - axs[2].legend() - plt.suptitle(f'Ramp up limit: {time_between_ramp_up} h') - plt.savefig(f'dispatch_limit_{time_between_ramp_up}.png') - -def print_solution(m): - """ - Debug printing for post-solve model setup - @ In, m, pyo.ConcreteModel, model containing problem - @ Out, None - """ - print('') - print('*'*80) - print('solution:') - print(' objective value:', m.OBJ()) - print('time | price | steam src | elec gen (s, e) | elec sink') - for t in m.T: - print(f'{m.Times[t]:1.2e} | ' + - f'{elec_price[t]: 3.3f} | ' + - f'{m.NPP_production[0, t].value: 1.3e} | ' + - f'({m.BOP_production[0, t].value: 1.3e}, {m.BOP_production[1, t].value: 1.3e}) | ' + - f'{m.Grid_production[0, t].value: 1.3e}' - ) - print('*'*80) - - -####### -# -# Solver. -# -def solve_model(m): - """ - Solves the model. - @ In, m, pyo.ConcreteModel, model containing problem - @ Out, m, pyo.ConcreteModel, results - """ - soln = pyo.SolverFactory(SOLVER).solve(m) - return soln - -if __name__ == '__main__': - m = make_concrete_model() - print_setup(m) - s = solve_model(m) - print_solution(m) - plot_solution(m) - plt.show() diff --git a/src/dispatch/twin_pyomo_test.py b/src/dispatch/twin_pyomo_test.py deleted file mode 100644 index 68d17e24..00000000 --- a/src/dispatch/twin_pyomo_test.py +++ /dev/null @@ -1,251 +0,0 @@ -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" -Enables quick testing of Pyomo problem solves, intending to set up a system -similar to how it is set up in the Pyomo dispatcher. This allows rapid testing -of different configurations for the rolling window optimization. -""" -import platform - -import numpy as np - -import pyomo.environ as pyo -from pyomo.opt import SolverStatus, TerminationCondition - -if platform.system() == 'Windows': - SOLVER = 'glpk' -else: - SOLVER = 'cbc' - -# setup stuff -components = ['steam_source', 'elec_generator', 'steam_storage', 'elec_sink'] -resources = ['steam', 'electricity'] -time = np.linspace(0, 10, 11) # from @1 to @2 in @3 steps -dt = time[1] - time[0] -resource_map = {'steam_source': {'steam': 0}, - 'elec_generator': {'steam': 0, 'electricity': 1}, - 'steam_storage': {'steam': 0}, - 'elec_sink': {'electricity': 0}, - } -activity = {} -for comp in components: - activity[comp] = np.zeros((len(resources), len(time)), dtype=float) - -# sizing specifications -storage_initial = 50 # kg of steam -storage_limit = 400 # kg of steam -steam_produced = 100 # kg/h of steam -gen_consume_limit = 110 # consumes at most 110 kg/h steam -sink_limit = 10000 # kWh/h = kW of electricity - -def make_concrete_model(): - """ - Test writing a simple concrete model with terms typical to the pyomo dispatcher. - @ In, None - @ Out, m, pyo.ConcreteModel, instance of the model to solve - """ - 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(0, len(time), dtype=int) # indexes time - # move onto model - m.C = pyo.Set(initialize=C) - m.R = pyo.Set(initialize=R) - m.T = pyo.Set(initialize=T) - # store some stuff for reference -> NOT NOTICED by Pyomo, we hope - m.Times = time - m.Components = components - m.resource_index_map = resource_map - m.Activity = activity - #******************* - # set up optimization variables - # -> for now we just do this manually - # steam_source - m.steam_source_index_map = pyo.Set(initialize=range(len(m.resource_index_map['steam_source']))) - m.steam_source_production = pyo.Var(m.steam_source_index_map, m.T, initialize=0) - # elec_generator - m.elec_generator_index_map = pyo.Set(initialize=range(len(m.resource_index_map['elec_generator']))) - m.elec_generator_production = pyo.Var(m.elec_generator_index_map, m.T, initialize=0) - # steam_storage - m.steam_storage_index_map = pyo.Set(initialize=range(len(m.resource_index_map['steam_storage']))) - m.steam_storage_production = pyo.Var(m.steam_storage_index_map, m.T, initialize=0) - # elec_sink - m.elec_sink_index_map = pyo.Set(initialize=range(len(m.resource_index_map['elec_sink']))) - m.elec_sink_production = pyo.Var(m.elec_sink_index_map, m.T, initialize=0) - #******************* - # set up lower, upper bounds - # -> for now we just do this manually - # -> consuming is negative sign by convention! - # -> producing is positive sign by convention! - # steam source produces exactly 100 steam - m.steam_source_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.steam_source_production[0, t] >= steam_produced) - m.steam_source_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.steam_source_production[0, t] <= steam_produced) - # elec generator can consume steam to produce electricity; 0 < consumed steam < 1000 - # -> this effectively limits electricity production, but we're defining capacity in steam terms for fun - # -> therefore signs are negative, -1000 < consumed steam < 0! - m.elec_generator_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.elec_generator_production[0, t] >= -gen_consume_limit) - m.elec_generator_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.elec_generator_production[0, t] <= 0) - # elec sink can take any amount of electricity - # -> consuming, so -10000 < consumed elec < 0 - m.elec_sink_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.elec_sink_production[0, t] >= -sink_limit) - m.elec_sink_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.elec_sink_production[0, t] <= 0) - # storage is in LEVEL not ACTIVITY (e.g. kg not kg/s) -> lets say it can store X kg - m.steam_storage_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.steam_storage_production[0, t] >= 0) - m.steam_storage_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.steam_storage_production[0, t] <= storage_limit) - #******************* - # create transfer function - # 2 steam make 1 electricity (sure, why not) - m.elec_generator_transfer = pyo.Constraint(m.T, rule=_generator_transfer) - #******************* - # create conservation rules - # steam - m.steam_conservation = pyo.Constraint(m.T, rule=_conserve_steam) - # electricity - m.elec_conservation = pyo.Constraint(m.T, rule=_conserve_electricity) - #******************* - # create objective function - m.OBJ = pyo.Objective(sense=pyo.maximize, rule=_economics) - ####### - # return - return m - -####### -# -# Callback Functions -# -def _generator_transfer(m, t): - """ - Constraint rule for electricity generation in generator - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - return - m.elec_generator_production[0, t] == 2.0 * m.elec_generator_production[1, t] - -def _conserve_steam(m, t): - """ - Constraint rule for conserving steam - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - # signs are tricky here, consumption is negative and production is positive - # a positive delta in steam storage level means it absorbed steam, so it's a negative term - storage_source = - (m.steam_storage_production[0, t] - (storage_initial if t == 0 else m.steam_storage_production[0, t-1])) / dt - sources = storage_source + m.steam_source_production[0, t] - sinks = m.elec_generator_production[0, t] - return sources + sinks == 0 - -def _conserve_electricity(m, t): - """ - Constraint rule for conserving electricity - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - sources = m.elec_generator_production[1, t] - sinks = m.elec_sink_production[0, t] - return sources + sinks == 0 - -def _economics(m): - """ - Constraint rule for optimization target - @ In, m, pyo.ConcreteModel, model containing problem - @ Out, objective, float, constraining evaluation - """ - opex = sum(m.elec_generator_production[0, t] for t in m.T) * 10 # will be negative b/c consumed - sales = - sum((m.elec_sink_production[0, t] * (100 if t < 5 else 1)) for t in m.T) # net positive because consumed - return opex + sales - -####### -# -# Debug printing functions -# -def print_setup(m): - """ - Debug printing for pre-solve model setup - @ In, m, pyo.ConcreteModel, model containing problem - @ 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 print_solution(m): - """ - Debug printing for post-solve model setup - @ In, m, pyo.ConcreteModel, model containing problem - @ Out, None - """ - print('') - print('*'*80) - print('solution:') - print(' objective value:', m.OBJ()) - print('time | steam source | steam storage | elec gen (s, e) | elec sink') - for t in m.T: - print(f'{m.Times[t]:1.2e} | ' + - f'{m.steam_source_production[0, t].value: 1.3e} | ' + - f'{m.steam_storage_production[0, t].value: 1.3e} | ' + - f'({m.elec_generator_production[0, t].value: 1.3e}, {m.elec_generator_production[1, t].value: 1.3e}) | ' + - f'{m.elec_sink_production[0, t].value: 1.3e}' - ) - print('*'*80) - -####### -# -# Solver. -# -def solve_model(m): - """ - Solves the model. - @ In, m, pyo.ConcreteModel, model containing problem - @ Out, m, pyo.ConcreteModel, results - """ - soln = pyo.SolverFactory(SOLVER).solve(m) - return soln - -if __name__ == '__main__': - m = make_concrete_model() - print_setup(m) - s = solve_model(m) - print_solution(m) - -# solution using setup: -# time = np.linspace(0, 10, 11) -# storage_initial = 50 # kg of steam -# storage_limit = 400 # kg of steam -# steam_produced = 100 # kg/h of steam -# gen_consume_limit = 110 # consumes at most 110 kg/h steam -# sink_limit = 10000 # kWh/h = kW of electricity -# 1 steam = 2 * electricity -# cost for generator = 10 * kg/h steam consumed -# profit = 100 * electricity consumed at sink, t < 5, -# 1 * electricity consumed at sink, t >= 5, -# -# should look like: -# ******************************************************************************** -# solution: -# objective value: 20100.0 -# time | steam source | steam storage | elec gen (s, e) | elec sink -# 0.00e+00 | 1.000e+02 | 4.000e+01 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 1.00e+00 | 1.000e+02 | 3.000e+01 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 2.00e+00 | 1.000e+02 | 2.000e+01 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 3.00e+00 | 1.000e+02 | 1.000e+01 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 4.00e+00 | 1.000e+02 | 0.000e+00 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 5.00e+00 | 1.000e+02 | 0.000e+00 | (-1.000e+02, 5.000e+01) | -5.000e+01 -# 6.00e+00 | 1.000e+02 | 0.000e+00 | (-1.000e+02, 5.000e+01) | -5.000e+01 -# 7.00e+00 | 1.000e+02 | 1.000e+02 | ( 0.000e+00, 0.000e+00) | 0.000e+00 -# 8.00e+00 | 1.000e+02 | 2.000e+02 | ( 0.000e+00, 0.000e+00) | 0.000e+00 -# 9.00e+00 | 1.000e+02 | 3.000e+02 | ( 0.000e+00, 0.000e+00) | 0.000e+00 -# 1.00e+01 | 1.000e+02 | 4.000e+02 | ( 0.000e+00, 0.000e+00) | 0.000e+00 -# ******************************************************************************** diff --git a/src/dispatch/twin_pyomo_test_rte.py b/src/dispatch/twin_pyomo_test_rte.py deleted file mode 100644 index 1b9cb147..00000000 --- a/src/dispatch/twin_pyomo_test_rte.py +++ /dev/null @@ -1,285 +0,0 @@ -# Copyright 2020, Battelle Energy Alliance, LLC -# ALL RIGHTS RESERVED -""" -Enables quick testing of Pyomo problem solves, intending to set up a system -similar to how it is set up in the Pyomo dispatcher. This allows rapid testing -of different configurations for the rolling window optimization. -""" -import platform - -import numpy as np - -import pyomo.environ as pyo -from pyomo.opt import SolverStatus, TerminationCondition - -if platform.system() == 'Windows': - SOLVER = 'glpk' -else: - SOLVER = 'cbc' - -# setup stuff -components = ['steam_source', 'elec_generator', 'steam_storage', 'elec_sink'] -resources = ['steam', 'electricity'] -time = np.linspace(0, 10, 11) # from @1 to @2 in @3 steps -dt = time[1] - time[0] -resource_map = {'steam_source': {'steam': 0}, - 'elec_generator': {'steam': 0, 'electricity': 1}, - 'steam_storage': {'steam': 0}, - 'elec_sink': {'electricity': 0}, - } -activity = {} -for comp in components: - activity[comp] = np.zeros((len(resources), len(time)), dtype=float) - -# sizing specifications -storage_initial = 50 # kg of steam -storage_limit = 400 # kg of steam -steam_produced = 100 # kg/h of steam -gen_consume_limit = 110 # consumes at most 110 kg/h steam -sink_limit = 10000 # kWh/h = kW of electricity -rte = 0.90 # storage round trip efficiency -sq_rte = np.sqrt(rte) # square root of rte - -def make_concrete_model(): - """ - Test writing a simple concrete model with terms typical to the pyomo dispatcher. - @ In, None - @ Out, m, pyo.ConcreteModel, instance of the model to solve - """ - global storage_initial - 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(0, len(time), dtype=int) # indexes time - # move onto model - m.C = pyo.Set(initialize=C) - m.R = pyo.Set(initialize=R) - m.T = pyo.Set(initialize=T) - # store some stuff for reference -> NOT NOTICED by Pyomo, we hope - m.Times = time - m.Components = components - m.resource_index_map = resource_map - m.Activity = activity - #******************* - # set up optimization variables - # -> for now we just do this manually - # steam_source - m.steam_source_index_map = pyo.Set(initialize=range(len(m.resource_index_map['steam_source']))) - m.steam_source_production = pyo.Var(m.steam_source_index_map, m.T, initialize=steam_produced, bounds=(steam_produced, steam_produced)) - # elec_generator - m.elec_generator_index_map = pyo.Set(initialize=range(len(m.resource_index_map['elec_generator']))) - bounds = lambda m, r, t: ({0: -gen_consume_limit}.get(r, None), {0: 0}.get(r, None)) - m.elec_generator_production = pyo.Var(m.elec_generator_index_map, m.T, initialize=0, bounds=bounds) - # steam_storage - m.steam_storage_index_map = pyo.Set(initialize=range(len(m.resource_index_map['steam_storage']))) - m.steam_storage_level = pyo.Var(m.steam_storage_index_map, m.T, initialize=0, bounds=(0, storage_limit)) - m.steam_storage_charge = pyo.Var(m.steam_storage_index_map, m.T, initialize=0, bounds=(-storage_limit/dt, 0)) # TODO should be ramp rate limit - m.steam_storage_discharge = pyo.Var(m.steam_storage_index_map, m.T, initialize=0, bounds=(0, storage_limit/dt)) # TODO ramp rate limit - m.steam_oneway_bigM = pyo.Var(m.T, initialize=0, within=pyo.Binary) # 0 is charging, discharging is 1 - large_eps = storage_limit / dt * 1.1 - m.steam_oneway_charge = pyo.Constraint(m.T, rule=lambda m, t: -m.steam_storage_charge[0, t] <= (1 - m.steam_oneway_bigM[t]) * large_eps) - m.steam_oneway_discharge = pyo.Constraint(m.T, rule=lambda m, t: m.steam_storage_discharge[0, t] <= m.steam_oneway_bigM[t] * large_eps) - # OLD m.steam_storage_production = pyo.Var(m.steam_storage_index_map, m.T, initialize=0, bounds=(0, storage_limit)) - # elec_sink - m.elec_sink_index_map = pyo.Set(initialize=range(len(m.resource_index_map['elec_sink']))) - m.elec_sink_production = pyo.Var(m.elec_sink_index_map, m.T, initialize=0, bounds=(-sink_limit, 0)) - #******************* - # apply a periodic condition - storage_initial = m.steam_storage_level[(0, m.T[-1])] - #******************* - # set up lower, upper bounds - # -> for now we just do this manually - # -> consuming is negative sign by convention! - # -> producing is positive sign by convention! - # steam source produces exactly 100 steam - # OLD m.steam_source_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.steam_source_production[0, t] >= steam_produced) - # OLD m.steam_source_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.steam_source_production[0, t] <= steam_produced) - # elec generator can consume steam to produce electricity; 0 < consumed steam < 1000 - # -> this effectively limits electricity production, but we're defining capacity in steam terms for fun - # -> therefore signs are negative, -1000 < consumed steam < 0! - # OLD m.elec_generator_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.elec_generator_production[0, t] >= -gen_consume_limit) - # OLD m.elec_generator_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.elec_generator_production[0, t] <= 0) - # elec sink can take any amount of electricity - # -> consuming, so -10000 < consumed elec < 0 - # OLD m.elec_sink_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.elec_sink_production[0, t] >= -sink_limit) - # OLD m.elec_sink_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.elec_sink_production[0, t] <= 0) - # storage is in LEVEL not ACTIVITY (e.g. kg not kg/s) -> lets say it can store X kg - # OLD m.steam_storage_lower_limit = pyo.Constraint(m.T, rule=lambda m, t: m.steam_storage_production[0, t] >= 0) - # OLD m.steam_storage_upper_limit = pyo.Constraint(m.T, rule=lambda m, t: m.steam_storage_production[0, t] <= storage_limit) - #******************* - # create transfer function - # 2 steam make 1 electricity (sure, why not) - m.elec_generator_transfer = pyo.Constraint(m.T, rule=_generator_transfer) - #******************* - # create conservation rules - # steam - m.steam_conservation = pyo.Constraint(m.T, rule=_conserve_steam) - # electricity - m.elec_conservation = pyo.Constraint(m.T, rule=_conserve_electricity) - # track storage level by charge/discharge - m.steam_storage_management = pyo.Constraint(m.T, rule=_manage_storage) - #******************* - # create objective function - m.OBJ = pyo.Objective(sense=pyo.maximize, rule=_economics) - ####### - # return - return m - -####### -# -# Callback Functions -# -def _generator_transfer(m, t): - """ - Constraint rule for electricity generation in generator - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - return - m.elec_generator_production[0, t] == 2.0 * m.elec_generator_production[1, t] - -def _conserve_steam(m, t): - """ - Constraint rule for conserving steam - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - # signs are tricky here, consumption is negative and production is positive - # a positive delta in steam storage level means it absorbed steam, so it's a negative term - storage_source = m.steam_storage_charge[0, t] + m.steam_storage_discharge[0, t] #(storage_initial if t == 0 else m.steam_storage_level[0, t-1])) / dt - sources = storage_source + m.steam_source_production[0, t] - sinks = m.elec_generator_production[0, t] - return sources + sinks == 0 - -def _conserve_electricity(m, t): - """ - Constraint rule for conserving electricity - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - sources = m.elec_generator_production[1, t] - sinks = m.elec_sink_production[0, t] - return sources + sinks == 0 - -def _manage_storage(m, t): - """ - Balance level of the storage unit with charge, discharge activities - @ In, m, pyo.ConcreteModel, model containing problem - @ In, t, int, time indexer - @ Out, constraint, bool, constraining evaluation - """ - if t > 0: - previous = m.steam_storage_level[0, t-1] - else: - previous = storage_initial - # recall charge is negative (absorbing) and discharge is positive (emitting) - production = -sq_rte * m.steam_storage_charge[0, t] - m.steam_storage_discharge[0, t] / sq_rte - return m.steam_storage_level[0, t] == previous + production * dt - -def _economics(m): - """ - Constraint rule for optimization target - @ In, m, pyo.ConcreteModel, model containing problem - @ Out, objective, float, constraining evaluation - """ - opex = sum(m.elec_generator_production[0, t] for t in m.T) * 10 # will be negative b/c consumed - sales = - sum((m.elec_sink_production[0, t] * (100 if t < 5 else 1)) for t in m.T) # net positive because consumed - return opex + sales - -####### -# -# Debug printing functions -# -def print_setup(m): - """ - Debug printing for pre-solve model setup - @ In, m, pyo.ConcreteModel, model containing problem - @ 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 print_solution(m): - """ - Debug printing for post-solve model setup - @ In, m, pyo.ConcreteModel, model containing problem - @ Out, None - """ - print('') - print('*'*80) - print('solution:') - print(' objective value:', m.OBJ()) - print(' storage initial:', storage_initial) - print(' dt:', dt) - print(f'{"time":^8} | {"steam src":^10} | {"storage c, d, l":^34} | {"elec gen (s, e)":^22} | elec sink | bigM') - for t in m.T: - print(f'{m.Times[t]:1.2e} | ' + - f'{m.steam_source_production[0, t].value: 1.3e} | ' + - f'{m.steam_storage_charge[0, t].value: 1.3e}, {m.steam_storage_discharge[0, t].value: 1.3e}, {m.steam_storage_level[0, t].value: 1.3e} | ' + - f'{m.elec_generator_production[0, t].value: 1.3e}, {m.elec_generator_production[1, t].value: 1.3e} | ' + - f'{m.elec_sink_production[0, t].value: 1.3e} | ' + - f'{m.steam_oneway_bigM[t].value} ' - ) - print('*'*80) - -####### -# -# Solver. -# -def solve_model(m): - """ - Solves the model. - @ In, m, pyo.ConcreteModel, model containing problem - @ Out, m, pyo.ConcreteModel, results - """ - soln = pyo.SolverFactory(SOLVER).solve(m) - return soln - -if __name__ == '__main__': - m = make_concrete_model() - print_setup(m) - s = solve_model(m) - print_solution(m) - -# solution using setup: -# time = np.linspace(0, 10, 11) -# storage_initial = 50 # kg of steam -# storage_limit = 400 # kg of steam -# steam_produced = 100 # kg/h of steam -# gen_consume_limit = 110 # consumes at most 110 kg/h steam -# sink_limit = 10000 # kWh/h = kW of electricity -# 1 steam = 2 * electricity -# cost for generator = 10 * kg/h steam consumed -# profit = 100 * electricity consumed at sink, t < 5, -# 1 * electricity consumed at sink, t >= 5, -# -# should look like: -# ******************************************************************************** -# solution: -# objective value: 20100.0 -# time | steam source | steam storage | elec gen (s, e) | elec sink -# 0.00e+00 | 1.000e+02 | 4.000e+01 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 1.00e+00 | 1.000e+02 | 3.000e+01 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 2.00e+00 | 1.000e+02 | 2.000e+01 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 3.00e+00 | 1.000e+02 | 1.000e+01 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 4.00e+00 | 1.000e+02 | 0.000e+00 | (-1.100e+02, 5.500e+01) | -5.500e+01 -# 5.00e+00 | 1.000e+02 | 0.000e+00 | (-1.000e+02, 5.000e+01) | -5.000e+01 -# 6.00e+00 | 1.000e+02 | 0.000e+00 | (-1.000e+02, 5.000e+01) | -5.000e+01 -# 7.00e+00 | 1.000e+02 | 1.000e+02 | ( 0.000e+00, 0.000e+00) | 0.000e+00 -# 8.00e+00 | 1.000e+02 | 2.000e+02 | ( 0.000e+00, 0.000e+00) | 0.000e+00 -# 9.00e+00 | 1.000e+02 | 3.000e+02 | ( 0.000e+00, 0.000e+00) | 0.000e+00 -# 1.00e+01 | 1.000e+02 | 4.000e+02 | ( 0.000e+00, 0.000e+00) | 0.000e+00 -# ******************************************************************************** diff --git a/src/input_loader.py b/src/input_loader.py index 7e295c18..12dbd191 100644 --- a/src/input_loader.py +++ b/src/input_loader.py @@ -56,9 +56,9 @@ def parse(xml, loc, messageHandler): if components_node is None: raise IOError(' node is missing from HERON input file!') for comp_xml in components_node: - comp = Components.Component(messageHandler=messageHandler) + comp = Components.HeronComponent(messageHandler=messageHandler) # check parsing - comp.read_input(comp_xml, case.get_mode()) + comp.read_input(comp_xml, )#case.get_mode()) components.append(comp) if not components: raise IOError('No nodes were found in the section!') diff --git a/templates/template_driver.py b/templates/template_driver.py index 54046104..b7fc0235 100644 --- a/templates/template_driver.py +++ b/templates/template_driver.py @@ -1196,7 +1196,7 @@ def _modify_inner_components(self, template, case, components): ## For each interaction of each component, that means making sure the Function, ARMA, or constant makes it. ## Constants from outer (namely sweep/opt capacities) are set in the MC Sampler from the outer ## The Dispatch needs info from the Outer to know which capacity to use, so we can't pass it from here. - capacity = component.get_capacity(None, raw=True) + capacity = component.get_capacity()#None, raw=True) interaction = component.get_interaction() parametric = capacity.is_parametric()