From 16cbea2eea43fb888275a5b980921ffa55b8b0b9 Mon Sep 17 00:00:00 2001 From: Manish Khanra <43202860+Manish-Khanra@users.noreply.github.com> Date: Mon, 2 Dec 2024 12:31:08 +0100 Subject: [PATCH] Hydrogen plant and Season hydrogen storage (#468) Add Hydrogen Plant DSM unit and Seasonal Hydrogen Storage DST component --- assume/strategies/__init__.py | 8 +- assume/strategies/naive_strategies.py | 14 +- assume/units/__init__.py | 2 + assume/units/dsm_load_shift.py | 201 +++++++-- assume/units/dst_components.py | 113 ++++- assume/units/hydrogen_plant.py | 319 ++++++++++++++ assume/units/steel_plant.py | 142 +------ docs/source/release_notes.rst | 14 +- .../example_03/industrial_dsm_units.csv | 2 +- .../notebooks/10_DSU_and_flexibility.ipynb | 6 +- tests/test_hydrogen_plant.py | 397 ++++++++++++++++++ tests/test_seasonal_hydrogen_storage.py | 197 +++++++++ tests/test_steel_plant.py | 100 ++++- 13 files changed, 1305 insertions(+), 210 deletions(-) create mode 100644 assume/units/hydrogen_plant.py create mode 100644 tests/test_hydrogen_plant.py create mode 100644 tests/test_seasonal_hydrogen_storage.py diff --git a/assume/strategies/__init__.py b/assume/strategies/__init__.py index de440a8ce..3bddb5062 100644 --- a/assume/strategies/__init__.py +++ b/assume/strategies/__init__.py @@ -12,9 +12,9 @@ flexablePosCRMStorage, ) from assume.strategies.naive_strategies import ( - NaiveDASteelplantStrategy, + NaiveDADSMStrategy, NaiveProfileStrategy, - NaiveRedispatchSteelplantStrategy, + NaiveRedispatchDSMStrategy, NaiveRedispatchStrategy, NaiveSingleBidStrategy, ) @@ -38,8 +38,8 @@ "flexable_neg_crm_storage": flexableNegCRMStorage, "flexable_pos_crm_storage": flexablePosCRMStorage, "naive_redispatch": NaiveRedispatchStrategy, - "naive_da_steelplant": NaiveDASteelplantStrategy, - "naive_steel_redispatch": NaiveRedispatchSteelplantStrategy, + "naive_da_dsm": NaiveDADSMStrategy, + "naive_redispatch_dsm": NaiveRedispatchDSMStrategy, "manual_strategy": SimpleManualTerminalStrategy, "dmas_powerplant": DmasPowerplantStrategy, "dmas_storage": DmasStorageStrategy, diff --git a/assume/strategies/naive_strategies.py b/assume/strategies/naive_strategies.py index 95944c649..647e6aeec 100644 --- a/assume/strategies/naive_strategies.py +++ b/assume/strategies/naive_strategies.py @@ -142,7 +142,12 @@ def calculate_bids( return bids -class NaiveDASteelplantStrategy(BaseStrategy): +class NaiveDADSMStrategy(BaseStrategy): + """ + A naive strategy of a Demand Side Management (DSM) unit. The bid volume is the optimal power requirement of + the unit at the start time of the product. The bid price is the marginal cost of the unit at the start time of the product. + """ + def calculate_bids( self, unit: SupportsMinMax, @@ -176,7 +181,12 @@ def calculate_bids( return bids -class NaiveRedispatchSteelplantStrategy(BaseStrategy): +class NaiveRedispatchDSMStrategy(BaseStrategy): + """ + A naive strategy of a Demand Side Management (DSM) unit that bids the available flexibility of the unit on the redispatch market. + The bid volume is the flexible power requirement of the unit at the start time of the product. The bid price is the marginal cost of the unit at the start time of the product. + """ + def calculate_bids( self, unit: SupportsMinMax, diff --git a/assume/units/__init__.py b/assume/units/__init__.py index 5f5b9bc43..0d924d6a5 100644 --- a/assume/units/__init__.py +++ b/assume/units/__init__.py @@ -7,6 +7,7 @@ from assume.units.powerplant import PowerPlant from assume.units.storage import Storage from assume.units.steel_plant import SteelPlant +from assume.units.hydrogen_plant import HydrogenPlant from assume.units.dst_components import demand_side_technologies unit_types: dict[str, BaseUnit] = { @@ -14,4 +15,5 @@ "demand": Demand, "storage": Storage, "steel_plant": SteelPlant, + "hydrogen_plant": HydrogenPlant, } diff --git a/assume/units/dsm_load_shift.py b/assume/units/dsm_load_shift.py index f475d2204..856a93728 100644 --- a/assume/units/dsm_load_shift.py +++ b/assume/units/dsm_load_shift.py @@ -2,10 +2,16 @@ # # SPDX-License-Identifier: AGPL-3.0-or-later +import logging + import pyomo.environ as pyo +from pyomo.opt import SolverStatus, TerminationCondition +from assume.common.fast_pandas import FastSeries from assume.units.dst_components import demand_side_technologies +logger = logging.getLogger(__name__) + class DSMFlex: def __init__(self, components, **kwargs): @@ -50,6 +56,137 @@ def initialize_components(self): self.model, self.model.dsm_blocks[technology] ) + def switch_to_opt(self, instance): + """ + Switches the instance to solve a cost based optimisation problem by deactivating the flexibility constraints and objective. + + Args: + instance (pyomo.ConcreteModel): The instance of the Pyomo model. + + Returns: + pyomo.ConcreteModel: The modified instance with flexibility constraints and objective deactivated. + """ + # deactivate the flexibility constraints and objective + instance.obj_rule_flex.deactivate() + + instance.total_cost_upper_limit.deactivate() + instance.total_power_input_constraint_with_flex.deactivate() + + return instance + + def switch_to_flex(self, instance): + """ + Switches the instance to flexibility mode by deactivating few constraints and objective function. + + Args: + instance (pyomo.ConcreteModel): The instance of the Pyomo model. + + Returns: + pyomo.ConcreteModel: The modified instance with optimal constraints and objective deactivated. + """ + # deactivate the optimal constraints and objective + instance.obj_rule_opt.deactivate() + instance.total_power_input_constraint.deactivate() + + # fix values of model.total_power_input + for t in instance.time_steps: + instance.total_power_input[t].fix(self.opt_power_requirement.iloc[t]) + instance.total_cost = self.total_cost + + return instance + + def determine_optimal_operation_without_flex(self): + """ + Determines the optimal operation of the steel plant without considering flexibility. + """ + # create an instance of the model + instance = self.model.create_instance() + # switch the instance to the optimal mode by deactivating the flexibility constraints and objective + instance = self.switch_to_opt(instance) + # solve the instance + results = self.solver.solve(instance, options=self.solver_options) + + # Check solver status and termination condition + if (results.solver.status == SolverStatus.ok) and ( + results.solver.termination_condition == TerminationCondition.optimal + ): + logger.debug("The model was solved optimally.") + + # Display the Objective Function Value + objective_value = instance.obj_rule_opt() + logger.debug("The value of the objective function is %s.", objective_value) + + elif results.solver.termination_condition == TerminationCondition.infeasible: + logger.debug("The model is infeasible.") + + else: + logger.debug("Solver Status: ", results.solver.status) + logger.debug( + "Termination Condition: ", results.solver.termination_condition + ) + + opt_power_requirement = [ + pyo.value(instance.total_power_input[t]) for t in instance.time_steps + ] + self.opt_power_requirement = FastSeries( + index=self.index, value=opt_power_requirement + ) + + self.total_cost = sum( + instance.variable_cost[t].value for t in instance.time_steps + ) + + # Variable cost series + variable_cost = [ + pyo.value(instance.variable_cost[t]) for t in instance.time_steps + ] + self.variable_cost_series = FastSeries(index=self.index, value=variable_cost) + + def determine_optimal_operation_with_flex(self): + """ + Determines the optimal operation of the steel plant without considering flexibility. + """ + # create an instance of the model + instance = self.model.create_instance() + # switch the instance to the flexibility mode by deactivating the optimal constraints and objective + instance = self.switch_to_flex(instance) + # solve the instance + results = self.solver.solve(instance, options=self.solver_options) + + # Check solver status and termination condition + if (results.solver.status == SolverStatus.ok) and ( + results.solver.termination_condition == TerminationCondition.optimal + ): + logger.debug("The model was solved optimally.") + + # Display the Objective Function Value + objective_value = instance.obj_rule_flex() + logger.debug("The value of the objective function is %s.", objective_value) + + elif results.solver.termination_condition == TerminationCondition.infeasible: + logger.debug("The model is infeasible.") + + else: + logger.debug("Solver Status: ", results.solver.status) + logger.debug( + "Termination Condition: ", results.solver.termination_condition + ) + + flex_power_requirement = [ + pyo.value(instance.total_power_input[t]) for t in instance.time_steps + ] + self.flex_power_requirement = FastSeries( + index=self.index, value=flex_power_requirement + ) + + # Variable cost series + flex_variable_cost = [ + instance.variable_cost[t].value for t in instance.time_steps + ] + self.flex_variable_cost_series = FastSeries( + index=self.index, value=flex_variable_cost + ) + def flexibility_cost_tolerance(self, model): """ Modify the optimization model to include constraints for flexibility within cost tolerance. @@ -69,53 +206,25 @@ def total_cost_upper_limit(m, t): @model.Constraint(model.time_steps) def total_power_input_constraint_with_flex(m, t): - if self.has_electrolyser: + # Apply constraints based on the technology type + if self.technology == "hydrogen_plant": + # Hydrogen plant constraint return ( m.total_power_input[t] - m.load_shift[t] == self.model.dsm_blocks["electrolyser"].power_in[t] - + self.model.dsm_blocks["eaf"].power_in[t] - + self.model.dsm_blocks["dri_plant"].power_in[t] - ) - else: - return ( - m.total_power_input[t] - m.load_shift[t] - == self.model.dsm_blocks["eaf"].power_in[t] - + self.model.dsm_blocks["dri_plant"].power_in[t] ) - - def recalculate_with_accepted_offers(self, model): - self.reference_power = self.forecaster[f"{self.id}_power"] - self.accepted_pos_capacity = self.forecaster[ - f"{self.id}_power_steelneg" - ] # _accepted_pos_res - - # Parameters - model.reference_power = pyo.Param( - model.time_steps, - initialize={t: value for t, value in enumerate(self.reference_power)}, - ) - - model.accepted_pos_capacity = pyo.Param( - model.time_steps, - initialize={t: value for t, value in enumerate(self.accepted_pos_capacity)}, - ) - - # Variables - model.capacity_upper_bound = pyo.Var( - model.time_steps, within=pyo.NonNegativeReals - ) - - # Constraints - @model.Constraint(model.time_steps) - def capacity_upper_bound_constraint(m, t): - return ( - m.capacity_upper_bound[t] - == m.reference_power[t] - m.accepted_pos_capacity[t] - ) - - @model.Constraint(model.time_steps) - def total_power_upper_limit(m, t): - if m.accepted_pos_capacity[t] > 0: - return m.total_power_input[t] <= m.capacity_upper_bound[t] - else: - return pyo.Constraint.Skip + elif self.technology == "steel_plant": + # Steel plant constraint with conditional electrolyser inclusion + if self.has_electrolyser: + return ( + m.total_power_input[t] - m.load_shift[t] + == self.model.dsm_blocks["electrolyser"].power_in[t] + + self.model.dsm_blocks["eaf"].power_in[t] + + self.model.dsm_blocks["dri_plant"].power_in[t] + ) + else: + return ( + m.total_power_input[t] - m.load_shift[t] + == self.model.dsm_blocks["eaf"].power_in[t] + + self.model.dsm_blocks["dri_plant"].power_in[t] + ) diff --git a/assume/units/dst_components.py b/assume/units/dst_components.py index b23bdfdc0..ae2aff45e 100644 --- a/assume/units/dst_components.py +++ b/assume/units/dst_components.py @@ -1397,7 +1397,7 @@ def charging_profile_constraint(b, t): return model_block -class HydrogenStorage(GenericStorage): +class HydrogenBufferStorage(GenericStorage): """ A class to represent a hydrogen storage unit in an energy system model. @@ -1473,6 +1473,114 @@ def max_discharge_power_constraint(b, t): return model_block +class SeasonalHydrogenStorage(GenericStorage): + """ + A class to represent a seasonal hydrogen storage unit with specific attributes for + seasonal operation. This class internally handles conversion of `time_steps`. + """ + + def __init__( + self, + max_capacity: float, + time_steps: list[int], + horizon: int = 0, + min_capacity: float = 0.0, + max_power_charge: float | None = None, + max_power_discharge: float | None = None, + efficiency_charge: float = 1.0, + efficiency_discharge: float = 1.0, + initial_soc: float = 1.0, + final_soc_target: float = 0.5, + ramp_up: float | None = None, + ramp_down: float | None = None, + storage_loss_rate: float = 0.0, + off_season_start: int = 0, + off_season_end: int = 0, + on_season_start: int = 0, + on_season_end: int = 0, + **kwargs, + ): + horizon = int(horizon) + + super().__init__( + max_capacity=max_capacity, + time_steps=time_steps, + min_capacity=min_capacity, + max_power_charge=max_power_charge, + max_power_discharge=max_power_discharge, + efficiency_charge=efficiency_charge, + efficiency_discharge=efficiency_discharge, + initial_soc=initial_soc, + final_soc_target=final_soc_target, + ramp_up=ramp_up, + ramp_down=ramp_down, + storage_loss_rate=storage_loss_rate, + **kwargs, + ) + + # Convert `time_steps` to a list of integers representing each time step + self.time_steps = time_steps + self.final_soc_target = final_soc_target + + # Check if initial SOC is within the bounds [0, 1] + if initial_soc > 1: + initial_soc /= max_capacity + logger.warning( + f"Initial SOC is greater than 1.0. Setting it to {initial_soc}." + ) + + # Parse comma-separated season start and end values into lists of integers + off_season_start_list = [int(x) for x in off_season_start.split(",")] + off_season_end_list = [int(x) for x in off_season_end.split(",")] + on_season_start_list = [int(x) for x in on_season_start.split(",")] + on_season_end_list = [int(x) for x in on_season_end.split(",")] + + # Generate `off_season` and `on_season` lists based on parsed start and end values + self.off_season = [] + for start, end in zip(off_season_start_list, off_season_end_list): + self.off_season.extend(list(range(start, end + 1))) + + self.on_season = [] + for start, end in zip(on_season_start_list, on_season_end_list): + self.on_season.extend(list(range(start, end + 1))) + + def add_to_model( + self, model: pyo.ConcreteModel, model_block: pyo.Block + ) -> pyo.Block: + """ + Adds seasonal storage block to the Pyomo model, defining constraints based on seasonal operation. + """ + model_block = super().add_to_model(model, model_block) + + # Add final_soc_target as a parameter to model_block + model_block.final_soc_target = pyo.Param(initialize=self.final_soc_target) + + # Seasonal Constraints + @model_block.Constraint(self.off_season) + def off_season_no_discharge(b, t): + """ + Prevent discharging during the off-season. + """ + return b.discharge[t] == 0 + + @model_block.Constraint(self.on_season) + def on_season_no_charge(b, t): + """ + Prevent charging during the on-season. + """ + return b.charge[t] == 0 + + # Final SOC Constraint + @model_block.Constraint() + def final_soc_constraint(b): + """ + Ensure SOC at the end of the time steps meets the target. + """ + return b.soc[self.time_steps[-1]] >= b.final_soc_target * b.max_capacity + + return model_block + + class DRIStorage(GenericStorage): """ A class to represent a Direct Reduced Iron (DRI) storage unit in an energy system model. @@ -1552,7 +1660,8 @@ def max_discharge_power_constraint(b, t): # Mapping of component type identifiers to their respective classes demand_side_technologies: dict = { "electrolyser": Electrolyser, - "hydrogen_storage": HydrogenStorage, + "hydrogen_buffer_storage": HydrogenBufferStorage, + "hydrogen_seasonal_storage": SeasonalHydrogenStorage, "dri_plant": DRIPlant, "dri_storage": DRIStorage, "eaf": ElectricArcFurnace, diff --git a/assume/units/hydrogen_plant.py b/assume/units/hydrogen_plant.py new file mode 100644 index 000000000..8c8e51b16 --- /dev/null +++ b/assume/units/hydrogen_plant.py @@ -0,0 +1,319 @@ +# SPDX-FileCopyrightText: ASSUME Developers +# +# SPDX-License-Identifier: AGPL-3.0-or-later + +import logging + +import pandas as pd +import pyomo.environ as pyo +from pyomo.opt import ( + SolverFactory, + check_available_solvers, +) + +from assume.common.base import SupportsMinMax +from assume.units.dsm_load_shift import DSMFlex + +SOLVERS = ["appsi_highs", "gurobi", "glpk", "cbc", "cplex"] + +logger = logging.getLogger(__name__) + +# Set the log level to ERROR +logging.getLogger("pyomo").setLevel(logging.WARNING) + + +class HydrogenPlant(DSMFlex, SupportsMinMax): + """ + Represents a hydrogen plant in the energy system, including electrolyser and optional seasonal hydrogen storage. + + Args: + id (str): Unique identifier of the plant. + unit_operator (str): Operator of the plant. + bidding_strategies (dict): Bidding strategies. + node (str): Node location of the plant. + index (pd.DatetimeIndex): Time index for plant data. + location (tuple): Plant's geographical location. + components (dict): Components including electrolyser and hydrogen storage. + objective (str): Optimization objective. + flexibility_measure (str): Flexibility measure for load shifting. + demand (float): Hydrogen demand. + cost_tolerance (float): Maximum allowable cost increase. + """ + + required_technologies = ["electrolyser"] + optional_technologies = ["hydrogen_seasonal_storage"] + + def __init__( + self, + id: str, + unit_operator: str, + bidding_strategies: dict, + technology: str = "hydrogen_plant", + node: str = "node0", + index: pd.DatetimeIndex = None, + location: tuple[float, float] = (0.0, 0.0), + components: dict[str, dict] = None, + objective: str = None, + flexibility_measure: str = "max_load_shift", + demand: float = 0, + cost_tolerance: float = 10, + **kwargs, + ): + super().__init__( + id=id, + unit_operator=unit_operator, + technology=technology, + components=components, + bidding_strategies=bidding_strategies, + index=index, + node=node, + location=location, + **kwargs, + ) + + # check if the required components are present in the components dictionary + for component in self.required_technologies: + if component not in components.keys(): + raise ValueError( + f"Component {component} is required for the hydrogen plant unit." + ) + + # check if the provided components are valid and do not contain any unknown components + for component in components.keys(): + if ( + component not in self.required_technologies + and component not in self.optional_technologies + ): + raise ValueError( + f"Components {component} is not a valid component for the hydrogen plant unit." + ) + + # Initialize parameters + self.electricity_price = self.forecaster["price_EOM"] + self.demand = demand + + self.objective = objective + self.flexibility_measure = flexibility_measure + self.cost_tolerance = cost_tolerance + + # Check for the presence of components + self.has_h2seasonal_storage = ( + "hydrogen_seasonal_storage" in self.components.keys() + ) + self.has_electrolyser = "electrolyser" in self.components.keys() + + # Define the Pyomo model + self.model = pyo.ConcreteModel() + self.define_sets() + self.define_parameters() + self.define_variables() + + self.initialize_components() + self.initialize_process_sequence() + + self.define_constraints() + self.define_objective_opt() + + if self.flexibility_measure == "max_load_shift": + self.flexibility_cost_tolerance(self.model) + self.define_objective_flex() + + solvers = check_available_solvers(*SOLVERS) + if len(solvers) < 1: + raise Exception(f"None of {SOLVERS} are available") + + self.solver = SolverFactory(solvers[0]) + self.solver_options = { + "output_flag": False, + "log_to_console": False, + "LogToConsole": 0, + } + + self.opt_power_requirement = None + self.flex_power_requirement = None + + self.variable_cost_series = None + + def define_sets(self): + self.model.time_steps = pyo.Set( + initialize=[idx for idx, _ in enumerate(self.index)] + ) + + def define_parameters(self): + self.model.electricity_price = pyo.Param( + self.model.time_steps, + initialize={t: value for t, value in enumerate(self.electricity_price)}, + ) + self.model.absolute_hydrogen_demand = pyo.Param(initialize=self.demand) + + def define_variables(self): + self.model.total_power_input = pyo.Var( + self.model.time_steps, within=pyo.NonNegativeReals + ) + self.model.variable_cost = pyo.Var( + self.model.time_steps, within=pyo.NonNegativeReals + ) + self.model.hydrogen_demand = pyo.Var( + self.model.time_steps, within=pyo.NonNegativeReals + ) + + def initialize_process_sequence(self): + """ + Initializes the process sequence and constraints for the hydrogen plant. + Distributes hydrogen produced by the electrolyser between hydrogen demand + and optional hydrogen storage. + """ + + @self.model.Constraint(self.model.time_steps) + def hydrogen_production_distribution(m, t): + """ + Balances hydrogen produced by the electrolyser to either satisfy the hydrogen demand + directly, be stored in hydrogen storage, or both if storage is available. + """ + electrolyser_output = self.model.dsm_blocks["electrolyser"].hydrogen_out[t] + + if self.has_h2seasonal_storage: + # With storage: demand can be fulfilled by electrolyser, storage discharge, or both + storage_discharge = self.model.dsm_blocks[ + "hydrogen_seasonal_storage" + ].discharge[t] + storage_charge = self.model.dsm_blocks[ + "hydrogen_seasonal_storage" + ].charge[t] + + # Hydrogen can be supplied to demand and/or storage, and storage can also discharge to meet demand + return ( + electrolyser_output + storage_discharge + == self.model.hydrogen_demand[t] + storage_charge + ) + else: + # Without storage: demand is met solely by electrolyser output + return electrolyser_output == self.model.hydrogen_demand[t] + + def define_constraints(self): + """ + Defines the constraints for the hydrogen plant model, ensuring that the total hydrogen output + over all time steps meets the absolute hydrogen demand. Hydrogen can be sourced from the + electrolyser alone or combined with storage discharge if storage is available. + """ + + @self.model.Constraint() + def total_hydrogen_demand_constraint(m): + """ + Ensures that the total hydrogen output over all time steps meets the absolute hydrogen demand. + If storage is available, the total demand can be fulfilled by both electrolyser output and storage discharge. + If storage is unavailable, the electrolyser output alone must meet the demand. + """ + if self.has_h2seasonal_storage: + # With storage: sum of electrolyser output and storage discharge must meet the total hydrogen demand + return ( + sum( + self.model.dsm_blocks["electrolyser"].hydrogen_out[t] + + self.model.dsm_blocks["hydrogen_seasonal_storage"].discharge[ + t + ] + for t in self.model.time_steps + ) + == self.model.absolute_hydrogen_demand + ) + else: + # Without storage: sum of electrolyser output alone must meet the total hydrogen demand + return ( + sum( + self.model.dsm_blocks["electrolyser"].hydrogen_out[t] + for t in self.model.time_steps + ) + == self.model.absolute_hydrogen_demand + ) + + # Constraint for total power input + @self.model.Constraint(self.model.time_steps) + def total_power_input_constraint(m, t): + """ + Ensures the total power input is the sum of power inputs of all components. + """ + return ( + m.total_power_input[t] + == self.model.dsm_blocks["electrolyser"].power_in[t] + ) + + # Constraint for variable cost per time step + @self.model.Constraint(self.model.time_steps) + def cost_per_time_step(m, t): + """ + Calculates the variable cost per time step. + """ + + return ( + self.model.variable_cost[t] + == self.model.dsm_blocks["electrolyser"].operating_cost[t] + ) + + def define_objective_opt(self): + if self.objective == "min_variable_cost": + + @self.model.Objective(sense=pyo.minimize) + def obj_rule_opt(m): + return sum(self.model.variable_cost[t] for t in self.model.time_steps) + + def define_objective_flex(self): + if self.flexibility_measure == "max_load_shift": + + @self.model.Objective(sense=pyo.maximize) + def obj_rule_flex(m): + return sum(m.load_shift[t] for t in self.model.time_steps) + + def calculate_optimal_operation_if_needed(self): + if ( + self.opt_power_requirement is not None + and self.flex_power_requirement is None + and self.flexibility_measure == "max_load_shift" + ): + self.determine_optimal_operation_with_flex() + + if self.opt_power_requirement is None and self.objective == "min_variable_cost": + self.determine_optimal_operation_without_flex() + + def calculate_marginal_cost(self, start: pd.Timestamp, power: float) -> float: + """ + Calculate the marginal cost of the unit based on the provided time and power. + + Args: + start (pandas.Timestamp): The start time of the dispatch. + power (float): The power output of the unit. + + Returns: + float: the marginal cost of the unit for the given power. + """ + # Initialize marginal cost + marginal_cost = 0 + + if self.opt_power_requirement[start] > 0: + marginal_cost = ( + self.variable_cost_series[start] / self.opt_power_requirement[start] + ) + + return marginal_cost + + def as_dict(self) -> dict: + """ + Returns the attributes of the unit as a dictionary, including specific attributes. + + Returns: + dict: The attributes of the unit as a dictionary. + """ + # Assuming unit_dict is a dictionary that you want to save to the database + components_list = [component for component in self.model.dsm_blocks.keys()] + + # Convert the list to a delimited string + components_string = ",".join(components_list) + + unit_dict = super().as_dict() + unit_dict.update( + { + "unit_type": "demand", + "components": components_string, + } + ) + + return unit_dict diff --git a/assume/units/steel_plant.py b/assume/units/steel_plant.py index d8e100281..fcb41f9ae 100644 --- a/assume/units/steel_plant.py +++ b/assume/units/steel_plant.py @@ -8,13 +8,10 @@ import pyomo.environ as pyo from pyomo.opt import ( SolverFactory, - SolverStatus, - TerminationCondition, check_available_solvers, ) from assume.common.base import SupportsMinMax -from assume.common.fast_pandas import FastSeries from assume.common.forecasts import Forecaster from assume.units.dsm_load_shift import DSMFlex @@ -45,7 +42,7 @@ class SteelPlant(DSMFlex, SupportsMinMax): """ required_technologies = ["dri_plant", "eaf"] - optional_technologies = ["electrolyser", "hydrogen_storage", "dri_storage"] + optional_technologies = ["electrolyser", "hydrogen_buffer_storage", "dri_storage"] def __init__( self, @@ -106,7 +103,7 @@ def __init__( self.cost_tolerance = cost_tolerance # Check for the presence of components - self.has_h2storage = "hydrogen_storage" in self.components.keys() + self.has_h2storage = "hydrogen_buffer_storage" in self.components.keys() self.has_dristorage = "dri_storage" in self.components.keys() self.has_electrolyser = "electrolyser" in self.components.keys() @@ -239,9 +236,9 @@ def direct_hydrogen_flow_constraint(m, t): if self.has_h2storage: return ( self.model.dsm_blocks["electrolyser"].hydrogen_out[t] - + self.model.dsm_blocks["hydrogen_storage"].discharge[t] + + self.model.dsm_blocks["hydrogen_buffer_storage"].discharge[t] == self.model.dsm_blocks["dri_plant"].hydrogen_in[t] - + self.model.dsm_blocks["hydrogen_storage"].charge[t] + + self.model.dsm_blocks["hydrogen_buffer_storage"].charge[t] ) else: return ( @@ -388,137 +385,6 @@ def calculate_optimal_operation_if_needed(self): if self.opt_power_requirement is None and self.objective == "min_variable_cost": self.determine_optimal_operation_without_flex() - def determine_optimal_operation_without_flex(self): - """ - Determines the optimal operation of the steel plant without considering flexibility. - """ - # create an instance of the model - instance = self.model.create_instance() - # switch the instance to the optimal mode by deactivating the flexibility constraints and objective - instance = self.switch_to_opt(instance) - # solve the instance - results = self.solver.solve(instance, options=self.solver_options) - - # Check solver status and termination condition - if (results.solver.status == SolverStatus.ok) and ( - results.solver.termination_condition == TerminationCondition.optimal - ): - logger.debug("The model was solved optimally.") - - # Display the Objective Function Value - objective_value = instance.obj_rule_opt() - logger.debug("The value of the objective function is %s.", objective_value) - - elif results.solver.termination_condition == TerminationCondition.infeasible: - logger.debug("The model is infeasible.") - - else: - logger.debug("Solver Status: ", results.solver.status) - logger.debug( - "Termination Condition: ", results.solver.termination_condition - ) - - opt_power_requirement = [ - pyo.value(instance.total_power_input[t]) for t in instance.time_steps - ] - self.opt_power_requirement = FastSeries( - index=self.index, value=opt_power_requirement - ) - - self.total_cost = sum( - instance.variable_cost[t].value for t in instance.time_steps - ) - - # Variable cost series - variable_cost = [ - pyo.value(instance.variable_cost[t]) for t in instance.time_steps - ] - self.variable_cost_series = FastSeries(index=self.index, value=variable_cost) - - def determine_optimal_operation_with_flex(self): - """ - Determines the optimal operation of the steel plant without considering flexibility. - """ - # create an instance of the model - instance = self.model.create_instance() - # switch the instance to the flexibility mode by deactivating the optimal constraints and objective - instance = self.switch_to_flex(instance) - # solve the instance - results = self.solver.solve(instance, options=self.solver_options) - - # Check solver status and termination condition - if (results.solver.status == SolverStatus.ok) and ( - results.solver.termination_condition == TerminationCondition.optimal - ): - logger.debug("The model was solved optimally.") - - # Display the Objective Function Value - objective_value = instance.obj_rule_flex() - logger.debug("The value of the objective function is %s.", objective_value) - - elif results.solver.termination_condition == TerminationCondition.infeasible: - logger.debug("The model is infeasible.") - - else: - logger.debug("Solver Status: ", results.solver.status) - logger.debug( - "Termination Condition: ", results.solver.termination_condition - ) - - flex_power_requirement = [ - pyo.value(instance.total_power_input[t]) for t in instance.time_steps - ] - self.flex_power_requirement = FastSeries( - index=self.index, value=flex_power_requirement - ) - - # Variable cost series - flex_variable_cost = [ - instance.variable_cost[t].value for t in instance.time_steps - ] - self.flex_variable_cost_series = FastSeries( - index=self.index, value=flex_variable_cost - ) - - def switch_to_opt(self, instance): - """ - Switches the instance to solve a cost based optimisation problem by deactivating the flexibility constraints and objective. - - Args: - instance (pyomo.ConcreteModel): The instance of the Pyomo model. - - Returns: - pyomo.ConcreteModel: The modified instance with flexibility constraints and objective deactivated. - """ - # deactivate the flexibility constraints and objective - instance.obj_rule_flex.deactivate() - - instance.total_cost_upper_limit.deactivate() - instance.total_power_input_constraint_with_flex.deactivate() - - return instance - - def switch_to_flex(self, instance): - """ - Switches the instance to flexibility mode by deactivating few constraints and objective function. - - Args: - instance (pyomo.ConcreteModel): The instance of the Pyomo model. - - Returns: - pyomo.ConcreteModel: The modified instance with optimal constraints and objective deactivated. - """ - # deactivate the optimal constraints and objective - instance.obj_rule_opt.deactivate() - instance.total_power_input_constraint.deactivate() - - # fix values of model.total_power_input - for t in instance.time_steps: - instance.total_power_input[t].fix(self.opt_power_requirement.iloc[t]) - instance.total_cost = self.total_cost - - return instance - def calculate_marginal_cost(self, start: datetime, power: float) -> float: """ Calculate the marginal cost of the unit based on the provided time and power. diff --git a/docs/source/release_notes.rst b/docs/source/release_notes.rst index 101da5c26..e8669c4d3 100644 --- a/docs/source/release_notes.rst +++ b/docs/source/release_notes.rst @@ -14,13 +14,23 @@ Upcoming Release e.g. ``pip install git+https://github.com/assume-framework/assume`` **Improvements:** -- **Performance Optimization:** Switched to a custom `FastSeries` class, which is based on the pandas Series +- **Timeseries Performance Optimization:** Switched to a custom `FastIndex` and `FastSeries` class, which is based on the pandas Series but utilizes NumPy arrays for internal data storage and indexing. This change significantly improves the performance of read and write operations, achieving an average speedup of **2x to 3x** compared to standard pandas Series. The `FastSeries` class retains a close resemblance to the pandas Series, including core functionalities like indexing, slicing, and arithmetic operations. This ensures seamless integration, allowing users to work with the new class without requiring significant code adaptation. -- **Performance Optimization:** Output role handles dict data directly and only converts to DataFrame on Database write. +- **Outputs Role Performance Optimization:** Output role handles dict data directly and only converts to DataFrame on Database write. +- **Overall Performance Optimization:** The overall performance of the framework has been improved by a factor of 5x to 12x + depending on the size of the simulation (number of units, markets, and time steps). + +**New Features:** +- **Hydrogen Plant:** A new demand side unit representing a hydrogen plant has been added. The hydrogen plant consists of an + electrolyzer and a seasonal hydrogen storage unit. The electrolyzer converts electricity into hydrogen, which can be + stored in the hydrogen storage unit and later used. +- **Seasonal Hydrogen Storage:** A new storage unit representing a seasonal hydrogen storage has been added. The seasonal hydrogen + storage unit can store hydrogen over long periods and release it when needed. It has specific constraints to avoid charging or + discharging during off-season or on-season time as well as a target level to be reached at the end of the season. **Bugfixes:** - **Tutorials**: General fixes of the tutorials, to align with updated functionalitites of Assume diff --git a/examples/inputs/example_03/industrial_dsm_units.csv b/examples/inputs/example_03/industrial_dsm_units.csv index 6f4522ce5..2cea5e67e 100644 --- a/examples/inputs/example_03/industrial_dsm_units.csv +++ b/examples/inputs/example_03/industrial_dsm_units.csv @@ -1,4 +1,4 @@ name,unit_type,technology,node,bidding_EOM,bidding_redispatch,unit_operator,objective,flexibility_measure,cost_tolerance,demand,fuel_type,max_power,min_power,ramp_up,ramp_down,min_operating_time,min_down_time,efficiency,specific_dri_demand,specific_electricity_consumption,specific_natural_gas_consumption,specific_hydrogen_consumption,specific_iron_ore_consumption,specific_lime_demand,lime_co2_factor,start_price,max_capacity,min_capacity,initial_soc,storage_loss_rate -A360,steel_plant,electrolyser,north,naive_da_steelplant,,dsm_operator_1,min_variable_cost,max_load_shift,20,223693.1507,hydrogen,922,0,922,922,0,0,0.8,,,,,,,,5,,,, +A360,steel_plant,electrolyser,north,naive_da_dsm,,dsm_operator_1,min_variable_cost,max_load_shift,20,223693.1507,hydrogen,922,0,922,922,0,0,0.8,,,,,,,,5,,,, A360,steel_plant,dri_plant,north,,,,,,,,hydrogen,120,36,120,120,0,0,,,0.3,0.9,1.83,1.43,,,,,,, A360,steel_plant,eaf,north,,,,,,,,hydrogen,162,0,162,162,0,0,,1.09,0.44,,,,0.046,0.1,,,,, diff --git a/examples/notebooks/10_DSU_and_flexibility.ipynb b/examples/notebooks/10_DSU_and_flexibility.ipynb index 9761f6e2a..1afdf1b88 100644 --- a/examples/notebooks/10_DSU_and_flexibility.ipynb +++ b/examples/notebooks/10_DSU_and_flexibility.ipynb @@ -826,7 +826,7 @@ " \"unit_type\": [\"steel_plant\", \"steel_plant\", \"steel_plant\"],\n", " \"technology\": [\"electrolyser\", \"dri_plant\", \"eaf\"],\n", " \"node\": [\"south\", \"south\", \"south\"],\n", - " \"bidding_EOM\": [\"naive_da_steelplant\", \"\", \"\"],\n", + " \"bidding_EOM\": [\"naive_da_dsm\", \"\", \"\"],\n", " \"unit_operator\": [\"dsm_operator_1\", \"dsm_operator_1\", \"dsm_operator_1\"],\n", " \"objective\": [\"min_variable_cost\", \"\", \"\"],\n", " \"flexibility_measure\": [None, \"\", \"\"], # max_load_shift\n", @@ -2014,7 +2014,7 @@ "# Update the 'bidding strategy'\n", "industrial_dsm_units.loc[\n", " industrial_dsm_units[\"technology\"] == \"electrolyser\", \"bidding_EOM\"\n", - "] = \"naive_steel_redispatch\"\n", + "] = \"naive_redispatch_dsm\"\n", "\n", "\n", "# Save the updated industrial_dsm_units file\n", @@ -2719,7 +2719,7 @@ "# Update the 'bidding strategy'\n", "industrial_dsm_units.loc[\n", " industrial_dsm_units[\"technology\"] == \"electrolyser\", \"bidding_EOM\"\n", - "] = \"naive_da_steelplant\"\n", + "] = \"naive_da_dsm\"\n", "\n", "# Save the updated industrial_dsm_units CSV to the new folder\n", "industrial_dsm_units.to_csv(\n", diff --git a/tests/test_hydrogen_plant.py b/tests/test_hydrogen_plant.py new file mode 100644 index 000000000..929c1c1bb --- /dev/null +++ b/tests/test_hydrogen_plant.py @@ -0,0 +1,397 @@ +# SPDX-FileCopyrightText: ASSUME Developers +# +# SPDX-License-Identifier: AGPL-3.0-or-later + + +import pandas as pd +import pyomo.environ as pyo +import pytest + +from assume.common.fast_pandas import FastSeries +from assume.common.forecasts import NaiveForecast +from assume.strategies.naive_strategies import NaiveDADSMStrategy +from assume.units.hydrogen_plant import HydrogenPlant + + +@pytest.fixture +def hydrogen_components(): + # Define full component configuration for hydrogen plant components including electrolyser and seasonal storage + return { + "electrolyser": { + "max_power": 100, # Maximum power input in MW + "min_power": 0, # Minimum power input in MW + "ramp_up": 100, # Ramp-up rate in MW per time step + "ramp_down": 100, # Ramp-down rate in MW per time step + "efficiency": 1, # Efficiency of the electrolyser + "min_operating_time": 0, # Minimum number of operating steps + "min_down_time": 0, # Minimum downtime steps + }, + "hydrogen_seasonal_storage": { + "max_capacity": 500, # Maximum storage capacity in MWh + "min_capacity": 50, # Minimum storage capacity in MWh + "max_power_charge": 30, # Maximum charging power in MW + "max_power_discharge": 30, # Maximum discharging power in MW + "efficiency_charge": 0.9, # Efficiency for charging + "efficiency_discharge": 0.9, # Efficiency for discharging + "initial_soc": 1.0, # Initial state of charge (SOC) as fraction of max capacity + "final_soc_target": 0.5, # Target SOC at end of horizon, as a fraction of max capacity + "ramp_up": 10, # Maximum increase in charge/discharge per step + "ramp_down": 10, + "horizon": 23, # Maximum decrease in charge/discharge per step + "storage_loss_rate": 0.01, # 1% storage loss per time step + "off_season_start": "0", # Off-season start times, as comma-separated values + "off_season_end": "12", # Off-season end times, as comma-separated values + "on_season_start": "13", # On-season start times, as comma-separated values + "on_season_end": "23", # On-season end times, as comma-separated values + }, + } + + +@pytest.fixture +def hydrogen_plant(hydrogen_components) -> HydrogenPlant: + # Define the time index and forecast data + index = pd.date_range("2023-01-01", periods=24, freq="h") + forecast = NaiveForecast( + index, + price_EOM=[60] * 24, # Forecast electricity prices + ) + + # Define a bidding strategy for the hydrogen plant + bidding_strategy = { + "EOM": NaiveDADSMStrategy(), + } + + # Initialize the HydrogenPlant with the specified components, forecast, and strategy + return HydrogenPlant( + id="test_hydrogen_plant", + unit_operator="test_operator", + objective="min_variable_cost", + flexibility_measure="max_load_shift", + cost_tolerance=50, + bidding_strategies=bidding_strategy, + components=hydrogen_components, + forecaster=forecast, + demand=500, # Total hydrogen demand over the horizon + ) + + +def test_optimal_operation_without_flex_initialization(hydrogen_plant): + # Run the initial without-flexibility operation to populate opt_power_requirement + hydrogen_plant.determine_optimal_operation_without_flex() + + # Check that opt_power_requirement is populated and of the correct type + assert ( + hydrogen_plant.opt_power_requirement is not None + ), "opt_power_requirement should be populated" + assert isinstance(hydrogen_plant.opt_power_requirement, FastSeries) + + # Create an instance of the model and switch to optimization mode + instance = hydrogen_plant.model.create_instance() + instance = hydrogen_plant.switch_to_opt(instance) + hydrogen_plant.solver.solve(instance, tee=False) + + # Check that the total power input is greater than zero, indicating operation + total_power_input = sum( + instance.total_power_input[t].value + for t in instance.time_steps + if instance.total_power_input[t].value is not None + ) + assert total_power_input > 0, "Total power input should be greater than zero" + + # Check hydrogen balance after solving + for t in instance.time_steps: + electrolyser_out = instance.dsm_blocks["electrolyser"].hydrogen_out[t].value + hydrogen_demand = ( + instance.hydrogen_demand[t].value + if hasattr(instance.hydrogen_demand[t], "value") + else instance.hydrogen_demand[t] + ) + storage_charge = ( + instance.dsm_blocks["hydrogen_seasonal_storage"].charge[t].value + ) + storage_discharge = ( + instance.dsm_blocks["hydrogen_seasonal_storage"].discharge[t].value + ) + + # Check that the balance holds within a small tolerance + balance_difference = abs( + electrolyser_out - (hydrogen_demand + storage_charge - storage_discharge) + ) + assert ( + balance_difference < 1e-3 + ), f"Hydrogen balance mismatch at time {t}, difference: {balance_difference}" + + # Now test the with-flexibility operation to verify interaction with opt_power_requirement + hydrogen_plant.determine_optimal_operation_with_flex() + + # Re-solve the instance to verify the operation with flexibility + instance = hydrogen_plant.model.create_instance() + instance = hydrogen_plant.switch_to_flex(instance) + hydrogen_plant.solver.solve(instance, tee=False) + + # Verify that opt_power_requirement values are used correctly + for t in instance.time_steps: + assert ( + instance.total_power_input[t].value is not None + ), f"Total power input at time {t} should be set" + expected_power = hydrogen_plant.opt_power_requirement.iloc[t] + actual_power = instance.total_power_input[t].value + assert ( + abs(expected_power - actual_power) < 1e-3 + ), f"Mismatch in power input at time {t}: expected {expected_power}, got {actual_power}" + + +def test_ramping_constraints_without_flex(hydrogen_plant): + # Test that ramping constraints are respected in operation without flexibility + hydrogen_plant.determine_optimal_operation_without_flex() + instance = hydrogen_plant.model.create_instance() + instance = hydrogen_plant.switch_to_opt(instance) + hydrogen_plant.solver.solve(instance, tee=False) + + # Access ramp_up and ramp_down as attributes + ramp_up = hydrogen_plant.components["hydrogen_seasonal_storage"].ramp_up + ramp_down = hydrogen_plant.components["hydrogen_seasonal_storage"].ramp_down + + for t in list(instance.time_steps)[1:]: + charge_prev = ( + instance.dsm_blocks["hydrogen_seasonal_storage"].charge[t - 1].value + ) + charge_curr = instance.dsm_blocks["hydrogen_seasonal_storage"].charge[t].value + discharge_prev = ( + instance.dsm_blocks["hydrogen_seasonal_storage"].discharge[t - 1].value + ) + discharge_curr = ( + instance.dsm_blocks["hydrogen_seasonal_storage"].discharge[t].value + ) + + # Check charge ramping + if charge_prev is not None and charge_curr is not None: + change_in_charge = abs(charge_curr - charge_prev) + assert ( + change_in_charge <= ramp_up + 1e3 + ), f"Charge ramp-up at time {t} exceeds limit" + if discharge_prev is not None and discharge_curr is not None: + change_in_charge = abs(discharge_curr - discharge_prev) + assert ( + change_in_charge <= ramp_down + 1e3 + ), f"Charge ramp-down at time {t} exceeds limit" + + +def test_final_soc_target_without_flex(hydrogen_plant): + hydrogen_plant.determine_optimal_operation_without_flex() + instance = hydrogen_plant.model.create_instance() + instance = hydrogen_plant.switch_to_opt(instance) + hydrogen_plant.solver.solve(instance, tee=False) + + # Access final SOC using integer index + final_step_index = instance.time_steps[-1] + final_soc = ( + instance.dsm_blocks["hydrogen_seasonal_storage"].soc[final_step_index].value + ) + final_soc_target = ( + hydrogen_plant.components["hydrogen_seasonal_storage"].final_soc_target + * hydrogen_plant.components["hydrogen_seasonal_storage"].max_capacity + ) + + assert final_soc is not None, "Final SOC should not be None" + assert final_soc >= final_soc_target, "Final SOC does not meet target" + + +def test_initial_soc_greater_than_capacity(hydrogen_plant): + # After initialization, check if initial SOC > max capacity was adjusted + storage = hydrogen_plant.components["hydrogen_seasonal_storage"] + adjusted_soc = ( + storage.initial_soc * storage.max_capacity + if storage.initial_soc > 1 + else storage.initial_soc + ) + assert ( + adjusted_soc <= storage.max_capacity + ), "Initial SOC should be adjusted if set above max capacity" + + +def test_optimal_operation_with_flex_initialization(hydrogen_plant): + # Set synthetic values for opt_power_requirement and total_cost + hydrogen_plant.opt_power_requirement = FastSeries( + value=30, index=hydrogen_plant.index + ) + hydrogen_plant.total_cost = 100000 # Assign a synthetic cost value for testing + + # Trigger flexibility operation + hydrogen_plant.determine_optimal_operation_with_flex() + + # Verify that flex_power_requirement is populated and is of the correct type + assert ( + hydrogen_plant.flex_power_requirement is not None + ), "flex_power_requirement should be populated" + assert isinstance(hydrogen_plant.flex_power_requirement, FastSeries) + + # Create an instance of the model and switch to flexibility mode + instance = hydrogen_plant.model.create_instance() + instance = hydrogen_plant.switch_to_flex(instance) + hydrogen_plant.solver.solve(instance, tee=False) + + # Check that the total power input reflects flexible operation + total_power_input = sum( + instance.total_power_input[t].value + for t in instance.time_steps + if instance.total_power_input[t].value is not None + ) + assert ( + total_power_input > 0 + ), "Total power input should be greater than zero under flexibility" + + # Verify that the synthetic opt_power_requirement values are being used + for t in instance.time_steps: + expected_power = hydrogen_plant.opt_power_requirement.iloc[t] + actual_power = instance.total_power_input[t].value + assert ( + abs(expected_power - actual_power) < 1e-3 + ), f"Mismatch in power input at time {t}: expected {expected_power}, got {actual_power}" + + # Additional check to ensure the hydrogen demand is met + for t in instance.time_steps: + electrolyser_out = instance.dsm_blocks["electrolyser"].hydrogen_out[t].value + hydrogen_demand = ( + instance.hydrogen_demand[t].value + if hasattr(instance.hydrogen_demand[t], "value") + else instance.hydrogen_demand[t] + ) + storage_charge = ( + instance.dsm_blocks["hydrogen_seasonal_storage"].charge[t].value + ) + storage_discharge = ( + instance.dsm_blocks["hydrogen_seasonal_storage"].discharge[t].value + ) + + # Check that the hydrogen demand balance holds within a small tolerance + balance_difference = abs( + electrolyser_out - (hydrogen_demand + storage_charge - storage_discharge) + ) + assert ( + balance_difference < 1e-3 + ), f"Hydrogen balance mismatch at time {t}, difference: {balance_difference}" + + +def test_unknown_technology_error(): + # Test for unknown technology error with required components present + hydrogen_components = { + "electrolyser": { # Required component + "max_power": 50, + "min_power": 0, + "ramp_up": 10, + "ramp_down": 10, + "efficiency": 0.8, + }, + "unknown_tech": { # Unknown component to trigger the error + "max_power": 20, + "min_power": 0, + "efficiency": 0.5, + }, + } + + with pytest.raises(ValueError, match=r"unknown_tech"): + HydrogenPlant( + id="test_unknown_technology", + unit_operator="test_operator", + objective="min_variable_cost", + flexibility_measure="max_load_shift", + bidding_strategies={"EOM": NaiveDADSMStrategy()}, + components=hydrogen_components, + forecaster=NaiveForecast( + index=pd.date_range("2023-01-01", periods=24, freq="h"), + price_EOM=[60] * 24, + ), + demand=500, + ) + + +@pytest.fixture +def hydrogen_components_no_storage(): + # Define component configuration for hydrogen plant with only an electrolyser (no storage) + return { + "electrolyser": { + "max_power": 100, + "min_power": 0, + "ramp_up": 50, + "ramp_down": 50, + "efficiency": 0.7, + "min_operating_time": 1, + "min_down_time": 1, + } + } + + +@pytest.fixture +def hydrogen_plant_no_storage(hydrogen_components_no_storage) -> HydrogenPlant: + # Define the time index and forecast data + index = pd.date_range("2023-01-01", periods=24, freq="h") + forecast = NaiveForecast( + index, + price_EOM=[60] * 24, # Forecast electricity prices + ) + + # Define a bidding strategy for the hydrogen plant + bidding_strategy = { + "EOM": NaiveDADSMStrategy(), + } + + # Initialize the HydrogenPlant with only an electrolyser component (no storage) + return HydrogenPlant( + id="test_hydrogen_plant_no_storage", + unit_operator="test_operator", + objective="min_variable_cost", + flexibility_measure="max_load_shift", + bidding_strategies=bidding_strategy, + components=hydrogen_components_no_storage, + forecaster=forecast, + demand=800, # Total hydrogen demand over the horizon + ) + + +def test_electrolyser_only_operation(hydrogen_plant_no_storage): + # Test that hydrogen demand is met solely by electrolyser output when there's no storage + hydrogen_plant_no_storage.determine_optimal_operation_without_flex() + instance = hydrogen_plant_no_storage.model.create_instance() + instance = hydrogen_plant_no_storage.switch_to_opt(instance) + hydrogen_plant_no_storage.solver.solve(instance, tee=False) + + for t in instance.time_steps: + electrolyser_out = instance.dsm_blocks["electrolyser"].hydrogen_out[t].value + hydrogen_demand = ( + instance.hydrogen_demand[t].value + if hasattr(instance.hydrogen_demand[t], "value") + else instance.hydrogen_demand[t] + ) + + # Verify that electrolyser output meets hydrogen demand without storage + assert ( + electrolyser_out == hydrogen_demand + ), f"Hydrogen demand not met solely by electrolyser at time {t}" + + +def test_electrolyser_meets_total_hydrogen_demand(hydrogen_plant_no_storage): + # Test that the sum of electrolyser output alone meets total hydrogen demand when there is no storage + hydrogen_plant_no_storage.determine_optimal_operation_without_flex() + instance = hydrogen_plant_no_storage.model.create_instance() + instance = hydrogen_plant_no_storage.switch_to_opt(instance) + hydrogen_plant_no_storage.solver.solve(instance, tee=False) + + # Calculate the total electrolyser output across all time steps + total_electrolyser_output = sum( + instance.dsm_blocks["electrolyser"].hydrogen_out[t].value + for t in instance.time_steps + ) + + # Access absolute hydrogen demand from the model + absolute_hydrogen_demand = pyo.value(instance.absolute_hydrogen_demand) + + # Assert that the total electrolyser output matches the absolute hydrogen demand + assert total_electrolyser_output == absolute_hydrogen_demand, ( + f"Total electrolyser output {total_electrolyser_output} does not match " + f"absolute hydrogen demand {absolute_hydrogen_demand}" + ) + + +if __name__ == "__main__": + pytest.main(["-s", __file__]) diff --git a/tests/test_seasonal_hydrogen_storage.py b/tests/test_seasonal_hydrogen_storage.py new file mode 100644 index 000000000..9eec0090d --- /dev/null +++ b/tests/test_seasonal_hydrogen_storage.py @@ -0,0 +1,197 @@ +# SPDX-FileCopyrightText: ASSUME Developers +# +# SPDX-License-Identifier: AGPL-3.0-or-later + +import pandas as pd +import pyomo.environ as pyo +import pytest + +from assume.units.dst_components import SeasonalHydrogenStorage + +# Define solver name to be used for testing +USE_SOLVER = "appsi_highs" # Update to your specific solver if necessary + + +# Fixture to create a price profile over a 100-step horizon +@pytest.fixture +def price_profile(): + # A price profile with varying values to simulate different cost conditions + return pd.Series([50, 45, 55, 40, 60, 70, 20, 65, 50, 80] * 10, index=range(100)) + + +# Fixture to define seasonal storage configuration for tests +@pytest.fixture +def seasonal_storage_config(): + return { + "max_capacity": 500, # Maximum energy capacity in MWh + "min_capacity": 50, # Minimum energy capacity in MWh + "max_power_charge": 30, # Maximum charging power in MW + "max_power_discharge": 30, # Maximum discharging power in MW + "efficiency_charge": 0.9, # Charging efficiency + "efficiency_discharge": 0.9, # Discharging efficiency + "initial_soc": 0.5, # Initial SOC as a fraction of max capacity + "final_soc_target": 0.6, # SOC target as a fraction of max capacity at the end of the horizon + "ramp_up": 10, # Maximum increase in charge/discharge power per step + "ramp_down": 10, # Maximum decrease in charge/discharge power per step + "storage_loss_rate": 0.01, # 1% storage loss rate per step + # Multiple seasonal periods specified as comma-separated ranges + "off_season_start": "0,50", # Off-season start steps + "off_season_end": "10,60", # Off-season end steps + "on_season_start": "11,61", # On-season start steps + "on_season_end": "49,99", # On-season end steps + } + + +# Fixture to create and solve a model using SeasonalHydrogenStorage +@pytest.fixture +def seasonal_storage_model(seasonal_storage_config, price_profile): + # Define the Pyomo model + model = pyo.ConcreteModel() + model.time_steps = pyo.Set(initialize=range(100)) # 100 time steps + + # Add electricity price parameter to the model + model.electricity_price = pyo.Param( + model.time_steps, initialize=price_profile.to_dict() + ) + + # Initialize the SeasonalHydrogenStorage unit + storage = SeasonalHydrogenStorage( + **seasonal_storage_config, time_steps=model.time_steps, horizon=100 + ) + model.storage = pyo.Block() + storage.add_to_model(model, model.storage) + + # Define an objective to minimize the cost of charging and maximize the benefit from discharging + model.total_cost = pyo.Objective( + expr=sum( + model.storage.charge[t] * model.electricity_price[t] + - model.storage.discharge[t] * model.electricity_price[t] + for t in model.time_steps + ), + sense=pyo.minimize, + ) + + # Solve the model + solver = pyo.SolverFactory(USE_SOLVER) + results = solver.solve(model, tee=False) + + return model, results + + +# Test to check if the model solves successfully +def test_model_solves_successfully(seasonal_storage_model): + model, results = seasonal_storage_model + assert (results.solver.status == pyo.SolverStatus.ok) and ( + results.solver.termination_condition == pyo.TerminationCondition.optimal + ), "Solver did not find an optimal solution." + + +# Test to verify that SOC stays within the defined bounds +def test_state_of_charge_constraints(seasonal_storage_model, seasonal_storage_config): + model, _ = seasonal_storage_model + min_soc = seasonal_storage_config["min_capacity"] + max_soc = seasonal_storage_config["max_capacity"] + + for t in model.time_steps: + soc = pyo.value(model.storage.soc[t]) + assert ( + min_soc <= soc <= max_soc + ), f"SOC at time {t} is {soc}, outside bounds {min_soc} - {max_soc}." + + +# Test to verify that charging and discharging powers do not exceed defined limits +def test_power_bounds(seasonal_storage_model, seasonal_storage_config): + model, _ = seasonal_storage_model + max_charge = seasonal_storage_config["max_power_charge"] + max_discharge = seasonal_storage_config["max_power_discharge"] + + for t in model.time_steps: + charge = pyo.value(model.storage.charge[t]) + discharge = pyo.value(model.storage.discharge[t]) + assert ( + 0 <= charge <= max_charge + ), f"Charge at time {t} exceeds max charge {max_charge}." + assert ( + 0 <= discharge <= max_discharge + ), f"Discharge at time {t} exceeds max discharge {max_discharge}." + + +# Test to ensure no discharging occurs during off-season +def test_off_season_no_discharge(seasonal_storage_model): + model, _ = seasonal_storage_model + off_season = [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 50, + 51, + 52, + 53, + 54, + 55, + 56, + 57, + 58, + 59, + 60, + ] + + for t in off_season: + discharge = pyo.value(model.storage.discharge[t]) + assert discharge == 0, f"Discharge at time {t} during off-season is not zero." + + +# Test to ensure no charging occurs during on-season +def test_on_season_no_charge(seasonal_storage_model): + model, _ = seasonal_storage_model + on_season = list(range(11, 50)) + list(range(61, 99)) + + for t in on_season: + charge = pyo.value(model.storage.charge[t]) + assert charge == 0, f"Charge at time {t} during on-season is not zero." + + +# Test to ensure the SOC at the end of the horizon meets the final SOC target +def test_final_soc_target(seasonal_storage_model, seasonal_storage_config): + model, _ = seasonal_storage_model + final_soc_target = seasonal_storage_config["final_soc_target"] + max_capacity = seasonal_storage_config["max_capacity"] + final_soc = pyo.value(model.storage.soc[99]) + + expected_soc = final_soc_target * max_capacity + assert ( + final_soc >= expected_soc + ), f"Final SOC is {final_soc}, below target {expected_soc}." + + +# Test to verify that the storage loss rate affects the SOC over time +def test_storage_loss_rate(seasonal_storage_model, seasonal_storage_config): + model, _ = seasonal_storage_model + storage_loss_rate = seasonal_storage_config["storage_loss_rate"] + + # Check that SOC reflects the loss rate over time + previous_soc = ( + seasonal_storage_config["initial_soc"] * seasonal_storage_config["max_capacity"] + ) + for t in model.time_steps: + soc = pyo.value(model.storage.soc[t]) + charge = pyo.value(model.storage.charge[t]) + discharge = pyo.value(model.storage.discharge[t]) + + expected_soc = ( + previous_soc + + seasonal_storage_config["efficiency_charge"] * charge + - discharge / seasonal_storage_config["efficiency_discharge"] + - storage_loss_rate * previous_soc + ) + assert ( + abs(soc - expected_soc) < 1e-4 + ), f"SOC at time {t} deviates from expected due to storage loss." + previous_soc = soc diff --git a/tests/test_steel_plant.py b/tests/test_steel_plant.py index 6c8bbecf6..6c0f644ee 100644 --- a/tests/test_steel_plant.py +++ b/tests/test_steel_plant.py @@ -8,8 +8,8 @@ from assume.common.fast_pandas import FastSeries from assume.common.forecasts import NaiveForecast from assume.strategies.naive_strategies import ( - NaiveDASteelplantStrategy, - NaiveRedispatchSteelplantStrategy, + NaiveDADSMStrategy, + NaiveRedispatchDSMStrategy, ) from assume.units.steel_plant import SteelPlant @@ -65,8 +65,8 @@ def steel_plant(dsm_components) -> SteelPlant: ) bidding_strategies = { - "EOM": NaiveDASteelplantStrategy(), - "RD": NaiveRedispatchSteelplantStrategy(), + "EOM": NaiveDADSMStrategy(), + "RD": NaiveRedispatchDSMStrategy(), } return SteelPlant( id="test_steel_plant", @@ -74,10 +74,10 @@ def steel_plant(dsm_components) -> SteelPlant: objective="min_variable_cost", flexibility_measure="max_load_shift", bidding_strategies=bidding_strategies, - index=forecast.index, components=dsm_components, forecaster=forecast, demand=1000, + technology="steel_plant", ) @@ -116,13 +116,6 @@ def test_determine_optimal_operation_without_flex(steel_plant): dri_output == dri_input ), f"DRI output at time {t} does not match DRI input" - for t in instance.time_steps: - dri_output = instance.dsm_blocks["dri_plant"].dri_output[t].value - dri_input = instance.dsm_blocks["eaf"].dri_input[t].value - assert ( - dri_output == dri_input - ), f"Material flow from DRI plant to EAF at time {t} is inconsistent" - total_steel_output = sum( instance.dsm_blocks["eaf"].steel_output[t].value for t in instance.time_steps ) @@ -131,6 +124,31 @@ def test_determine_optimal_operation_without_flex(steel_plant): ), f"Total steel output {total_steel_output} does not match steel demand {instance.steel_demand}" +def test_ramping_constraints(steel_plant): + steel_plant.determine_optimal_operation_without_flex() + instance = steel_plant.model.create_instance() + instance = steel_plant.switch_to_opt(instance) + steel_plant.solver.solve(instance, tee=False) + + # Loop through time steps to check that the electrolyser ramp constraints hold + for t in list(instance.time_steps)[1:]: + # Current and previous power values for electrolyser + power_prev = instance.dsm_blocks["electrolyser"].power_in[t - 1].value + power_curr = instance.dsm_blocks["electrolyser"].power_in[t].value + + # Access ramp constraints using dot notation + ramp_up = steel_plant.components["electrolyser"].ramp_up + ramp_down = steel_plant.components["electrolyser"].ramp_down + + # Verify ramping constraints + assert ( + power_curr - power_prev <= ramp_up + ), f"Electrolyser ramp-up constraint violated at time {t}" + assert ( + power_prev - power_curr <= ramp_down + ), f"Electrolyser ramp-down constraint violated at time {t}" + + def test_determine_optimal_operation_with_flex(steel_plant): # Ensure that the optimal operation without flexibility is determined first steel_plant.determine_optimal_operation_without_flex() @@ -152,5 +170,63 @@ def test_determine_optimal_operation_with_flex(steel_plant): assert total_power_input == 3000 +# New fixture without electrolyser to test line 89 +@pytest.fixture +def steel_plant_without_electrolyser(dsm_components) -> SteelPlant: + index = pd.date_range("2023-01-01", periods=24, freq="h") + forecast = NaiveForecast( + index, + price_EOM=[50] * 24, + fuel_price_natural_gas=[30] * 24, + co2_price=[20] * 24, + ) + + # Remove 'electrolyser' from the components to trigger line 89 + dsm_components_without_electrolyser = dsm_components.copy() + dsm_components_without_electrolyser.pop("electrolyser", None) + + bidding_strategies = { + "EOM": NaiveDADSMStrategy(), + "RD": NaiveRedispatchDSMStrategy(), + } + return SteelPlant( + id="test_steel_plant_without_electrolyser", + unit_operator="test_operator", + objective="min_variable_cost", + flexibility_measure="max_load_shift", + bidding_strategies=bidding_strategies, + components=dsm_components_without_electrolyser, + forecaster=forecast, + demand=1000, + technology="steel_plant", + ) + + +def test_handle_missing_electrolyser(steel_plant_without_electrolyser): + # This should raise an error because 'electrolyser' is required for steel plant + steel_plant_without_electrolyser.determine_optimal_operation_without_flex() + instance = steel_plant_without_electrolyser.model.create_instance() + instance = steel_plant_without_electrolyser.switch_to_opt(instance) + steel_plant_without_electrolyser.solver.solve(instance, tee=False) + + # Loop through time steps to check that the electrolyser ramp constraints hold + for t in list(instance.time_steps)[1:]: + # Current and previous power values for electrolyser + power_prev = instance.dsm_blocks["dri_plant"].power_in[t - 1].value + power_curr = instance.dsm_blocks["dri_plant"].power_in[t].value + + # Access ramp constraints using dot notation + ramp_up = steel_plant_without_electrolyser.components["dri_plant"].ramp_up + ramp_down = steel_plant_without_electrolyser.components["dri_plant"].ramp_down + + # Verify ramping constraints + assert ( + power_curr - power_prev <= ramp_up + ), f"Dri plant ramp-up constraint violated at time {t}" + assert ( + power_prev - power_curr <= ramp_down + ), f"Dri plant ramp-down constraint violated at time {t}" + + if __name__ == "__main__": pytest.main(["-s", __file__])