From 2c2aaaad2662edb10156ccbf50353733408eb1f8 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Fri, 6 Mar 2026 16:47:33 -0700 Subject: [PATCH 01/21] Breakout simulation code --- .../heuristic_pyomo_controller.py | 166 ++++++++++++++++++ .../optimized_pyomo_controller.py | 157 +++++++++++++++++ 2 files changed, 323 insertions(+) diff --git a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py index 85729da98..8039788ee 100644 --- a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py +++ b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py @@ -1,5 +1,7 @@ from typing import TYPE_CHECKING +import numpy as np +import pyomo.environ as pyomo from attrs import field, define from h2integrate.core.utilities import merge_shared_inputs @@ -41,6 +43,21 @@ def setup(self): super().setup() + for connection in self.dispatch_connections: + # get connection definition + source_tech, intended_dispatch_tech = connection + if any(intended_dispatch_tech in name for name in self.tech_group_name): + if source_tech == intended_dispatch_tech: + # When getting rules for the same tech, the tech name is not used in order to + # allow for automatic connections rather than complicating the h2i model set up + self.add_discrete_input("dispatch_block_rule_function", val=self.dummy_method) + else: + self.add_discrete_input( + f"{'dispatch_block_rule_function'}_{source_tech}", val=self.dummy_method + ) + else: + continue + self.add_input( "storage_capacity", val=self.config.max_capacity, @@ -57,6 +74,155 @@ def setup(self): if self.config.discharge_efficiency is not None: self.discharge_efficiency = self.config.discharge_efficiency + def pyomo_setup(self, discrete_inputs): + """Create the Pyomo model, attach per-tech Blocks, and return dispatch solver. + + Returns: + callable: Function(performance_model, performance_model_kwargs, inputs, commodity) + executing rolling-window heuristic dispatch or optimization and returning: + (total_out, storage_out, unmet_demand, unused_commodity, soc) + """ + # initialize the pyomo model + self.pyomo_model = pyomo.ConcreteModel() + + index_set = pyomo.Set(initialize=range(self.config.n_control_window)) + + # run each pyomo rule set up function for each technology + 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): + # names are specified differently if connecting within the tech group vs + # connecting from an external tech group. This facilitates OM connections + if source_tech == intended_dispatch_tech: + dispatch_block_rule_function = discrete_inputs["dispatch_block_rule_function"] + self.dispatch_tech.append(source_tech) + else: + dispatch_block_rule_function = discrete_inputs[ + f"{'dispatch_block_rule_function'}_{source_tech}" + ] + # create pyomo block and set attr + blocks = pyomo.Block(index_set, rule=dispatch_block_rule_function) + setattr(self.pyomo_model, source_tech, blocks) + 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, invokes a heuristic control strategy to set fixed + 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, np.ndarray, np.ndarray, np.ndarray]: + total_commodity_out : + Net commodity supplied to demand each timestep (min(demand, storage + gen)). + storage_commodity_out : + Commodity supplied (positive) by the storage asset each timestep. + unmet_demand : + Positive shortfall = demand - total_out (0 if fully met). + unused_commodity : + Surplus generation + storage discharge not used to meet demand. + 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 + unmet_demand = np.zeros(self.n_timesteps) + storage_commodity_out = np.zeros(self.n_timesteps) + total_commodity_out = np.zeros(self.n_timesteps) + unused_commodity = np.zeros(self.n_timesteps) + soc = 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)) + + self.initialize_parameters() + + # 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[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] + + # Update time series parameters for the heuristic method + self.update_time_series_parameters() + # determine dispatch commands for the current control window + # using the heuristic method + self.set_fixed_dispatch( + commodity_in, + self.config.system_commodity_interface_limit, + demand_in, + ) + + # 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, + ) + + # 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] + total_commodity_out[j] = np.minimum( + demand_in[j - t], storage_commodity_out[j] + commodity_in[j - t] + ) + unmet_demand[j] = np.maximum(0, demand_in[j - t] - total_commodity_out[j]) + unused_commodity[j] = np.maximum( + 0, storage_commodity_out[j] + commodity_in[j - t] - demand_in[j - t] + ) + + return total_commodity_out, storage_commodity_out, unmet_demand, unused_commodity, soc + + return pyomo_dispatch_solver + def initialize_parameters(self, inputs): """Initializes parameters.""" diff --git a/h2integrate/control/control_strategies/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index b9f0da125..d8a5eeb9b 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -1,5 +1,6 @@ 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 @@ -137,6 +138,162 @@ def setup(self): self.dispatch_inputs = self.config.make_dispatch_inputs() + def pyomo_setup(self, discrete_inputs): + """Create the Pyomo model, attach per-tech Blocks, and return dispatch solver. + + Returns: + callable: Function(performance_model, performance_model_kwargs, inputs, commodity) + executing rolling-window heuristic dispatch or optimization 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 = [] + + # run each pyomo rule set up function for each technology + 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) + # create pyomo block and set attr + 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, invokes a heuristic control strategy to set fixed + 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, np.ndarray, np.ndarray, np.ndarray]: + total_commodity_out : + Net commodity supplied to demand each timestep (min(demand, storage + gen)). + storage_commodity_out : + Commodity supplied (positive) by the storage asset each timestep. + unmet_demand : + Positive shortfall = demand - total_out (0 if fully met). + unused_commodity : + Surplus generation + storage discharge not used to meet demand. + 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 + unmet_demand = np.zeros(self.n_timesteps) + storage_commodity_out = np.zeros(self.n_timesteps) + total_commodity_out = np.zeros(self.n_timesteps) + unused_commodity = np.zeros(self.n_timesteps) + soc = 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[f"{commodity_name}_in"], inputs[f"{commodity_name}_demand"] + ) + + # 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] + + # 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( + commodity_in=commodity_in, + commodity_demand=demand_in, + 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] + total_commodity_out[j] = np.minimum( + demand_in[j - t], storage_commodity_out[j] + commodity_in[j - t] + ) + unmet_demand[j] = np.maximum(0, demand_in[j - t] - total_commodity_out[j]) + unused_commodity[j] = np.maximum( + 0, storage_commodity_out[j] + commodity_in[j - t] - demand_in[j - t] + ) + + return total_commodity_out, storage_commodity_out, unmet_demand, unused_commodity, soc + + return pyomo_dispatch_solver + def initialize_parameters(self, commodity_in, commodity_demand): """Initialize parameters for optimization model From d6ee5f03f0c46eb8d24192a14c6a9f76b1805d65 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Sat, 7 Mar 2026 08:59:06 -0700 Subject: [PATCH 02/21] Remove duplicate code --- .../pyomo_controller_baseclass.py | 80 +++++++++---------- 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py index b2693d207..bf6109bf9 100644 --- a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py +++ b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py @@ -123,20 +123,20 @@ def setup(self): # create inputs for all pyomo object creation functions from all connected technologies self.dispatch_connections = self.options["plant_config"]["tech_to_dispatch_connections"] - for connection in self.dispatch_connections: - # get connection definition - source_tech, intended_dispatch_tech = connection - if any(intended_dispatch_tech in name for name in self.tech_group_name): - if source_tech == intended_dispatch_tech: - # When getting rules for the same tech, the tech name is not used in order to - # allow for automatic connections rather than complicating the h2i model set up - self.add_discrete_input("dispatch_block_rule_function", val=self.dummy_method) - else: - self.add_discrete_input( - f"{'dispatch_block_rule_function'}_{source_tech}", val=self.dummy_method - ) - else: - continue + # for connection in self.dispatch_connections: + # # get connection definition + # source_tech, intended_dispatch_tech = connection + # if any(intended_dispatch_tech in name for name in self.tech_group_name): + # if source_tech == intended_dispatch_tech: + # # When getting rules for the same tech, the tech name is not used in order to + # # allow for automatic connections rather than complicating the h2i model setup + # self.add_discrete_input("dispatch_block_rule_function", val=self.dummy_method) + # else: + # self.add_discrete_input( + # f"{'dispatch_block_rule_function'}_{source_tech}", val=self.dummy_method + # ) + # else: + # continue # create output for the pyomo control model self.add_discrete_output( @@ -161,32 +161,32 @@ def pyomo_setup(self, discrete_inputs): # initialize the pyomo model self.pyomo_model = pyomo.ConcreteModel() - index_set = pyomo.Set(initialize=range(self.config.n_control_window)) - - self.source_techs = [] - self.dispatch_tech = [] - - # run each pyomo rule set up function for each technology - 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): - # names are specified differently if connecting within the tech group vs - # connecting from an external tech group. This facilitates OM connections - if source_tech == intended_dispatch_tech: - dispatch_block_rule_function = discrete_inputs["dispatch_block_rule_function"] - self.dispatch_tech.append(source_tech) - else: - dispatch_block_rule_function = discrete_inputs[ - f"{'dispatch_block_rule_function'}_{source_tech}" - ] - # create pyomo block and set attr - blocks = pyomo.Block(index_set, rule=dispatch_block_rule_function) - setattr(self.pyomo_model, source_tech, blocks) - self.source_techs.append(source_tech) - else: - continue + # index_set = pyomo.Set(initialize=range(self.config.n_control_window)) + + # self.source_techs = [] + # self.dispatch_tech = [] + + # # run each pyomo rule set up function for each technology + # 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): + # # names are specified differently if connecting within the tech group vs + # # connecting from an external tech group. This facilitates OM connections + # if source_tech == intended_dispatch_tech: + # dispatch_block_rule_function = discrete_inputs["dispatch_block_rule_function"] + # self.dispatch_tech.append(source_tech) + # else: + # dispatch_block_rule_function = discrete_inputs[ + # f"{'dispatch_block_rule_function'}_{source_tech}" + # ] + # # create pyomo block and set attr + # blocks = pyomo.Block(index_set, rule=dispatch_block_rule_function) + # setattr(self.pyomo_model, source_tech, blocks) + # self.source_techs.append(source_tech) + # else: + # continue # define dispatch solver def pyomo_dispatch_solver( From 14be9ddd0660aeb1e4cb16139a6dd5ffe845b2ec Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 9 Mar 2026 11:58:23 -0400 Subject: [PATCH 03/21] Fix testing --- .../heuristic_pyomo_controller.py | 3 +- .../pyomo_controller_baseclass.py | 142 +----------------- .../control/test/test_optimal_controllers.py | 6 - h2integrate/core/h2integrate_model.py | 13 +- 4 files changed, 18 insertions(+), 146 deletions(-) diff --git a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py index 8039788ee..278228545 100644 --- a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py +++ b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py @@ -97,7 +97,6 @@ def pyomo_setup(self, discrete_inputs): # connecting from an external tech group. This facilitates OM connections if source_tech == intended_dispatch_tech: dispatch_block_rule_function = discrete_inputs["dispatch_block_rule_function"] - self.dispatch_tech.append(source_tech) else: dispatch_block_rule_function = discrete_inputs[ f"{'dispatch_block_rule_function'}_{source_tech}" @@ -174,7 +173,7 @@ def pyomo_dispatch_solver( # get the starting index for each control window window_start_indices = list(range(0, self.n_timesteps, self.config.n_control_window)) - self.initialize_parameters() + self.initialize_parameters(inputs) # loop over all control windows, where t is the starting index of each window for t in window_start_indices: diff --git a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py index bf6109bf9..4c66173a7 100644 --- a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py +++ b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py @@ -123,20 +123,6 @@ def setup(self): # create inputs for all pyomo object creation functions from all connected technologies self.dispatch_connections = self.options["plant_config"]["tech_to_dispatch_connections"] - # for connection in self.dispatch_connections: - # # get connection definition - # source_tech, intended_dispatch_tech = connection - # if any(intended_dispatch_tech in name for name in self.tech_group_name): - # if source_tech == intended_dispatch_tech: - # # When getting rules for the same tech, the tech name is not used in order to - # # allow for automatic connections rather than complicating the h2i model setup - # self.add_discrete_input("dispatch_block_rule_function", val=self.dummy_method) - # else: - # self.add_discrete_input( - # f"{'dispatch_block_rule_function'}_{source_tech}", val=self.dummy_method - # ) - # else: - # continue # create output for the pyomo control model self.add_discrete_output( @@ -161,33 +147,6 @@ def pyomo_setup(self, discrete_inputs): # initialize the pyomo model self.pyomo_model = pyomo.ConcreteModel() - # index_set = pyomo.Set(initialize=range(self.config.n_control_window)) - - # self.source_techs = [] - # self.dispatch_tech = [] - - # # run each pyomo rule set up function for each technology - # 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): - # # names are specified differently if connecting within the tech group vs - # # connecting from an external tech group. This facilitates OM connections - # if source_tech == intended_dispatch_tech: - # dispatch_block_rule_function = discrete_inputs["dispatch_block_rule_function"] - # self.dispatch_tech.append(source_tech) - # else: - # dispatch_block_rule_function = discrete_inputs[ - # f"{'dispatch_block_rule_function'}_{source_tech}" - # ] - # # create pyomo block and set attr - # blocks = pyomo.Block(index_set, rule=dispatch_block_rule_function) - # setattr(self.pyomo_model, source_tech, blocks) - # self.source_techs.append(source_tech) - # else: - # continue - # define dispatch solver def pyomo_dispatch_solver( performance_model: callable, @@ -251,99 +210,14 @@ def pyomo_dispatch_solver( unused_commodity = np.zeros(self.n_timesteps) soc = 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)) - - control_strategy = self.options["tech_config"]["control_strategy"]["model"] - - # TODO: implement optional kwargs for this method: maybe this will remove if statement here - if "Heuristic" in control_strategy: - # Initialize parameters for heuristic dispatch strategy - self.initialize_parameters(inputs) - elif "Optimized" in control_strategy: - # Initialize parameters for optimized dispatch strategy - self.initialize_parameters( - inputs[f"{commodity_name}_in"], inputs[f"{commodity_name}_demand"] - ) - - else: - raise ( - NotImplementedError( - f"Control strategy '{control_strategy}' was given, \ - but has not been implemented yet." - ) - ) - - # 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] - - if "Heuristic" in control_strategy: - # Update time series parameters for the heuristic method - self.update_time_series_parameters() - # determine dispatch commands for the current control window - # using the heuristic method - self.set_fixed_dispatch( - commodity_in, - self.config.system_commodity_interface_limit, - demand_in, - ) - - elif "Optimized" in control_strategy: - # 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( - commodity_in=commodity_in, - commodity_demand=demand_in, - 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, - ) - - else: - raise ( - NotImplementedError( - f"Control strategy '{control_strategy}' was given, \ - but has not been implemented yet." - ) - ) - - # 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] - total_commodity_out[j] = np.minimum( - demand_in[j - t], storage_commodity_out[j] + commodity_in[j - t] - ) - unmet_demand[j] = np.maximum(0, demand_in[j - t] - total_commodity_out[j]) - unused_commodity[j] = np.maximum( - 0, storage_commodity_out[j] + commodity_in[j - t] - demand_in[j - t] - ) + # This is where the pyomo controller interface is defined in individual + # pyomo controllers + + # The structure should be as follows: + # 1. initialize_parameters() + # 2. update_time_series_parameters() + # 3. solve dispatch model or set fixed dispatch + # 4. update outputs return total_commodity_out, storage_commodity_out, unmet_demand, unused_commodity, soc diff --git a/h2integrate/control/test/test_optimal_controllers.py b/h2integrate/control/test/test_optimal_controllers.py index 3dd9b1018..f6b78b850 100644 --- a/h2integrate/control/test/test_optimal_controllers.py +++ b/h2integrate/control/test/test_optimal_controllers.py @@ -38,7 +38,6 @@ "description": "...", "technologies": { "battery": { - "dispatch_rule_set": {"model": "PyomoRuleStorageBaseclass"}, "control_strategy": {"model": "OptimizedDispatchController"}, "performance_model": {"model": "PySAMBatteryPerformanceModel"}, "model_inputs": { @@ -73,17 +72,12 @@ }, "combiner": { "performance_model": {"model": "GenericCombinerPerformanceModel"}, - "dispatch_rule_set": {"model": "PyomoDispatchGenericConverter"}, "model_inputs": { "performance_parameters": { "commodity": "electricity", "commodity_rate_units": "kW", "in_streams": 1, }, - "dispatch_rule_parameters": { - "commodity": "electricity", - "commodity_rate_units": "kW", - }, }, }, }, diff --git a/h2integrate/core/h2integrate_model.py b/h2integrate/core/h2integrate_model.py index b0c9c2ece..1fbe63e43 100644 --- a/h2integrate/core/h2integrate_model.py +++ b/h2integrate/core/h2integrate_model.py @@ -1204,11 +1204,16 @@ def connect_technologies(self): if tech_name == dispatching_tech_name: continue else: - # Connect the dispatch rules output to the dispatching_tech_name input - self.model.connect( - f"{tech_name}.dispatch_block_rule_function", - f"{dispatching_tech_name}.dispatch_block_rule_function_{tech_name}", + # Only connect dispatch rules if they are defined in the tech_config + tech_dispatch_rule = self.technology_config.get(tech_name, {}).get( + "dispatch_rule_set", False ) + if tech_dispatch_rule: + # Connect the dispatch rules output to the dispatching_tech_name input + self.model.connect( + f"{tech_name}.dispatch_block_rule_function", + f"{dispatching_tech_name}.dispatch_block_rule_function_{tech_name}", + ) if (pyxdsm is not None) and (len(technology_interconnections) > 0): try: From 8db833251adb37de875852d5dc6a986541fa646a Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 9 Mar 2026 14:07:20 -0400 Subject: [PATCH 04/21] Update doc strings and docs --- CHANGELOG.md | 1 + docs/control/pyomo_controllers.md | 2 +- .../30_pyomo_optimized_dispatch/tech_config.yaml | 4 ---- .../heuristic_pyomo_controller.py | 2 +- .../optimized_pyomo_controller.py | 10 ++++------ .../pyomo_controller_baseclass.py | 14 ++++++-------- 6 files changed, 13 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bdbc7f06e..757030172 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ `min_charge_fraction`, `max_charge_fraction`, and `init_charge_fraction` across all configuration classes, YAML configs, tests, and examples. These values are fractions between 0 and 1, so the previous "percent" naming was misleading. [PR 581](https://github.com/NatLabRockies/H2Integrate/pull/581) +- Breaks out pyomo controller simulation code from base class to individual controllers. [PR 587](https://github.com/NatLabRockies/H2Integrate/pull/587) ## 0.7 [March 3, 2026] diff --git a/docs/control/pyomo_controllers.md b/docs/control/pyomo_controllers.md index 1f6c517fc..8f7a35708 100644 --- a/docs/control/pyomo_controllers.md +++ b/docs/control/pyomo_controllers.md @@ -21,7 +21,7 @@ For an example of how to use the heuristic Pyomo control framework with the `Heu (optimized-load-following-controller)= ## Optimized Load Following Controller -The optimized dispatch method is specified by setting the storage control to `OptimizedDispatchController`. The same `dispatch_rule_set` for each technology connected to the storage technology is followed as in the heuristic case, where each technology must also have a `dispatch_rule_set` defined in the `tech_config`. The optimized dispatch currently uses the same dispatch rules as the heuristic dispatch. This method maximizes the load met while minimizing the cost of the system (operating cost) over each specified time window. +The optimized dispatch method is specified by setting the storage control to `OptimizedDispatchController`. Unlike the heuristic method, the optimized dispatch method does not use `dispatch_rule_set` as in input in the `tech_config`. This method maximizes the load met while minimizing the cost of the system (operating cost) over each specified time window. The optimized dispatch using Pyomo is implemented differently than the heuristic dispatch in order to be able to properly aggregate the individual Pyomo technology models into a cohesive Pyomo plant model for the optimization solver. Practically, this means that the Pyomo elements of the dispatch (including the individual technology models and the plant model) are not exposed to the main H2I code flow, and do not appear in the N2 diagram. The figure below shows a flow diagram of how the dispatch is implemented. The green blocks below represent what is represented in the N2 diagram of the system. The dispatch routine is currently self-contained within the storage technology of the system, though it includes solving an aggregated plant model in the optimization diff --git a/examples/30_pyomo_optimized_dispatch/tech_config.yaml b/examples/30_pyomo_optimized_dispatch/tech_config.yaml index 947562759..b34b9ddd1 100644 --- a/examples/30_pyomo_optimized_dispatch/tech_config.yaml +++ b/examples/30_pyomo_optimized_dispatch/tech_config.yaml @@ -6,8 +6,6 @@ technologies: model: PYSAMWindPlantPerformanceModel cost_model: model: ATBWindPlantCostModel - dispatch_rule_set: - model: PyomoDispatchGenericConverter resource: type: pysam_wind wind_speed: 9. @@ -37,8 +35,6 @@ technologies: commodity: electricity commodity_rate_units: kW battery: - dispatch_rule_set: - model: PyomoRuleStorageBaseclass control_strategy: model: OptimizedDispatchController performance_model: diff --git a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py index 278228545..f3f02d7b0 100644 --- a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py +++ b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py @@ -79,7 +79,7 @@ def pyomo_setup(self, discrete_inputs): Returns: callable: Function(performance_model, performance_model_kwargs, inputs, commodity) - executing rolling-window heuristic dispatch or optimization and returning: + executing rolling-window heuristic dispatch and returning: (total_out, storage_out, unmet_demand, unused_commodity, soc) """ # initialize the pyomo model diff --git a/h2integrate/control/control_strategies/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index d8a5eeb9b..3cd829dfd 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -139,11 +139,11 @@ def setup(self): self.dispatch_inputs = self.config.make_dispatch_inputs() def pyomo_setup(self, discrete_inputs): - """Create the Pyomo model, attach per-tech Blocks, and return dispatch solver. + """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 heuristic dispatch or optimization and returning: + executing rolling-window optimization to determine dispatch and returning: (total_out, storage_out, unmet_demand, unused_commodity, soc) """ # initialize the pyomo model @@ -154,16 +154,14 @@ def pyomo_setup(self, discrete_inputs): self.source_techs = [] self.dispatch_tech = [] - # run each pyomo rule set up function for each technology 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 + # record source and dispatch techs if source_tech == intended_dispatch_tech: self.dispatch_tech.append(source_tech) - # create pyomo block and set attr self.source_techs.append(source_tech) else: continue @@ -181,7 +179,7 @@ def pyomo_dispatch_solver( Iterates over the full simulation period in chunks of size `self.config.n_control_window`, (re)configures per-window dispatch - parameters, invokes a heuristic control strategy to set fixed + 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. diff --git a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py index 4c66173a7..41f8de74d 100644 --- a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py +++ b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py @@ -109,9 +109,7 @@ def dummy_method(self, in1, in2): def setup(self): """Register per-technology dispatch rule inputs and expose the solver callable. - Adds discrete inputs named 'dispatch_block_rule_function' (and variants - suffixed with source tech names for cross-tech connections) plus a - discrete output 'pyomo_dispatch_solver' that will hold the assembled + Adds discrete a discrete output 'pyomo_dispatch_solver' that will hold the assembled callable after compute(). """ @@ -137,11 +135,11 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): discrete_outputs["pyomo_dispatch_solver"] = self.pyomo_setup(discrete_inputs) def pyomo_setup(self, discrete_inputs): - """Create the Pyomo model, attach per-tech Blocks, and return dispatch solver. + """Create the Pyomo model and return dispatch solver. Returns: callable: Function(performance_model, performance_model_kwargs, inputs, commodity) - executing rolling-window heuristic dispatch or optimization and returning: + executing rolling-window dispatch and returning: (total_out, storage_out, unmet_demand, unused_commodity, soc) """ # initialize the pyomo model @@ -160,9 +158,9 @@ def pyomo_dispatch_solver( Iterates over the full simulation period in chunks of size `self.config.n_control_window`, (re)configures per-window dispatch - parameters, invokes a heuristic control strategy to set fixed - dispatch decisions, and then calls the provided performance_model - over each window to obtain storage output and SOC trajectories. + parameters, applies the chosen control strategy, and then calls the + provided performance_model over each window to obtain storage output and + SOC trajectories. Args: performance_model (callable): From 4ecac615b2924e669dcff25fde65facb4b6f599f Mon Sep 17 00:00:00 2001 From: elenya-grant <116225007+elenya-grant@users.noreply.github.com> Date: Mon, 9 Mar 2026 15:12:50 -0600 Subject: [PATCH 05/21] added docstring for heuristic controller config --- .../heuristic_pyomo_controller.py | 46 +++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py index f3f02d7b0..699d4b133 100644 --- a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py +++ b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py @@ -18,8 +18,54 @@ @define(kw_only=True) class HeuristicLoadFollowingControllerConfig(PyomoControllerBaseConfig): + """Configuration class for the HeuristicLoadFollowingController. + + Attributes: + charge_efficiency (float | None, optional): Efficiency of charging the storage, represented + as a decimal between 0 and 1 (e.g., 0.9 for 90% efficiency). Optional if + `round_trip_efficiency` is provided. + discharge_efficiency (float | None, optional): Efficiency of discharging the storage, + represented as a decimal between 0 and 1 (e.g., 0.9 for 90% efficiency). Optional if + `round_trip_efficiency` is provided. + round_trip_efficiency (float | None, optional): Combined efficiency of charging and + discharging the storage, represented as a decimal between 0 and 1 (e.g., 0.81 for + 81% efficiency). Optional if `charge_efficiency` and `discharge_efficiency` are + provided. + """ + charge_efficiency: float = field(validator=range_val_or_none(0, 1), default=None) discharge_efficiency: float = field(validator=range_val_or_none(0, 1), default=None) + round_trip_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) + + def __attrs_post_init__(self): + """ + Post-initialization logic to validate and calculate efficiencies. + + Ensures that either `charge_efficiency` and `discharge_efficiency` are provided, + or `round_trip_efficiency` is provided. If `round_trip_efficiency` is provided, + it calculates `charge_efficiency` and `discharge_efficiency` as the square root + of `round_trip_efficiency`. + """ + if self.round_trip_efficiency is not None: + if self.charge_efficiency is not None or self.discharge_efficiency is not None: + raise ValueError( + "Exactly one of the following sets of parameters must be set: (a) " + "`round_trip_efficiency`, or (b) both `charge_efficiency` " + "and `discharge_efficiency`." + ) + + # Calculate charge and discharge efficiencies from round-trip efficiency + self.charge_efficiency = np.sqrt(self.round_trip_efficiency) + self.discharge_efficiency = np.sqrt(self.round_trip_efficiency) + elif self.charge_efficiency is not None and self.discharge_efficiency is not None: + # Ensure both charge and discharge efficiencies are provided + pass + else: + raise ValueError( + "Exactly one of the following sets of parameters must be set: (a) " + "`round_trip_efficiency`, or (b) both `charge_efficiency` " + "and `discharge_efficiency`." + ) class HeuristicLoadFollowingController(PyomoControllerBaseClass): From a58c135b53217e115fa37409f61fc8915a908cac Mon Sep 17 00:00:00 2001 From: elenya-grant <116225007+elenya-grant@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:04:09 -0600 Subject: [PATCH 06/21] fixed issue that I made with the HeuristicLoadFollowingControllerConfig --- .../control_strategies/heuristic_pyomo_controller.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py index 699d4b133..fb2460ee8 100644 --- a/h2integrate/control/control_strategies/heuristic_pyomo_controller.py +++ b/h2integrate/control/control_strategies/heuristic_pyomo_controller.py @@ -46,6 +46,8 @@ def __attrs_post_init__(self): it calculates `charge_efficiency` and `discharge_efficiency` as the square root of `round_trip_efficiency`. """ + + super().__attrs_post_init__() if self.round_trip_efficiency is not None: if self.charge_efficiency is not None or self.discharge_efficiency is not None: raise ValueError( @@ -57,15 +59,6 @@ def __attrs_post_init__(self): # Calculate charge and discharge efficiencies from round-trip efficiency self.charge_efficiency = np.sqrt(self.round_trip_efficiency) self.discharge_efficiency = np.sqrt(self.round_trip_efficiency) - elif self.charge_efficiency is not None and self.discharge_efficiency is not None: - # Ensure both charge and discharge efficiencies are provided - pass - else: - raise ValueError( - "Exactly one of the following sets of parameters must be set: (a) " - "`round_trip_efficiency`, or (b) both `charge_efficiency` " - "and `discharge_efficiency`." - ) class HeuristicLoadFollowingController(PyomoControllerBaseClass): From 81743ec832a7ce26e044a174b0f470be5b79acfd Mon Sep 17 00:00:00 2001 From: genevievestarke <103534902+genevievestarke@users.noreply.github.com> Date: Wed, 11 Mar 2026 14:14:14 -0400 Subject: [PATCH 07/21] Update docs/control/pyomo_controllers.md Co-authored-by: Jared Thomas --- docs/control/pyomo_controllers.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/control/pyomo_controllers.md b/docs/control/pyomo_controllers.md index 8f7a35708..228926e74 100644 --- a/docs/control/pyomo_controllers.md +++ b/docs/control/pyomo_controllers.md @@ -21,7 +21,7 @@ For an example of how to use the heuristic Pyomo control framework with the `Heu (optimized-load-following-controller)= ## Optimized Load Following Controller -The optimized dispatch method is specified by setting the storage control to `OptimizedDispatchController`. Unlike the heuristic method, the optimized dispatch method does not use `dispatch_rule_set` as in input in the `tech_config`. This method maximizes the load met while minimizing the cost of the system (operating cost) over each specified time window. +The optimized dispatch method is specified by setting the storage control to `OptimizedDispatchController`. Unlike the heuristic method, the optimized dispatch method does not use `dispatch_rule_set` as an input in the `tech_config`. The `OptimizedDispatchController` method maximizes the load met while minimizing the cost of the system (operating cost) over each specified time window. The optimized dispatch using Pyomo is implemented differently than the heuristic dispatch in order to be able to properly aggregate the individual Pyomo technology models into a cohesive Pyomo plant model for the optimization solver. Practically, this means that the Pyomo elements of the dispatch (including the individual technology models and the plant model) are not exposed to the main H2I code flow, and do not appear in the N2 diagram. The figure below shows a flow diagram of how the dispatch is implemented. The green blocks below represent what is represented in the N2 diagram of the system. The dispatch routine is currently self-contained within the storage technology of the system, though it includes solving an aggregated plant model in the optimization From c7817fa02e588c626765d9426ee03c8736db5b92 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 11 Mar 2026 14:38:17 -0400 Subject: [PATCH 08/21] Make initialize_parameters() functions the same for both controllers --- .../generic_converter_min_operating_cost.py | 10 +++++----- .../control_rules/plant_dispatch_model.py | 15 +++++++-------- .../pyomo_storage_rule_min_operating_cost.py | 12 +++++++----- .../optimized_pyomo_controller.py | 16 +++++++--------- 4 files changed, 26 insertions(+), 27 deletions(-) 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 1dd66d987..6deaa9d10 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 @@ -51,14 +51,14 @@ def __init__( # The units of this are in hours, so half an hour would be 0.5, etc. self.time_duration = [time_duration] * len(self.blocks.index_set()) - def initialize_parameters( - self, commodity_in: list, commodity_demand: list, dispatch_inputs: dict - ): + def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): """Initialize parameters for optimization model Args: - commodity_in (list): List of generated commodity in for this time slice. - commodity_demand (list): The demanded commodity for this time slice. + 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. dispatch_inputs (dict): Dictionary of the dispatch input parameters from config """ diff --git a/h2integrate/control/control_rules/plant_dispatch_model.py b/h2integrate/control/control_rules/plant_dispatch_model.py index 298ff6ccd..a28f29007 100644 --- a/h2integrate/control/control_rules/plant_dispatch_model.py +++ b/h2integrate/control/control_rules/plant_dispatch_model.py @@ -98,21 +98,20 @@ def dispatch_block_rule(self, hybrid, t): ################################## self._create_hybrid_constraints(hybrid, t) - def initialize_parameters( - self, commodity_in: list, commodity_demand: list, dispatch_params: dict - ): + def initialize_parameters(self, inputs: dict, dispatch_params: dict): """Initialize parameters for optimization model - Args: - commodity_in (list): List of generated commodity in for this time slice. - commodity_demand (list): The demanded commodity for this time slice. - dispatch_inputs (dict): Dictionary of the dispatch input parameters from config + 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. + dispatch_inputs (dict): Dictionary of the dispatch input parameters from config """ self.time_weighting_factor = self.time_weighting_factor_input # Discount factor for tech in self.source_techs: pyomo_block = self.tech_dispatch_models.__getattribute__(f"{tech}_rule") - pyomo_block.initialize_parameters(commodity_in, commodity_demand, dispatch_params) + pyomo_block.initialize_parameters(inputs, dispatch_params) def _create_variables_and_ports(self, hybrid, t): """Connect variables and ports from individual technology model 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 7950984aa..0c69ffdc9 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 @@ -53,17 +53,19 @@ def __init__( # The units of this are in hours, so half an hour would be 0.5, etc. self.time_duration = [time_duration] * len(self.blocks.index_set()) - def initialize_parameters( - self, commodity_in: list, commodity_demand: list, dispatch_inputs: dict - ): + def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): """Initialize parameters for optimization model Args: - commodity_in (list): List of generated commodity in for this time slice. - commodity_demand (list): The demanded commodity for this time slice. + 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. dispatch_inputs (dict): Dictionary of the dispatch input parameters from config """ + commodity_demand = inputs[f"{self.commodity_name}_demand"] + # 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"]) diff --git a/h2integrate/control/control_strategies/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index 3cd829dfd..540a28549 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -233,9 +233,7 @@ def pyomo_dispatch_solver( window_start_indices = list(range(0, self.n_timesteps, self.config.n_control_window)) # Initialize parameters for optimized dispatch strategy - self.initialize_parameters( - inputs[f"{commodity_name}_in"], inputs[f"{commodity_name}_demand"] - ) + self.initialize_parameters(inputs) # loop over all control windows, where t is the starting index of each window for t in window_start_indices: @@ -292,12 +290,14 @@ def pyomo_dispatch_solver( return pyomo_dispatch_solver - def initialize_parameters(self, commodity_in, commodity_demand): + def initialize_parameters(self, inputs): """Initialize parameters for optimization model Args: - commodity_in (list): List of generated commodity in for this time slice. - commodity_demand (list): The demanded commodity for this time slice. + 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 @@ -314,9 +314,7 @@ def initialize_parameters(self, commodity_in, commodity_demand): # hybrid_dispatch_rule is the thing where you can access variables and hybrid_rule \ # functions from - self.hybrid_dispatch_rule.initialize_parameters( - commodity_in, commodity_demand, self.dispatch_inputs - ) + 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 From dd7f0787efc3fc7a432114cb3a7d97c96865570b Mon Sep 17 00:00:00 2001 From: genevievestarke <103534902+genevievestarke@users.noreply.github.com> Date: Wed, 11 Mar 2026 14:42:06 -0400 Subject: [PATCH 09/21] Apply suggestion from @jaredthomas68 Co-authored-by: Jared Thomas --- .../control/control_strategies/pyomo_controller_baseclass.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py index 41f8de74d..30dcacb3d 100644 --- a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py +++ b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py @@ -109,7 +109,7 @@ def dummy_method(self, in1, in2): def setup(self): """Register per-technology dispatch rule inputs and expose the solver callable. - Adds discrete a discrete output 'pyomo_dispatch_solver' that will hold the assembled + Adds discrete output 'pyomo_dispatch_solver' that will hold the assembled callable after compute(). """ From 39643da77d024ddc1eb67069dfcff28f9af29c5a Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 11 Mar 2026 14:44:11 -0400 Subject: [PATCH 10/21] Fix doc string format --- .../control/control_strategies/pyomo_controller_baseclass.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py index 30dcacb3d..582fd623b 100644 --- a/h2integrate/control/control_strategies/pyomo_controller_baseclass.py +++ b/h2integrate/control/control_strategies/pyomo_controller_baseclass.py @@ -109,7 +109,7 @@ def dummy_method(self, in1, in2): def setup(self): """Register per-technology dispatch rule inputs and expose the solver callable. - Adds discrete output 'pyomo_dispatch_solver' that will hold the assembled + Adds discrete output 'pyomo_dispatch_solver' that will hold the assembled callable after compute(). """ From e1c226eb20454ff42e29162b2345e29a72a6770d Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Fri, 13 Mar 2026 18:55:55 -0400 Subject: [PATCH 11/21] Initial changes for enabling grid charging for load following --- .../run_pyomo_optimized_dispatch.py | 30 ++++ .../tech_config.yaml | 2 + .../control_rules/plant_dispatch_model.py | 8 +- .../pyomo_storage_rule_min_operating_cost.py | 138 ++++++++++++++++-- .../optimized_pyomo_controller.py | 22 ++- 5 files changed, 178 insertions(+), 22 deletions(-) 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..9aeaad329 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,42 @@ # Create an H2Integrate model model = H2IntegrateModel("pyomo_optimized_dispatch.yaml") +# --- Parameters --- +amplitude = 1.0 # 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 +noise = np.random.Generator(noise_mean, 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.electricity_met_value_in", commodity_met_value_profile, units="MW") +model.prob.set_val("battery.electricity_buy_price_in", commodity_buy_price_profile, units="MW") # Run the model model.run() diff --git a/examples/30_pyomo_optimized_dispatch/tech_config.yaml b/examples/30_pyomo_optimized_dispatch/tech_config.yaml index b34b9ddd1..da7fa0cac 100644 --- a/examples/30_pyomo_optimized_dispatch/tech_config.yaml +++ b/examples/30_pyomo_optimized_dispatch/tech_config.yaml @@ -69,3 +69,5 @@ technologies: 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 + max_system_capacity: 207500 # rating of the wind farm + allow_grid_charging: true diff --git a/h2integrate/control/control_rules/plant_dispatch_model.py b/h2integrate/control/control_rules/plant_dispatch_model.py index a28f29007..e1d587e47 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): """ 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 0c69ffdc9..2e930e39f 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 @@ -61,15 +61,23 @@ def initialize_parameters(self, inputs: dict, dispatch_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. + if allow_grid_charging: + f"{commodity}_met_value_in" : Variable weight for meeting the load + (e.g. could be a grid price) + f"{commodity}_buy_price_in" : Variable cost of energy from the grid dispatch_inputs (dict): Dictionary of the dispatch input parameters from config """ commodity_demand = inputs[f"{self.commodity_name}_demand"] + self.allow_grid_charging = dispatch_inputs["allow_grid_charging"] + + if self.allow_grid_charging: + commodity_met_value_in = inputs[f"{self.commodity_name}_met_value_in"] + commodity_buy_price_in = inputs[f"{self.commodity_name}_buy_price_in"] # 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) @@ -87,7 +95,19 @@ 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.commodity_load_demand = [commodity_demand[t] for t in self.blocks.index_set()] + if self.allow_grid_charging: + # This preserves the possibility of a variable interconnect limit + load_production_value = dispatch_inputs.get("max_system_capacity", None) + self.load_production_limit = [load_production_value for t in self.blocks.index_set()] + self.commodity_met_value = [commodity_met_value_in[t] for t in self.blocks.index_set()] + self.commodity_buy_price = [commodity_buy_price_in[t] for t in self.blocks.index_set()] + + else: + self.set_timeseries_parameter( + "commodity_met_value", dispatch_inputs["commodity_met_value"] + ) self._set_initial_soc_constraint() @@ -246,6 +266,22 @@ def _create_parameters(self, pyomo_model: pyo.ConcreteModel, t): mutable=True, units=pyo_commodity_storage_unit, ) + # Add variables that allow charging from the grid + if self.allow_grid_charging: + 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=pyo_usd_per_commodity_storage_unit_hrs, + ) + pyomo_model.load_production_limit = pyo.Param( + doc=f"Grid inport limit for load [$/{self.commodity_storage_units}]", + default=1000.0, + within=pyo.NonNegativeReals, + mutable=True, + units=pyo_commodity_storage_unit, + ) def _create_variables(self, pyomo_model: pyo.ConcreteModel, t): """Create storage-related decision variables in the Pyomo model. @@ -319,6 +355,14 @@ def _create_variables(self, pyomo_model: pyo.ConcreteModel, t): domain=pyo.Binary, units=pyo.units.dimensionless, ) + # Add variables that allow charging from the grid + if self.allow_grid_charging: + pyomo_model.commodity_bought = pyo.Var( + doc=f"Commodity bought for the system [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + bounds=(0, pyomo_model.load_production_limit), + units=eval("pyo.units." + self.commodity_storage_units), + ) def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): """Create operational and state-of-charge constraints for storage and the system. @@ -336,7 +380,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( @@ -378,6 +422,17 @@ def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): expr=pyomo_model.commodity_out <= pyomo_model.commodity_load_demand * pyomo_model.is_generating, ) + # Add variables that allow charging from the grid + # NOTE: This constraint prevents buying and selling to the grid at the same time + if self.allow_grid_charging: + pyomo_model.purchases_transmission_limit = pyo.Constraint( + doc="Transmission limit on commodity purchases", + expr=( + pyomo_model.commodity_bought + # <= pyomo_model.load_production_limit + <= pyomo_model.load_production_limit * (1 - pyomo_model.is_generating) + ), + ) ################################## # SOC Inventory Constraints # @@ -461,10 +516,15 @@ 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 charging from the grid + if self.allow_grid_charging: + 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 @@ -473,11 +533,27 @@ def update_time_series_parameters( 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. + if allow_grid_charging: + commodity_met_value_in (list): List of variable value of meeting the provided load + commodity_buy_price_in (list): List of variable electricity price from the grid. """ + # 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.commodity_load_demand = [ + time_update_inputs[f"{self.commodity_name}_demand"][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 charging from the grid + if self.allow_grid_charging: + self.commodity_met_value = [ + time_update_inputs[f"{self.commodity_name}_met_value_in"][t] + for t in self.blocks.index_set() + ] + self.commodity_buy_price = [ + time_update_inputs[f"{self.commodity_name}_buy_price_in"][t] + for t in self.blocks.index_set() + ] # Objective functions def min_operating_cost_objective(self, hybrid_blocks, tech_name: str): @@ -504,6 +580,13 @@ def min_operating_cost_objective(self, hybrid_blocks, tech_name: str): ) for t in self.blocks.index_set() ) + if self.allow_grid_charging: + 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 @@ -515,15 +598,27 @@ 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, - } - ) + if self.allow_grid_charging: + 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, + "commodity_bought": hybrid_model.commodity_bought, + } + ) + else: + 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, + } + ) hybrid_model.__setattr__(f"{tech_name}_port", tech_port) return hybrid_model.__getattribute__(f"{tech_name}_port") @@ -556,6 +651,13 @@ def _create_hybrid_variables(self, hybrid_model: pyo.ConcreteModel, tech_name: s domain=pyo.NonNegativeReals, units=pyo_commodity_units, ) + # Add variables that allow charging from the grid + if self.allow_grid_charging: + hybrid_model.commodity_bought = pyo.Var( + doc=f"{self.commodity_name} bought [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=pyo_commodity_units, + ) ################################## # Storage Variables # ################################## @@ -609,6 +711,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/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index 540a28549..263ee8af8 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -61,6 +61,10 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): 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). + max_system_capacity (float): + Maximum amount the storage can charge (i.e. transmission limit for grid charging) + allow_grid_charging (bool): + This sets whether the storage can buy commodity to charge or not """ max_charge_rate: int | float = field() @@ -72,6 +76,9 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): commodity_met_value: float = 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? + max_system_capacity: float = field(default=None) + allow_grid_charging: bool = field(default=False) def make_dispatch_inputs(self): dispatch_keys = [ @@ -85,6 +92,8 @@ def make_dispatch_inputs(self): "charge_efficiency", "discharge_efficiency", "max_charge_rate", + "max_system_capacity", + "allow_grid_charging", ] dispatch_inputs = {k: self.as_dict()[k] for k in dispatch_keys} @@ -237,20 +246,23 @@ def pyomo_dispatch_solver( # loop over all control windows, where t is the starting index of each window for t in window_start_indices: + time_update_inputs = self.create_time_update_dictionary(inputs, t) # 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] + if self.config.allow_grid_charging: + 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( - 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 @@ -316,6 +328,12 @@ def initialize_parameters(self, inputs): # functions from self.hybrid_dispatch_rule.initialize_parameters(inputs, self.dispatch_inputs) + def create_time_update_dictionary(self, inputs, t): + time_update_inputs = { + k: self.as_dict()[k][t : t + self.config.n_control_window] for k in inputs.keys() + } + return time_update_inputs + def update_time_series_parameters( self, commodity_in=None, commodity_demand=None, updated_initial_soc=None ): From b657a44cdec91fc6cd315a808a5a2d4a573003e7 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Fri, 20 Mar 2026 17:09:40 -0400 Subject: [PATCH 12/21] Updates to parameter names and logic --- .../run_pyomo_optimized_dispatch.py | 7 +- .../tech_config.yaml | 2 +- .../generic_converter_min_operating_cost.py | 19 ++++- .../pyomo_storage_rule_min_operating_cost.py | 74 +++++++++-------- .../optimized_pyomo_controller.py | 82 +++++++++++++------ ...optimal_controller_with_generic_storage.py | 1 + .../control/test/test_optimal_controllers.py | 1 + 7 files changed, 119 insertions(+), 67 deletions(-) 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 9aeaad329..a62383149 100644 --- a/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py +++ b/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py @@ -27,7 +27,8 @@ # --- Generate the Random Gaussian Noise --- # Create noise with the same shape as the time vector -noise = np.random.Generator(noise_mean, noise_std_dev, size=t.shape) +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)) @@ -41,8 +42,8 @@ # 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.electricity_met_value_in", commodity_met_value_profile, units="MW") -model.prob.set_val("battery.electricity_buy_price_in", commodity_buy_price_profile, units="MW") +model.prob.set_val("battery.electricity_met_value_in", commodity_met_value_profile, units="USD/kW") +model.prob.set_val("battery.electricity_buy_price_in", commodity_buy_price_profile, units="USD/kW") # Run the model model.run() diff --git a/examples/30_pyomo_optimized_dispatch/tech_config.yaml b/examples/30_pyomo_optimized_dispatch/tech_config.yaml index 0ed6b6c9e..6f1555e6a 100644 --- a/examples/30_pyomo_optimized_dispatch/tech_config.yaml +++ b/examples/30_pyomo_optimized_dispatch/tech_config.yaml @@ -71,4 +71,4 @@ technologies: 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 max_system_capacity: 207500 # rating of the wind farm - allow_grid_charging: true + allow_commodity_buying: true 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 6edf299e3..79600815f 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. + f"{commodity}_met_value_in" : Variable weight for meeting the load + if allow_grid_charging: + f"{commodity}_buy_price_in" : Variable cost of energy from the grid + (e.g. could be a grid price) dispatch_inputs (dict): Dictionary of the dispatch input parameters from config """ @@ -183,18 +187,25 @@ 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. + commodity_met_value_in (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_grid_charging: + commodity_buy_price_in (list): List of variable electricity price from the grid. """ 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/storage/pyomo_storage_rule_min_operating_cost.py b/h2integrate/control/control_rules/storage/pyomo_storage_rule_min_operating_cost.py index 0292f6a0c..bd9d7046f 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,22 +73,22 @@ 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. - if allow_grid_charging: - f"{commodity}_met_value_in" : Variable weight for meeting the load - (e.g. could be a grid price) + f"{commodity}_in" : Available generated commodity profile. + f"{commodity}_demand" : Demanded commodity output profile. + f"{commodity}_met_value_in" : Variable weight for meeting the load + if allow_commodity_buying: f"{commodity}_buy_price_in" : Variable cost of energy from the grid + (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"] - self.allow_grid_charging = dispatch_inputs["allow_grid_charging"] - if self.allow_grid_charging: - commodity_met_value_in = inputs[f"{self.commodity_name}_met_value_in"] + if self.allow_commodity_buying: commodity_buy_price_in = inputs[f"{self.commodity_name}_buy_price_in"] + commodity_met_value_in = inputs[f"{self.commodity_name}_met_value_in"] + # 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"]) @@ -105,18 +110,14 @@ def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): # System parameters self.set_timeseries_from_list("commodity_load_demand", commodity_demand) + self.set_timeseries_from_list("commodity_met_value", commodity_met_value_in) # self.commodity_load_demand = [commodity_demand[t] for t in self.blocks.index_set()] - if self.allow_grid_charging: + if self.allow_commodity_buying: # This preserves the possibility of a variable interconnect limit - load_production_value = dispatch_inputs.get("max_system_capacity", None) - self.load_production_limit = [load_production_value for t in self.blocks.index_set()] - self.commodity_met_value = [commodity_met_value_in[t] for t in self.blocks.index_set()] - self.commodity_buy_price = [commodity_buy_price_in[t] for t in self.blocks.index_set()] - - else: self.set_timeseries_parameter( - "commodity_met_value", dispatch_inputs["commodity_met_value"] + "load_production_limit", dispatch_inputs["max_system_capacity"] ) + self.set_timeseries_from_list("commodity_buy_price", commodity_buy_price_in) self._set_initial_soc_constraint() @@ -271,20 +272,20 @@ def _create_parameters(self, pyomo_model: pyo.ConcreteModel, t): units=self.rate_units_pyo, ) # Add variables that allow charging from the grid - if self.allow_grid_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=pyo_usd_per_commodity_storage_unit_hrs, + units=pyo.unitless, ) pyomo_model.load_production_limit = pyo.Param( doc=f"Grid inport limit for load [$/{self.commodity_storage_units}]", default=1000.0, within=pyo.NonNegativeReals, mutable=True, - units=pyo_commodity_storage_unit, + units=self.cost_units_per_amount_pyo, ) def _create_variables(self, pyomo_model: pyo.ConcreteModel, t): @@ -361,7 +362,7 @@ def _create_variables(self, pyomo_model: pyo.ConcreteModel, t): units=pyo.units.dimensionless, ) # Add variables that allow charging from the grid - if self.allow_grid_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, @@ -429,7 +430,7 @@ def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): ) # Add variables that allow charging from the grid # NOTE: This constraint prevents buying and selling to the grid at the same time - if self.allow_grid_charging: + if self.allow_commodity_buying: pyomo_model.purchases_transmission_limit = pyo.Constraint( doc="Transmission limit on commodity purchases", expr=( @@ -522,7 +523,7 @@ def _create_ports(self, pyomo_model: pyo.ConcreteModel, t): pyomo_model.port.add(pyomo_model.system_load) pyomo_model.port.add(pyomo_model.commodity_out) # Add variables that allow charging from the grid - if self.allow_grid_charging: + if self.allow_commodity_buying: pyomo_model.port.add(pyomo_model.commodity_bought) # Update time series parameters for next optimization window @@ -536,25 +537,28 @@ def update_time_series_parameters( Args: commodity_in (list): List of generated commodity in for this time slice. commodity_demand (list): The demanded commodity for this time slice. + commodity_met_value_in (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_grid_charging: - commodity_met_value_in (list): List of variable value of meeting the provided load + if allow_commodity_buying: commodity_buy_price_in (list): List of variable electricity price from the grid. """ # TODO: Standardize the inputs for this function self.time_duration = [1.0] * len(self.blocks.index_set()) - self.commodity_load_demand = [ - time_update_inputs[f"{self.commodity_name}_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.commodity_met_value = [ + time_update_inputs[f"{self.commodity_name}_met_value_in"][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 charging from the grid - if self.allow_grid_charging: - self.commodity_met_value = [ - time_update_inputs[f"{self.commodity_name}_met_value_in"][t] - for t in self.blocks.index_set() - ] + if self.allow_commodity_buying: self.commodity_buy_price = [ time_update_inputs[f"{self.commodity_name}_buy_price_in"][t] for t in self.blocks.index_set() @@ -585,7 +589,7 @@ def min_operating_cost_objective(self, hybrid_blocks, tech_name: str): ) for t in self.blocks.index_set() ) - if self.allow_grid_charging: + if self.allow_commodity_buying: self.obj += sum( hybrid_blocks[t].time_weighting_factor * self.blocks[t].time_duration @@ -603,7 +607,7 @@ 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. """ - if self.allow_grid_charging: + if self.allow_commodity_buying: tech_port = Port( initialize={ "system_production": hybrid_model.system_production, @@ -656,11 +660,11 @@ def _create_hybrid_variables(self, hybrid_model: pyo.ConcreteModel, tech_name: s units=self.rate_units_pyo, ) # Add variables that allow charging from the grid - if self.allow_grid_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=pyo_commodity_units, + units=self.rate_units_pyo, ) ################################## # Storage Variables # diff --git a/h2integrate/control/control_strategies/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index 61f89cf39..08dd2dfaf 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -63,8 +63,10 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): The default of this parameter is 1.0 (i.e., 1 hour time steps). max_system_capacity (float): Maximum amount the storage can charge (i.e. transmission limit for grid charging) - allow_grid_charging (bool): + 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() @@ -73,12 +75,13 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): 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) + commodity_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? - max_system_capacity: float = field(default=None) - allow_grid_charging: bool = field(default=False) + max_system_capacity: float = field(default=0) + allow_commodity_buying: bool = field(default=None) + commodity_buy_price: int | float | list = field(default=0.0) def make_dispatch_inputs(self): dispatch_keys = [ @@ -93,7 +96,8 @@ def make_dispatch_inputs(self): "discharge_efficiency", "max_charge_rate", "max_system_capacity", - "allow_grid_charging", + "allow_commodity_buying", + "commodity_buy_price", ] dispatch_inputs = {k: self.as_dict()[k] for k in dispatch_keys} @@ -147,6 +151,23 @@ def setup(self): self.dispatch_inputs = self.config.make_dispatch_inputs() + self.add_input( + f"{self.config.commodity}_met_value_in", + val=self.config.commodity_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_in", + val=self.config.commodity_buy_price, + shape=self.n_timesteps, + units="USD/" + self.config.commodity_rate_units, + desc="Value of meeting the demand", + ) + def pyomo_setup(self, discrete_inputs): """Create the Pyomo model, extract dispatch technology names, and return dispatch solver. @@ -232,21 +253,35 @@ def pyomo_dispatch_solver( # get the starting index for each control window window_start_indices = list(range(0, self.n_timesteps, self.config.n_control_window)) + # Check inputs for grid parameters + if self.config.allow_commodity_buying: + # Check grid buy price + if all(inputs[f"{self.config.commodity}_buy_price_in"]) == 0: + raise ValueError( + "commodity_buy_price must be an input and" " >0 if using grid charging" + ) + + # Check max system capacity + if self.config.max_system_capacity == 0: + raise ValueError( + "max_system_capacity must be an input and " ">0 if using grid charging" + ) + # 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: - time_update_inputs = self.create_time_update_dictionary(inputs, t) # 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] - - if self.config.allow_grid_charging: - 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] + 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) @@ -313,14 +348,14 @@ def initialize_parameters(self, inputs): self.hybrid_dispatch_rule.initialize_parameters(inputs, self.dispatch_inputs) def create_time_update_dictionary(self, inputs, t): - time_update_inputs = { - k: self.as_dict()[k][t : t + self.config.n_control_window] for k in inputs.keys() - } + 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] + print("time update inputs:", time_update_inputs.keys(), time_update_inputs) return time_update_inputs - def update_time_series_parameters( - self, commodity_in=None, commodity_demand=None, updated_initial_soc=None - ): + def update_time_series_parameters(self, inputs, updated_initial_soc=None): """Updates the pyomo optimization problem with parameters that change with time Args: @@ -329,9 +364,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, @@ -374,6 +407,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) diff --git a/h2integrate/control/test/test_optimal_controller_with_generic_storage.py b/h2integrate/control/test/test_optimal_controller_with_generic_storage.py index ba8ffd9cd..4bbca6d80 100644 --- a/h2integrate/control/test/test_optimal_controller_with_generic_storage.py +++ b/h2integrate/control/test/test_optimal_controller_with_generic_storage.py @@ -60,6 +60,7 @@ def tech_config_generic(): "time_weighting_factor": 0.995, "system_commodity_interface_limit": 10.0, "n_control_window": 24, + "allow_commodity_buying": False, }, }, } diff --git a/h2integrate/control/test/test_optimal_controllers.py b/h2integrate/control/test/test_optimal_controllers.py index 3a40d161b..6dc155cb3 100644 --- a/h2integrate/control/test/test_optimal_controllers.py +++ b/h2integrate/control/test/test_optimal_controllers.py @@ -62,6 +62,7 @@ def tech_config_generic(): "round_digits": 4, "time_weighting_factor": 0.995, "n_control_window": 24, + "allow_commodity_buying": False, }, }, }, From c6b51a719b9abace5df4355048366a7eba72c92e Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 23 Mar 2026 18:32:11 -0400 Subject: [PATCH 13/21] Update grid inputs with openmdao structure --- .../run_pyomo_optimized_dispatch.py | 6 +- .../tech_config.yaml | 3 +- .../generic_converter_min_operating_cost.py | 4 +- .../pyomo_storage_rule_min_operating_cost.py | 31 ++++---- .../optimized_pyomo_controller.py | 70 +++++++++++++------ ...optimal_controller_with_generic_storage.py | 2 +- .../control/test/test_optimal_controllers.py | 2 +- 7 files changed, 72 insertions(+), 46 deletions(-) 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 a62383149..056b79a0f 100644 --- a/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py +++ b/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py @@ -42,8 +42,8 @@ # 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.electricity_met_value_in", commodity_met_value_profile, units="USD/kW") -model.prob.set_val("battery.electricity_buy_price_in", commodity_buy_price_profile, units="USD/kW") +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") # Run the model model.run() @@ -92,7 +92,7 @@ ) 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.battery_electricity_out", units="MW")[start_hour:end_hour], linestyle="-.", label="Battery Electricity Out (MW)", ) diff --git a/examples/30_pyomo_optimized_dispatch/tech_config.yaml b/examples/30_pyomo_optimized_dispatch/tech_config.yaml index 6f1555e6a..35c81f61f 100644 --- a/examples/30_pyomo_optimized_dispatch/tech_config.yaml +++ b/examples/30_pyomo_optimized_dispatch/tech_config.yaml @@ -66,9 +66,10 @@ technologies: control_parameters: 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 + 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 max_system_capacity: 207500 # rating of the wind farm allow_commodity_buying: true + commodity_buy_price: 0.4 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 79600815f..1a46c9ba2 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 @@ -196,11 +196,11 @@ def update_time_series_parameters( Args: commodity_in (list): List of generated commodity in for this time slice. commodity_demand (list): The demanded commodity for this time slice. - commodity_met_value_in (list): List of variable value of meeting the provided load + commodity_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_grid_charging: - commodity_buy_price_in (list): List of variable electricity price from the grid. + commodity_buy_price (list): List of variable electricity price from the grid. """ self.time_duration = [1.0] * len(self.blocks.index_set()) self.available_production = [ 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 bd9d7046f..a310421c4 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 @@ -75,9 +75,9 @@ def initialize_parameters(self, inputs: dict, dispatch_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}_met_value_in" : Variable weight for meeting the load + "demand_met_value" : Variable weight for meeting the load if allow_commodity_buying: - f"{commodity}_buy_price_in" : Variable cost of energy from the grid + f"{commodity}_buy_price" : Variable cost of energy from the grid (e.g. could be a grid price) dispatch_inputs (dict): Dictionary of the dispatch input parameters from config @@ -85,9 +85,9 @@ def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): commodity_demand = inputs[f"{self.commodity_name}_demand"] if self.allow_commodity_buying: - commodity_buy_price_in = inputs[f"{self.commodity_name}_buy_price_in"] + commodity_buy_price_in = inputs[f"{self.commodity_name}_buy_price"] - commodity_met_value_in = inputs[f"{self.commodity_name}_met_value_in"] + demand_met_value_in = inputs["demand_met_value"] # Dispatch Parameters self.set_timeseries_parameter("cost_per_charge", dispatch_inputs["cost_per_charge"]) @@ -110,7 +110,7 @@ def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): # System parameters self.set_timeseries_from_list("commodity_load_demand", commodity_demand) - self.set_timeseries_from_list("commodity_met_value", commodity_met_value_in) + 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 @@ -257,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, @@ -278,14 +278,14 @@ def _create_parameters(self, pyomo_model: pyo.ConcreteModel, t): default=0.0, within=pyo.Reals, mutable=True, - units=pyo.unitless, + units=self.cost_units_per_amount_pyo, ) pyomo_model.load_production_limit = pyo.Param( doc=f"Grid inport limit for load [$/{self.commodity_storage_units}]", default=1000.0, within=pyo.NonNegativeReals, mutable=True, - units=self.cost_units_per_amount_pyo, + units=self.rate_units_pyo, ) def _create_variables(self, pyomo_model: pyo.ConcreteModel, t): @@ -367,7 +367,7 @@ def _create_variables(self, pyomo_model: pyo.ConcreteModel, t): doc=f"Commodity bought for the system [{self.commodity_storage_units}]", domain=pyo.NonNegativeReals, bounds=(0, pyomo_model.load_production_limit), - units=eval("pyo.units." + self.commodity_storage_units), + units=self.rate_units_pyo, ) def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): @@ -537,11 +537,11 @@ def update_time_series_parameters( Args: commodity_in (list): List of generated commodity in for this time slice. commodity_demand (list): The demanded commodity for this time slice. - commodity_met_value_in (list): List of variable value of meeting the provided load + 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_in (list): List of variable electricity price from the grid. + commodity_buy_price (list): List of variable electricity price from the grid. """ # TODO: Standardize the inputs for this function self.time_duration = [1.0] * len(self.blocks.index_set()) @@ -551,16 +551,15 @@ def update_time_series_parameters( # self.commodity_load_demand = [ # time_update_inputs[f"{self.commodity_name}_demand"][t] for t in self.blocks.index_set() # ] - self.commodity_met_value = [ - time_update_inputs[f"{self.commodity_name}_met_value_in"][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 charging from the grid if self.allow_commodity_buying: self.commodity_buy_price = [ - time_update_inputs[f"{self.commodity_name}_buy_price_in"][t] + time_update_inputs[f"{self.commodity_name}_buy_price"][t] for t in self.blocks.index_set() ] @@ -585,7 +584,7 @@ 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() ) diff --git a/h2integrate/control/control_strategies/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index 08dd2dfaf..037d7f1a1 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -53,7 +53,7 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): 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 (float): 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 @@ -70,25 +70,54 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): """ 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) cost_per_production: float = field(default=None) cost_per_charge: float = field(default=None) cost_per_discharge: float = field(default=None) - commodity_met_value: int | float | list = 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? max_system_capacity: float = field(default=0) - allow_commodity_buying: bool = field(default=None) commodity_buy_price: int | float | list = field(default=0.0) + def __attrs_post_init__(self): + # # Handle demand_met_value list/float here + # if isinstance(self.system_commodity_interface_limit, float | int): + # self.system_commodity_interface_limit = [ + # self.system_commodity_interface_limit + # ] * self.n_control_window + + # Check inputs for grid parameters + if self.allow_commodity_buying: + # Check grid 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 grid charging" + ) + if isinstance(self.commodity_buy_price, list): + if all(self.commodity_buy_price) == 0: + raise ValueError( + "commodity_buy_price must be defined as an input and >0 \ + if using grid charging" + ) + + # Check max system capacity + if self.max_system_capacity == 0: + raise ValueError( + "max_system_capacity must be defined as an input and >0 if using grid charging" + ) + 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", @@ -152,8 +181,8 @@ def setup(self): self.dispatch_inputs = self.config.make_dispatch_inputs() self.add_input( - f"{self.config.commodity}_met_value_in", - val=self.config.commodity_met_value, + "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", @@ -161,7 +190,7 @@ def setup(self): if self.config.allow_commodity_buying: self.add_input( - f"{self.config.commodity}_buy_price_in", + f"{self.config.commodity}_buy_price", val=self.config.commodity_buy_price, shape=self.n_timesteps, units="USD/" + self.config.commodity_rate_units, @@ -253,20 +282,6 @@ def pyomo_dispatch_solver( # get the starting index for each control window window_start_indices = list(range(0, self.n_timesteps, self.config.n_control_window)) - # Check inputs for grid parameters - if self.config.allow_commodity_buying: - # Check grid buy price - if all(inputs[f"{self.config.commodity}_buy_price_in"]) == 0: - raise ValueError( - "commodity_buy_price must be an input and" " >0 if using grid charging" - ) - - # Check max system capacity - if self.config.max_system_capacity == 0: - raise ValueError( - "max_system_capacity must be an input and " ">0 if using grid charging" - ) - # 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) @@ -352,7 +367,13 @@ def create_time_update_dictionary(self, inputs, t): input_keys = inputs.keys() for i in input_keys: time_update_inputs[i] = inputs[i][t : t + self.config.n_control_window] - print("time update inputs:", time_update_inputs.keys(), time_update_inputs) + + 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): @@ -437,6 +458,11 @@ 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) diff --git a/h2integrate/control/test/test_optimal_controller_with_generic_storage.py b/h2integrate/control/test/test_optimal_controller_with_generic_storage.py index 4bbca6d80..6fc517fd5 100644 --- a/h2integrate/control/test/test_optimal_controller_with_generic_storage.py +++ b/h2integrate/control/test/test_optimal_controller_with_generic_storage.py @@ -55,7 +55,7 @@ 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, diff --git a/h2integrate/control/test/test_optimal_controllers.py b/h2integrate/control/test/test_optimal_controllers.py index 6dc155cb3..f5d5ef88a 100644 --- a/h2integrate/control/test/test_optimal_controllers.py +++ b/h2integrate/control/test/test_optimal_controllers.py @@ -58,7 +58,7 @@ def tech_config_generic(): "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, From 21810aa69da20787d7d567848954971087e41974 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 25 Mar 2026 16:26:41 -0400 Subject: [PATCH 14/21] Update constraints --- .../run_pyomo_optimized_dispatch.py | 17 +++++++++++++---- .../tech_config.yaml | 4 ++-- .../pyomo_storage_rule_min_operating_cost.py | 3 ++- 3 files changed, 17 insertions(+), 7 deletions(-) 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 056b79a0f..7bcfdd2e9 100644 --- a/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py +++ b/examples/30_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py @@ -8,7 +8,7 @@ model = H2IntegrateModel("pyomo_optimized_dispatch.yaml") # --- Parameters --- -amplitude = 1.0 # Amplitude of the sine wave +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) @@ -42,7 +42,7 @@ # 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.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") # Run the model @@ -51,7 +51,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 @@ -96,6 +96,8 @@ 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], @@ -104,7 +106,14 @@ ) ax[1].set_ylim([-1e2, 2.5e2]) ax[1].set_ylabel("Electricity Hourly (MW)") -ax[1].set_xlabel("Timestep (hr)") + +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 35c81f61f..7aa591bb3 100644 --- a/examples/30_pyomo_optimized_dispatch/tech_config.yaml +++ b/examples/30_pyomo_optimized_dispatch/tech_config.yaml @@ -11,7 +11,7 @@ technologies: wind_speed: 9. model_inputs: performance_parameters: - num_turbines: 25 + num_turbines: 15 turbine_rating_kw: 8300 rotor_diameter: 196. hub_height: 130. @@ -65,7 +65,7 @@ technologies: opex_fraction: 0.25 # 0.25% of capex per year from 2024 ATB control_parameters: 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 + 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 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 a310421c4..2f0d97f25 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 @@ -420,7 +420,8 @@ def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): pyomo_model.balance = pyo.Constraint( doc="Transmission energy balance", expr=( - pyomo_model.commodity_out == pyomo_model.system_production - pyomo_model.system_load + pyomo_model.commodity_out - pyomo_model.commodity_bought + == pyomo_model.system_production - pyomo_model.system_load ), ) pyomo_model.production_limit = pyo.Constraint( From c4af28067c613de4a365f272d47404ba609eae50 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Wed, 25 Mar 2026 16:27:47 -0400 Subject: [PATCH 15/21] Update constraints --- .../pyomo_storage_rule_min_operating_cost.py | 23 +++++++++++++------ 1 file changed, 16 insertions(+), 7 deletions(-) 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 2f0d97f25..d3996bec1 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 @@ -417,13 +417,22 @@ 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.commodity_bought - == pyomo_model.system_production - pyomo_model.system_load - ), - ) + if self.allow_commodity_buying: + pyomo_model.balance = pyo.Constraint( + doc="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="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", expr=pyomo_model.commodity_out From 6d5b5176967324306209da53dcb76c1e204fd12f Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Fri, 27 Mar 2026 17:36:49 -0400 Subject: [PATCH 16/21] Add tests for commodity buying --- .../optimized_pyomo_controller.py | 50 ++-- .../control/test/test_optimal_controllers.py | 240 +++++++++++++++++- 2 files changed, 265 insertions(+), 25 deletions(-) diff --git a/h2integrate/control/control_strategies/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index f5937cbd4..76cef25c1 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -80,34 +80,34 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): 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? - max_system_capacity: float = field(default=0) - commodity_buy_price: int | float | list = field(default=0.0) + max_system_capacity: float = field(default=None) + commodity_buy_price: int | float | list = field(default=None) def __attrs_post_init__(self): - # # Handle demand_met_value list/float here - # if isinstance(self.system_commodity_interface_limit, float | int): - # self.system_commodity_interface_limit = [ - # self.system_commodity_interface_limit - # ] * self.n_control_window - # Check inputs for grid parameters if self.allow_commodity_buying: - # Check grid 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 grid charging" - ) - if isinstance(self.commodity_buy_price, list): - if all(self.commodity_buy_price) == 0: - raise ValueError( - "commodity_buy_price must be defined as an input and >0 \ - if using grid charging" - ) + if self.commodity_buy_price: + # Check grid 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 grid charging" + ) + 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 grid charging" + ) + else: + raise ValueError( + "commodity_buy_price must be defined as an input and >0 \ + if using grid charging" + ) # Check max system capacity - if self.max_system_capacity == 0: + if self.max_system_capacity == 0 or self.max_system_capacity is None: raise ValueError( "max_system_capacity must be defined as an input and >0 if using grid charging" ) @@ -124,11 +124,13 @@ def make_dispatch_inputs(self): "charge_efficiency", "discharge_efficiency", "max_charge_rate", - "max_system_capacity", "allow_commodity_buying", - "commodity_buy_price", ] + if self.allow_commodity_buying: + dispatch_keys.append("max_system_capacity") + 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 diff --git a/h2integrate/control/test/test_optimal_controllers.py b/h2integrate/control/test/test_optimal_controllers.py index c09ac3f63..8a174e225 100644 --- a/h2integrate/control/test/test_optimal_controllers.py +++ b/h2integrate/control/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.optimized_pyomo_controller import ( OptimizedDispatchController, + OptimizedDispatchControllerConfig, ) @@ -245,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() @@ -627,3 +629,239 @@ 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 = OptimizedDispatchControllerConfig.from_dict(config_data) + + with subtests.test("check commodity_buy_price is None"): + assert config.commodity_buy_price is None + with subtests.test("check max_system_capacity is None"): + assert config.max_system_capacity 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 + OptimizedDispatchControllerConfig.from_dict(data) + + with pytest.raises(ValueError): + data = deepcopy(config_data) + data["allow_commodity_buying"] = True + data["commodity_buy_price"] = 0.4 + OptimizedDispatchControllerConfig.from_dict(data) + + with pytest.raises(ValueError): + data = deepcopy(config_data) + data["allow_commodity_buying"] = True + data["commodity_buy_price"] = 0.4 + data["max_system_capacity"] = 0.0 + OptimizedDispatchControllerConfig.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) + max_system_capacity = 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, + "max_system_capacity": max_system_capacity, + } + + # Setup the OpenMDAO problem and add subsystems + prob = om.Problem() + + prob.model.add_subsystem( + "h2_storage_optimized_load_following_controller", + OptimizedDispatchController( + 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("storage_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")) + + # 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 == storage_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("storage_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("storage_hydrogen_out", units="kg/h")).max() <= capacity + assert np.cumsum(prob.get_val("storage_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(), + -max_system_capacity, + rtol=1e-6, + ) From d78036a7cee61556f143d38e516c93324fe73178 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 30 Mar 2026 16:19:52 -0400 Subject: [PATCH 17/21] Update doc strings --- .../generic_converter_min_operating_cost.py | 13 +++++++------ .../pyomo_storage_rule_min_operating_cost.py | 7 ++++--- 2 files changed, 11 insertions(+), 9 deletions(-) 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 1a46c9ba2..77aa61476 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 @@ -67,9 +67,9 @@ def initialize_parameters(self, inputs: dict, dispatch_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}_met_value_in" : Variable weight for meeting the load - if allow_grid_charging: - f"{commodity}_buy_price_in" : Variable cost of energy from the grid + "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 @@ -196,11 +196,12 @@ def update_time_series_parameters( Args: commodity_in (list): List of generated commodity in for this time slice. commodity_demand (list): The demanded commodity for this time slice. - commodity_met_value (list): List of variable value of meeting the provided load + 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_grid_charging: - commodity_buy_price (list): List of variable electricity price from the grid. + 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 = [ 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 d3996bec1..fcbe73f0b 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 @@ -75,9 +75,9 @@ def initialize_parameters(self, inputs: dict, dispatch_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. - "demand_met_value" : Variable weight for meeting the load + "demand_met_value" : Variable weight for meeting the load if allow_commodity_buying: - f"{commodity}_buy_price" : Variable cost of energy from the grid + 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 @@ -551,7 +551,8 @@ def update_time_series_parameters( 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 electricity price from the grid. + 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()) From 65f8383b6f25c68b2b91834daf644a4b0c9f849d Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 30 Mar 2026 16:31:00 -0400 Subject: [PATCH 18/21] Update pyomo rules for min operating cost --- .../pyomo_storage_rule_min_operating_cost.py | 67 ++++++++----------- 1 file changed, 29 insertions(+), 38 deletions(-) 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 fcbe73f0b..28df5680c 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 @@ -115,7 +115,7 @@ def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): if self.allow_commodity_buying: # This preserves the possibility of a variable interconnect limit self.set_timeseries_parameter( - "load_production_limit", dispatch_inputs["max_system_capacity"] + "commodity_import_limit", dispatch_inputs["max_system_capacity"] ) self.set_timeseries_from_list("commodity_buy_price", commodity_buy_price_in) @@ -271,7 +271,7 @@ def _create_parameters(self, pyomo_model: pyo.ConcreteModel, t): mutable=True, units=self.rate_units_pyo, ) - # Add variables that allow charging from the grid + # 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}]", @@ -280,8 +280,9 @@ def _create_parameters(self, pyomo_model: pyo.ConcreteModel, t): mutable=True, units=self.cost_units_per_amount_pyo, ) - pyomo_model.load_production_limit = pyo.Param( - doc=f"Grid inport limit for load [$/{self.commodity_storage_units}]", + 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, @@ -361,12 +362,12 @@ def _create_variables(self, pyomo_model: pyo.ConcreteModel, t): domain=pyo.Binary, units=pyo.units.dimensionless, ) - # Add variables that allow charging from the grid + # 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.load_production_limit), + bounds=(0, pyomo_model.commodity_import_limit), units=self.rate_units_pyo, ) @@ -419,7 +420,7 @@ def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): ################################## if self.allow_commodity_buying: pyomo_model.balance = pyo.Constraint( - doc="System balance", + doc="Commodity system balance", expr=( pyomo_model.commodity_out - pyomo_model.commodity_bought == pyomo_model.system_production - pyomo_model.system_load @@ -427,26 +428,25 @@ def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): ) else: pyomo_model.balance = pyo.Constraint( - doc="System balance", + 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 charging from the grid - # NOTE: This constraint prevents buying and selling to the grid at the same time + # 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.purchases_transmission_limit = pyo.Constraint( - doc="Transmission limit on commodity purchases", + pyomo_model.upper_purchase_limit = pyo.Constraint( + doc="Upper limit on commodity purchases", expr=( pyomo_model.commodity_bought - # <= pyomo_model.load_production_limit - <= pyomo_model.load_production_limit * (1 - pyomo_model.is_generating) + <= pyomo_model.commodity_import_limit * (1 - pyomo_model.is_generating) ), ) @@ -532,7 +532,7 @@ 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 charging from the grid + # Add variables that allow for buying commodity for storage charging if self.allow_commodity_buying: pyomo_model.port.add(pyomo_model.commodity_bought) @@ -567,7 +567,7 @@ def update_time_series_parameters( ] self.model.initial_soc = updated_initial_soc self.initial_soc = updated_initial_soc - # Add variables that allow charging from the grid + # 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] @@ -617,27 +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. """ + 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: - 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, - "commodity_bought": hybrid_model.commodity_bought, - } - ) - else: - 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["commodity_bought"] = hybrid_model.commodity_bought + + tech_port = Port(initialize_dict) + hybrid_model.__setattr__(f"{tech_name}_port", tech_port) return hybrid_model.__getattribute__(f"{tech_name}_port") @@ -669,7 +660,7 @@ def _create_hybrid_variables(self, hybrid_model: pyo.ConcreteModel, tech_name: s domain=pyo.NonNegativeReals, units=self.rate_units_pyo, ) - # Add variables that allow charging from the grid + # 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}]", From 369eca9fbc4248649ba6c0564eb30b41a6ef84f6 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 30 Mar 2026 17:20:31 -0400 Subject: [PATCH 19/21] Address PR comments --- .../tech_config.yaml | 4 +- .../pyomo_storage_rule_min_operating_cost.py | 4 +- .../optimized_pyomo_controller.py | 47 ++++++++++--------- .../control/test/test_optimal_controllers.py | 12 ++--- 4 files changed, 34 insertions(+), 33 deletions(-) diff --git a/examples/30_pyomo_optimized_dispatch/tech_config.yaml b/examples/30_pyomo_optimized_dispatch/tech_config.yaml index 7aa591bb3..d22e2f5b5 100644 --- a/examples/30_pyomo_optimized_dispatch/tech_config.yaml +++ b/examples/30_pyomo_optimized_dispatch/tech_config.yaml @@ -11,7 +11,7 @@ technologies: wind_speed: 9. model_inputs: performance_parameters: - num_turbines: 15 + num_turbines: 25 turbine_rating_kw: 8300 rotor_diameter: 196. hub_height: 130. @@ -70,6 +70,6 @@ technologies: 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 - max_system_capacity: 207500 # rating of the wind farm + commodity_import_limit: 207500 # rating of the wind farm allow_commodity_buying: true commodity_buy_price: 0.4 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 28df5680c..550617a2f 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 @@ -115,7 +115,7 @@ def initialize_parameters(self, inputs: dict, dispatch_inputs: dict): if self.allow_commodity_buying: # This preserves the possibility of a variable interconnect limit self.set_timeseries_parameter( - "commodity_import_limit", dispatch_inputs["max_system_capacity"] + "commodity_import_limit", dispatch_inputs["commodity_import_limit"] ) self.set_timeseries_from_list("commodity_buy_price", commodity_buy_price_in) @@ -627,7 +627,7 @@ def _create_hybrid_port(self, hybrid_model: pyo.ConcreteModel, tech_name: str): if self.allow_commodity_buying: initialize_dict["commodity_bought"] = hybrid_model.commodity_bought - tech_port = Port(initialize_dict) + tech_port = Port(initialize=initialize_dict) hybrid_model.__setattr__(f"{tech_name}_port", tech_port) diff --git a/h2integrate/control/control_strategies/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index 76cef25c1..1766e4538 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -53,7 +53,7 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): 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 (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,8 +61,9 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): 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). - max_system_capacity (float): - Maximum amount the storage can charge (i.e. transmission limit for grid charging) + 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): @@ -73,6 +74,9 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): 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) @@ -80,7 +84,7 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): 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? - max_system_capacity: float = field(default=None) + commodity_import_limit: float = field(default=None) commodity_buy_price: int | float | list = field(default=None) def __attrs_post_init__(self): @@ -107,9 +111,10 @@ def __attrs_post_init__(self): ) # Check max system capacity - if self.max_system_capacity == 0 or self.max_system_capacity is None: + if self.commodity_import_limit == 0 or self.commodity_import_limit is None: raise ValueError( - "max_system_capacity must be defined as an input and >0 if using grid charging" + "commodity_import_limit must be defined as an input\ + and >0 if using grid charging" ) def make_dispatch_inputs(self): @@ -128,7 +133,7 @@ def make_dispatch_inputs(self): ] if self.allow_commodity_buying: - dispatch_keys.append("max_system_capacity") + dispatch_keys.append("commodity_import_limit") dispatch_keys.append("commodity_buy_price") dispatch_inputs = {k: self.as_dict()[k] for k in dispatch_keys} @@ -166,22 +171,6 @@ def setup(self): self.n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"] - 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() - self.add_input( "demand_met_value", val=self.config.demand_met_value, @@ -199,6 +188,18 @@ def setup(self): desc="Value of meeting the demand", ) + 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): """Create the Pyomo model, extract dispatch technology names, and return dispatch solver. diff --git a/h2integrate/control/test/test_optimal_controllers.py b/h2integrate/control/test/test_optimal_controllers.py index 8a174e225..1a618e838 100644 --- a/h2integrate/control/test/test_optimal_controllers.py +++ b/h2integrate/control/test/test_optimal_controllers.py @@ -658,8 +658,8 @@ def test_optimal_control_config_with_commodity_buying(subtests): with subtests.test("check commodity_buy_price is None"): assert config.commodity_buy_price is None - with subtests.test("check max_system_capacity is None"): - assert config.max_system_capacity is None + with subtests.test("check commodity_import_limit is None"): + assert config.commodity_import_limit is None config_data["allow_commodity_buying"] = True @@ -679,7 +679,7 @@ def test_optimal_control_config_with_commodity_buying(subtests): data = deepcopy(config_data) data["allow_commodity_buying"] = True data["commodity_buy_price"] = 0.4 - data["max_system_capacity"] = 0.0 + data["commodity_import_limit"] = 0.0 OptimizedDispatchControllerConfig.from_dict(data) @@ -690,7 +690,7 @@ def test_optimal_control_with_commodity_buying_generic_storage( 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) - max_system_capacity = 7 + commodity_import_limit = 7 # Set grid charging parameters tech_config_generic["technologies"]["h2_storage"]["model_inputs"]["control_parameters"] = { @@ -704,7 +704,7 @@ def test_optimal_control_with_commodity_buying_generic_storage( "n_control_window": 24, "allow_commodity_buying": True, "commodity_buy_price": 1, - "max_system_capacity": max_system_capacity, + "commodity_import_limit": commodity_import_limit, } # Setup the OpenMDAO problem and add subsystems @@ -862,6 +862,6 @@ def test_optimal_control_with_commodity_buying_generic_storage( 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(), - -max_system_capacity, + -commodity_import_limit, rtol=1e-6, ) From f56aaf8f99b3d2fa0f6f4e00f682cf8a87e20a67 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 30 Mar 2026 17:20:47 -0400 Subject: [PATCH 20/21] Address PR comments --- .../control/control_strategies/optimized_pyomo_controller.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/h2integrate/control/control_strategies/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index 1766e4538..6bd159620 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -113,8 +113,8 @@ def __attrs_post_init__(self): # 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 grid charging" + "commodity_import_limit must be defined as an input and \ + >0 if using grid charging" ) def make_dispatch_inputs(self): From 2152fd016f389c5644468936096b2d82913dc0f3 Mon Sep 17 00:00:00 2001 From: Genevieve Starke Date: Mon, 30 Mar 2026 17:23:13 -0400 Subject: [PATCH 21/21] Fix mentions of grid --- .../control_strategies/optimized_pyomo_controller.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/h2integrate/control/control_strategies/optimized_pyomo_controller.py b/h2integrate/control/control_strategies/optimized_pyomo_controller.py index 6bd159620..6a7b2f8cc 100644 --- a/h2integrate/control/control_strategies/optimized_pyomo_controller.py +++ b/h2integrate/control/control_strategies/optimized_pyomo_controller.py @@ -88,33 +88,33 @@ class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): commodity_buy_price: int | float | list = field(default=None) def __attrs_post_init__(self): - # Check inputs for grid parameters + # Check inputs for commodity buying parameters if self.allow_commodity_buying: if self.commodity_buy_price: - # Check grid 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 grid charging" + 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 grid charging" + if using commodity buying" ) else: raise ValueError( "commodity_buy_price must be defined as an input and >0 \ - if using grid charging" + 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 grid charging" + >0 if using commodity buying" ) def make_dispatch_inputs(self):