diff --git a/CHANGELOG.md b/CHANGELOG.md index 1b8f207d5..485b79a53 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ - Added standalone iron DRI and steel EAF performance and cost models - Added capability to have transport models that require user input parameters - Add geologic hydrogen surface processing converter +- Add optimal dispatch of storage for load following - Add baseclass for caching functionality - Added postprocessing function to save timeseries - Minor reorg for profast tools diff --git a/docs/control/pyomo_controllers.md b/docs/control/pyomo_controllers.md index 18de1d99d..bba16e924 100644 --- a/docs/control/pyomo_controllers.md +++ b/docs/control/pyomo_controllers.md @@ -2,15 +2,25 @@ # Pyomo control framework [Pyomo](https://www.pyomo.org/about) is an open-source optimization software package. It is used in H2Integrate to facilitate modeling and solving control problems, specifically to determine optimal dispatch strategies for dispatchable technologies. -Pyomo control, allows for the possibility of feedback control at specified intervals, but can also be used for open-loop control if desired. In the pyomo control framework in H2Integrate, each technology can have control rules associated with them that are in turn passed to the pyomo control component, which is owned by the storage technology. The pyomo control component combines the technology rules into a single pyomo model, which is then passed to the storage technology performance model inside a callable dispatch function. The dispatch function also accepts a simulation method from the performance model and iterates between the pyomo model for dispatch commands and the performance simulation function to simulated performance with the specified commands. The dispatch function runs in specified time windows for dispatch and performance until the whole simulation time has been run. +Pyomo control allows for the possibility of feedback control at specified intervals, but can also be used for open-loop control if desired. In the pyomo control framework in H2Integrate, each technology can have control rules associated with them that are in turn passed to the pyomo control component, which is owned by the storage technology. The pyomo control component combines the technology rules into a single pyomo model, which is then passed to the storage technology performance model inside a callable dispatch function. The dispatch function also accepts a simulation method from the performance model and iterates between the pyomo model for dispatch commands and the performance simulation function to simulate performance with the specified commands. The dispatch function runs in specified time windows for dispatch and performance until the whole simulation time has been run. An example of an N2 diagram for a system using the pyomo control framework for hydrogen storage and dispatch is shown below ([click here for an interactive version](./figures/pyomo-n2.html)). Note the control rules being passed to the dispatch component and the dispatch function, containing the full pyomo model, being passed to the performance model for the battery/storage technology. Another important thing to recognize, in contrast to the open-loop control framework, is that the storage technology outputs (commodity out, SOC, unused commodity, etc) are passed out of the performance model when using the Pyomo control framework rather than from the control component. ![](./figures/pyomo-n2.png) +The pyomo control framework currently supports both a simple heuristic method and an optimized dispatch method for load following control. + (heuristic-load-following-controller)= ## Heuristic Load Following Controller -The pyomo control framework currently supports only a simple heuristic method, `HeuristicLoadFollowingController`, but we plan to extend the framework to be able to run a full dispatch optimization using a pyomo solver. When using the pyomo framework, a `dispatch_rule_set` for each technology connected to the storage technology must also be specified. These will typically be `PyomoDispatchGenericConverter` for generating technologies, and `PyomoRuleStorageBaseclass` for storage technologies. More complex rule sets may be developed as needed. -For an example of how to use the pyomo control framework with the `HeuristicLoadFollowingController`, see +The simple heuristic method is specified by setting the storage control to `HeuristicLoadFollowingController`. When using the pyomo framework, a `dispatch_rule_set` for each technology connected to the storage technology must also be specified. These will typically be `PyomoDispatchGenericConverter` for generating technologies, and `PyomoRuleStorageBaseclass` for storage technologies. More complex rule sets may be developed as needed. + +For an example of how to use the heuristic pyomo control framework with the `HeuristicLoadFollowingController`, see - `examples/18_pyomo_heuristic_wind_battery_dispatch` + +(optimized-load-following-controller)= +## Optimized Load Following Controller +The optimized dispatch method is specified by setting the storage control to `optimized_dispatch_controller`. 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`. This method maximizes the load met while minimizing the cost of the system (operating cost) over each specified time window. + +For an example of how to use the optimized pyomo control framework with the `optimized_dispatch_controller`, see +- `examples/27_pyomo_optimized_dispatch` diff --git a/examples/01_onshore_steel_mn/tech_config.yaml b/examples/01_onshore_steel_mn/tech_config.yaml index 22811d68a..b70c886bc 100644 --- a/examples/01_onshore_steel_mn/tech_config.yaml +++ b/examples/01_onshore_steel_mn/tech_config.yaml @@ -82,6 +82,7 @@ technologies: model_inputs: shared_parameters: commodity_name: "electricity" + commodity_storage_units: "kW" max_charge_rate: 375740.4 #kW max_capacity: 375745.2 #kWh n_control_window: 24 @@ -100,11 +101,7 @@ technologies: power_capex: 311 # $/kW from 2024 ATB year 2025 opex_fraction: 0.024999840573439444 control_parameters: - commodity_storage_units: "kW" tech_name: "battery" - dispatch_rule_parameters: - commodity_name: "electricity" - commodity_storage_units: "kW" electrolyzer: performance_model: diff --git a/examples/02_texas_ammonia/tech_config.yaml b/examples/02_texas_ammonia/tech_config.yaml index d2cc1d20d..c523a8f6c 100644 --- a/examples/02_texas_ammonia/tech_config.yaml +++ b/examples/02_texas_ammonia/tech_config.yaml @@ -80,6 +80,7 @@ technologies: model_inputs: shared_parameters: commodity_name: "electricity" + commodity_storage_units: "kW" max_charge_rate: 96.0 #kW max_capacity: 96.0 #kWh n_control_window: 24 @@ -98,11 +99,8 @@ technologies: power_capex: 311 # $/kW from 2024 ATB year 2025 opex_fraction: 0.025 control_parameters: - commodity_storage_units: "kW" tech_name: "battery" - dispatch_rule_parameters: - commodity_name: "electricity" - commodity_storage_units: "kW" + electrolyzer: performance_model: model: "ECOElectrolyzerPerformanceModel" diff --git a/examples/09_co2/direct_ocean_capture/tech_config.yaml b/examples/09_co2/direct_ocean_capture/tech_config.yaml index dad6870c1..d5f6a637b 100644 --- a/examples/09_co2/direct_ocean_capture/tech_config.yaml +++ b/examples/09_co2/direct_ocean_capture/tech_config.yaml @@ -73,6 +73,7 @@ technologies: model_inputs: shared_parameters: commodity_name: "electricity" + commodity_storage_units: "kW" max_charge_rate: 50000 #kW max_capacity: 200000 #kWh n_control_window: 24 @@ -91,11 +92,8 @@ technologies: power_capex: 317 # $/kW from 2024 ATB year 2025 opex_fraction: 0.02536510376633359 control_parameters: - commodity_storage_units: "kW" tech_name: "battery" - dispatch_rule_parameters: - commodity_name: "electricity" - commodity_storage_units: "kW" + doc: performance_model: model: "DOCPerformanceModel" diff --git a/examples/09_co2/ocean_alkalinity_enhancement/tech_config.yaml b/examples/09_co2/ocean_alkalinity_enhancement/tech_config.yaml index 7670f0f98..59e71b81c 100644 --- a/examples/09_co2/ocean_alkalinity_enhancement/tech_config.yaml +++ b/examples/09_co2/ocean_alkalinity_enhancement/tech_config.yaml @@ -50,6 +50,7 @@ technologies: model_inputs: shared_parameters: commodity_name: "electricity" + commodity_storage_units: "kW" max_charge_rate: 50000 #kW max_capacity: 200000 #kWh n_control_window: 24 @@ -68,11 +69,8 @@ technologies: power_capex: 317 # $/kW from 2024 ATB year 2025 opex_fraction: 0.02536510376633359 control_parameters: - commodity_storage_units: "kW" tech_name: "battery" - dispatch_rule_parameters: - commodity_name: "electricity" - commodity_storage_units: "kW" + oae: performance_model: model: "OAEPerformanceModel" diff --git a/examples/12_ammonia_synloop/tech_config.yaml b/examples/12_ammonia_synloop/tech_config.yaml index ea41e075a..2138ad394 100644 --- a/examples/12_ammonia_synloop/tech_config.yaml +++ b/examples/12_ammonia_synloop/tech_config.yaml @@ -80,6 +80,7 @@ technologies: model_inputs: shared_parameters: commodity_name: "electricity" + commodity_storage_units: "kW" max_charge_rate: 96.0 #kW max_capacity: 96.0 #kWh n_control_window: 24 @@ -98,11 +99,8 @@ technologies: power_capex: 311 # $/kW from 2024 ATB year 2025 opex_fraction: 0.025 control_parameters: - commodity_storage_units: "kW" tech_name: "battery" - dispatch_rule_parameters: - commodity_name: "electricity" - commodity_storage_units: "kW" + electrolyzer: performance_model: model: "ECOElectrolyzerPerformanceModel" diff --git a/examples/18_pyomo_heuristic_dispatch/run_pyomo_heuristic_dispatch.py b/examples/18_pyomo_heuristic_dispatch/run_pyomo_heuristic_dispatch.py index 642dd03a7..8c3aeef24 100644 --- a/examples/18_pyomo_heuristic_dispatch/run_pyomo_heuristic_dispatch.py +++ b/examples/18_pyomo_heuristic_dispatch/run_pyomo_heuristic_dispatch.py @@ -67,7 +67,7 @@ range(start_hour, end_hour), demand_profile[start_hour:end_hour], linestyle="--", - label="Eletrical Demand (MW)", + label="Electrical Demand (MW)", ) ax[1].set_ylim([-7e2, 7e2]) ax[1].set_ylabel("Electricity Hourly (MW)") diff --git a/examples/18_pyomo_heuristic_dispatch/tech_config.yaml b/examples/18_pyomo_heuristic_dispatch/tech_config.yaml index 70f5029ae..d7765b996 100644 --- a/examples/18_pyomo_heuristic_dispatch/tech_config.yaml +++ b/examples/18_pyomo_heuristic_dispatch/tech_config.yaml @@ -46,6 +46,7 @@ technologies: model_inputs: shared_parameters: commodity_name: "electricity" + commodity_storage_units: "kW" max_charge_rate: 100000 max_capacity: 500000 n_control_window: 24 @@ -64,8 +65,4 @@ technologies: power_capex: 311 # $/kW from 2024 ATB year 2025 opex_fraction: 0.25 # 0.25% of capex per year from 2024 ATB control_parameters: - commodity_storage_units: "kW" tech_name: "battery" - dispatch_rule_parameters: - commodity_name: "electricity" - commodity_storage_units: "kW" diff --git a/examples/18_pyomo_heuristic_dispatch/tech_config_error_for_testing.yaml b/examples/18_pyomo_heuristic_dispatch/tech_config_error_for_testing.yaml index 634bd1e7a..05ba61664 100644 --- a/examples/18_pyomo_heuristic_dispatch/tech_config_error_for_testing.yaml +++ b/examples/18_pyomo_heuristic_dispatch/tech_config_error_for_testing.yaml @@ -48,6 +48,8 @@ technologies: model: "ATBBatteryCostModel" model_inputs: shared_parameters: + commodity_name: "electricity" + commodity_storage_units: "kW" max_charge_rate: 100000 max_capacity: 500000 n_control_window: 24 @@ -65,9 +67,4 @@ technologies: power_capex: 311 # $/kW from 2024 ATB year 2025 opex_fraction: 0.25 # 0.25% of capex per year from 2024 ATB control_parameters: - commodity_name: "electricity" - commodity_storage_units: "kW" tech_name: "wrong_tech_name" - dispatch_rule_parameters: - commodity_name: "electricity" - commodity_storage_units: "kW" diff --git a/examples/27_pyomo_optimized_dispatch/driver_config.yaml b/examples/27_pyomo_optimized_dispatch/driver_config.yaml new file mode 100644 index 000000000..6589a3934 --- /dev/null +++ b/examples/27_pyomo_optimized_dispatch/driver_config.yaml @@ -0,0 +1,5 @@ +name: "driver_config" +description: "This analysis runs a hybrid plant to dispatch storage optimally to meet an electrical load." + +general: + folder_output: outputs diff --git a/examples/27_pyomo_optimized_dispatch/plant_config.yaml b/examples/27_pyomo_optimized_dispatch/plant_config.yaml new file mode 100644 index 000000000..76fa25a40 --- /dev/null +++ b/examples/27_pyomo_optimized_dispatch/plant_config.yaml @@ -0,0 +1,76 @@ +name: "plant_config" +description: "This plant is located in TX, USA..." + +sites: + site: + latitude: 35.2018863 + longitude: -101.945027 + + resources: + wind_resource: + resource_model: "WTKNRELDeveloperAPIWindResource" + resource_parameters: + resource_year: 2012 + +plant: + plant_life: 30 + +# array of arrays containing left-to-right technology +# interconnections; can support bidirectional connections +# with the reverse definition. +# this will naturally grow as we mature the interconnected tech +technology_interconnections: [ + ["wind", "battery", "electricity", "cable"], +] + +# array of arrays containing left-to-right technology, technology doing the dispatching +# in this case, battery is connected to battery because there are controls rules for +# the battery and battery is controlling the dispatching +tech_to_dispatch_connections: [ + ["wind", "battery"], + ["battery", "battery"], +] + +resource_to_tech_connections: [ + # connect the wind resource to the wind technology + ['site.wind_resource', 'wind', 'wind_resource_data'], +] + +finance_parameters: + finance_groups: + profast_model: + commodity: "electricity" + finance_model: "ProFastLCO" + model_inputs: + params: + analysis_start_year: 2032 + installation_time: 36 # months + inflation_rate: 0.0 # 0 for nominal analysis + discount_rate: 0.09 # nominal return based on 2024 ATB baseline workbook for land-based wind + debt_equity_ratio: 2.62 # 2024 ATB uses 72.4% debt for land-based wind + property_tax_and_insurance: 0.03 # p-tax https://www.house.mn.gov/hrd/issinfo/clsrates.aspx # insurance percent of CAPEX estimated based on https://www.nrel.gov/docs/fy25osti/91775.pdf + total_income_tax_rate: 0.257 # 0.257 tax rate in 2024 atb baseline workbook, value here is based on federal (21%) and state in MN (9.8) + capital_gains_tax_rate: 0.15 # H2FAST default + sales_tax_rate: 0.07375 # total state and local sales tax in St. Louis County https://taxmaps.state.mn.us/salestax/ + debt_interest_rate: 0.07 # based on 2024 ATB nominal interest rate for land-based wind + debt_type: "Revolving debt" # can be "Revolving debt" or "One time loan". Revolving debt is H2FAST default and leads to much lower LCOH + loan_period_if_used: 0 # H2FAST default, not used for revolving debt + cash_onhand_months: 1 # H2FAST default + admin_expense: 0.00 # percent of sales H2FAST default + capital_items: + depr_type: "MACRS" # can be "MACRS" or "Straight line" - MACRS may be better and can reduce LCOH by more than $1/kg and is spec'd in the IRS MACRS schedule https://www.irs.gov/publications/p946#en_US_2020_publink1000107507 + depr_period: 5 # years - for clean energy facilities as specified by the IRS MACRS schedule https://www.irs.gov/publications/p946#en_US_2020_publink1000107507 + cost_adjustment_parameters: + cost_year_adjustment_inflation: 0.025 # used to adjust modeled costs to target_dollar_year + target_dollar_year: 2022 + finance_subgroups: + all_electricity: + commodity: "electricity" + finance_groups: ["profast_model"] + technologies: ["wind", "battery"] + + + # dispatched_electricity: + # commodity: "electricity" + # commodity_stream: "battery" #use only dispatched electricity from battery in finance calc + # technologies: ["wind", "battery"] diff --git a/examples/27_pyomo_optimized_dispatch/pyomo_optimized_dispatch.yaml b/examples/27_pyomo_optimized_dispatch/pyomo_optimized_dispatch.yaml new file mode 100644 index 000000000..b68a2d872 --- /dev/null +++ b/examples/27_pyomo_optimized_dispatch/pyomo_optimized_dispatch.yaml @@ -0,0 +1,7 @@ +name: "H2Integrate_config" + +system_summary: "This hybrid plant contains wind and battery storage technologies. The system is designed to dispatch storage optimally meet a specific electrical load." + +driver_config: "driver_config.yaml" +technology_config: "tech_config.yaml" +plant_config: "plant_config.yaml" diff --git a/examples/27_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py b/examples/27_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py new file mode 100644 index 000000000..dfe84bf1c --- /dev/null +++ b/examples/27_pyomo_optimized_dispatch/run_pyomo_optimized_dispatch.py @@ -0,0 +1,78 @@ +import numpy as np +from matplotlib import pyplot as plt + +from h2integrate.core.h2integrate_model import H2IntegrateModel + + +# Create an H2Integrate model +model = H2IntegrateModel("pyomo_optimized_dispatch.yaml") + +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") + +# Run the model +model.run() + +# Plot the results +fig, ax = plt.subplots(2, 1, sharex=True, figsize=(8, 6)) + +start_hour = 0 +end_hour = 200 + +ax[0].plot( + range(start_hour, end_hour), + model.prob.get_val("battery.SOC", units="percent")[start_hour:end_hour], + label="SOC", +) +ax[0].set_ylabel("SOC (%)") +ax[0].set_ylim([0, 110]) +ax[0].axhline(y=90.0, linestyle=":", color="k", alpha=0.5, label="Max Charge") +ax[0].legend() + +ax[1].plot( + range(start_hour, end_hour), + model.prob.get_val("battery.electricity_in", units="MW")[start_hour:end_hour], + linestyle="-", + label="Electricity In (MW)", +) +ax[1].plot( + range(start_hour, end_hour), + model.prob.get_val("battery.unused_electricity_out", units="MW")[start_hour:end_hour], + linestyle=":", + label="Unused Electricity (MW)", +) +ax[1].plot( + range(start_hour, end_hour), + model.prob.get_val("battery.unmet_electricity_demand_out", units="MW")[start_hour:end_hour], + linestyle=":", + label="Unmet Electrical Demand (MW)", +) +ax[1].plot( + range(start_hour, end_hour), + model.prob.get_val("battery.electricity_out", units="MW")[start_hour:end_hour], + linestyle="-", + label="Electricity Out (MW)", +) +ax[1].plot( + range(start_hour, end_hour), + model.prob.get_val("battery.battery_electricity_discharge", units="MW")[start_hour:end_hour], + linestyle="-.", + label="Battery Electricity Out (MW)", +) +ax[1].plot( + range(start_hour, end_hour), + demand_profile[start_hour:end_hour], + linestyle="--", + label="Electrical Demand (MW)", +) +ax[1].set_ylim([-1e2, 2.5e2]) +ax[1].set_ylabel("Electricity Hourly (MW)") +ax[1].set_xlabel("Timestep (hr)") + +plt.legend(ncol=2, frameon=False) +plt.tight_layout() +plt.savefig("optimized_dispatch_plot.png", dpi=300) diff --git a/examples/27_pyomo_optimized_dispatch/tech_config.yaml b/examples/27_pyomo_optimized_dispatch/tech_config.yaml new file mode 100644 index 000000000..80787e0bf --- /dev/null +++ b/examples/27_pyomo_optimized_dispatch/tech_config.yaml @@ -0,0 +1,80 @@ +name: "technology_config" +description: "This hybrid plant produces electricity from wind and battery storage." + + +technologies: + wind: + performance_model: + model: "PYSAMWindPlantPerformanceModel" + cost_model: + model: "ATBWindPlantCostModel" + dispatch_rule_set: + model: "PyomoDispatchGenericConverter" + resource: + type: "pysam_wind" + wind_speed: 9. + model_inputs: + performance_parameters: + num_turbines: 25 + turbine_rating_kw: 8300 + rotor_diameter: 196. + hub_height: 130. + create_model_from: "default" + config_name: "WindPowerSingleOwner" + pysam_options: !include pysam_options_8.3MW.yaml + run_recalculate_power_curve: False + layout: + layout_mode: "basicgrid" + layout_options: + row_D_spacing: 10.0 + turbine_D_spacing: 10.0 + rotation_angle_deg: 0.0 + row_phase_offset: 0.0 + layout_shape: "square" + cost_parameters: + capex_per_kW: 1500.0 + opex_per_kW_per_year: 45 + cost_year: 2019 + dispatch_rule_parameters: + commodity_name: "electricity" + commodity_storage_units: "kW" + battery: + dispatch_rule_set: + model: "PyomoRuleStorageBaseclass" + control_strategy: + model: "OptimizedDispatchController" + performance_model: + model: "PySAMBatteryPerformanceModel" + cost_model: + model: "ATBBatteryCostModel" + model_inputs: + shared_parameters: + commodity_name: "electricity" + max_charge_rate: 100000 + max_capacity: 400000 + init_charge_percent: 0.5 # Initial SOC for the storage + max_charge_percent: 0.9 # Maximum SOC allowable for the storage technology + min_charge_percent: 0.1 # Minimum SOC allowable for the storage technology + system_commodity_interface_limit: 1e12 + charge_efficiency: 0.95 # Charge efficiency of the storage technology + discharge_efficiency: 0.95 # Discharge efficiency of the storage technology + commodity_storage_units: "kW" + + performance_parameters: + system_model_source: "pysam" + chemistry: "LFPGraphite" + cost_parameters: + cost_year: 2022 + commodity_units: "kW" + energy_capex: 310 # $/kWh from 2024 ATB year 2025 + power_capex: 311 # $/kW from 2024 ATB year 2025 + opex_fraction: 0.25 # 0.25% of capex per year from 2024 ATB + control_parameters: + tech_name: "battery" + cost_per_charge: 0.03 # in $/kW, cost to charge the storage (note that charging is incentivized) + cost_per_discharge: 0.05 # in $/kW, cost to discharge the storage + commodity_met_value: 0.1 # in $/kW, penalty for not meeting the desired load demand + cost_per_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 + n_horizon_window: 48 # in timesteps (currently hours), The horizon window the optimization is run over 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 new file mode 100644 index 000000000..1dd66d987 --- /dev/null +++ b/h2integrate/control/control_rules/converters/generic_converter_min_operating_cost.py @@ -0,0 +1,296 @@ +import pyomo.environ as pyo +from pyomo.network import Port + + +class PyomoDispatchGenericConverterMinOperatingCosts: + """Class defining Pyomo rules for the optimized dispatch for load following + for generic commodity production components. + + Args: + commodity_info (dict): Dictionary of commodity information. This must contain the keys + "commodity_name" and "commodity_storage_units". + pyomo_model (pyo.ConcreteModel): Externally defined Pyomo model that works as the base + model that this class builds off of. + index_set (pyo.Set): Externally defined Pyomo index set for time steps. This should be + consistent with the forecast horizon of the optimization problem. + round_digits (int): Number of digits to round to in the Pyomo model. + block_set_name (str, optional): Name of the block set (model variables). + Defaults to "converter". + """ + + def __init__( + self, + commodity_info: dict, + pyomo_model: pyo.ConcreteModel, + index_set: pyo.Set, + round_digits: int, + time_duration: float, + block_set_name: str = "converter", + ): + # Set the number of digits to round to in the Pyomo model + self.round_digits = round_digits + # Set the block set name and commodity information + self.block_set_name = block_set_name + # Commodity information, this will be used to define variable and parameter + # names and units in the Pyomo model + self.commodity_name = commodity_info["commodity_name"] + self.commodity_storage_units = commodity_info["commodity_storage_units"] + + # 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 + # Set of time steps for the optimization problem, this will be used to define the Pyomo + # blocks for the dispatch model. This is where the internal variables, parameters, + # constraints, and ports are defined for the storage dispatch model in the + # dispatch_block_rule_function. + self.blocks = pyo.Block(index_set, rule=self.dispatch_block_rule_function) + + # Add the blocks to the Pyomo model with the specified block set name. + self.model.__setattr__(self.block_set_name, self.blocks) + # Set time steps for pyomo model. 1.0 here means that the time step is 1 hour. + # 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 + ): + """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 + + """ + self.cost_per_production = dispatch_inputs["cost_per_production"] + + def dispatch_block_rule_function(self, pyomo_model: pyo.ConcreteModel): + """ + Creates and initializes pyomo dispatch model components for a specific technology. + + This method sets up all model elements (parameters, variables, constraints, + and ports) associated with a technology block within the dispatch model. + It is typically called in the setup_pyomo() method of the PyomoControllerBaseClass. + + Args: + pyomo_model (pyo.ConcreteModel): The Pyomo model to which the technology + components will be added. + tech_name (str): The name or key identifying the technology (e.g., "battery", + "electrolyzer") for which model components are created. + """ + # Parameters + self._create_parameters(pyomo_model) + # Variables + self._create_variables(pyomo_model) + # Constraints + self._create_constraints(pyomo_model) + # Ports + self._create_ports(pyomo_model) + + # Base model setup + def _create_variables(self, pyomo_model: pyo.ConcreteModel): + """Create generic converter variables to add to Pyomo model instance. + + Args: + pyomo_model (pyo.ConcreteModel): pyomo_model the variables should be added to. + + """ + tech_var = pyo.Var( + doc=f"{self.commodity_name} production \ + from {self.block_set_name} [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + bounds=(0, pyomo_model.available_production), + units=eval("pyo.units." + self.commodity_storage_units), + initialize=0.0, + ) + + pyomo_model.__setattr__( + f"{self.block_set_name}_{self.commodity_name}", + tech_var, + ) + + def _create_ports(self, pyomo_model: pyo.ConcreteModel): + """Create generic converter ports to add to pyomo model instance. + + Args: + pyomo_model (pyo.ConcreteModel): pyomo_model the ports should be added to. + + """ + # create port + pyomo_model.port = Port() + # get port attribute from generic converter pyomo model + tech_port = pyomo_model.__getattribute__(f"{self.block_set_name}_{self.commodity_name}") + # add port to pyomo_model + pyomo_model.port.add(tech_port) + + def _create_parameters(self, pyomo_model: pyo.ConcreteModel): + """Create generic converter Pyomo parameters to add to the Pyomo model instance. + + This method defines converter parameters such as available production and the + cost per generation for the technology + + Args: + pyomo_model (pyo.ConcreteModel): pyomo_model that parameters are added to. + + """ + ################################## + # Parameters # + ################################## + pyomo_model.time_duration = pyo.Param( + doc=f"{pyomo_model.name} time step [hour]", + default=1.0, + within=pyo.NonNegativeReals, + mutable=True, + units=pyo.units.hr, + ) + pyomo_model.cost_per_production = pyo.Param( + doc=f"Production cost for generator [$/{self.commodity_storage_units}]", + default=0.0, + within=pyo.NonNegativeReals, + mutable=True, + units=eval(f"pyo.units.USD / pyo.units.{self.commodity_storage_units}h"), + ) + pyomo_model.available_production = pyo.Param( + doc=f"Available production for the generator [{self.commodity_storage_units}]", + default=0.0, + within=pyo.Reals, + mutable=True, + units=eval(f"pyo.units.{self.commodity_storage_units}"), + ) + + def _create_constraints(self, pyomo_model: pyo.ConcreteModel): + """Create generic converter Pyomo constraints to add to the Pyomo model instance. + + Method is currently empty but this serves as a placeholder to add constraints to the Pyomo + model instance if this class is inherited. + + Args: + pyomo_model (pyo.ConcreteModel): pyomo_model that constraints are added to. + + """ + + pass + + # Update time series parameters for next optimization window + def update_time_series_parameters( + self, commodity_in: list, commodity_demand: list, 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. + updated_initial_soc (float): The updated initial state of charge for storage + technologies for the current time slice. + """ + self.time_duration = [1.0] * len(self.blocks.index_set()) + self.available_production = [commodity_in[t] for t in self.blocks.index_set()] + + # Objective functions + def min_operating_cost_objective(self, hybrid_blocks, tech_name: str): + """Generic converter instance of minimum operating cost objective. + + Args: + hybrid_blocks (Pyomo.block): A generalized container for defining hierarchical + models by adding modeling components as attributes. + tech_name (str): The name or key identifying the technology for which + ports are created. + """ + + self.obj = sum( + hybrid_blocks[t].time_weighting_factor + * self.blocks[t].time_duration + * self.blocks[t].cost_per_production + * hybrid_blocks[t].__getattribute__(f"{tech_name}_{self.commodity_name}") + for t in hybrid_blocks.index_set() + ) + return self.obj + + # System-level functions + def _create_hybrid_port(self, hybrid_model: pyo.ConcreteModel, tech_name: str): + """Create generic converter ports to add to system-level pyomo model instance. + + Args: + hybrid_model (pyo.ConcreteModel): hybrid_model the ports should be added to. + tech_name (str): The name or key identifying the technology for which + ports are created. + """ + hybrid_model_tech = hybrid_model.__getattribute__(f"{tech_name}_{self.commodity_name}") + tech_port = Port(initialize={f"{tech_name}_{self.commodity_name}": hybrid_model_tech}) + hybrid_model.__setattr__(f"{tech_name}_port", tech_port) + + return hybrid_model.__getattribute__(f"{tech_name}_port") + + def _create_hybrid_variables(self, hybrid_model: pyo.ConcreteModel, tech_name: str): + """Create generic converter variables to add to system-level pyomo model instance. + + Args: + hybrid_model (pyo.ConcreteModel): hybrid_model the variables should be added to. + tech_name (str): The name or key identifying the technology for which + variables are created. + """ + tech_var = pyo.Var( + doc=f"{self.commodity_name} production \ + from {tech_name} [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=eval("pyo.units." + self.commodity_storage_units), + initialize=0.0, + ) + + hybrid_model.__setattr__(f"{tech_name}_{self.commodity_name}", tech_var) + + # Returns to power_source_gen_vars and load_vars in hybrid_rule + # load var is zero for converters + power_source_gen_var = hybrid_model.__getattribute__(f"{tech_name}_{self.commodity_name}") + load_var = 0 + return power_source_gen_var, load_var + + # Property getters and setters for time series parameters + @property + def available_production(self) -> list: + """Available production. + + Returns: + list: List of available production. + + """ + return [self.blocks[t].available_production.value for t in self.blocks.index_set()] + + @available_production.setter + def available_production(self, resource: list): + if len(resource) == len(self.blocks): + for t, gen in zip(self.blocks, resource): + self.blocks[t].available_production.set_value(round(gen, self.round_digits)) + else: + raise ValueError( + f"'resource' list ({len(resource)}) must be the same length as\ + time horizon ({len(self.blocks)})" + ) + + @property + def cost_per_production(self) -> float: + """Cost per generation [$/commodity_storage_units].""" + for t in self.blocks.index_set(): + return self.blocks[t].cost_per_production.value + + @cost_per_production.setter + def cost_per_production(self, om_dollar_per_kwh: float): + for t in self.blocks.index_set(): + self.blocks[t].cost_per_production.set_value( + round(om_dollar_per_kwh, self.round_digits) + ) + + @property + def time_duration(self) -> list: + """Time duration.""" + return [self.blocks[t].time_duration.value for t in self.blocks.index_set()] + + @time_duration.setter + def time_duration(self, time_duration: list): + if len(time_duration) == len(self.blocks): + for t, delta in zip(self.blocks, time_duration): + self.blocks[t].time_duration = round(delta, self.round_digits) + else: + raise ValueError( + self.time_duration.__name__ + " list must be the same length as time horizon" + ) diff --git a/h2integrate/control/control_rules/plant_dispatch_model.py b/h2integrate/control/control_rules/plant_dispatch_model.py new file mode 100644 index 000000000..60e009e81 --- /dev/null +++ b/h2integrate/control/control_rules/plant_dispatch_model.py @@ -0,0 +1,255 @@ +import pyomo.environ as pyo +from pyomo.network import Arc + + +class PyomoDispatchPlantModel: + """Class defining Pyomo model and rule for the optimized dispatch for load following + for the overall optimization problem describing the system. + + Args: + pyomo_model (pyo.ConcreteModel): Externally defined Pyomo model that works as the base + model that this class builds off of. + index_set (pyo.Set): Externally defined Pyomo index set for time steps. This should be + consistent with the forecast horizon of the optimization problem. + source_techs (list): List of technology names that are being dispatched in the system. + tech_dispatch_models (pyo.ConcreteModel): Externally defined Pyomo model that contains + the technology-specific dispatch rules and components for each technology in the system. + time_weighting_factor (float): Exponential time weighting factor for the + optimization problem that defines if/how future time steps are discounted relative to + the current time step in the optimization problem. + round_digits (int): Number of digits to round to in the Pyomo model. + block_set_name (str, optional): Name of the block set (model variables). + Defaults to "plant". + """ + + def __init__( + self, + pyomo_model: pyo.ConcreteModel, + index_set: pyo.Set, + source_techs: list, + tech_dispatch_models: pyo.ConcreteModel, + time_weighting_factor: float, + round_digits: int, + block_set_name: str = "plant", + ): + # Set the source technologies that are being dispatched in the system, this will be used to + # pull the generation and load variables from the technology-specific dispatch models. + # This includes all technologies in the dispatch model, e.g., converters and storage. + # These are pulled from self.pyomo_model in pyomo_controller.py + self.source_techs = source_techs + # Set the technology-specific dispatch models. These contain the technology-specific + # Pyomo models for each technology in the system. These are created in + # pyomo_controller.py and contain the technology-specific variables, parameters, + # constraints, and ports for each technology. + self.tech_dispatch_models = tech_dispatch_models + # Set the time weighting factor for the optimization problem, this defines if/how future + # time steps are discounted relative to the current time step in the optimization problem. + self.time_weighting_factor_input = time_weighting_factor + # Initialize lists to hold the generation and load variables from each technology + self.load_vars = {key: [] for key in index_set} + self.power_source_gen_vars = {key: [] for key in index_set} + # Initialize lists to hold the ports from each technology that will be connected to the + # system-level ports in the dispatch model. These are used to define the arcs between + # the technology-specific ports and the system-level ports. + self.ports = {key: [] for key in index_set} + # Initialize list to hold the arcs that connect the technology-specific ports to the + # system-level ports in the dispatch model. + self.arcs = [] + + # Set the number of digits to round to in the Pyomo model + self.round_digits = round_digits + + # The Pyomo model that this class builds off of, where all of the variables, parameters, + # constraints, and ports of the plant model will be added to. The plant model collects + # the generation and load variables from the technology-specific dispatch models and + # defines the system-level constraints that connect the technologies together in the + # overall system. + self.model = pyomo_model + # Set of time steps for the optimization problem, this will be used to define the Pyomo + # blocks for the dispatch model. This is where the internal variables, parameters, + # constraints, and ports are defined for the storage dispatch model in the + # dispatch_block_rule_function. + self.blocks = pyo.Block(index_set, rule=self.dispatch_block_rule) + # Add the blocks to the Pyomo model with the specified block set name. + self.model.__setattr__(block_set_name, self.blocks) + + def dispatch_block_rule(self, hybrid, t): + """ + Creates and initializes pyomo dispatch model components for a the system-level dispatch + + This method sets up all model elements (parameters, variables, constraints, + and ports) associated with a pyomo block within the dispatch model. + + Args: + hybrid (pyo.ConcreteModel): The Pyomo model to which the technology + components will be added. + t (int): integer location of variables in the control time window + """ + ################################## + # Parameters # + ################################## + self._create_parameters(hybrid) + ################################## + # Variables / Ports # + ################################## + self._create_variables_and_ports(hybrid, t) + ################################## + # Constraints # + ################################## + self._create_hybrid_constraints(hybrid, t) + + def initialize_parameters( + self, commodity_in: list, commodity_demand: list, 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 + + """ + 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) + + def _create_variables_and_ports(self, hybrid, t): + """Connect variables and ports from individual technology model + to system-level pyomo model instance. + + Args: + hybrid (pyo.ConcreteModel): The Pyomo model to which the technology + components will be added. + t (int): integer location of variables in the control time window + """ + + for tech in self.source_techs: + pyomo_block = self.tech_dispatch_models.__getattribute__(f"{tech}_rule") + gen_var, load_var = pyomo_block._create_hybrid_variables(hybrid, f"{tech}_rule") + + # Add production and load variables to system-level list + self.power_source_gen_vars[t].append(gen_var) + self.load_vars[t].append(load_var) + self.ports[t].append(pyomo_block._create_hybrid_port(hybrid, f"{tech}_rule")) + + @staticmethod + def _create_parameters(hybrid): + """Create system-level pyomo model parameters + + Args: + hybrid (pyo.ConcreteModel): The Pyomo model to which the technology + components will be added. + """ + hybrid.time_weighting_factor = pyo.Param( + doc="Exponential time weighting factor [-]", + initialize=1.0, + within=pyo.PercentFraction, + mutable=True, + units=pyo.units.dimensionless, + ) + + def _create_hybrid_constraints(self, hybrid, t): + """Define system-level constraints for pyomo model. + + Args: + hybrid (pyo.ConcreteModel): The Pyomo model to which the technology + components will be added. + t (int): integer location of variables in the control time window + """ + hybrid.production_total = pyo.Constraint( + doc="hybrid system generation total", + rule=hybrid.system_production == sum(self.power_source_gen_vars[t]), + ) + + hybrid.load_total = pyo.Constraint( + doc="hybrid system load total", + rule=hybrid.system_load == sum(self.load_vars[t]), + ) + + def create_arcs(self): + """ + Defines the mapping between individual technology variables to system level + """ + ################################## + # Arcs # + ################################## + for tech in self.source_techs: + pyomo_block = self.tech_dispatch_models.__getattribute__(f"{tech}_rule") + + def arc_rule(m, t): + source_port = pyomo_block.blocks[t].port + destination_port = self.blocks[t].__getattribute__(f"{tech}_rule_port") + return {"source": source_port, "destination": destination_port} + + tech_hybrid_arc = Arc(self.blocks.index_set(), rule=arc_rule) + self.model.__setattr__(f"{tech}_hybrid_arc", tech_hybrid_arc) + + tech_arc = self.model.__getattribute__(f"{tech}_hybrid_arc") + self.arcs.append(tech_arc) + + 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 + ): + """ + Updates the pyomo optimization problem with parameters that change with time + + Args: + commodity_in (list): List of generated commodity in for this time slice. + commodity_demand (list): The demanded commodity for this time slice. + updated_initial_soc (float): The updated initial state of charge for storage + technologies for the current time slice. + """ + # Note: currently, storage techs use commodity_demand and converter techs use commodity_in + # Better way to do this? + 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 + ) + + def create_min_operating_cost_expression(self): + """ + Creates system-level instance of minimum operating cost objective for pyomo solver. + """ + + self._delete_objective() + + def operating_cost_objective_rule(m) -> float: + obj = 0.0 + for tech in self.source_techs: + name = tech + "_rule" + # Create the min_operating_cost expression for each technology + pyomo_block = self.tech_dispatch_models.__getattribute__(name) + # Add to the overall hybrid operating cost expression + obj += pyomo_block.min_operating_cost_objective(self.blocks, name) + return obj + + # Set operating cost rule in Pyomo problem objective + self.model.objective = pyo.Objective(rule=operating_cost_objective_rule, sense=pyo.minimize) + + def _delete_objective(self): + if hasattr(self.model, "objective"): + self.model.del_component(self.model.objective) + + @property + def time_weighting_factor(self) -> float: + for t in self.blocks.index_set(): + return self.blocks[t + 1].time_weighting_factor.value + + @time_weighting_factor.setter + def time_weighting_factor(self, weighting: float): + for t in self.blocks.index_set(): + self.blocks[t].time_weighting_factor = round(weighting**t, self.round_digits) + + @property + def storage_commodity_out(self) -> list: + # This is used in the storage_dispatch_commands method of the control strategy + """Storage commodity out.""" + return [ + self.blocks[t].discharge_commodity.value - self.blocks[t].charge_commodity.value + for t in self.blocks.index_set() + ] diff --git a/h2integrate/control/control_rules/pyomo_rule_baseclass.py b/h2integrate/control/control_rules/pyomo_rule_baseclass.py index 4eed272e6..bb5bc621e 100644 --- a/h2integrate/control/control_rules/pyomo_rule_baseclass.py +++ b/h2integrate/control/control_rules/pyomo_rule_baseclass.py @@ -2,7 +2,7 @@ import pyomo.environ as pyo from attrs import field, define -from h2integrate.core.utilities import BaseConfig +from h2integrate.core.utilities import BaseConfig, merge_shared_inputs @define(kw_only=True) @@ -29,7 +29,8 @@ def initialize(self): def setup(self): self.config = PyomoRuleBaseConfig.from_dict( - self.options["tech_config"]["model_inputs"]["dispatch_rule_parameters"], + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "dispatch_rule"), + strict=False, additional_cls_name=self.__class__.__name__, ) diff --git a/h2integrate/control/control_rules/storage/pyomo_storage_rule_baseclass.py b/h2integrate/control/control_rules/storage/pyomo_storage_rule_baseclass.py index 29e2adbd0..053e6426c 100644 --- a/h2integrate/control/control_rules/storage/pyomo_storage_rule_baseclass.py +++ b/h2integrate/control/control_rules/storage/pyomo_storage_rule_baseclass.py @@ -5,7 +5,7 @@ class PyomoRuleStorageBaseclass(PyomoRuleBaseClass): - """Base class defining PYomo rules for generic commodity storage components.""" + """Base class defining Pyomo rules for generic commodity storage components.""" def _create_parameters(self, pyomo_model: pyo.ConcreteModel, t): """Create storage-related parameters in the Pyomo 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 new file mode 100644 index 000000000..55df5510e --- /dev/null +++ b/h2integrate/control/control_rules/storage/pyomo_storage_rule_min_operating_cost.py @@ -0,0 +1,645 @@ +import pyomo.environ as pyo +from pyomo.network import Port + + +class PyomoRuleStorageMinOperatingCosts: + """Class defining Pyomo rules for the optimized dispatch for load following + for generic commodity storage components. + + Args: + commodity_info (dict): Dictionary of commodity information. This must contain the keys + "commodity_name" and "commodity_storage_units". + pyomo_model (pyo.ConcreteModel): Externally defined Pyomo model that works as the base + model that this class builds off of. + index_set (pyo.Set): Externally defined Pyomo index set for time steps. This should be + consistent with the forecast horizon of the optimization problem. + round_digits (int): Number of digits to round to in the Pyomo model. + block_set_name (str, optional): Name of the block set (model variables). + Defaults to "storage". + """ + + def __init__( + self, + commodity_info: dict, + pyomo_model: pyo.ConcreteModel, + index_set: pyo.Set, + round_digits: int, + time_duration: float, + block_set_name: str = "storage", + ): + # Set the number of digits to round to in the Pyomo model + self.round_digits = round_digits + # Set the block set name and commodity information + self.block_set_name = block_set_name + # Commodity information, this will be used to define variable and parameter + # 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]"]) + + # 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 + # Set of time steps for the optimization problem, this will be used to define the Pyomo + # blocks for the dispatch model. This is where the internal variables, parameters, + # constraints, and ports are defined for the storage dispatch model in the + # dispatch_block_rule_function. + self.blocks = pyo.Block(index_set, rule=self.dispatch_block_rule_function) + + # Add the blocks to the Pyomo model with the specified block set name. + self.model.__setattr__(self.block_set_name, self.blocks) + # Set time steps for pyomo model. 1.0 here means that the time step is 1 hour. + # 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 + ): + """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 + + """ + # 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) + self.set_timeseries_parameter("maximum_storage", dispatch_inputs["max_capacity"]) + + self.set_timeseries_parameter("minimum_soc", dispatch_inputs["min_charge_percent"]) + self.set_timeseries_parameter("maximum_soc", dispatch_inputs["max_charge_percent"]) + + self.initial_soc = dispatch_inputs["initial_soc_percent"] + self.charge_efficiency = dispatch_inputs.get("charge_efficiency", 0.94) + self.discharge_efficiency = dispatch_inputs.get("discharge_efficiency", 0.94) + + # Set charge and discharge rate equal to each other for now + self.set_timeseries_parameter("max_charge", dispatch_inputs["max_charge_rate"]) + 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_initial_soc_constraint() + + def dispatch_block_rule_function(self, pyomo_model: pyo.ConcreteModel, tech_name: str): + """ + Creates and initializes pyomo dispatch model components for a specific technology. + + This method sets up all model elements (parameters, variables, constraints, + and ports) associated with a technology block within the dispatch model. + + Args: + pyomo_model (pyo.ConcreteModel): The Pyomo model to which the technology + components will be added. + tech_name (str): The name or key identifying the technology (e.g., "battery", + "electrolyzer") for which model components are created. + """ + # Parameters + self._create_parameters(pyomo_model, tech_name) + # Variables + self._create_variables(pyomo_model, tech_name) + # Constraints + self._create_constraints(pyomo_model, tech_name) + # Ports + self._create_ports(pyomo_model, tech_name) + + # Base model setup + def _create_parameters(self, pyomo_model: pyo.ConcreteModel, t): + """Create storage-related parameters in the Pyomo model. + + This method defines key storage parameters such as capacity limits, + state-of-charge (SOC) bounds, efficiencies, and time duration for each + time step. This method also defined system parameters such as the value of + load the load met and the production limit of the system. + + Args: + pyomo_model (pyo.ConcreteModel): Pyomo model instance representing + the storage system. + t: Time index or iterable representing time steps (unused in this method). + """ + ################################## + # Storage Parameters # + ################################## + + pyo_commodity_storage_unit = eval(f"pyo.units.{self.commodity_storage_units}") + pyo_commodity_storage_unit_hrs = eval(f"pyo.units.{self.commodity_storage_units}h") + pyo_usd_per_commodity_storage_unit_hrs = eval( + f"pyo.units.USD / pyo.units.{self.commodity_storage_units}h" + ) + usd_pr_units_str = f"[$/{self.commodity_storage_units}]" + + pyomo_model.time_duration = pyo.Param( + doc=f"{pyomo_model.name} time step [hour]", + default=1.0, + within=pyo.NonNegativeReals, + mutable=True, + units=pyo.units.hr, + ) + + pyomo_model.cost_per_charge = pyo.Param( + doc=f"Operating cost of {pyomo_model.name} charging {usd_pr_units_str}", + default=0.0, + within=pyo.NonNegativeReals, + mutable=True, + units=pyo_usd_per_commodity_storage_unit_hrs, + ) + pyomo_model.cost_per_discharge = pyo.Param( + doc=f"Operating cost of {pyomo_model.name} discharging {usd_pr_units_str}", + default=0.0, + within=pyo.NonNegativeReals, + mutable=True, + units=pyo_usd_per_commodity_storage_unit_hrs, + ) + pyomo_model.minimum_storage = pyo.Param( + doc=f"{pyomo_model.name} minimum storage rating [{self.commodity_storage_units}]", + default=0.0, + within=pyo.NonNegativeReals, + mutable=True, + units=pyo_commodity_storage_unit, + ) + pyomo_model.maximum_storage = pyo.Param( + doc=f"{pyomo_model.name} maximum storage rating [{self.commodity_storage_units}]", + default=1000.0, + within=pyo.NonNegativeReals, + mutable=True, + units=pyo_commodity_storage_unit_hrs, + ) + pyomo_model.minimum_soc = pyo.Param( + doc=f"{pyomo_model.name} minimum state-of-charge [-]", + default=0.1, + within=pyo.PercentFraction, + mutable=True, + units=pyo.units.dimensionless, + ) + pyomo_model.maximum_soc = pyo.Param( + doc=f"{pyomo_model.name} maximum state-of-charge [-]", + default=0.9, + within=pyo.PercentFraction, + mutable=True, + units=pyo.units.dimensionless, + ) + + ################################## + # Efficiency Parameters # + ################################## + pyomo_model.charge_efficiency = pyo.Param( + doc=f"{pyomo_model.name} Charging efficiency [-]", + default=0.938, + within=pyo.PercentFraction, + mutable=True, + units=pyo.units.dimensionless, + ) + pyomo_model.discharge_efficiency = pyo.Param( + doc=f"{pyomo_model.name} discharging efficiency [-]", + default=0.938, + within=pyo.PercentFraction, + mutable=True, + units=pyo.units.dimensionless, + ) + ################################## + # Capacity Parameters # + ################################## + pyomo_model.max_charge = pyo.Param( + doc=f"{pyomo_model.name} maximum charge [{self.commodity_storage_units}]", + within=pyo.NonNegativeReals, + mutable=True, + units=pyo_commodity_storage_unit, + ) + pyomo_model.max_discharge = pyo.Param( + doc=f"{pyomo_model.name} maximum discharge [{self.commodity_storage_units}]", + within=pyo.NonNegativeReals, + mutable=True, + units=pyo_commodity_storage_unit, + ) + ################################## + # System Parameters # + ################################## + pyomo_model.epsilon = pyo.Param( + doc="A small value used in objective for binary logic", + default=1e-3, + within=pyo.NonNegativeReals, + mutable=True, + units=pyo.units.USD, + ) + pyomo_model.commodity_met_value = pyo.Param( + doc=f"Commodity demand met value per generation [$/{self.commodity_storage_units}]", + default=0.0, + within=pyo.Reals, + mutable=True, + units=pyo_usd_per_commodity_storage_unit_hrs, + ) + pyomo_model.commodity_load_demand = pyo.Param( + doc=f"Load demand for the commodity [{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. + + This method defines binary and continuous variables representing + charging/discharging modes, energy flows, and state-of-charge, as well + as system variables such as system load, system production, and commodity produced. + + Args: + pyomo_model (pyo.ConcreteModel): Pyomo model instance representing + the storage system. + t: Time index or iterable representing time steps (unused in this method). + """ + ################################## + # Variables # + ################################## + pyomo_model.is_charging = pyo.Var( + doc=f"1 if {pyomo_model.name} is charging; 0 Otherwise [-]", + domain=pyo.Binary, + units=pyo.units.dimensionless, + ) + pyomo_model.is_discharging = pyo.Var( + doc=f"1 if {pyomo_model.name} is discharging; 0 Otherwise [-]", + domain=pyo.Binary, + units=pyo.units.dimensionless, + ) + pyomo_model.soc0 = pyo.Var( + doc=f"{pyomo_model.name} initial state-of-charge at beginning of period[-]", + domain=pyo.PercentFraction, + bounds=(pyomo_model.minimum_soc, pyomo_model.maximum_soc), + units=pyo.units.dimensionless, + ) + pyomo_model.soc = pyo.Var( + doc=f"{pyomo_model.name} state-of-charge at end of period [-]", + domain=pyo.PercentFraction, + bounds=(pyomo_model.minimum_soc, pyomo_model.maximum_soc), + units=pyo.units.dimensionless, + ) + + pyomo_model.charge_commodity = pyo.Var( + doc=f"{self.commodity_name} into {pyomo_model.name} [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=eval(f"pyo.units.{self.commodity_storage_units}"), + ) + pyomo_model.discharge_commodity = pyo.Var( + doc=f"{self.commodity_name} out of {pyomo_model.name} [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=eval(f"pyo.units.{self.commodity_storage_units}"), + ) + ################################## + # System Variables # + ################################## + pyomo_model.system_production = pyo.Var( + doc=f"System generation [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=eval(f"pyo.units.{self.commodity_storage_units}"), + ) + pyomo_model.system_load = pyo.Var( + doc=f"System load [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=eval(f"pyo.units.{self.commodity_storage_units}"), + ) + pyomo_model.commodity_out = pyo.Var( + doc=f"Commodity out of the system [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + bounds=(0, pyomo_model.commodity_load_demand), + units=eval(f"pyo.units.{self.commodity_storage_units}"), + ) + pyomo_model.is_generating = pyo.Var( + doc="System is producing commodity binary [-]", + domain=pyo.Binary, + units=pyo.units.dimensionless, + ) + + def _create_constraints(self, pyomo_model: pyo.ConcreteModel, t): + """Create operational and state-of-charge constraints for storage and the system. + + This method defines constraints that enforce: + - Mutual exclusivity between charging and discharging. + - Upper and lower bounds on charge/discharge flows. + - The state-of-charge balance over time. + - The system balance of output with system production and load + - The system output is less than or equal to the load (because of linear optimization) + + Args: + pyomo_model (pyo.ConcreteModel): Pyomo model instance representing + the storage system. + t: Time index or iterable representing time steps (unused in this method). + """ + ################################## + # Charging Constraints # + ################################## + # Charge commodity bounds + pyomo_model.charge_commodity_ub = pyo.Constraint( + doc=f"{pyomo_model.name} charging storage upper bound", + expr=pyomo_model.charge_commodity <= pyomo_model.max_charge * pyomo_model.is_charging, + ) + pyomo_model.charge_commodity_lb = pyo.Constraint( + doc=f"{pyomo_model.name} charging storage lower bound", + expr=pyomo_model.charge_commodity + >= pyomo_model.minimum_storage * pyomo_model.is_charging, + ) + # Discharge commodity bounds + pyomo_model.discharge_commodity_lb = pyo.Constraint( + doc=f"{pyomo_model.name} Discharging storage lower bound", + expr=pyomo_model.discharge_commodity + >= pyomo_model.minimum_storage * pyomo_model.is_discharging, + ) + pyomo_model.discharge_commodity_ub = pyo.Constraint( + doc=f"{pyomo_model.name} Discharging storage upper bound", + expr=pyomo_model.discharge_commodity + <= pyomo_model.max_discharge * pyomo_model.is_discharging, + ) + # Storage packing constraint + pyomo_model.charge_discharge_packing = pyo.Constraint( + doc=f"{pyomo_model.name} packing constraint for charging and discharging binaries", + expr=pyomo_model.is_charging + pyomo_model.is_discharging <= 1, + ) + ################################## + # System constraints # + ################################## + pyomo_model.balance = pyo.Constraint( + doc="Transmission energy 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 + <= pyomo_model.commodity_load_demand * pyomo_model.is_generating, + ) + + ################################## + # SOC Inventory Constraints # + ################################## + + def soc_inventory_rule(m): + return m.soc == ( + m.soc0 + + m.time_duration + * ( + m.charge_efficiency * m.charge_commodity + - (1 / m.discharge_efficiency) * m.discharge_commodity + ) + / m.maximum_storage + ) + + # Storage State-of-charge balance + pyomo_model.soc_inventory = pyo.Constraint( + doc=f"{pyomo_model.name} state-of-charge inventory balance", + rule=soc_inventory_rule, + ) + + def _set_initial_soc_constraint(self): + """ + This method links the SOC between the end of one control period and the beginning + of the next control period. + """ + ################################## + # SOC Linking # + ################################## + self.model.initial_soc = pyo.Param( + doc=f"{self.commodity_name} initial state-of-charge at beginning of the horizon[-]", + within=pyo.PercentFraction, + default=0.5, + mutable=True, + units=pyo.units.dimensionless, + ) + + ################################## + # SOC Constraints # + ################################## + # Linking time periods together + def storage_soc_linking_rule(m, t): + if t == self.blocks.index_set().first(): + return self.blocks[t].soc0 == m.initial_soc + return self.blocks[t].soc0 == self.blocks[t - 1].soc + + self.model.soc_linking = pyo.Constraint( + self.blocks.index_set(), + doc=self.block_set_name + " state-of-charge block linking constraint", + rule=storage_soc_linking_rule, + ) + + def _create_ports(self, pyomo_model: pyo.ConcreteModel, t): + """Create Pyomo ports for connecting the storage component. + + Ports are used to connect inflows and outflows of the storage system + (e.g., charging and discharging commodities) to the overall Pyomo model. + + Args: + pyomo_model (pyo.ConcreteModel): Pyomo model instance representing + the storage system. + t: Time index or iterable representing time steps (unused in this method). + """ + ################################## + # Ports # + ################################## + pyomo_model.port = Port() + pyomo_model.port.add(pyomo_model.charge_commodity) + pyomo_model.port.add(pyomo_model.discharge_commodity) + pyomo_model.port.add(pyomo_model.system_production) + pyomo_model.port.add(pyomo_model.system_load) + pyomo_model.port.add(pyomo_model.commodity_out) + + # Update time series parameters for next optimization window + def update_time_series_parameters( + self, commodity_in: list, commodity_demand: list, updated_initial_soc: None + ): + """Updates the pyomo optimization problem with parameters that change with time + + Args: + commodity_in (list): List of generated commodity in for this time slice. + commodity_demand (list): The demanded commodity for this time slice. + updated_initial_soc (float): The updated initial state of charge for storage + technologies for the current time slice. + """ + self.time_duration = [1.0] * len(self.blocks.index_set()) + self.commodity_load_demand = [commodity_demand[t] for t in self.blocks.index_set()] + self.model.initial_soc = updated_initial_soc + self.initial_soc = updated_initial_soc + + # Objective functions + def min_operating_cost_objective(self, hybrid_blocks, tech_name: str): + """Storage instance of minimum operating cost objective. + + Args: + hybrid_blocks (Pyomo.block): A generalized container for defining hierarchical + models by adding modeling components as attributes. + + """ + # Note that this objective function incentivizes charging the storage and penalizes + # discharging the storage. This is to help the storage model maintain a state of charge. + # This is also why cost_per_discharge should not equal cost_per_charge, which can lead + # to battery oscillation behavior. + self.obj = sum( + hybrid_blocks[t].time_weighting_factor + * self.blocks[t].time_duration + * ( + 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 + ) + for t in self.blocks.index_set() + ) + return self.obj + + # System-level functions + def _create_hybrid_port(self, hybrid_model: pyo.ConcreteModel, tech_name: str): + """Create generic storage ports to add to system-level pyomo model instance. + + Args: + hybrid_model (pyo.ConcreteModel): hybrid_model the ports should be added to. + 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, + } + ) + hybrid_model.__setattr__(f"{tech_name}_port", tech_port) + + return hybrid_model.__getattribute__(f"{tech_name}_port") + + def _create_hybrid_variables(self, hybrid_model: pyo.ConcreteModel, tech_name: str): + """Create generic storage variables to add to system-level pyomo model instance. + + Args: + hybrid_model (pyo.ConcreteModel): hybrid_model the variables should be added to. + tech_name (str): The name or key identifying the technology for which + variables are created. + """ + ################################## + # System Variables # + ################################## + pyo_commodity_units = eval("pyo.units." + self.commodity_storage_units) + + hybrid_model.system_production = pyo.Var( + doc=f"System generation [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=pyo_commodity_units, + ) + hybrid_model.system_load = pyo.Var( + doc=f"System load [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=pyo_commodity_units, + ) + hybrid_model.commodity_out = pyo.Var( + doc=f"{self.commodity_name} sold [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=pyo_commodity_units, + ) + ################################## + # Storage Variables # + ################################## + + hybrid_model.charge_commodity = pyo.Var( + doc=f"{self.commodity_name} into {tech_name} [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=pyo_commodity_units, + ) + hybrid_model.discharge_commodity = pyo.Var( + doc=f"{self.commodity_name} out of {tech_name} [{self.commodity_storage_units}]", + domain=pyo.NonNegativeReals, + units=pyo_commodity_units, + ) + # Returns to power_source_gen_vars and load_vars in hybrid_rule + power_source_gen_var = hybrid_model.discharge_commodity + load_var = hybrid_model.charge_commodity + return power_source_gen_var, load_var + + @staticmethod + def _check_efficiency_value(efficiency): + """Checks efficiency is between 0 and 1 or 0 and 100. Returns fractional value""" + if efficiency < 0: + raise ValueError("Efficiency value must greater than 0") + elif efficiency > 1: + efficiency /= 100 + if efficiency > 1: + raise ValueError("Efficiency value must between 0 and 1 or 0 and 100") + return efficiency + + # INPUTS + @property + def time_duration(self) -> list: + """Time duration.""" + return [self.blocks[t].time_duration.value for t in self.blocks.index_set()] + + @time_duration.setter + def time_duration(self, time_duration: list): + if len(time_duration) == len(self.blocks): + for t, delta in zip(self.blocks, time_duration): + self.blocks[t].time_duration = round(delta, self.round_digits) + else: + raise ValueError( + self.time_duration.__name__ + " list must be the same length as time horizon" + ) + + # Property getters and setters for time series parameters + + def set_timeseries_parameter(self, param_name: str, param_val: float): + for t in self.blocks.index_set(): + val_rounded = round(param_val, self.round_digits) + self.blocks[t].__setattr__(param_name, val_rounded) + + @property + def charge_efficiency(self) -> float: + """Charge efficiency.""" + for t in self.blocks.index_set(): + return self.blocks[t].charge_efficiency.value + + @charge_efficiency.setter + def charge_efficiency(self, efficiency: float): + efficiency = self._check_efficiency_value(efficiency) + for t in self.blocks.index_set(): + self.blocks[t].charge_efficiency = round(efficiency, self.round_digits) + + @property + def discharge_efficiency(self) -> float: + """Discharge efficiency.""" + for t in self.blocks.index_set(): + return self.blocks[t].discharge_efficiency.value + + @discharge_efficiency.setter + def discharge_efficiency(self, efficiency: float): + efficiency = self._check_efficiency_value(efficiency) + for t in self.blocks.index_set(): + self.blocks[t].discharge_efficiency = round(efficiency, self.round_digits) + + @property + def round_trip_efficiency(self) -> float: + """Round trip efficiency.""" + return self.charge_efficiency * self.discharge_efficiency + + @round_trip_efficiency.setter + def round_trip_efficiency(self, round_trip_efficiency: float): + round_trip_efficiency = self._check_efficiency_value(round_trip_efficiency) + # Assumes equal charge and discharge efficiencies + efficiency = round_trip_efficiency ** (1 / 2) + self.charge_efficiency = efficiency + self.discharge_efficiency = efficiency + + @property + def commodity_load_demand(self) -> list: + return [self.blocks[t].commodity_load_demand.value for t in self.blocks.index_set()] + + @commodity_load_demand.setter + def commodity_load_demand(self, commodity_demand: list): + if len(commodity_demand) == len(self.blocks): + for t, limit in zip(self.blocks, commodity_demand): + self.blocks[t].commodity_load_demand.set_value(round(limit, self.round_digits)) + else: + raise ValueError("'commodity_demand' list must be the same length as time horizon") diff --git a/h2integrate/control/control_strategies/controller_opt_problem_state.py b/h2integrate/control/control_strategies/controller_opt_problem_state.py new file mode 100644 index 000000000..a69e78138 --- /dev/null +++ b/h2integrate/control/control_strategies/controller_opt_problem_state.py @@ -0,0 +1,101 @@ +from pyomo.opt import TerminationCondition + + +class DispatchProblemState: + """Class for tracking dispatch problem solve state and metrics""" + + def __init__(self): + self._start_time = () + self._n_days = () + self._termination_condition = () + self._solve_time = () + self._objective = () + self._upper_bound = () + self._lower_bound = () + self._constraints = () + self._variables = () + self._non_zeros = () + self._gap = () + self._n_non_optimal_solves = 0 + + def store_problem_metrics(self, solver_results, start_time, n_days, objective_value): + """ + This method takes the solver results and formats them for debugging. + The outputs of this method are not actively used in the H2I simulation, but they are useful + for debugging and tracking solver performance over time. + + NOTE: that this method was brought over from HOPP. The link to the original method is here: + https://github.com/NatLabRockies/HOPP/blob/dde63faf6ea804828b2a7054cd6ec2c0a2f19614/hopp/simulation/technologies/dispatch/dispatch_problem_state.py + + + Args: + solver_results: The results object returned by the optimization solver. + start_time: The starting time of the optimization problem. + n_days: The number of days in the optimization horizon. + objective_value: The value of the objective function from the optimization problem. + """ + self.value("start_time", start_time) + self.value("n_days", n_days) + + # Unpack the solver results dictionary + solver_results_dict = { + k.lower().replace(" ", "_"): v.value + for k, v in solver_results.solver._list[0].items() + if k != "Statistics" + } + # Reformat variables names to take out spaces for programmatic simplicity + solver_problem_dict = { + k.lower().replace(" ", "_"): v.value for k, v in solver_results.problem._list[0].items() + } + # Simplify names of quantities of interest + prob_to_attr_map = { + "number_of_nonzeros": "non_zeros", + "number_of_variables": "variables", + "number_of_constraints": "constraints", + "lower_bound": "lower_bound", + "upper_bound": "upper_bound", + } + + # Save termination condition, solve time, and objective values + self.termination_condition = str(solver_results_dict["termination_condition"]) + if "time" in solver_results_dict: + # if optimally solved, the results will have a time variable + self.value("solve_time", solver_results_dict["time"]) + else: + # if the solve timed-out (not optimal), the results will have a wallclock_time variable + self.value("solve_time", solver_results_dict["wallclock_time"]) + + self.value("objective", objective_value) + + # Map values into this class structure + for solver_prob_key, attribute_name in prob_to_attr_map.items(): + self.value(attribute_name, solver_problem_dict[solver_prob_key]) + + # solver_results.solution.Gap not defined + upper_bound = solver_problem_dict["upper_bound"] + lower_bound = solver_problem_dict["lower_bound"] + if upper_bound != 0.0: + # if the upper bound is not equal to zero, calculate the gap between the lower and upper + # bounds (i.e. the feasible space for the solution) + gap = abs(upper_bound - lower_bound) / abs(upper_bound) + elif lower_bound == 0.0: + # If upper bound = 0 from the previous if, and the lower bound also equals 0, then the + # gap is zero + gap = 0.0 + else: + # Otherwise, the upper bound is infinite, and thus so is the solution space + gap = float("inf") + self.value("gap", gap) + + # This keeps track of the number of non-optimal solves (i.e. satisfactory) + if not solver_results_dict["termination_condition"] == TerminationCondition.optimal: + self._n_non_optimal_solves += 1 + + def value(self, metric_name: str, set_value=None): + if set_value is not None: + data = list(self.__getattribute__(f"_{metric_name}")) + data.append(set_value) + self.__setattr__(f"_{metric_name}", tuple(data)) + + else: + return self.__getattribute__(f"_{metric_name}") diff --git a/h2integrate/control/control_strategies/pyomo_controllers.py b/h2integrate/control/control_strategies/pyomo_controllers.py index 21ac733b0..632774636 100644 --- a/h2integrate/control/control_strategies/pyomo_controllers.py +++ b/h2integrate/control/control_strategies/pyomo_controllers.py @@ -3,10 +3,19 @@ import numpy as np import pyomo.environ as pyomo from attrs import field, define +from pyomo.util.check_units import assert_units_consistent from h2integrate.core.utilities import BaseConfig, merge_shared_inputs from h2integrate.core.validators import range_val +from h2integrate.control.control_rules.plant_dispatch_model import PyomoDispatchPlantModel from h2integrate.control.control_strategies.controller_baseclass import ControllerBaseClass +from h2integrate.control.control_strategies.controller_opt_problem_state import DispatchProblemState +from h2integrate.control.control_rules.storage.pyomo_storage_rule_min_operating_cost import ( + PyomoRuleStorageMinOperatingCosts, +) +from h2integrate.control.control_rules.converters.generic_converter_min_operating_cost import ( + PyomoDispatchGenericConverterMinOperatingCosts, +) if TYPE_CHECKING: # to avoid circular imports @@ -108,6 +117,9 @@ def setup(self): # get technology group name self.tech_group_name = self.pathname.split(".") + # initialize dispatch inputs to None + self.dispatch_options = None + # 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: @@ -142,7 +154,7 @@ def pyomo_setup(self, discrete_inputs): Returns: callable: Function(performance_model, performance_model_kwargs, inputs, commodity_name) - executing rolling-window heuristic dispatch and returning: + executing rolling-window heuristic dispatch or optimization and returning: (total_out, storage_out, unmet_demand, unused_commodity, soc) """ # initialize the pyomo model @@ -150,6 +162,9 @@ def pyomo_setup(self, discrete_inputs): 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 @@ -160,6 +175,7 @@ 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}" @@ -167,6 +183,7 @@ def pyomo_setup(self, discrete_inputs): # 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 @@ -175,13 +192,14 @@ def pyomo_dispatch_solver( performance_model: callable, performance_model_kwargs, inputs, + pyomo_model=self.pyomo_model, commodity_name: str = self.config.commodity_name, ): - r""" + """ 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 + `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. @@ -224,7 +242,6 @@ def pyomo_dispatch_solver( Notes: 1. Arrays returned have length self.n_timesteps (full simulation period). """ - self.initialize_parameters() # initialize outputs unmet_demand = np.zeros(self.n_timesteps) @@ -238,16 +255,35 @@ def pyomo_dispatch_solver( 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() + 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: - self.update_time_series_parameters() # get the inputs over the current control window commodity_in = inputs[self.config.commodity_name + "_in"][ t : t + self.config.n_control_window ] demand_in = inputs[f"{commodity_name}_demand"][t : t + self.config.n_control_window] - if control_strategy == "HeuristicLoadFollowingController": + 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( @@ -256,6 +292,23 @@ def pyomo_dispatch_solver( 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( @@ -271,6 +324,8 @@ def pyomo_dispatch_solver( **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)) @@ -316,24 +371,21 @@ def _check_efficiency_value(efficiency): def blocks(self) -> pyomo.Block: return getattr(self.pyomo_model, self.config.tech_name) - @property - def model(self) -> pyomo.ConcreteModel: - return self._model - class SimpleBatteryControllerHeuristic(PyomoControllerBaseClass): - """Fixes battery dispatch operations based on user input. + """Fixes storage dispatch operations based on user input. - Currently, enforces available generation and grid limit assuming no battery charging from grid. + Currently, enforces available generation and system interface limit assuming no + storage charging from external sources. Enforces: - Available generation cannot be exceeded for charging. - - Interface (grid / export) limit bounds discharge. - - No grid charging (unless logic extended elsewhere). + - Interface (e.g., grid / export) limit bounds discharge. + - No external charging (unless logic extended elsewhere). """ def setup(self): - """Initialize SimpleBatteryControllerHeuristic.""" + """Initialize the heuristic storage controller.""" super().setup() self.round_digits = 4 @@ -344,7 +396,6 @@ def setup(self): def initialize_parameters(self): """Initializes parameters.""" - self.minimum_storage = 0.0 self.maximum_storage = self.config.max_capacity self.minimum_soc = self.config.min_charge_percent @@ -449,13 +500,15 @@ def enforce_power_fraction_simple_bounds( minimum_soc: float, maximum_soc: float, ) -> float: - """Enforces simple bounds (0, .9) for battery power fractions. + """Enforces simple bounds for storage power fractions. Args: storage_fraction (float): Storage fraction from heuristic method. + minimum_soc (float): Minimum state of charge fraction. + maximum_soc (float): Maximum state of charge fraction. Returns: - storage_fraction (float): Bounded storage fraction. + float: Bounded storage fraction within [minimum_soc, maximum_soc]. """ if storage_fraction > maximum_soc: @@ -586,7 +639,7 @@ def user_fixed_dispatch(self, fixed_dispatch: list): def storage_dispatch_commands(self) -> list: """ Commanded dispatch including available commodity at current time step that has not - been used to charge the battery. + been used to charge storage. """ return [ (self.blocks[t].discharge_commodity.value - self.blocks[t].charge_commodity.value) @@ -640,6 +693,7 @@ def maximum_soc(self, maximum_soc: float): for t in self.blocks.index_set(): self.blocks[t].maximum_soc = round(maximum_soc, self.round_digits) + # Need these properties to define these values for methods in this class @property def charge_efficiency(self) -> float: """Charge efficiency.""" @@ -664,38 +718,25 @@ def discharge_efficiency(self, efficiency: float): for t in self.blocks.index_set(): self.blocks[t].discharge_efficiency = round(efficiency, self.round_digits) - @property - def round_trip_efficiency(self) -> float: - """Round trip efficiency.""" - return self.charge_efficiency * self.discharge_efficiency - - @round_trip_efficiency.setter - def round_trip_efficiency(self, round_trip_efficiency: float): - round_trip_efficiency = self._check_efficiency_value(round_trip_efficiency) - # Assumes equal charge and discharge efficiencies - efficiency = round_trip_efficiency ** (1 / 2) - self.charge_efficiency = efficiency - self.discharge_efficiency = efficiency - @define(kw_only=True) class HeuristicLoadFollowingControllerConfig(PyomoControllerBaseConfig): max_charge_rate: int | float = field() charge_efficiency: float = field(default=None) discharge_efficiency: float = field(default=None) - include_lifecycle_count: bool = field(default=False) class HeuristicLoadFollowingController(SimpleBatteryControllerHeuristic): - """Operates the battery based on heuristic rules to meet the demand profile based power - available from power generation profiles and power demand profile. + """Operates storage based on heuristic rules to meet the demand profile based on + available commodity from generation profiles and demand profile. - Currently, enforces available generation and grid limit assuming no battery charging from grid. + Currently, enforces available generation and system interface limit assuming no + storage charging from external sources. """ def setup(self): - """Initialize HeuristicLoadFollowingController.""" + """Initialize the heuristic load-following controller.""" self.config = HeuristicLoadFollowingControllerConfig.from_dict( merge_shared_inputs(self.options["tech_config"]["model_inputs"], "control"), additional_cls_name=self.__class__.__name__, @@ -716,7 +757,7 @@ def set_fixed_dispatch( system_commodity_interface_limit: list, commodity_demand: list, ): - """Sets charge and discharge power of battery dispatch using fixed_dispatch attribute + """Sets charge and discharge amount of storage dispatch using fixed_dispatch attribute and enforces available generation and charge/discharge limits. Args: @@ -750,3 +791,293 @@ def _heuristic_method(self, commodity_in, commodity_demand): if -fd > self.max_charge_fraction[t]: fd = -self.max_charge_fraction[t] self._fixed_dispatch[t] = fd + + +@define +class OptimizedDispatchControllerConfig(PyomoControllerBaseConfig): + """ + Configuration data container for Pyomo-based optimal dispatch. + + This class groups the parameters needed by the optimized dispatch controller. + Values are typically populated from the technology + `tech_config.yaml` (merged under the "control" section). + + Attributes: + max_charge_rate (float): + The maximum charge that the storage can accept + (in units of the commodity per time step). + charge_efficiency (float): + The efficiency of charging the storage (between 0 and 1). + discharge_efficiency (float): + The efficiency of discharging the storage (between 0 and 1). + commodity_name (str): + The name of the commodity being stored (e.g., "electricity", "hydrogen"). + commodity_storage_units (str): + The units of the commodity being stored (e.g., "kW", "kg"). + cost_per_production (float): + The cost to use the incoming produced commodity (in $/commodity_storage_units). + cost_per_charge (float): + The cost per unit of charging the storage (in $/commodity_storage_units). + cost_per_discharge (float): + The cost per unit of discharging the storage (in $/commodity_storage_units). + commodity_met_value (float): + The penalty for not meeting the desired load demand (in $/commodity_storage_units). + time_weighting_factor (float): + The weighting factor applied to future time steps in the optimization objective + (between 0 and 1). + round_digits (int): + The number of digits to round to in the Pyomo model for numerical stability. + The default of this parameter is 4. + 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_charge_rate: int | float = field() + charge_efficiency: float = field(default=None) + discharge_efficiency: float = field(default=None) + commodity_name: str = field(default=None) + commodity_storage_units: str = field(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: float = field(default=None) + time_weighting_factor: float = field(default=0.995) + round_digits: int = field(default=4) + time_duration: float = field(default=1.0) # hours + + def make_dispatch_inputs(self): + dispatch_keys = [ + "cost_per_production", + "cost_per_charge", + "cost_per_discharge", + "commodity_met_value", + "max_capacity", + "max_charge_percent", + "min_charge_percent", + "charge_efficiency", + "discharge_efficiency", + "max_charge_rate", + ] + + dispatch_inputs = {k: self.as_dict()[k] for k in dispatch_keys} + dispatch_inputs.update({"initial_soc_percent": self.init_charge_percent}) + return dispatch_inputs + + +class OptimizedDispatchController(PyomoControllerBaseClass): + """Operates storage based on optimization to meet the demand profile based on + available commodity from generation profiles and demand profile while minimizing costs. + + Uses a rolling-window optimization approach with configurable horizon and control windows. + + """ + + def setup(self): + """Initialize the optimized dispatch controller.""" + self.config = OptimizedDispatchControllerConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "control") + ) + + self.add_input( + "max_charge_rate", + val=self.config.max_charge_rate, + units=self.config.commodity_storage_units, + desc="Storage charge rate", + ) + + self.add_input( + "storage_capacity", + val=self.config.max_capacity, + units=f"{self.config.commodity_storage_units}*h", + desc="Storage capacity", + ) + + 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_charge_percent + + # Is this the best place to put this??? + self.commodity_info = { + "commodity_name": self.config.commodity_name, + "commodity_storage_units": self.config.commodity_storage_units, + } + # TODO: note that this definition of cost_per_production is not generalizable to multiple + # production technologies. Would need a name adjustment to connect it to + # production tech + + self.dispatch_inputs = self.config.make_dispatch_inputs() + + def initialize_parameters(self, commodity_in, commodity_demand): + """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. + + """ + # Where pyomo model communicates with the rest of the controller + # self.hybrid_dispatch_model is the pyomo model, this is the thing in hybrid_rule + self.hybrid_dispatch_model = self._create_dispatch_optimization_model() + self.hybrid_dispatch_rule.create_min_operating_cost_expression() + self.hybrid_dispatch_rule.create_arcs() + assert_units_consistent(self.hybrid_dispatch_model) + + # This calls a class that stores problem state information such as solver metrics and + # the objective function. This is directly used in the H2I simulation, but is + # useful for tracking solver performance and debugging. + self.problem_state = DispatchProblemState() + + # hybrid_dispatch_rule is the thing where you can access variables and hybrid_rule \ + # functions from + self.hybrid_dispatch_rule.initialize_parameters( + commodity_in, commodity_demand, self.dispatch_inputs + ) + + def update_time_series_parameters( + self, commodity_in=None, commodity_demand=None, updated_initial_soc=None + ): + """Updates the pyomo optimization problem with parameters that change with time + + Args: + commodity_in (list): List of generated commodity in for this time slice. + commodity_demand (list): The demanded commodity for this time slice. + updated_initial_soc (float): The updated initial state of charge for storage + technologies for the current time slice. + """ + self.hybrid_dispatch_rule.update_time_series_parameters( + commodity_in, commodity_demand, updated_initial_soc + ) + + def solve_dispatch_model( + self, + start_time: int = 0, + n_days: int = 0, + ): + """Solves the dispatch optimization model and stores problem metrics. + + Args: + start_time (int): Starting timestep index for the current solve window. + n_days (int): Total number of days in the simulation. + + """ + + solver_results = self.glpk_solve_call(self.hybrid_dispatch_model) + # The outputs of the store_problem_metrics method are not actively used in the H2I + # simulation, but they are useful for debugging and tracking solver performance over time. + self.problem_state.store_problem_metrics( + solver_results, start_time, n_days, pyomo.value(self.hybrid_dispatch_model.objective) + ) + + def _create_dispatch_optimization_model(self): + """ + Creates monolith dispatch model by creating pyomo models for each technology, then + aggregating them into hybrid_rule + """ + model = pyomo.ConcreteModel(name="hybrid_dispatch") + ################################# + # Sets # + ################################# + model.forecast_horizon = pyomo.Set( + doc="Set of time periods in time horizon", + initialize=range(self.config.n_control_window), + ) + for tech in self.source_techs: + if tech == self.dispatch_tech[0]: + dispatch = PyomoRuleStorageMinOperatingCosts( + self.commodity_info, + model, + model.forecast_horizon, + self.config.round_digits, + self.config.time_duration, + block_set_name=f"{tech}_rule", + ) + self.pyomo_model.__setattr__(f"{tech}_rule", dispatch) + else: + dispatch = PyomoDispatchGenericConverterMinOperatingCosts( + self.commodity_info, + model, + model.forecast_horizon, + self.config.round_digits, + self.config.time_duration, + block_set_name=f"{tech}_rule", + ) + self.pyomo_model.__setattr__(f"{tech}_rule", dispatch) + + # Create hybrid pyomo model, inputting individual technology models + self.hybrid_dispatch_rule = PyomoDispatchPlantModel( + model, + model.forecast_horizon, + self.source_techs, + self.pyomo_model, + self.config.time_weighting_factor, + self.config.round_digits, + ) + return model + + def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): + """Build Pyomo model blocks and assign the dispatch solver.""" + self.dispatch_inputs["max_charge_rate"] = inputs["max_charge_rate"][0] + self.dispatch_inputs["max_capacity"] = inputs["storage_capacity"][0] + self.config.max_capacity = inputs["storage_capacity"][0] + self.config.max_charge_rate = inputs["max_charge_rate"][0] + + discrete_outputs["pyomo_dispatch_solver"] = self.pyomo_setup(discrete_inputs) + + @staticmethod + def glpk_solve_call( + pyomo_model: pyomo.ConcreteModel, + log_name: str = "", + user_solver_options: dict | None = None, + ): + """ + This method takes in the dispatch system-level pyomo model that we have built, + gives it to the solver, and gives back solver results. + """ + + # log_name = "annual_solve_GLPK.log" # For debugging MILP solver + # Ref. on solver options: https://en.wikibooks.org/wiki/GLPK/Using_GLPSOL + glpk_solver_options = { + "cuts": None, + "presol": None, + # 'mostf': None, + # 'mipgap': 0.001, + "tmlim": 30, + } + solver_options = SolverOptions(glpk_solver_options, log_name, user_solver_options, "log") + with pyomo.SolverFactory("glpk") as solver: + results = solver.solve(pyomo_model, options=solver_options.constructed) + + return results + + @property + def storage_dispatch_commands(self) -> list: + """ + Commanded dispatch including available commodity at current time step that has not + been used to charge storage. + """ + return self.hybrid_dispatch_rule.storage_commodity_out + + +class SolverOptions: + """Class for housing solver options""" + + def __init__( + self, + solver_spec_options: dict, + log_name: str = "", + user_solver_options: dict | None = None, + solver_spec_log_key: str = "logfile", + ): + self.instance_log = "dispatch_solver.log" + self.solver_spec_options = solver_spec_options + self.user_solver_options = user_solver_options + + self.constructed = solver_spec_options + if log_name != "": + self.constructed[solver_spec_log_key] = self.instance_log + if user_solver_options is not None: + self.constructed.update(user_solver_options) diff --git a/h2integrate/control/test/test_pyomo_controllers.py b/h2integrate/control/test/test_heuristic_controllers.py similarity index 100% rename from h2integrate/control/test/test_pyomo_controllers.py rename to h2integrate/control/test/test_heuristic_controllers.py diff --git a/h2integrate/control/test/test_optimal_controllers.py b/h2integrate/control/test/test_optimal_controllers.py new file mode 100644 index 000000000..9a3f33b43 --- /dev/null +++ b/h2integrate/control/test/test_optimal_controllers.py @@ -0,0 +1,348 @@ +import numpy as np +import pytest + +from h2integrate.core.h2integrate_model import H2IntegrateModel + + +plant_config = { + "name": "plant_config", + "description": "...", + "sites": { + "site": {"latitude": 35.2018863, "longitude": -101.945027}, + }, + "plant": { + "plant_life": 1, + "grid_connection": False, + "ppa_price": 0.025, + "hybrid_electricity_estimated_cf": 0.492, + "simulation": { + "dt": 3600, + "n_timesteps": 8760, + }, + }, + "technology_interconnections": [["combiner", "battery", "electricity", "cable"]], + "tech_to_dispatch_connections": [ + ["combiner", "battery"], + ["battery", "battery"], + ], +} + +driver_config = { + "name": "driver_config", + "description": "Pyomo optimal min operating cost test", + "general": {}, +} + +tech_config = { + "name": "technology_config", + "description": "...", + "technologies": { + "battery": { + "dispatch_rule_set": {"model": "PyomoRuleStorageBaseclass"}, + "control_strategy": {"model": "OptimizedDispatchController"}, + "performance_model": {"model": "PySAMBatteryPerformanceModel"}, + "model_inputs": { + "shared_parameters": { + "max_charge_rate": 50000, + "max_capacity": 200000, + "n_control_window": 24, + "n_horizon_window": 48, + "init_charge_percent": 0.5, + "max_charge_percent": 0.9, + "min_charge_percent": 0.1, + "commodity_name": "electricity", + "commodity_storage_units": "kW", + "time_weighting_factor": 0.995, + "charge_efficiency": 0.95, + "discharge_efficiency": 0.95, + "cost_per_charge": 0.004, + "cost_per_discharge": 0.005, + "cost_per_production": 0.0, + "commodity_met_value": 0.1, + "round_digits": 4, + }, + "performance_parameters": { + "system_model_source": "pysam", + "chemistry": "LFPGraphite", + "control_variable": "input_power", + }, + "control_parameters": { + "tech_name": "battery", + "system_commodity_interface_limit": 1e12, + }, + }, + }, + "combiner": { + "performance_model": {"model": "GenericCombinerPerformanceModel"}, + "dispatch_rule_set": {"model": "PyomoDispatchGenericConverter"}, + "model_inputs": { + "performance_parameters": { + "commodity": "electricity", + "commodity_units": "kW", + "in_streams": 1, + }, + "dispatch_rule_parameters": { + "commodity_name": "electricity", + "commodity_storage_units": "kW", + }, + }, + }, + }, +} + + +def test_min_operating_cost_load_following_battery_dispatch(subtests): + # Fabricate some oscillating power generation data: 1000 kW for the first 12 hours, 10000 kW for + # the second twelve hours, and repeat that daily cycle over a year. + n_look_ahead_half = int(24 / 2) + + electricity_in = np.concatenate( + (np.ones(n_look_ahead_half) * 1000, np.ones(n_look_ahead_half) * 10000) + ) + electricity_in = np.tile(electricity_in, 365) + + demand_in = np.ones(8760) * 6000.0 + + # Create an H2Integrate model + model = H2IntegrateModel( + { + "driver_config": driver_config, + "technology_config": tech_config, + "plant_config": plant_config, + } + ) + + # Setup the system and required values + model.setup() + model.prob.set_val("combiner.electricity_in1", electricity_in) + model.prob.set_val("battery.electricity_demand", demand_in) + + # Run the model + model.prob.run_model() + + # Test the case where the charging/discharging cycle remains within the max and min SOC limits + # Check the expected outputs to actual outputs + expected_electricity_out = [ + 5999.99997732, + 5992.25494845, + 5991.96052468, + 5991.63342842, + 5991.26824325, + 5990.86174194, + 5990.40961477, + 5989.90607785, + 5989.34362595, + 5988.71271658, + 5988.00134229, + 5987.19448473, + 6000.0, + 6000.0, + 6000.0, + 6000.0, + 6000.0, + 6000.0, + 6000.0, + 6000.0, + 6000.0, + 6000.0, + 6000.0, + 6000.0, + ] + + expected_battery_electricity_discharge = [ + 4999.99997732, + 4992.25494845, + 4991.96052468, + 4991.63342842, + 4991.26824325, + 4990.86174194, + 4990.40961477, + 4989.90607785, + 4989.34362595, + 4988.71271658, + 4988.00134229, + 4987.19448473, + -3990.28117686, + -3990.74350731, + -3991.15657455, + -3991.53821244, + -3991.89075932, + -3992.21669822, + -3992.51846124, + -3992.79833559, + -3993.05841677, + -3993.30060036, + -3993.52658465, + -3993.73788536, + ] + + expected_battery_soc = [ + 49.87479765, + 47.50390223, + 45.12932011, + 42.75108798, + 40.3682277, + 37.9797332, + 35.58459541, + 33.18174418, + 30.76997461, + 28.34786178, + 25.91365422, + 23.4651282, + 25.4168656, + 27.36180226, + 29.29938411, + 31.23051435, + 33.15594392, + 35.07630244, + 36.99212221, + 38.90385706, + 40.81189705, + 42.71658006, + 44.61820098, + 46.517019, + 44.14026872, + 41.75708657, + 39.3692475, + 36.97562686, + 34.57512794, + 32.16658699, + 29.7486839, + 27.31984419, + 24.87811568, + 22.42099596, + 19.94517118, + 17.44609219, + ] + + expected_unmet_demand = np.array( + [ + 2.26821512e-05, + 7.74505155, + 8.03947532, + 8.36657158, + 8.73175675, + 9.13825806, + 9.59038523, + 1.00939222e01, + 1.06563740e01, + 1.12872834e01, + 1.19986577e01, + 1.28055153e01, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + ] + ) + + expected_unused_electricity = np.array( + [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 9.71882314, + 9.25649269, + 8.84342545, + 8.46178756, + 8.10924068, + 7.78330178, + 7.48153876, + 7.20166441, + 6.94158323, + 6.69939964, + 6.47341535, + 6.26211464, + ] + ) + + with subtests.test("Check battery.electricity_out"): + assert ( + pytest.approx(expected_electricity_out) + == model.prob.get_val("battery.electricity_out")[0:24] + ) + + with subtests.test("Check battery_electricity_discharge"): + assert ( + pytest.approx(expected_battery_electricity_discharge) + == model.prob.get_val("battery.battery_electricity_discharge")[0:24] + ) + + # Check a longer portion of SOC to make sure SOC is getting linked between optimization periods + with subtests.test("Check SOC"): + assert pytest.approx(expected_battery_soc) == model.prob.get_val("battery.SOC")[0:36] + + with subtests.test("Check unmet_demand"): + assert ( + pytest.approx(expected_unmet_demand, abs=1e-4) + == model.prob.get_val("battery.unmet_electricity_demand_out")[0:24] + ) + + with subtests.test("Check unused_electricity_out"): + assert ( + pytest.approx(expected_unused_electricity) + == model.prob.get_val("battery.unused_electricity_out")[0:24] + ) + + # Test the case where the battery efficiency is lower + tech_config["technologies"]["battery"]["model_inputs"]["shared_parameters"][ + "charge_efficiency" + ] = 0.632 + tech_config["technologies"]["battery"]["model_inputs"]["shared_parameters"][ + "discharge_efficiency" + ] = 0.632 + + model = H2IntegrateModel( + { + "driver_config": driver_config, + "technology_config": tech_config, + "plant_config": plant_config, + } + ) + + # Setup the system and required values + model.setup() + model.prob.set_val("combiner.electricity_in1", electricity_in) + model.prob.set_val("battery.electricity_demand", demand_in) + + # Run the model + model.prob.run_model() + + expected_electricity_out = [ + 5999.99997732, + 5992.25494845, + 5991.96052468, + 5991.63342842, + 5991.26824325, + 5990.86174194, + 5990.40961477, + 5989.90607785, + 5989.34362595, + 5988.71271658, + 1558.72773849, + 1000.0, + ] + + # Make sure output changes if efficiency is changed + with subtests.test("Check electricity_out for different efficiency"): + assert ( + pytest.approx(expected_electricity_out) + == model.prob.get_val("battery.electricity_out")[:12] + ) diff --git a/h2integrate/core/supported_models.py b/h2integrate/core/supported_models.py index 060289fb9..d42433b56 100644 --- a/h2integrate/core/supported_models.py +++ b/h2integrate/core/supported_models.py @@ -92,6 +92,7 @@ from h2integrate.converters.hydrogen.singlitico_cost_model import SingliticoCostModel from h2integrate.converters.co2.marine.direct_ocean_capture import DOCCostModel, DOCPerformanceModel from h2integrate.control.control_strategies.pyomo_controllers import ( + OptimizedDispatchController, HeuristicLoadFollowingController, ) from h2integrate.converters.hydrogen.geologic.mathur_modified import GeoH2SubsurfaceCostModel @@ -147,6 +148,12 @@ from h2integrate.control.control_strategies.converters.demand_openloop_controller import ( DemandOpenLoopConverterController, ) +from h2integrate.control.control_rules.storage.pyomo_storage_rule_min_operating_cost import ( + PyomoRuleStorageMinOperatingCosts, +) +from h2integrate.control.control_rules.converters.generic_converter_min_operating_cost import ( + PyomoDispatchGenericConverterMinOperatingCosts, +) from h2integrate.control.control_strategies.converters.flexible_demand_openloop_controller import ( FlexibleDemandOpenLoopConverterController, ) @@ -249,11 +256,14 @@ "PassThroughOpenLoopController": PassThroughOpenLoopController, "DemandOpenLoopStorageController": DemandOpenLoopStorageController, "HeuristicLoadFollowingController": HeuristicLoadFollowingController, + "OptimizedDispatchController": OptimizedDispatchController, "DemandOpenLoopConverterController": DemandOpenLoopConverterController, "FlexibleDemandOpenLoopConverterController": FlexibleDemandOpenLoopConverterController, # Dispatch "PyomoDispatchGenericConverter": PyomoDispatchGenericConverter, "PyomoRuleStorageBaseclass": PyomoRuleStorageBaseclass, + "PyomoRuleStorageMinOperatingCosts": PyomoRuleStorageMinOperatingCosts, + "PyomoDispatchGenericConverterMinOperatingCosts": PyomoDispatchGenericConverterMinOperatingCosts, # noqa: E501 # Feedstock "FeedstockPerformanceModel": FeedstockPerformanceModel, "FeedstockCostModel": FeedstockCostModel,