diff --git a/examples/30_pyomo_optimized_dispatch/plant_config.yaml b/examples/30_pyomo_optimized_dispatch/plant_config.yaml index ea111f4f8..dc72703a1 100644 --- a/examples/30_pyomo_optimized_dispatch/plant_config.yaml +++ b/examples/30_pyomo_optimized_dispatch/plant_config.yaml @@ -62,6 +62,7 @@ finance_parameters: all_electricity: commodity: electricity finance_groups: [profast_model] + # technologies: [wind, battery, grid_buy] technologies: [wind, battery] # dispatched_electricity: diff --git a/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py b/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py index accec7347..0de525d6f 100644 --- a/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py +++ b/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py @@ -7,12 +7,44 @@ # Create an H2Integrate model model = H2IntegrateModel("pyomo_optimized_dispatch.yaml") +# --- Parameters --- +amplitude = 0.9 # Amplitude of the sine wave +frequency = 0.05 # Frequency of the sine wave in Hz +duration = 8760.0 # Duration of the signal in seconds +sampling_rate = 1 # Number of samples per second (Fs) + +# Noise parameters +noise_mean = 0.0 +noise_std_dev = 0.1 # Standard deviation controls the noise intensity + +# --- Generate the Time Vector --- +# Create a time array from 0 to duration with a specific sampling rate +t = np.linspace(1.0, duration, int(sampling_rate * duration), endpoint=True) + +# --- Generate the Pure Sine Wave Signal --- +# Formula: y(t) = A * sin(2 * pi * f * t) +pure_signal = amplitude * np.sin(2.0 * np.pi * frequency * t) + +# --- Generate the Random Gaussian Noise --- +# Create noise with the same shape as the time vector +rng = np.random.default_rng() +noise = rng.normal(loc=noise_mean, scale=noise_std_dev, size=t.shape) + +# --- Create the Noisy Signal --- +noisy_signal = (pure_signal + noise) * 0.04 + 0.04 * np.ones(len(t)) + +commodity_met_value_profile = np.ones(8760) * 1 +commodity_buy_price_profile = noisy_signal + demand_profile = np.ones(8760) * 100.0 # TODO: Update with demand module once it is developed model.setup() model.prob.set_val("battery.electricity_demand", demand_profile, units="MW") +# model.prob.set_val("battery.demand_met_value", commodity_met_value_profile, units="USD/kW") +model.prob.set_val("battery.electricity_buy_price", commodity_buy_price_profile, units="USD/kW") +# model.prob.set_val("grid_buy.electricity_buy_price", commodity_buy_price_profile, units="USD/kW") # Run the model model.run() @@ -20,7 +52,7 @@ model.post_process() # Plot the results -fig, ax = plt.subplots(2, 1, sharex=True, figsize=(8, 6)) +fig, ax = plt.subplots(3, 1, sharex=True, figsize=(8, 6)) start_hour = 0 end_hour = 200 @@ -61,10 +93,19 @@ ) ax[1].plot( range(start_hour, end_hour), - model.prob.get_val("battery.battery_electricity", units="MW")[start_hour:end_hour], + model.prob.get_val("battery.electricity_bought_for_storage", units="MW")[start_hour:end_hour], + linestyle="-", + label="Electricity Bought (MW)", +) +ax[1].plot( + range(start_hour, end_hour), + model.prob.get_val("battery.battery_electricity_out", units="MW")[start_hour:end_hour], linestyle="-.", label="Battery Electricity Out (MW)", ) + +print(min(model.prob.get_val("battery.electricity_out", units="MW"))) +print(min(model.prob.get_val("battery.battery_electricity_out", units="MW"))) ax[1].plot( range(start_hour, end_hour), demand_profile[start_hour:end_hour], @@ -73,7 +114,15 @@ ) ax[1].set_ylim([-1e2, 2.5e2]) ax[1].set_ylabel("Electricity Hourly (MW)") -ax[1].set_xlabel("Timestep (hr)") +ax[1].legend() + +ax[2].plot( + range(start_hour, end_hour), + model.prob.get_val("battery.electricity_buy_price", units="USD/MW")[start_hour:end_hour], + label="Grid Purchase Price ($/MW)", +) +ax[2].set_ylabel("Grid Purchase Price ($/MW)") +ax[2].set_xlabel("Timestep (hr)") plt.legend(ncol=2, frameon=False) plt.tight_layout() diff --git a/examples/30_pyomo_optimized_dispatch/tech_config.yaml b/examples/30_pyomo_optimized_dispatch/tech_config.yaml index abe1b397a..c96565d61 100644 --- a/examples/30_pyomo_optimized_dispatch/tech_config.yaml +++ b/examples/30_pyomo_optimized_dispatch/tech_config.yaml @@ -62,12 +62,14 @@ technologies: charge_efficiency: 0.95 # Charge efficiency of the storage technology discharge_efficiency: 0.95 # Discharge efficiency of the storage technology cost_per_charge: 0.03 # in $/kW, cost to charge the storage (note that charging is incentivized) - cost_per_discharge: 0.05 # in $/kW, cost to discharge the storage - commodity_met_value: 0.1 # in $/kW, penalty for not meeting the desired load demand + cost_per_discharge: 0.07 # in $/kW, cost to discharge the storage + demand_met_value: 0.1 # in $/kW, penalty for not meeting the desired load demand cost_per_production: 0.0 # in $/kW, cost to use the incoming produced commodity (i.e. electricity from wind) time_weighting_factor: 0.995 # This parameter discounts each subsequent time step incrementally in the future in the horizon window by this amount n_control_window: 24 # in timesteps (currently hours), The length of time that the control is applied to in the rolling window optimization - system_commodity_interface_limit: 1e12 + system_commodity_interface_limit: 207500 # rating of the wind farm + allow_commodity_buying: true + commodity_buy_price: 0.4 elec_combiner: performance_model: model: GenericCombinerPerformanceModel diff --git a/h2integrate/control/control_rules/converters/generic_converter_min_operating_cost.py b/h2integrate/control/control_rules/converters/generic_converter_min_operating_cost.py index 56cd25d47..83645b06f 100644 --- a/h2integrate/control/control_rules/converters/generic_converter_min_operating_cost.py +++ b/h2integrate/control/control_rules/converters/generic_converter_min_operating_cost.py @@ -65,8 +65,12 @@ def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): Args: inputs (dict): Dictionary of numpy arrays (length = self.n_timesteps) containing at least: - f"{commodity}_in" : Available generated commodity profile. - f"{commodity}_demand" : Demanded commodity output profile. + f"{commodity}_in" : Available generated commodity profile. + f"{commodity}_demand" : Demanded commodity output profile. + "demand_met_value" : Variable weight for meeting the load + if allow_commodity_buying: + f"{commodity}_buy_price" : Variable cost of buying commodity for charging + (e.g. could be a grid price) dispatch_inputs (dict): Dictionary of the dispatch input parameters from config """ @@ -183,18 +187,26 @@ def _create_constraints(self, pyomo_model: pyo.ConcreteModel): # Update time series parameters for next optimization window def update_time_series_parameters( - self, commodity_in: list, commodity_demand: list, updated_initial_soc: float + self, + time_update_inputs: dict, + updated_initial_soc: float, ): """Updates the pyomo optimization problem with parameters that change with time Args: commodity_in (list): List of generated commodity in for this time slice. commodity_demand (list): The demanded commodity for this time slice. + demand_met_value (list): List of variable value of meeting the provided load updated_initial_soc (float): The updated initial state of charge for storage technologies for the current time slice. + if allow_commodity_buying: + commodity_buy_price (list): List of variable cost of buying commodity for charging + (e.g. could be a grid price) """ self.time_duration = [1.0] * len(self.blocks.index_set()) - self.available_production = [commodity_in[t] for t in self.blocks.index_set()] + self.available_production = [ + time_update_inputs[f"{self.commodity_name}_in"][t] for t in self.blocks.index_set() + ] # Objective functions def min_operating_cost_objective(self, hybrid_blocks, tech_name: str): diff --git a/h2integrate/control/control_rules/plant_dispatch_model.py b/h2integrate/control/control_rules/plant_dispatch_model.py index a28f29007..3b7caf5e8 100644 --- a/h2integrate/control/control_rules/plant_dispatch_model.py +++ b/h2integrate/control/control_rules/plant_dispatch_model.py @@ -189,9 +189,7 @@ def arc_rule(m, t): pyo.TransformationFactory("network.expand_arcs").apply_to(self.model) - def update_time_series_parameters( - self, commodity_in=list, commodity_demand=list, updated_initial_soc=float - ): + def update_time_series_parameters(self, time_update_inputs=dict, updated_initial_soc=float): """ Updates the pyomo optimization problem with parameters that change with time @@ -206,9 +204,7 @@ def update_time_series_parameters( for tech in self.source_techs: name = tech + "_rule" pyomo_block = self.tech_dispatch_models.__getattribute__(name) - pyomo_block.update_time_series_parameters( - commodity_in, commodity_demand, updated_initial_soc - ) + pyomo_block.update_time_series_parameters(time_update_inputs, updated_initial_soc) def create_min_operating_cost_expression(self): """ @@ -252,3 +248,9 @@ def storage_commodity_out(self) -> list: self.blocks[t].discharge_commodity.value - self.blocks[t].charge_commodity.value for t in self.blocks.index_set() ] + + @property + def storage_commodity_bought(self) -> list: + # This is used in the storage_dispatch_commands method of the control strategy + """Storage commodity bought.""" + return [self.blocks[t].commodity_bought.value for t in self.blocks.index_set()] diff --git a/h2integrate/control/control_rules/storage/pyomo_storage_rule_min_operating_cost.py b/h2integrate/control/control_rules/storage/pyomo_storage_rule_min_operating_cost.py index 6faf5b513..42a29257f 100644 --- a/h2integrate/control/control_rules/storage/pyomo_storage_rule_min_operating_cost.py +++ b/h2integrate/control/control_rules/storage/pyomo_storage_rule_min_operating_cost.py @@ -25,6 +25,7 @@ def __init__( index_set: pyo.Set, round_digits: int, time_duration: float, + allow_commodity_buying: bool, block_set_name: str = "storage", ): # Set the number of digits to round to in the Pyomo model @@ -35,6 +36,7 @@ def __init__( # names and units in the Pyomo model self.commodity_name = commodity_info["commodity_name"] self.commodity_storage_units = commodity_info["commodity_storage_units"] + # This loads the currency unit definition into pyomo pyo.units.load_definitions_from_strings(["USD = [currency]"]) @@ -47,6 +49,9 @@ def __init__( self.amount_units_pyo = eval(amount_units_pyo_str) self.cost_units_per_amount_pyo = eval(f"pyo.units.USD / ({amount_units_pyo_str})") + # Define whether storage is allowed to buy commodity for charging + self.allow_commodity_buying = allow_commodity_buying + # The Pyomo model that this class builds off of, where all of the variables, parameters, # constraints, and ports will be added to. self.model = pyomo_model @@ -68,17 +73,25 @@ def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): Args: inputs (dict): Dictionary of numpy arrays (length = self.n_timesteps) containing at least: - f"{commodity}_in" : Available generated commodity profile. - f"{commodity}_demand" : Demanded commodity output profile. + f"{commodity}_in" : Available generated commodity profile. + f"{commodity}_demand" : Demanded commodity output profile. + "demand_met_value" : Variable weight for meeting the load + if allow_commodity_buying: + f"{commodity}_buy_price" : Variable cost of buying commodity for charging + (e.g. could be a grid price) dispatch_inputs (dict): Dictionary of the dispatch input parameters from config """ commodity_demand = inputs[f"{self.commodity_name}_demand"] + if self.allow_commodity_buying: + commodity_buy_price_in = inputs[f"{self.commodity_name}_buy_price"] + + demand_met_value_in = inputs["demand_met_value"] + # Dispatch Parameters self.set_timeseries_parameter("cost_per_charge", dispatch_inputs["cost_per_charge"]) self.set_timeseries_parameter("cost_per_discharge", dispatch_inputs["cost_per_discharge"]) - self.set_timeseries_parameter("commodity_met_value", dispatch_inputs["commodity_met_value"]) # Storage parameters self.set_timeseries_parameter("minimum_storage", 0.0) @@ -96,7 +109,15 @@ def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): self.set_timeseries_parameter("max_discharge", dispatch_inputs["max_charge_rate"]) # System parameters - self.commodity_load_demand = [commodity_demand[t] for t in self.blocks.index_set()] + self.set_timeseries_from_list("commodity_load_demand", commodity_demand) + self.set_timeseries_from_list("demand_met_value", demand_met_value_in) + # self.commodity_load_demand = [commodity_demand[t] for t in self.blocks.index_set()] + if self.allow_commodity_buying: + # This preserves the possibility of a variable interconnect limit + self.set_timeseries_parameter( + "commodity_import_limit", dispatch_inputs["commodity_import_limit"] + ) + self.set_timeseries_from_list("commodity_buy_price", commodity_buy_price_in) self._set_initial_soc_constraint() @@ -236,7 +257,7 @@ def _create_parameters(self, pyomo_model: pyo.ConcreteModel, t): mutable=True, units=pyo.units.USD, ) - pyomo_model.commodity_met_value = pyo.Param( + pyomo_model.demand_met_value = pyo.Param( doc=f"Commodity demand met value per generation [$/{self.commodity_storage_units}]", default=0.0, within=pyo.Reals, @@ -250,6 +271,23 @@ def _create_parameters(self, pyomo_model: pyo.ConcreteModel, t): mutable=True, units=self.rate_units_pyo, ) + # Add parameters that allow for buying commodity for storage charging + if self.allow_commodity_buying: + pyomo_model.commodity_buy_price = pyo.Param( + doc=f"Commodity buy price per generation [$/{self.commodity_storage_units}]", + default=0.0, + within=pyo.Reals, + mutable=True, + units=self.cost_units_per_amount_pyo, + ) + pyomo_model.commodity_import_limit = pyo.Param( + doc=f"External commodity source import limit \ + [$/{self.commodity_storage_units}]", + default=1000.0, + within=pyo.NonNegativeReals, + mutable=True, + units=self.rate_units_pyo, + ) def _create_variables(self, pyomo_model: pyo.ConcreteModel, t): """Create storage-related decision variables in the Pyomo model. @@ -324,6 +362,14 @@ def _create_variables(self, pyomo_model: pyo.ConcreteModel, t): domain=pyo.Binary, units=pyo.units.dimensionless, ) + # Add variables that allow for buying commodity for storage charging + if self.allow_commodity_buying: + pyomo_model.commodity_bought = pyo.Var( + doc=f"Commodity bought for the system [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + bounds=(0, pyomo_model.commodity_import_limit), + units=self.rate_units_pyo, + ) def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): """Create operational and state-of-charge constraints for storage and the system. @@ -341,7 +387,7 @@ def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): t: Time index or iterable representing time steps (unused in this method). """ ################################## - # Charging Constraints # + # Storage Constraints # ################################## # Charge commodity bounds pyomo_model.charge_commodity_ub = pyo.Constraint( @@ -372,17 +418,37 @@ def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): ################################## # System constraints # ################################## - pyomo_model.balance = pyo.Constraint( - doc="Transmission energy balance", - expr=( - pyomo_model.commodity_out == pyomo_model.system_production - pyomo_model.system_load - ), - ) + if self.allow_commodity_buying: + pyomo_model.balance = pyo.Constraint( + doc="Commodity system balance", + expr=( + pyomo_model.commodity_out - pyomo_model.commodity_bought + == pyomo_model.system_production - pyomo_model.system_load + ), + ) + else: + pyomo_model.balance = pyo.Constraint( + doc="Commodity system balance", + expr=( + pyomo_model.commodity_out + == pyomo_model.system_production - pyomo_model.system_load + ), + ) pyomo_model.production_limit = pyo.Constraint( - doc="Transmission limit on electricity sales", + doc="Upper limit on system production", expr=pyomo_model.commodity_out <= pyomo_model.commodity_load_demand * pyomo_model.is_generating, ) + # Add variables that allow for buying commodity for storage charging + # NOTE: This constraint prevents buying and selling the commodity at the same time + if self.allow_commodity_buying: + pyomo_model.purchase_limit = pyo.Constraint( + doc="Upper limit on commodity purchases", + expr=( + pyomo_model.commodity_bought + <= pyomo_model.commodity_import_limit * (1 - pyomo_model.is_generating) + ), + ) ################################## # SOC Inventory Constraints # @@ -466,23 +532,47 @@ def _create_ports(self, pyomo_model: pyo.ConcreteModel, t): pyomo_model.port.add(pyomo_model.system_production) pyomo_model.port.add(pyomo_model.system_load) pyomo_model.port.add(pyomo_model.commodity_out) + # Add variables that allow for buying commodity for storage charging + if self.allow_commodity_buying: + pyomo_model.port.add(pyomo_model.commodity_bought) # Update time series parameters for next optimization window def update_time_series_parameters( - self, commodity_in: list, commodity_demand: list, updated_initial_soc: float + self, + time_update_inputs: dict, + updated_initial_soc: float, ): """Updates the pyomo optimization problem with parameters that change with time Args: commodity_in (list): List of generated commodity in for this time slice. commodity_demand (list): The demanded commodity for this time slice. + demand_met_value (list): List of variable value of meeting the provided load updated_initial_soc (float): The updated initial state of charge for storage technologies for the current time slice. + if allow_commodity_buying: + commodity_buy_price (list): List of variable cost of buying commodity for charging + (e.g. could be a grid price) """ + # TODO: Standardize the inputs for this function self.time_duration = [1.0] * len(self.blocks.index_set()) - self.commodity_load_demand = [commodity_demand[t] for t in self.blocks.index_set()] + self.set_timeseries_from_list( + "commodity_load_demand", time_update_inputs[f"{self.commodity_name}_demand"] + ) + # self.commodity_load_demand = [ + # time_update_inputs[f"{self.commodity_name}_demand"][t] for t in self.blocks.index_set() + # ] + self.demand_met_value = [ + time_update_inputs["demand_met_value"][t] for t in self.blocks.index_set() + ] self.model.initial_soc = updated_initial_soc self.initial_soc = updated_initial_soc + # Add variables that allow for buying commodity for storage charging + if self.allow_commodity_buying: + self.commodity_buy_price = [ + time_update_inputs[f"{self.commodity_name}_buy_price"][t] + for t in self.blocks.index_set() + ] # Objective functions def min_operating_cost_objective(self, hybrid_blocks, tech_name: str): @@ -505,10 +595,17 @@ def min_operating_cost_objective(self, hybrid_blocks, tech_name: str): self.blocks[t].cost_per_discharge * hybrid_blocks[t].discharge_commodity - self.blocks[t].cost_per_charge * hybrid_blocks[t].charge_commodity + (self.blocks[t].commodity_load_demand - hybrid_blocks[t].commodity_out) - * self.blocks[t].commodity_met_value + * self.blocks[t].demand_met_value ) for t in self.blocks.index_set() ) + if self.allow_commodity_buying: + self.obj += sum( + hybrid_blocks[t].time_weighting_factor + * self.blocks[t].time_duration + * (hybrid_blocks[t].commodity_bought * self.blocks[t].commodity_buy_price) + for t in self.blocks.index_set() + ) return self.obj # System-level functions @@ -520,15 +617,18 @@ def _create_hybrid_port(self, hybrid_model: pyo.ConcreteModel, tech_name: str): tech_name (str): The name or key identifying the technology for which ports are created. """ - tech_port = Port( - initialize={ - "system_production": hybrid_model.system_production, - "system_load": hybrid_model.system_load, - "commodity_out": hybrid_model.commodity_out, - "charge_commodity": hybrid_model.charge_commodity, - "discharge_commodity": hybrid_model.discharge_commodity, - } - ) + initialize_dict = { + "system_production": hybrid_model.system_production, + "system_load": hybrid_model.system_load, + "commodity_out": hybrid_model.commodity_out, + "charge_commodity": hybrid_model.charge_commodity, + "discharge_commodity": hybrid_model.discharge_commodity, + } + if self.allow_commodity_buying: + initialize_dict["commodity_bought"] = hybrid_model.commodity_bought + + tech_port = Port(initialize=initialize_dict) + hybrid_model.__setattr__(f"{tech_name}_port", tech_port) return hybrid_model.__getattribute__(f"{tech_name}_port") @@ -560,6 +660,13 @@ def _create_hybrid_variables(self, hybrid_model: pyo.ConcreteModel, tech_name: s domain=pyo.NonNegativeReals, units=self.rate_units_pyo, ) + # Add variables that allow for buying commodity for storage charging + if self.allow_commodity_buying: + hybrid_model.commodity_bought = pyo.Var( + doc=f"{self.commodity_name} bought [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=self.rate_units_pyo, + ) ################################## # Storage Variables # ################################## @@ -613,6 +720,14 @@ def set_timeseries_parameter(self, param_name: str, param_val: float): val_rounded = round(param_val, self.round_digits) self.blocks[t].__setattr__(param_name, val_rounded) + def set_timeseries_from_list(self, param_name: str, param_list: list): + if len(param_list) == len(self.blocks): + for t, param_val in zip(self.blocks, param_list): + val_rounded = round(param_val, self.round_digits) + self.blocks[t].__setattr__(param_name, val_rounded) + else: + raise ValueError("param_list must be the same length as time horizon") + @property def charge_efficiency(self) -> float: """Charge efficiency.""" diff --git a/h2integrate/control/control_strategies/storage/optimized_load_following_with_grid.py b/h2integrate/control/control_strategies/storage/optimized_load_following_with_grid.py new file mode 100644 index 000000000..9808ceaa2 --- /dev/null +++ b/h2integrate/control/control_strategies/storage/optimized_load_following_with_grid.py @@ -0,0 +1,549 @@ +from typing import TYPE_CHECKING + +import numpy as np +import pyomo.environ as pyomo +from attrs import field, define +from pyomo.util.check_units import assert_units_consistent + +from h2integrate.core.utilities import merge_shared_inputs +from h2integrate.core.validators import range_val +from h2integrate.control.control_rules.plant_dispatch_model import PyomoDispatchPlantModel +from h2integrate.control.control_strategies.pyomo_controller_baseclass import ( + SolverOptions, + PyomoControllerBaseClass, + PyomoControllerBaseConfig, +) +from h2integrate.control.control_strategies.controller_opt_problem_state import DispatchProblemState +from h2integrate.control.control_rules.storage.pyomo_storage_rule_min_operating_cost import ( + PyomoRuleStorageMinOperatingCosts, +) +from h2integrate.control.control_rules.converters.generic_converter_min_operating_cost import ( + PyomoDispatchGenericConverterMinOperatingCosts, +) + + +if TYPE_CHECKING: # to avoid circular imports + pass + + +@define +class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): + """ + Configuration data container for Pyomo-based optimal dispatch. + + This class groups the parameters needed by the optimized dispatch controller. + Values are typically populated from the technology + `tech_config.yaml` (merged under the "control" section). + + Attributes: + max_charge_rate (float): + The maximum charge that the storage can accept + (in units of the commodity per time step). + charge_efficiency (float): + The efficiency of charging the storage (between 0 and 1). + discharge_efficiency (float): + The efficiency of discharging the storage (between 0 and 1). + commodity (str): + The name of the commodity being stored (e.g., "electricity", "hydrogen"). + commodity_rate_units (str): + The rate units of the commodity being stored (e.g., "kW", "kg/h"). + cost_per_production (float): + The cost to use the incoming produced commodity (in $/commodity_rate_units). + cost_per_charge (float): + The cost per unit of charging the storage (in $/commodity_rate_units). + cost_per_discharge (float): + The cost per unit of discharging the storage (in $/commodity_rate_units). + demand_met_value (int, float, list): + The penalty for not meeting the desired load demand (in $/commodity_rate_units). + time_weighting_factor (float): + The weighting factor applied to future time steps in the optimization objective + (between 0 and 1). + time_duration (float): + The duration of each time step in the Pyomo model in hours. + The default of this parameter is 1.0 (i.e., 1 hour time steps). + commodity_import_limit (float): + Maximum amount of the commodity that the storage/system can buy + (i.e. transmission limit for grid charging) + allow_commodity_buying (bool): + This sets whether the storage can buy commodity to charge or not + commodity_buy_price (int, float, list): + Price of the commodity that can be bought for storage. + """ + + max_charge_rate: int | float = field() + allow_commodity_buying: bool = field() + charge_efficiency: float = field(validator=range_val(0, 1), default=None) + discharge_efficiency: float = field(validator=range_val(0, 1), default=None) + # TODO: note that this definition of cost_per_production is not generalizable to multiple + # production technologies. Would need a name adjustment to connect it to + # production tech + cost_per_production: float = field(default=None) + cost_per_charge: float = field(default=None) + cost_per_discharge: float = field(default=None) + demand_met_value: int | float | list = field(default=None) + time_weighting_factor: float = field(validator=range_val(0, 1), default=0.995) + time_duration: float = field(default=1.0) # hours + # Can we set this to interconnection? do we want to? + commodity_import_limit: float = field(default=None) + commodity_buy_price: int | float | list = field(default=None) + + def __attrs_post_init__(self): + # Check inputs for commodity buying parameters + if self.allow_commodity_buying: + if self.commodity_buy_price: + # Check commodity buy price + if isinstance(self.commodity_buy_price, float | int): + if self.commodity_buy_price == 0: + raise ValueError( + "commodity_buy_price must be defined as an input and >0 \ + if using commodity buying" + ) + if isinstance(self.commodity_buy_price, list) or self.commodity_buy_price is None: + if all(self.commodity_buy_price) == 0: + raise ValueError( + "commodity_buy_price must be defined as an input and >0 \ + if using commodity buying" + ) + else: + raise ValueError( + "commodity_buy_price must be defined as an input and >0 \ + if using commodity buying" + ) + + # Check max system capacity + if self.commodity_import_limit == 0 or self.commodity_import_limit is None: + raise ValueError( + "commodity_import_limit must be defined as an input and \ + >0 if using commodity buying" + ) + + def make_dispatch_inputs(self): + dispatch_keys = [ + "cost_per_production", + "cost_per_charge", + "cost_per_discharge", + "demand_met_value", + "max_capacity", + "max_soc_fraction", + "min_soc_fraction", + "charge_efficiency", + "discharge_efficiency", + "max_charge_rate", + "allow_commodity_buying", + ] + + if self.allow_commodity_buying: + dispatch_keys.append("commodity_import_limit") + dispatch_keys.append("commodity_buy_price") + + dispatch_inputs = {k: self.as_dict()[k] for k in dispatch_keys} + dispatch_inputs.update({"initial_soc_fraction": self.init_soc_fraction}) + return dispatch_inputs + + +class OptimizedDispatchController(PyomoControllerBaseClass): + """Operates storage based on optimization to meet the demand profile based on + available commodity from generation profiles and demand profile while minimizing costs. + + Uses a rolling-window optimization approach with configurable horizon and control windows. + + """ + + def setup(self): + """Initialize the optimized dispatch controller.""" + self.config = OptimizedDispatchControllerConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "control") + ) + + self.add_input( + "max_charge_rate", + val=self.config.max_charge_rate, + units=self.config.commodity_rate_units, + desc="Storage charge rate", + ) + + self.add_input( + "storage_capacity", + val=self.config.max_capacity, + units=f"{self.config.commodity_rate_units}*h", + desc="Storage capacity", + ) + + self.n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"] + + self.add_input( + "demand_met_value", + val=self.config.demand_met_value, + shape=self.n_timesteps, + units="USD/" + self.config.commodity_rate_units, + desc="Value of meeting the demand", + ) + + if self.config.allow_commodity_buying: + self.add_input( + f"{self.config.commodity}_buy_price", + val=self.config.commodity_buy_price, + shape=self.n_timesteps, + units="USD/" + self.config.commodity_rate_units, + desc="Value of meeting the demand", + ) + self.add_input( + f"{self.config.commodity}_available_to_buy", + val=self.config.commodity_import_limit, + shape=self.n_timesteps, + units=self.config.commodity_rate_units, + ) + + self.add_output( + f"{self.config.commodity}_bought_for_storage", + val=0, + shape=self.n_timesteps, + units=self.config.commodity_rate_units, + ) + + self.add_output( + "controller_estimated_SOC", + val=0, + shape=self.n_timesteps, + units=self.config.commodity_rate_units, + ) + + self.add_output( + "controller_dispatch_commands", + val=0, + shape=self.n_timesteps, + units=self.config.commodity_rate_units, + ) + + super().setup() + + self.n_control_window = self.config.n_control_window + self.updated_initial_soc = self.config.init_soc_fraction + + self.commodity_info = { + "commodity_name": self.config.commodity, + "commodity_storage_units": self.config.commodity_rate_units, + } + + self.dispatch_inputs = self.config.make_dispatch_inputs() + + def pyomo_setup(self, discrete_inputs, outputs): + """Create the Pyomo model, extract dispatch technology names, and return dispatch solver. + + Returns: + callable: Function(performance_model, performance_model_kwargs, inputs, commodity) + executing rolling-window optimization to determine dispatch and returning: + (total_out, storage_out, unmet_demand, unused_commodity, soc) + """ + # initialize the pyomo model + self.pyomo_model = pyomo.ConcreteModel() + + pyomo.Set(initialize=range(self.config.n_control_window)) + + self.source_techs = [] + self.dispatch_tech = [] + + for connection in self.dispatch_connections: + # get connection definition + source_tech, intended_dispatch_tech = connection + # only add connections to intended dispatch tech + if any(intended_dispatch_tech in name for name in self.tech_group_name): + # record source and dispatch techs + if source_tech == intended_dispatch_tech: + self.dispatch_tech.append(source_tech) + self.source_techs.append(source_tech) + else: + continue + + # define dispatch solver + def pyomo_dispatch_solver( + performance_model: callable, + performance_model_kwargs, + inputs, + pyomo_model=self.pyomo_model, + commodity_name: str = self.config.commodity, + ): + """ + Execute rolling-window dispatch for the controlled technology. + + Iterates over the full simulation period in chunks of size + `self.config.n_control_window`, (re)configures per-window dispatch + parameters, solves the Pyomo optimization model to determine + dispatch decisions, and then calls the provided performance_model + over each window to obtain storage output and SOC trajectories. + + Args: + performance_model (callable): + Function implementing the technology performance over a control + window. Signature must accept (storage_dispatch_commands, + **performance_model_kwargs, sim_start_index=) + and return (storage_out_window, soc_window) arrays of length + n_control_window. + performance_model_kwargs (dict): + Extra keyword arguments forwarded unchanged to performance_model + at window (e.g., efficiencies, timestep size). + inputs (dict): + Dictionary of numpy arrays (length = self.n_timesteps) containing at least: + f"{commodity}_in" : available generated commodity profile. + f"{commodity}_demand" : demanded commodity output profile. + commodity (str, optional): + Base commodity name (e.g. "electricity", "hydrogen"). Default: + self.config.commodity. + + Returns: + tuple[np.ndarray, np.ndarray]: + storage_commodity_out : + Commodity supplied (positive) by the storage asset each timestep. + soc : + State of charge trajectory (percent of capacity). + + Raises: + NotImplementedError: + If the configured control strategy is not implemented. + + Notes: + 1. Arrays returned have length self.n_timesteps (full simulation period). + """ + + # initialize outputs + storage_commodity_out = np.zeros(self.n_timesteps) + soc = np.zeros(self.n_timesteps) + controller_storage_commands = np.zeros(self.n_timesteps) + if self.config.allow_commodity_buying: + commodity_bought = np.zeros(self.n_timesteps) + + # get the starting index for each control window + window_start_indices = list(range(0, self.n_timesteps, self.config.n_control_window)) + + # Initialize parameters for optimized dispatch strategy + time_update_inputs = self.create_time_update_dictionary(inputs, window_start_indices[0]) + self.initialize_parameters(time_update_inputs) + + # loop over all control windows, where t is the starting index of each window + for t in window_start_indices: + # get the inputs over the current control window + time_update_inputs = self.create_time_update_dictionary(inputs, t) + # commodity_in = inputs[f"{self.config.commodity}_in"][ + # t : t + self.config.n_control_window + # ] + # if self.config.allow_commodity_buying: + # inputs[f"{commodity_name}_met_value_in"][t : t + self.config.n_control_window] + # inputs[f"{commodity_name}_buy_price_in"][t : t + self.config.n_control_window] + + # Progress report + if t % (self.n_timesteps // 4) < self.n_control_window: + percentage = round((t / self.n_timesteps) * 100) + print(f"{percentage}% done with optimal dispatch") + # Update time series parameters for the optimization method + self.update_time_series_parameters( + time_update_inputs, + updated_initial_soc=self.updated_initial_soc, + ) + # Run dispatch optimization to minimize costs while meeting demand + self.solve_dispatch_model( + start_time=t, + n_days=self.n_timesteps // 24, + ) + + # run the performance/simulation model for the current control window + # using the dispatch commands + storage_commodity_out_control_window, soc_control_window = performance_model( + self.storage_dispatch_commands, + **performance_model_kwargs, + sim_start_index=t, + ) + # update SOC for next time window + self.updated_initial_soc = soc_control_window[-1] / 100 # turn into ratio + + # get a list of all time indices belonging to the current control window + window_indices = list(range(t, t + self.config.n_control_window)) + + # loop over all time steps in the current control window + for j in window_indices: + # save the output for the control window to the output for the full + # simulation + storage_commodity_out[j] = storage_commodity_out_control_window[j - t] + soc[j] = soc_control_window[j - t] + controller_storage_commands[j] = self.storage_dispatch_commands[j - t] + if self.config.allow_commodity_buying: + commodity_bought[j] = self.hybrid_dispatch_rule.storage_commodity_bought[ + j - t + ] + + if self.config.allow_commodity_buying: + outputs[f"{self.config.commodity}_bought_for_storage"] = commodity_bought + # Note that this SOC is from the performance model and not the + # controller's internal SOC variable + outputs["controller_estimated_SOC"] = soc + outputs["controller_dispatch_commands"] = controller_storage_commands + return storage_commodity_out, soc + + return pyomo_dispatch_solver + + def initialize_parameters(self, inputs): + """Initialize parameters for optimization model + + Args: + inputs (dict): + Dictionary of numpy arrays (length = self.n_timesteps) containing at least: + f"{commodity}_in" : Available generated commodity profile. + f"{commodity}_demand" : Demanded commodity output profile. + + """ + # Where pyomo model communicates with the rest of the controller + # self.hybrid_dispatch_model is the pyomo model, this is the thing in hybrid_rule + if "max_charge_rate" in inputs: + self.dispatch_inputs["max_charge_rate"] = inputs["max_charge_rate"][0] + if "storage_capacity" in inputs: + self.dispatch_inputs["max_capacity"] = inputs["storage_capacity"][0] + self.hybrid_dispatch_model = self._create_dispatch_optimization_model() + self.hybrid_dispatch_rule.create_min_operating_cost_expression() + self.hybrid_dispatch_rule.create_arcs() + assert_units_consistent(self.hybrid_dispatch_model) + + # This calls a class that stores problem state information such as solver metrics and + # the objective function. This is directly used in the H2I simulation, but is + # useful for tracking solver performance and debugging. + self.problem_state = DispatchProblemState() + + # hybrid_dispatch_rule is the thing where you can access variables and hybrid_rule \ + # functions from + self.hybrid_dispatch_rule.initialize_parameters(inputs, self.dispatch_inputs) + + def create_time_update_dictionary(self, inputs, t): + time_update_inputs = {} + input_keys = inputs.keys() + for i in input_keys: + time_update_inputs[i] = inputs[i][t : t + self.config.n_control_window] + + additional_keys = ["demand_met_value"] + if self.config.allow_commodity_buying: + additional_keys.append(f"{self.config.commodity}_buy_price") + for i in additional_keys: + time_update_inputs[i] = self.dispatch_inputs[i][t : t + self.config.n_control_window] + + return time_update_inputs + + def update_time_series_parameters(self, inputs, updated_initial_soc=None): + """Updates the pyomo optimization problem with parameters that change with time + + Args: + commodity_in (list): List of generated commodity in for this time slice. + commodity_demand (list): The demanded commodity for this time slice. + updated_initial_soc (float): The updated initial state of charge for storage + technologies for the current time slice. + """ + self.hybrid_dispatch_rule.update_time_series_parameters(inputs, updated_initial_soc) + + def solve_dispatch_model( + self, + start_time: int = 0, + n_days: int = 0, + ): + """Solves the dispatch optimization model and stores problem metrics. + + Args: + start_time (int): Starting timestep index for the current solve window. + n_days (int): Total number of days in the simulation. + + """ + + solver_results = self.glpk_solve_call(self.hybrid_dispatch_model) + # The outputs of the store_problem_metrics method are not actively used in the H2I + # simulation, but they are useful for debugging and tracking solver performance over time. + self.problem_state.store_problem_metrics( + solver_results, start_time, n_days, pyomo.value(self.hybrid_dispatch_model.objective) + ) + + def _create_dispatch_optimization_model(self): + """ + Creates monolith dispatch model by creating pyomo models for each technology, then + aggregating them into hybrid_rule + """ + model = pyomo.ConcreteModel(name="hybrid_dispatch") + ################################# + # Sets # + ################################# + model.forecast_horizon = pyomo.Set( + doc="Set of time periods in time horizon", + initialize=range(self.config.n_control_window), + ) + for tech in self.source_techs: + if tech == self.dispatch_tech[0]: + dispatch = PyomoRuleStorageMinOperatingCosts( + self.commodity_info, + model, + model.forecast_horizon, + self.config.round_digits, + self.config.time_duration, + self.config.allow_commodity_buying, + block_set_name=f"{tech}_rule", + ) + self.pyomo_model.__setattr__(f"{tech}_rule", dispatch) + else: + dispatch = PyomoDispatchGenericConverterMinOperatingCosts( + self.commodity_info, + model, + model.forecast_horizon, + self.config.round_digits, + self.config.time_duration, + block_set_name=f"{tech}_rule", + ) + self.pyomo_model.__setattr__(f"{tech}_rule", dispatch) + + # Create hybrid pyomo model, inputting individual technology models + self.hybrid_dispatch_rule = PyomoDispatchPlantModel( + model, + model.forecast_horizon, + self.source_techs, + self.pyomo_model, + self.config.time_weighting_factor, + self.config.round_digits, + ) + return model + + def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): + """Build Pyomo model blocks and assign the dispatch solver.""" + self.dispatch_inputs["max_charge_rate"] = inputs["max_charge_rate"][0] + self.dispatch_inputs["max_capacity"] = inputs["storage_capacity"][0] + self.dispatch_inputs["demand_met_value"] = inputs["demand_met_value"][:] + if self.config.allow_commodity_buying: + self.dispatch_inputs[f"{self.config.commodity}_buy_price"] = inputs[ + f"{self.config.commodity}_buy_price" + ][:] + + discrete_outputs["pyomo_dispatch_solver"] = self.pyomo_setup(discrete_inputs, outputs) + + @staticmethod + def glpk_solve_call( + pyomo_model: pyomo.ConcreteModel, + log_name: str = "", + user_solver_options: dict | None = None, + ): + """ + This method takes in the dispatch system-level pyomo model that we have built, + gives it to the solver, and gives back solver results. + """ + + # log_name = "annual_solve_GLPK.log" # For debugging MILP solver + # Ref. on solver options: https://en.wikibooks.org/wiki/GLPK/Using_GLPSOL + glpk_solver_options = { + "cuts": None, + "presol": None, + # 'mostf': None, + # 'mipgap': 0.001, + "tmlim": 30, + } + solver_options = SolverOptions(glpk_solver_options, log_name, user_solver_options, "log") + with pyomo.SolverFactory("glpk") as solver: + results = solver.solve(pyomo_model, options=solver_options.constructed) + + return results + + @property + def storage_dispatch_commands(self) -> list: + """ + Commanded dispatch including available commodity at current time step that has not + been used to charge storage. + """ + return self.hybrid_dispatch_rule.storage_commodity_out diff --git a/h2integrate/control/control_strategies/storage/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/storage/optimized_pyomo_controller.py index 1b89510c5..f228faadf 100644 --- a/h2integrate/control/control_strategies/storage/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/storage/optimized_pyomo_controller.py @@ -53,7 +53,7 @@ class OptimizedDispatchStorageControllerConfig(PyomoStorageControllerBaseConfig) The cost per unit of charging the storage (in $/commodity_rate_units). cost_per_discharge (float): The cost per unit of discharging the storage (in $/commodity_rate_units). - commodity_met_value (float): + demand_met_value (int, float, list): The penalty for not meeting the desired load demand (in $/commodity_rate_units). time_weighting_factor (float): The weighting factor applied to future time steps in the optimization objective @@ -61,32 +61,81 @@ class OptimizedDispatchStorageControllerConfig(PyomoStorageControllerBaseConfig) time_duration (float): The duration of each time step in the Pyomo model in hours. The default of this parameter is 1.0 (i.e., 1 hour time steps). + commodity_import_limit (float): + Maximum amount of the commodity that the storage/system can buy + (i.e. transmission limit for grid charging) + allow_commodity_buying (bool): + This sets whether the storage can buy commodity to charge or not + commodity_buy_price (int, float, list): + Price of the commodity that can be bought for storage. """ max_charge_rate: int | float = field() + allow_commodity_buying: bool = field() charge_efficiency: float = field(validator=range_val(0, 1), default=None) discharge_efficiency: float = field(validator=range_val(0, 1), default=None) + # TODO: note that this definition of cost_per_production is not generalizable to multiple + # production technologies. Would need a name adjustment to connect it to + # production tech cost_per_production: float = field(default=None) cost_per_charge: float = field(default=None) cost_per_discharge: float = field(default=None) - commodity_met_value: float = field(default=None) + demand_met_value: int | float | list = field(default=None) time_weighting_factor: float = field(validator=range_val(0, 1), default=0.995) time_duration: float = field(default=1.0) # hours + # Can we set this to interconnection? do we want to? + commodity_import_limit: float = field(default=None) + commodity_buy_price: int | float | list = field(default=None) + + def __attrs_post_init__(self): + # Check inputs for commodity buying parameters + if self.allow_commodity_buying: + if self.commodity_buy_price: + # Check commodity buy price + if isinstance(self.commodity_buy_price, float | int): + if self.commodity_buy_price == 0: + raise ValueError( + "commodity_buy_price must be defined as an input and >0 \ + if using commodity buying" + ) + if isinstance(self.commodity_buy_price, list) or self.commodity_buy_price is None: + if all(self.commodity_buy_price) == 0: + raise ValueError( + "commodity_buy_price must be defined as an input and >0 \ + if using commodity buying" + ) + else: + raise ValueError( + "commodity_buy_price must be defined as an input and >0 \ + if using commodity buying" + ) + + # Check max system capacity + if self.commodity_import_limit == 0 or self.commodity_import_limit is None: + raise ValueError( + "commodity_import_limit must be defined as an input and \ + >0 if using commodity buying" + ) def make_dispatch_inputs(self): dispatch_keys = [ "cost_per_production", "cost_per_charge", "cost_per_discharge", - "commodity_met_value", + "demand_met_value", "max_capacity", "max_soc_fraction", "min_soc_fraction", "charge_efficiency", "discharge_efficiency", "max_charge_rate", + "allow_commodity_buying", ] + if self.allow_commodity_buying: + dispatch_keys.append("commodity_import_limit") + dispatch_keys.append("commodity_buy_price") + dispatch_inputs = {k: self.as_dict()[k] for k in dispatch_keys} dispatch_inputs.update({"initial_soc_fraction": self.init_soc_fraction}) return dispatch_inputs @@ -127,23 +176,63 @@ def setup(self): self.n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"] + self.add_input( + "demand_met_value", + val=self.config.demand_met_value, + shape=self.n_timesteps, + units="USD/" + self.config.commodity_rate_units, + desc="Value of meeting the demand", + ) + + if self.config.allow_commodity_buying: + self.add_input( + f"{self.config.commodity}_buy_price", + val=self.config.commodity_buy_price, + shape=self.n_timesteps, + units="USD/" + self.config.commodity_rate_units, + desc="Value of meeting the demand", + ) + self.add_input( + f"{self.config.commodity}_available_to_buy", + val=self.config.commodity_import_limit, + shape=self.n_timesteps, + units=self.config.commodity_rate_units, + ) + + self.add_output( + f"{self.config.commodity}_bought_for_storage", + val=0, + shape=self.n_timesteps, + units=self.config.commodity_rate_units, + ) + + self.add_output( + "controller_estimated_SOC", + val=0, + shape=self.n_timesteps, + units=self.config.commodity_rate_units, + ) + + self.add_output( + "controller_dispatch_commands", + val=0, + shape=self.n_timesteps, + units=self.config.commodity_rate_units, + ) + super().setup() self.n_control_window = self.config.n_control_window self.updated_initial_soc = self.config.init_soc_fraction - # Is this the best place to put this??? self.commodity_info = { "commodity_name": self.config.commodity, "commodity_storage_units": self.config.commodity_rate_units, } - # TODO: note that this definition of cost_per_production is not generalizable to multiple - # production technologies. Would need a name adjustment to connect it to - # production tech self.dispatch_inputs = self.config.make_dispatch_inputs() - def pyomo_setup(self, discrete_inputs): + def pyomo_setup(self, discrete_inputs, outputs): """Create the Pyomo model, extract dispatch technology names, and return dispatch solver. Returns: @@ -224,20 +313,27 @@ def pyomo_dispatch_solver( # initialize outputs storage_commodity_out = np.zeros(self.n_timesteps) soc = np.zeros(self.n_timesteps) + controller_storage_commands = np.zeros(self.n_timesteps) + if self.config.allow_commodity_buying: + commodity_bought = np.zeros(self.n_timesteps) # get the starting index for each control window window_start_indices = list(range(0, self.n_timesteps, self.config.n_control_window)) # Initialize parameters for optimized dispatch strategy - self.initialize_parameters(inputs) + time_update_inputs = self.create_time_update_dictionary(inputs, window_start_indices[0]) + self.initialize_parameters(time_update_inputs) # loop over all control windows, where t is the starting index of each window for t in window_start_indices: # get the inputs over the current control window - commodity_in = inputs[f"{self.config.commodity}_in"][ - t : t + self.config.n_control_window - ] - demand_in = inputs[f"{commodity_name}_demand"][t : t + self.config.n_control_window] + time_update_inputs = self.create_time_update_dictionary(inputs, t) + # commodity_in = inputs[f"{self.config.commodity}_in"][ + # t : t + self.config.n_control_window + # ] + # if self.config.allow_commodity_buying: + # inputs[f"{commodity_name}_met_value_in"][t : t + self.config.n_control_window] + # inputs[f"{commodity_name}_buy_price_in"][t : t + self.config.n_control_window] # Progress report if t % (self.n_timesteps // 4) < self.n_control_window: @@ -245,8 +341,7 @@ def pyomo_dispatch_solver( print(f"{percentage}% done with optimal dispatch") # Update time series parameters for the optimization method self.update_time_series_parameters( - commodity_in=commodity_in, - commodity_demand=demand_in, + time_update_inputs, updated_initial_soc=self.updated_initial_soc, ) # Run dispatch optimization to minimize costs while meeting demand @@ -274,7 +369,18 @@ def pyomo_dispatch_solver( # simulation storage_commodity_out[j] = storage_commodity_out_control_window[j - t] soc[j] = soc_control_window[j - t] - + controller_storage_commands[j] = self.storage_dispatch_commands[j - t] + if self.config.allow_commodity_buying: + commodity_bought[j] = self.hybrid_dispatch_rule.storage_commodity_bought[ + j - t + ] + + if self.config.allow_commodity_buying: + outputs[f"{self.config.commodity}_bought_for_storage"] = commodity_bought + # Note that this SOC is from the performance model and not the + # controller's internal SOC variable + outputs["controller_estimated_SOC"] = soc + outputs["controller_dispatch_commands"] = controller_storage_commands return storage_commodity_out, soc return pyomo_dispatch_solver @@ -309,9 +415,21 @@ def initialize_parameters(self, inputs): # functions from self.hybrid_dispatch_rule.initialize_parameters(inputs, self.dispatch_inputs) - def update_time_series_parameters( - self, commodity_in=None, commodity_demand=None, updated_initial_soc=None - ): + def create_time_update_dictionary(self, inputs, t): + time_update_inputs = {} + input_keys = inputs.keys() + for i in input_keys: + time_update_inputs[i] = inputs[i][t : t + self.config.n_control_window] + + additional_keys = ["demand_met_value"] + if self.config.allow_commodity_buying: + additional_keys.append(f"{self.config.commodity}_buy_price") + for i in additional_keys: + time_update_inputs[i] = self.dispatch_inputs[i][t : t + self.config.n_control_window] + + return time_update_inputs + + def update_time_series_parameters(self, inputs, updated_initial_soc=None): """Updates the pyomo optimization problem with parameters that change with time Args: @@ -320,9 +438,7 @@ def update_time_series_parameters( updated_initial_soc (float): The updated initial state of charge for storage technologies for the current time slice. """ - self.hybrid_dispatch_rule.update_time_series_parameters( - commodity_in, commodity_demand, updated_initial_soc - ) + self.hybrid_dispatch_rule.update_time_series_parameters(inputs, updated_initial_soc) def solve_dispatch_model( self, @@ -365,6 +481,7 @@ def _create_dispatch_optimization_model(self): model.forecast_horizon, self.config.round_digits, self.config.time_duration, + self.config.allow_commodity_buying, block_set_name=f"{tech}_rule", ) self.pyomo_model.__setattr__(f"{tech}_rule", dispatch) @@ -394,8 +511,13 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): """Build Pyomo model blocks and assign the dispatch solver.""" self.dispatch_inputs["max_charge_rate"] = inputs["max_charge_rate"][0] self.dispatch_inputs["max_capacity"] = inputs["storage_capacity"][0] + self.dispatch_inputs["demand_met_value"] = inputs["demand_met_value"][:] + if self.config.allow_commodity_buying: + self.dispatch_inputs[f"{self.config.commodity}_buy_price"] = inputs[ + f"{self.config.commodity}_buy_price" + ][:] - discrete_outputs["pyomo_dispatch_solver"] = self.pyomo_setup(discrete_inputs) + discrete_outputs["pyomo_dispatch_solver"] = self.pyomo_setup(discrete_inputs, outputs) @staticmethod def glpk_solve_call( diff --git a/h2integrate/control/control_strategies/storage/test/test_optimal_controllers.py b/h2integrate/control/control_strategies/storage/test/test_optimal_controllers.py index 0c0ec2633..3872d1b1f 100644 --- a/h2integrate/control/control_strategies/storage/test/test_optimal_controllers.py +++ b/h2integrate/control/control_strategies/storage/test/test_optimal_controllers.py @@ -1,3 +1,5 @@ +from copy import deepcopy + import numpy as np import pytest import openmdao.api as om @@ -8,6 +10,7 @@ from h2integrate.storage.simple_storage_auto_sizing import StorageAutoSizingModel from h2integrate.control.control_strategies.storage.optimized_pyomo_controller import ( OptimizedDispatchStorageController, + OptimizedDispatchStorageControllerConfig, ) @@ -60,10 +63,11 @@ def tech_config_battery(): "cost_per_charge": 0.004, "cost_per_discharge": 0.005, "cost_per_production": 0.0, - "commodity_met_value": 0.1, + "demand_met_value": 0.1, "round_digits": 4, "time_weighting_factor": 0.995, "n_control_window": 24, + "allow_commodity_buying": False, }, }, }, @@ -118,11 +122,12 @@ def tech_config_generic(): "tech_name": "h2_storage", "cost_per_charge": 0.03, # USD/kg "cost_per_discharge": 0.05, # USD/kg - "commodity_met_value": 0.1, # USD/kg + "demand_met_value": 0.1, # USD/kg "cost_per_production": 0.0, # USD/kg "time_weighting_factor": 0.995, "system_commodity_interface_limit": 10.0, "n_control_window": 24, + "allow_commodity_buying": False, }, }, } @@ -157,7 +162,7 @@ def tech_config_autosizing(): "tech_name": "h2_storage", "cost_per_charge": 0.03, # USD/kg "cost_per_discharge": 0.05, # USD/kg - "commodity_met_value": 0.1, # USD/kg + "demand_met_value": 0.1, # USD/kg "cost_per_production": 0.0, # USD/kg "time_weighting_factor": 0.995, "system_commodity_interface_limit": 10.0, @@ -165,6 +170,7 @@ def tech_config_autosizing(): "max_charge_rate": 5.0, "max_capacity": 5.0, "init_soc_fraction": 0.1, + "allow_commodity_buying": False, }, }, } @@ -242,7 +248,6 @@ def test_min_operating_cost_load_following_battery_dispatch( assert pytest.approx(prob.model.get_val("battery.SOC")[0], rel=1e-2) == 50 # Find where the signal increases, decreases, and stays at zero - print("SOC", prob.model.get_val("battery.SOC")) indx_soc_increase = np.argwhere( np.diff(prob.model.get_val("battery.SOC", units="unitless"), prepend=True) > 0 ).flatten() @@ -374,7 +379,7 @@ def test_optimal_control_with_generic_storage( assert np.all(prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h") >= 0) with subtests.test("Charge is always negative"): assert np.all(prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h") <= 0) - with subtests.test("Charge + Discharge == storage_hydrogen_out"): + with subtests.test("Charge + Discharge == hydrogen_out"): charge_plus_discharge = prob.get_val( "h2_storage.storage_hydrogen_charge", units="kg/h" ) + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h") @@ -624,3 +629,462 @@ def test_optimal_dispatch_with_autosizing_storage_demand_is_avg_in( expected_charge, rtol=1e-6, ) + + +@pytest.mark.regression +def test_optimal_control_config_with_commodity_buying(subtests): + config_data = { + "tech_name": "h2_storage", + "max_charge_rate": 10.0, + "charge_efficiency": 1.0, + "discharge_efficiency": 1.0, + "commodity": "hydrogen", + "commodity_rate_units": "kg/h", + "max_capacity": 40.0, + "init_soc_fraction": 0.2, + "max_soc_fraction": 1.0, + "min_soc_fraction": 0.1, + "cost_per_charge": 0.03, # USD/kg + "cost_per_discharge": 0.05, # USD/kg + "demand_met_value": 0.1, # USD/kg + "cost_per_production": 0.0, # USD/kg + "time_weighting_factor": 0.995, + "system_commodity_interface_limit": 10.0, + "n_control_window": 24, + "allow_commodity_buying": False, + } + + config = OptimizedDispatchStorageControllerConfig.from_dict(config_data) + + with subtests.test("check commodity_buy_price is None"): + assert config.commodity_buy_price is None + with subtests.test("check commodity_import_limit is None"): + assert config.commodity_import_limit is None + + config_data["allow_commodity_buying"] = True + + with subtests.test("with invalid commodity_buy_price"): + with pytest.raises(ValueError): + data = deepcopy(config_data) + data["allow_commodity_buying"] = True + OptimizedDispatchStorageControllerConfig.from_dict(data) + + with pytest.raises(ValueError): + data = deepcopy(config_data) + data["allow_commodity_buying"] = True + data["commodity_buy_price"] = 0.4 + OptimizedDispatchStorageControllerConfig.from_dict(data) + + with pytest.raises(ValueError): + data = deepcopy(config_data) + data["allow_commodity_buying"] = True + data["commodity_buy_price"] = 0.4 + data["commodity_import_limit"] = 0.0 + OptimizedDispatchStorageControllerConfig.from_dict(data) + + +@pytest.mark.regression +def test_optimal_control_with_commodity_buying_generic_storage( + plant_config_h2_storage, tech_config_generic, subtests +): + commodity_demand = np.full(48, 5.0) + commodity_in = np.tile(np.concat([np.zeros(3), np.cumsum(np.ones(15)), np.full(6, 4.0)]), 2) + commodity_buy_price = np.tile(np.concat([np.arange(-3, 9), np.arange(8, -4, -1)]), 2) + commodity_import_limit = 7 + + # Set grid charging parameters + tech_config_generic["technologies"]["h2_storage"]["model_inputs"]["control_parameters"] = { + "tech_name": "h2_storage", + "cost_per_charge": 0.03, # USD/kg + "cost_per_discharge": 0.05, # USD/kg + "demand_met_value": 0.1, # USD/kg + "cost_per_production": 0.0, # USD/kg + "time_weighting_factor": 0.995, + "system_commodity_interface_limit": 10.0, + "n_control_window": 24, + "allow_commodity_buying": True, + "commodity_buy_price": 1, + "commodity_import_limit": commodity_import_limit, + } + + # Setup the OpenMDAO problem and add subsystems + prob = om.Problem() + + prob.model.add_subsystem( + "h2_storage_optimized_load_following_controller", + OptimizedDispatchStorageController( + plant_config=plant_config_h2_storage, + tech_config=tech_config_generic["technologies"]["h2_storage"], + ), + promotes=["*"], + ) + + prob.model.add_subsystem( + "h2_storage", + StoragePerformanceModel( + plant_config=plant_config_h2_storage, + tech_config=tech_config_generic["technologies"]["h2_storage"], + ), + promotes=["*"], + ) + + # Setup the system and required values + prob.setup() + prob.set_val("h2_storage.hydrogen_in", commodity_in) + prob.set_val("h2_storage.hydrogen_demand", commodity_demand) + prob.set_val("hydrogen_buy_price", commodity_buy_price) + + # Run the model + prob.run_model() + + charge_rate = prob.get_val("h2_storage.max_charge_rate", units="kg/h")[0] + discharge_rate = prob.get_val("h2_storage.max_charge_rate", units="kg/h")[0] + capacity = prob.get_val("h2_storage.storage_capacity", units="kg")[0] + + print("outputs: ", prob.get_val("hydrogen_out")) + print("discharge: ", prob.get_val("h2_storage.storage_hydrogen_discharge")) + print("charge: ", prob.get_val("h2_storage.storage_hydrogen_charge")) + print("commodity in: ", prob.get_val("h2_storage.hydrogen_in")) + print("demand: ", prob.get_val("h2_storage.hydrogen_demand")) + print("commodity_buy: ", prob.get_val("hydrogen_buy_price")) + print("hydrogen_out: ", prob.get_val("hydrogen_out")) + print("hydrogen bought: ", prob.get_val("hydrogen_bought_for_storage")) + + # Test that discharge is always positive + with subtests.test("Discharge is always positive"): + assert np.all(prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h") >= 0) + with subtests.test("Charge is always negative"): + assert np.all(prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h") <= 0) + with subtests.test("Charge + Discharge == hydrogen_out"): + charge_plus_discharge = prob.get_val( + "h2_storage.storage_hydrogen_charge", units="kg/h" + ) + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h") + np.testing.assert_allclose( + charge_plus_discharge, prob.get_val("hydrogen_out", units="kg/h"), rtol=1e-6 + ) + with subtests.test("Initial SOC is correct"): + assert ( + pytest.approx(prob.model.get_val("h2_storage.SOC", units="unitless")[0], rel=1e-6) + == 0.375 + ) + + indx_soc_increase = np.argwhere( + np.diff(prob.model.get_val("h2_storage.SOC", units="unitless"), prepend=True) > 0 + ).flatten() + indx_soc_decrease = np.argwhere( + np.diff(prob.model.get_val("h2_storage.SOC", units="unitless"), prepend=False) < 0 + ).flatten() + indx_soc_same = np.argwhere( + np.diff(prob.model.get_val("h2_storage.SOC", units="unitless"), prepend=True) == 0.0 + ).flatten() + + with subtests.test("SOC increases when charging"): + assert np.all( + prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h")[indx_soc_increase] < 0 + ) + assert np.all( + prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h")[indx_soc_decrease] == 0 + ) + assert np.all( + prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h")[indx_soc_same] == 0 + ) + + with subtests.test("SOC decreases when discharging"): + assert np.all( + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h")[indx_soc_decrease] + > 0 + ) + assert np.all( + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h")[indx_soc_increase] + == 0 + ) + assert np.all( + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h")[indx_soc_same] == 0 + ) + + with subtests.test("Max SOC <= Max storage percent"): + assert prob.get_val("h2_storage.SOC", units="unitless").max() <= 1.0 + + with subtests.test("Min SOC >= Min storage percent"): + assert prob.get_val("h2_storage.SOC", units="unitless").min() >= 0.1 + + with subtests.test("Charge never exceeds charge rate"): + assert ( + prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h").min() + >= -1 * charge_rate + ) + + with subtests.test("Discharge never exceeds discharge rate"): + assert ( + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h").max() + <= discharge_rate + ) + + with subtests.test("Discharge never exceeds demand"): + np.testing.assert_allclose( + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h").max(), + commodity_demand, + rtol=1e-6, + ) + + with subtests.test("Sometimes discharges"): + assert any( + k > 1e-3 for k in prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h") + ) + + with subtests.test("Sometimes charges"): + assert any( + k < -1e-3 for k in prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h") + ) + + with subtests.test("Cumulative charge/discharge does not exceed storage capacity"): + assert np.cumsum(prob.get_val("hydrogen_out", units="kg/h")).max() <= capacity + assert np.cumsum(prob.get_val("hydrogen_out", units="kg/h")).min() >= -1 * capacity + + with subtests.test("Expected discharge from hour 10-30"): + expected_discharge = np.concat( + [np.zeros(8), np.ones(3), np.zeros(3), [5, 0], np.arange(5, 1, -1)] + ) + np.testing.assert_allclose( + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h")[10:30], + expected_discharge, + rtol=1e-6, + atol=1e-6, + ) + + with subtests.test("Expected charge hour 0-20"): + expected_charge = np.concat([np.ones(3) * -7, np.zeros(5), [-1, -2], np.zeros(10)]) + np.testing.assert_allclose( + prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h")[0:20], + expected_charge, + rtol=1e-6, + ) + + with subtests.test("Output never exceeds system commodity draw limit"): + np.testing.assert_allclose( + prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h").min(), + -commodity_import_limit, + rtol=1e-6, + ) + + with subtests.test("Expected bought commodity"): + expected_commodity_bought = np.concat( + [np.ones(3) * 7, np.zeros(18), np.ones(3) * 7, [0], [5], np.zeros(19), np.ones(3) * 7] + ) + np.testing.assert_allclose( + prob.get_val("hydrogen_bought_for_storage", units="kg/h"), + expected_commodity_bought, + rtol=1e-6, + atol=1e-6, + ) + + +def _setup_commodity_buying_problem( + plant_config, tech_config, commodity_buy_price_timeseries, commodity_import_limit=7 +): + """Helper to set up an OptimizedDispatchStorageController problem with commodity buying enabled. + + Args: + plant_config: Plant configuration dictionary. + tech_config: Technology configuration dictionary (will be modified in-place). + commodity_buy_price_timeseries: Numpy array of buy prices per timestep. + commodity_import_limit: Maximum commodity import rate (default: 7 kg/h). + + Returns: + tuple: (prob, commodity_in, commodity_demand) - the configured OpenMDAO problem + and the input arrays used. + """ + commodity_demand = np.full(48, 5.0) + commodity_in = np.tile(np.concat([np.zeros(3), np.cumsum(np.ones(15)), np.full(6, 4.0)]), 2) + + tech_config["technologies"]["h2_storage"]["model_inputs"]["control_parameters"] = { + "tech_name": "h2_storage", + "cost_per_charge": 0.03, + "cost_per_discharge": 0.05, + "demand_met_value": 0.1, + "cost_per_production": 0.0, + "time_weighting_factor": 0.995, + "system_commodity_interface_limit": 10.0, + "n_control_window": 24, + "allow_commodity_buying": True, + "commodity_buy_price": 1, + "commodity_import_limit": commodity_import_limit, + } + + prob = om.Problem() + + prob.model.add_subsystem( + "h2_storage_optimized_load_following_controller", + OptimizedDispatchStorageController( + plant_config=plant_config, + tech_config=tech_config["technologies"]["h2_storage"], + ), + promotes=["*"], + ) + + prob.model.add_subsystem( + "h2_storage", + StoragePerformanceModel( + plant_config=plant_config, + tech_config=tech_config["technologies"]["h2_storage"], + ), + promotes=["*"], + ) + + prob.setup() + prob.set_val("h2_storage.hydrogen_in", commodity_in) + prob.set_val("h2_storage.hydrogen_demand", commodity_demand) + prob.set_val("hydrogen_buy_price", commodity_buy_price_timeseries) + + prob.run_model() + + return prob, commodity_in, commodity_demand + + +def _run_standard_commodity_buying_assertions(prob, commodity_demand, subtests): + """Run standard physical constraint checks for commodity buying tests. + + Args: + prob: The solved OpenMDAO problem. + commodity_demand: The demand profile array. + subtests: pytest subtests fixture. + """ + charge_rate = prob.get_val("h2_storage.max_charge_rate", units="kg/h")[0] + discharge_rate = prob.get_val("h2_storage.max_charge_rate", units="kg/h")[0] + capacity = prob.get_val("h2_storage.storage_capacity", units="kg")[0] + + with subtests.test("Discharge is always positive"): + assert np.all(prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h") >= 0) + with subtests.test("Charge is always negative"): + assert np.all(prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h") <= 0) + with subtests.test("Charge + Discharge == hydrogen_out"): + charge_plus_discharge = prob.get_val( + "h2_storage.storage_hydrogen_charge", units="kg/h" + ) + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h") + np.testing.assert_allclose( + charge_plus_discharge, prob.get_val("hydrogen_out", units="kg/h"), rtol=1e-6 + ) + with subtests.test("Max SOC <= Max storage percent"): + assert prob.get_val("h2_storage.SOC", units="unitless").max() <= 1.0 + with subtests.test("Min SOC >= Min storage percent"): + assert prob.get_val("h2_storage.SOC", units="unitless").min() >= 0.1 + with subtests.test("Charge never exceeds charge rate"): + assert ( + prob.get_val("h2_storage.storage_hydrogen_charge", units="kg/h").min() + >= -1 * charge_rate + ) + with subtests.test("Discharge never exceeds discharge rate"): + assert ( + prob.get_val("h2_storage.storage_hydrogen_discharge", units="kg/h").max() + <= discharge_rate + ) + with subtests.test("Commodity bought is non-negative"): + assert np.all(prob.get_val("hydrogen_bought_for_storage", units="kg/h") >= -1e-6) + with subtests.test("Cumulative charge/discharge does not exceed storage capacity"): + assert np.cumsum(prob.get_val("hydrogen_out", units="kg/h")).max() <= capacity + assert np.cumsum(prob.get_val("hydrogen_out", units="kg/h")).min() >= -1 * capacity + + +@pytest.mark.regression +def test_optimal_control_commodity_buying_all_positive_prices( + plant_config_h2_storage, tech_config_generic, subtests +): + """Test commodity buying with all positive prices. + + When all prices are positive, buying commodity adds to the operating cost. + The optimizer should only buy when the cost of buying is less than the + penalty for not meeting demand (demand_met_value=0.1 $/kg). + With a high buy price (5.0 $/kg >> demand_met_value), buying is expensive + and the optimizer should prefer to not buy. + """ + commodity_buy_price = np.full(48, 5.0) + commodity_import_limit = 7 + + prob, commodity_in, commodity_demand = _setup_commodity_buying_problem( + plant_config_h2_storage, + tech_config_generic, + commodity_buy_price, + commodity_import_limit, + ) + + _run_standard_commodity_buying_assertions(prob, commodity_demand, subtests) + + commodity_bought = prob.get_val("hydrogen_bought_for_storage", units="kg/h") + + # With high positive prices (5.0 >> demand_met_value of 0.1), buying is too expensive + # relative to the penalty for unmet demand, so the optimizer should not buy any commodity. + with subtests.test("No commodity bought when price exceeds demand met value"): + np.testing.assert_allclose(commodity_bought, 0.0, atol=1e-6) + + +@pytest.mark.regression +def test_optimal_control_commodity_buying_all_negative_prices( + plant_config_h2_storage, tech_config_generic, subtests +): + """Test commodity buying with all negative prices. + + When all prices are negative, buying commodity *reduces* the operating cost + (the optimizer gets paid to buy). The optimizer should buy as much as + possible up to the import limit whenever it can. + """ + commodity_buy_price = np.full(48, -5.0) + commodity_import_limit = 7 + + prob, commodity_in, commodity_demand = _setup_commodity_buying_problem( + plant_config_h2_storage, + tech_config_generic, + commodity_buy_price, + commodity_import_limit, + ) + + _run_standard_commodity_buying_assertions(prob, commodity_demand, subtests) + + commodity_bought = prob.get_val("hydrogen_bought_for_storage", units="kg/h") + + # With negative prices, buying reduces cost. The optimizer should buy aggressively. + with subtests.test("Commodity is bought when prices are negative"): + assert np.sum(commodity_bought) > 0 + + # The purchase_limit constraint prevents buying when is_generating=1, but + # whenever the system is not generating, it should buy the maximum. + with subtests.test("Bought commodity never exceeds import limit"): + assert np.all(commodity_bought <= commodity_import_limit + 1e-6) + + # With strong negative prices, the optimizer should buy at the import limit + # whenever it is allowed (i.e. when not generating). Check that it buys at + # the maximum for at least some timesteps. + with subtests.test("Buys at import limit for some timesteps"): + assert np.any(np.isclose(commodity_bought, commodity_import_limit, atol=1e-4)) + + +@pytest.mark.regression +def test_optimal_control_commodity_buying_zero_prices( + plant_config_h2_storage, tech_config_generic, subtests +): + """Test commodity buying with zero prices. + + When the buy price is zero, buying commodity has no direct cost. The + optimizer should buy freely to help charge storage and meet demand since + there is no cost penalty for doing so. + """ + commodity_buy_price = np.full(48, 0.0) + commodity_import_limit = 7 + + prob, commodity_in, commodity_demand = _setup_commodity_buying_problem( + plant_config_h2_storage, + tech_config_generic, + commodity_buy_price, + commodity_import_limit, + ) + + _run_standard_commodity_buying_assertions(prob, commodity_demand, subtests) + + commodity_bought = prob.get_val("hydrogen_bought_for_storage", units="kg/h") + + # With zero prices, buying is free. The optimizer should buy commodity to help + # meet demand, especially during periods when commodity_in is insufficient. + with subtests.test("Commodity is bought when prices are zero"): + assert np.sum(commodity_bought) > 0 + + with subtests.test("Bought commodity never exceeds import limit"): + assert np.all(commodity_bought <= commodity_import_limit + 1e-6)