From 8aed31975792a2d4c43301caf4a3327dff444052 Mon Sep 17 00:00:00 2001 From: joda9 Date: Thu, 27 Nov 2025 09:19:13 +0100 Subject: [PATCH] =?UTF-8?q?feat:=20Implement=20=C2=A714a=20EnWG=20curtailm?= =?UTF-8?q?ent=20constraints=20for=20charging=20points=20and=20heat=20pump?= =?UTF-8?q?s.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Added support for §14a EnWG curtailment settings in the to_powermodels function. - Introduced new variables and constraints in the OPF to handle curtailment logic. - Enhanced logging for curtailment statistics in the from_powermodels function. - Updated the objective function to minimize costs associated with §14a curtailment. - Ensured all time-varying data in the PyPSA network is of float64 dtype to prevent errors. - Refactored various functions to improve readability and maintainability. - Added new Julia files for §14a constraints and integrated them into the OPF workflow. --- .../config/config_grid_expansion_default.cfg | 25 ++ edisgo/edisgo.py | 178 ++++++++ edisgo/flex_opt/__init__.py | 21 + edisgo/flex_opt/costs.py | 2 +- edisgo/flex_opt/curtailment_14a.py | 424 ++++++++++++++++++ edisgo/flex_opt/reinforce_grid.py | 35 ++ edisgo/io/electromobility_import.py | 6 +- edisgo/io/generators_import.py | 62 +-- edisgo/io/powermodels_io.py | 176 ++++++-- edisgo/io/pypsa_io.py | 12 + edisgo/network/timeseries.py | 26 +- edisgo/network/topology.py | 32 +- .../src/core/constraint_curtailment_14a.jl | 97 ++++ .../opf/eDisGo_OPF.jl/src/core/objective.jl | 20 + .../opf/eDisGo_OPF.jl/src/core/variables.jl | 31 ++ edisgo/opf/eDisGo_OPF.jl/src/eDisGo_OPF.jl | 1 + edisgo/opf/eDisGo_OPF.jl/src/prob/opf_bf.jl | 7 + edisgo/opf/powermodels_opf.py | 17 + edisgo/tools/tools.py | 1 + 19 files changed, 1068 insertions(+), 105 deletions(-) create mode 100644 edisgo/flex_opt/curtailment_14a.py create mode 100644 edisgo/opf/eDisGo_OPF.jl/src/core/constraint_curtailment_14a.jl diff --git a/edisgo/config/config_grid_expansion_default.cfg b/edisgo/config/config_grid_expansion_default.cfg index 74a612eb7..c4d511bde 100644 --- a/edisgo/config/config_grid_expansion_default.cfg +++ b/edisgo/config/config_grid_expansion_default.cfg @@ -96,6 +96,31 @@ lv_load_case_line = 1.0 lv_feed-in_case_transformer = 1.0 lv_feed-in_case_line = 1.0 +# §14a EnWG curtailment +# ===================== +# Settings for curtailment of controllable consumption devices according to §14a EnWG +# as an alternative to grid expansion + +[curtailment_14a_enwg] + +# enable_curtailment: whether to consider §14a curtailment as an option in grid +# reinforcement. If True, curtailment can be used instead of grid expansion +# to solve grid issues. +enable_curtailment = True + +# max_power_kw: maximum allowed power in kW after curtailment according to §14a EnWG. +# The legal default is 4.2 kW for controllable consumption devices. +max_power_kw = 4.2 + +# curtailment_priority: defines how to prioritize components for curtailment. +# Options: 'p_set' (largest nominal power first), 'random', 'grid_level' +# (start with LV, then MV) +curtailment_priority = p_set + +# components_type: which types of components can be curtailed. +# Options: 'heat_pump', 'charging_point', 'both' +components_type = both + # costs # ============ diff --git a/edisgo/edisgo.py b/edisgo/edisgo.py index c452dbc6e..77d949682 100755 --- a/edisgo/edisgo.py +++ b/edisgo/edisgo.py @@ -17,6 +17,11 @@ from edisgo.flex_opt.charging_strategies import charging_strategy from edisgo.flex_opt.check_tech_constraints import lines_relative_load +from edisgo.flex_opt.curtailment_14a import ( + apply_curtailment_14a, + check_curtailment_effect, + identify_components_for_curtailment, +) from edisgo.flex_opt.heat_pump_operation import ( operating_strategy as hp_operating_strategy, ) @@ -828,6 +833,7 @@ def to_powermodels( flexible_loads=None, flexible_storage_units=None, opf_version=1, + curtailment_14a=None, ): """ Convert eDisGo representation of the network topology and timeseries to @@ -854,6 +860,11 @@ def to_powermodels( Version of optimization models to choose from. Must be one of [1, 2, 3, 4]. For more information see :func:`edisgo.opf.powermodels_opf.pm_optimize`. Default: 1. + curtailment_14a : dict or None + Configuration for §14a EnWG curtailment of heat pumps and charging points. + Dictionary should contain keys 'apply_curtailment' (bool), 'max_power_mw' + (float), and 'components' (list of component names). + Default: None. Returns ------- @@ -870,6 +881,7 @@ def to_powermodels( flexible_loads=flexible_loads, flexible_storage_units=flexible_storage_units, opf_version=opf_version, + curtailment_14a=curtailment_14a, ) def pm_optimize( @@ -886,6 +898,7 @@ def pm_optimize( save_heat_storage=True, save_slack_gen=True, save_slacks=True, + curtailment_14a=None, ): """ Run OPF in julia subprocess and write results of OPF back to edisgo object. @@ -933,6 +946,11 @@ def pm_optimize( hence there will be no logging coming from julia subprocess in python process. Default: False. + curtailment_14a : dict or None + Configuration for §14a EnWG curtailment of heat pumps and charging points. + Dictionary should contain keys 'apply_curtailment' (bool), 'max_power_mw' + (float), and 'components' (list of component names). + Default: None. """ return powermodels_opf.pm_optimize( self, @@ -945,6 +963,7 @@ def pm_optimize( method=method, warm_start=warm_start, silence_moi=silence_moi, + curtailment_14a=curtailment_14a, ) def to_graph(self): @@ -2330,6 +2349,165 @@ def apply_heat_pump_operating_strategy( """ hp_operating_strategy(self, strategy=strategy, heat_pump_names=heat_pump_names) + def apply_curtailment_14a_enwg( + self, + components: list[str] | None = None, + max_power_kw: float = 4.2, + components_type: str | None = None, + ) -> dict[str, float]: + """ + Apply §14a EnWG curtailment to heat pumps and/or charging points. + + §14a EnWG (Netzausbaugebiet) allows network operators to temporarily curtail + controllable consumption devices (especially heat pumps and charging points) + to a maximum power of 4.2 kW instead of performing grid expansion. + This can be used as a cost-effective alternative to grid reinforcement measures. + + The curtailment is applied to the existing time series data in + :attr:`~.network.timeseries.TimeSeries.loads_active_power` and + :attr:`~.network.timeseries.TimeSeries.loads_reactive_power`. + + Parameters + ---------- + components : list of str or None + List of component names (heat pumps and/or charging points) to apply + curtailment to. If None, all heat pumps and charging points in the network + are curtailed. Default: None. + max_power_kw : float + Maximum allowed power in kW after curtailment. According to §14a EnWG, + this is typically 4.2 kW. Default: 4.2. + components_type : str or None + Type of components to curtail. Can be 'heat_pump', 'charging_point', + or None. If None, both heat pumps and charging points are considered. + Only used if `components` is None. Default: None. + + Returns + ------- + dict of str to float + Dictionary with component names as keys and the curtailed energy + (in MWh) as values. This shows how much energy was curtailed for + each component. + + Notes + ----- + The curtailment is applied by limiting the active power time series + to the specified maximum power. The reactive power is adjusted + proportionally based on the power factor. + + The function returns information about curtailed energy which can be + used for cost-benefit analysis when comparing grid expansion costs + with curtailment compensation. + + Examples + -------- + >>> # Apply curtailment to all heat pumps and charging points + >>> curtailed = edisgo.apply_curtailment_14a_enwg() + >>> + >>> # Apply curtailment only to heat pumps + >>> curtailed = edisgo.apply_curtailment_14a_enwg( + ... components_type='heat_pump' + ... ) + >>> + >>> # Apply curtailment to specific components with different limit + >>> curtailed = edisgo.apply_curtailment_14a_enwg( + ... components=['Heat_pump_LVGrid_1_1', 'ChargingPoint_MVGrid_1_home_1'], + ... max_power_kw=3.0 + ... ) + + """ + return apply_curtailment_14a( + self, + components=components, + max_power_kw=max_power_kw, + components_type=components_type, + ) + + def check_curtailment_14a_effect( + self, + components: list[str], + max_power_kw: float = 4.2, + ) -> dict[str, float]: + """ + Check the effect of §14a EnWG curtailment without applying it. + + This function analyzes what the effect of curtailing specified components + according to §14a EnWG would be, without modifying the time series data. + This is useful for planning and cost-benefit analysis. + + Parameters + ---------- + components : list of str + List of component names (heat pumps and/or charging points) to check + curtailment effect for. + max_power_kw : float + Maximum allowed power in kW after curtailment. Default: 4.2. + + Returns + ------- + dict of str to float + Dictionary with the following keys: + - 'total_curtailed_energy_mwh': Total energy that would be curtailed + - 'max_simultaneous_curtailment_mw': Maximum simultaneous curtailed power + - 'avg_curtailed_power_mw': Average curtailed power across all time steps + - 'hours_with_curtailment': Number of hours with active curtailment + + Examples + -------- + >>> # Check effect before applying + >>> effect = edisgo.check_curtailment_14a_effect( + ... components=['Heat_pump_LVGrid_1_1'] + ... ) + >>> print(f"Would curtail {effect['total_curtailed_energy_mwh']:.2f} MWh") + + """ + return check_curtailment_effect( + self, + components=components, + max_power_kw=max_power_kw, + ) + + def identify_components_for_curtailment_14a( + self, + critical_components: list[str] | None = None, + curtailment_priority: str = "p_set", + ) -> list[str]: + """ + Identify which components should be curtailed according to §14a EnWG. + + This function can be used as part of the grid reinforcement optimization + to identify which components should be curtailed instead of performing + grid expansion. + + Parameters + ---------- + critical_components : list of str or None + List of components that are causing grid issues (overloading or + voltage violations). If None, all curtailable components are considered. + Default: None. + curtailment_priority : str + Defines how to prioritize components for curtailment. Options: + 'p_set' (largest nominal power first), 'random', 'grid_level' + (start with LV, then MV). Default: 'p_set'. + + Returns + ------- + list of str + Sorted list of component names to curtail, ordered by priority. + + Examples + -------- + >>> # Identify components to curtail, prioritizing largest ones + >>> components = edisgo.identify_components_for_curtailment_14a( + ... curtailment_priority='p_set' + ... ) + + """ + return identify_components_for_curtailment( + self, + critical_components=critical_components, + curtailment_priority=curtailment_priority, + ) + def import_dsm(self, scenario: str, timeindex=None): """ Gets industrial and CTS DSM profiles from the diff --git a/edisgo/flex_opt/__init__.py b/edisgo/flex_opt/__init__.py index e69de29bb..f0d981009 100644 --- a/edisgo/flex_opt/__init__.py +++ b/edisgo/flex_opt/__init__.py @@ -0,0 +1,21 @@ +""" +Flexibility and optimization modules for eDisGo. + +This module contains functions for grid flexibility measures, including: +- Charging strategies for electric vehicles +- Heat pump operation strategies +- Battery storage operation +- Grid reinforcement and optimization +- §14a EnWG curtailment for controllable consumers +""" + +__all__ = [ + "charging_strategies", + "heat_pump_operation", + "battery_storage_operation", + "reinforce_grid", + "curtailment_14a", + "check_tech_constraints", + "costs", + "q_control", +] diff --git a/edisgo/flex_opt/costs.py b/edisgo/flex_opt/costs.py index 7cbe359b1..d36fc4893 100644 --- a/edisgo/flex_opt/costs.py +++ b/edisgo/flex_opt/costs.py @@ -162,7 +162,7 @@ def _get_line_costs(lines_added): ]["quantity"].to_frame() lines_added_unique = lines_added.index.unique() lines_added = ( - lines_added.groupby(axis=0, level=0) + lines_added.groupby(level=0) .sum() .loc[lines_added_unique, ["quantity"]] ) diff --git a/edisgo/flex_opt/curtailment_14a.py b/edisgo/flex_opt/curtailment_14a.py new file mode 100644 index 000000000..e0afaa284 --- /dev/null +++ b/edisgo/flex_opt/curtailment_14a.py @@ -0,0 +1,424 @@ +""" +This module implements curtailment according to §14a EnWG (Netzausbaugebiet). + +§14a EnWG allows network operators to temporarily curtail controllable consumption +devices (especially heat pumps and charging points) to a MINIMUM power of 4.2 kW. + +IMPORTANT: This means devices can operate ABOVE 4.2 kW normally, and the operator +can reduce them DOWN TO (but not below) 4.2 kW. + +The functions in this module are primarily for worst-case analysis and testing. +For realistic §14a implementation in grid optimization, use the OPF with the +curtailment_14a parameter in edisgo.optimize() or edisgo.reinforce(). + +""" + +from __future__ import annotations + +import logging + +from typing import TYPE_CHECKING + +import pandas as pd + +if TYPE_CHECKING: + from edisgo import EDisGo + +logger = logging.getLogger(__name__) + + +def apply_curtailment_14a( + edisgo_obj: EDisGo, + components: list[str] | None = None, + max_power_kw: float = 4.2, + components_type: str | None = None, +) -> dict[str, pd.Series]: + """ + Apply §14a EnWG curtailment to heat pumps and/or charging points. + + WARNING: This function applies PERMANENT curtailment to time series data, + limiting all power values to a maximum of 4.2 kW. This is a WORST-CASE + scenario and NOT how §14a is typically used in practice! + + For realistic §14a implementation, use the OPF with curtailment_14a parameter: + >>> edisgo.optimize(curtailment_14a={'apply_curtailment': True}) + + §14a EnWG allows network operators to curtail devices DOWN TO a minimum of + 4.2 kW (not TO a maximum of 4.2 kW). Devices normally operate at their + requested power, and curtailment to 4.2 kW is only applied when needed. + + This function simulates a permanent worst-case where ALL power values are + limited to 4.2 kW, which is useful for comparison but not realistic. + + Parameters + ---------- + edisgo_obj : :class:`~.EDisGo` + The eDisGo object containing the network and time series data. + components : list of str or None + List of component names to apply curtailment to. If None, all + heat pumps and charging points in the network are curtailed. + Default: None. + max_power_kw : float + Maximum allowed power in kW after curtailment. According to §14a EnWG, + this is typically 4.2 kW. Default: 4.2. + components_type : str or None + Type of components to curtail. Can be 'heat_pump', 'charging_point', + or None. If None, both heat pumps and charging points are considered. + Only used if `components` is None. Default: None. + + Returns + ------- + dict of str to :pandas:`pandas.Series` + Dictionary with component names as keys and the curtailed energy + (in MWh) as values. This shows how much energy was curtailed for + each component. + + Notes + ----- + The curtailment is applied by limiting ALL active power time series values + to the specified maximum power. The reactive power is adjusted + proportionally based on the power factor. + + This function PERMANENTLY modifies time series data. For a non-destructive + analysis, use check_curtailment_14a_effect() instead. + + For realistic §14a optimization (where curtailment only happens when needed + to avoid grid issues), use the OPF: + >>> edisgo.optimize(curtailment_14a={'apply_curtailment': True}) + + """ + max_power_mw = max_power_kw / 1000.0 # Convert kW to MW + + # Get components to curtail + if components is None: + components = _get_curtailable_components(edisgo_obj, components_type) + else: + # Validate that provided components exist + all_loads = edisgo_obj.topology.loads_df.index.tolist() + invalid_components = [c for c in components if c not in all_loads] + if invalid_components: + raise ValueError( + f"The following components do not exist in the network: " + f"{invalid_components}" + ) + + if not components: + logger.warning( + "No components found for §14a curtailment. No curtailment applied." + ) + return {} + + logger.info( + f"Applying §14a EnWG curtailment to {len(components)} components " + f"with max power {max_power_kw} kW." + ) + + curtailed_energy = {} + + # Get time series data + ts_active = edisgo_obj.timeseries.loads_active_power + ts_reactive = edisgo_obj.timeseries.loads_reactive_power + + for component in components: + if component not in ts_active.columns: + logger.warning( + f"Component {component} has no active power time series. Skipping." + ) + continue + + # Get original time series + original_active = ts_active[component].copy() + + # Apply curtailment + curtailed_active = original_active.clip(upper=max_power_mw) + + # Calculate curtailed energy (difference between original and curtailed) + time_delta = ( + edisgo_obj.timeseries.timeindex[1] - edisgo_obj.timeseries.timeindex[0] + ).total_seconds() / 3600.0 # in hours + energy_curtailed = ((original_active - curtailed_active) * time_delta).sum() + curtailed_energy[component] = energy_curtailed + + # Update active power time series + ts_active[component] = curtailed_active + + # Adjust reactive power proportionally if it exists + if component in ts_reactive.columns: + original_reactive = ts_reactive[component].copy() + # Calculate ratio of curtailed to original power (avoid division by zero) + ratio = pd.Series(1.0, index=original_active.index) + non_zero_mask = original_active != 0 + ratio[non_zero_mask] = ( + curtailed_active[non_zero_mask] / original_active[non_zero_mask] + ) + # Apply same ratio to reactive power + ts_reactive[component] = original_reactive * ratio + + logger.info( + f"§14a curtailment applied. Total curtailed energy: " + f"{sum(curtailed_energy.values()):.2f} MWh across {len(curtailed_energy)} " + f"components." + ) + + return curtailed_energy + + +def _get_curtailable_components( + edisgo_obj: EDisGo, components_type: str | None = None +) -> list[str]: + """ + Get list of components that can be curtailed according to §14a EnWG. + + Parameters + ---------- + edisgo_obj : :class:`~.EDisGo` + The eDisGo object containing the network topology. + components_type : str or None + Type of components to get. Can be 'heat_pump', 'charging_point', + or None (both types). Default: None. + + Returns + ------- + list of str + List of component names that can be curtailed. + + """ + loads_df = edisgo_obj.topology.loads_df + + if components_type == "heat_pump": + curtailable = loads_df[loads_df.type == "heat_pump"].index.tolist() + elif components_type == "charging_point": + curtailable = loads_df[loads_df.type == "charging_point"].index.tolist() + elif components_type is None or components_type == "both": + # Get both heat pumps and charging points + curtailable = loads_df[ + loads_df.type.isin(["heat_pump", "charging_point"]) + ].index.tolist() + else: + raise ValueError( + f"Invalid components_type '{components_type}'. Must be " + f"'heat_pump', 'charging_point', 'both', or None." + ) + + return curtailable + + +def identify_components_for_curtailment( + edisgo_obj: EDisGo, + critical_components: list[str] | None = None, + curtailment_priority: str = "p_set", +) -> list[str]: + """ + Identify which components should be curtailed to avoid grid reinforcement. + + This function can be used as part of the grid reinforcement optimization + to identify which components should be curtailed according to §14a EnWG + instead of performing grid expansion. + + Parameters + ---------- + edisgo_obj : :class:`~.EDisGo` + The eDisGo object containing the network and analysis results. + critical_components : list of str or None + List of components that are causing grid issues (overloading or + voltage violations). If None, all curtailable components are considered. + Default: None. + curtailment_priority : str + Defines how to prioritize components for curtailment. Options: + 'p_set' (largest nominal power first), 'random', 'grid_level' + (start with LV, then MV). Default: 'p_set'. + + Returns + ------- + list of str + Sorted list of component names to curtail, ordered by priority. + + """ + # Get all curtailable components + all_curtailable = _get_curtailable_components(edisgo_obj, components_type=None) + + if critical_components is not None: + # Filter to only include critical components that are curtailable + components = [c for c in critical_components if c in all_curtailable] + else: + components = all_curtailable + + if not components: + return [] + + # Apply prioritization + loads_df = edisgo_obj.topology.loads_df.loc[components] + + if curtailment_priority == "p_set": + # Prioritize components with highest nominal power + sorted_components = loads_df.sort_values( + "p_set", ascending=False + ).index.tolist() + elif curtailment_priority == "random": + # Random order + sorted_components = loads_df.sample(frac=1).index.tolist() + elif curtailment_priority == "grid_level": + # Prioritize based on grid level (LV first, then MV) + buses_df = edisgo_obj.topology.buses_df + loads_with_grid = loads_df.copy() + loads_with_grid["v_nom"] = loads_with_grid["bus"].map(buses_df["v_nom"]) + # Sort by voltage level (ascending = LV first) + sorted_components = loads_with_grid.sort_values("v_nom").index.tolist() + else: + raise ValueError( + f"Invalid curtailment_priority '{curtailment_priority}'. " + f"Must be 'p_set', 'random', or 'grid_level'." + ) + + return sorted_components + + +def check_curtailment_effect( + edisgo_obj: EDisGo, + components: list[str], + max_power_kw: float = 4.2, +) -> dict[str, float]: + """ + Check the effect of curtailing components without actually applying it. + + This function analyzes what the effect of curtailing specified components + would be, without modifying the time series data. + + Parameters + ---------- + edisgo_obj : :class:`~.EDisGo` + The eDisGo object containing the network and time series data. + components : list of str + List of component names to check curtailment effect for. + max_power_kw : float + Maximum allowed power in kW after curtailment. Default: 4.2. + + Returns + ------- + dict of str to float + Dictionary with the following keys: + - 'total_curtailed_energy_mwh': Total energy that would be curtailed + - 'max_simultaneous_curtailment_mw': Maximum simultaneous curtailed power + - 'avg_curtailed_power_mw': Average curtailed power across all time steps + - 'hours_with_curtailment': Number of hours with active curtailment + + """ + max_power_mw = max_power_kw / 1000.0 + + ts_active = edisgo_obj.timeseries.loads_active_power + + # Calculate curtailment for each component + total_curtailment = pd.Series(0.0, index=ts_active.index) + + for component in components: + if component not in ts_active.columns: + continue + + original = ts_active[component] + curtailed = original.clip(upper=max_power_mw) + curtailment = original - curtailed + total_curtailment += curtailment + + # Calculate statistics + time_delta = ( + edisgo_obj.timeseries.timeindex[1] - edisgo_obj.timeseries.timeindex[0] + ).total_seconds() / 3600.0 + + total_curtailed_energy = (total_curtailment * time_delta).sum() + max_simultaneous_curtailment = total_curtailment.max() + avg_curtailed_power = total_curtailment.mean() + hours_with_curtailment = (total_curtailment > 0).sum() * time_delta + + return { + "total_curtailed_energy_mwh": total_curtailed_energy, + "max_simultaneous_curtailment_mw": max_simultaneous_curtailment, + "avg_curtailed_power_mw": avg_curtailed_power, + "hours_with_curtailment": hours_with_curtailment, + } + + +def apply_curtailment_during_reinforcement( + edisgo_obj: EDisGo, + max_power_kw: float | None = None, + components_type: str | None = None, + curtailment_priority: str | None = None, +) -> dict[str, float]: + """ + Apply §14a EnWG curtailment as part of grid reinforcement optimization. + + This function integrates §14a curtailment into the grid reinforcement workflow. + It reads configuration from the config file and applies curtailment to + controllable consumers before grid expansion is considered. + + This function is intended to be called from within the grid reinforcement + process (e.g., from :func:`~.flex_opt.reinforce_grid.reinforce_grid`). + + Parameters + ---------- + edisgo_obj : :class:`~.EDisGo` + The eDisGo object containing the network and time series data. + max_power_kw : float or None + Maximum allowed power in kW after curtailment. If None, the value from + the config is used. Default: None. + components_type : str or None + Type of components to curtail ('heat_pump', 'charging_point', 'both'). + If None, the value from the config is used. Default: None. + curtailment_priority : str or None + Curtailment priority strategy. If None, the value from the config is used. + Default: None. + + Returns + ------- + dict of str to float + Dictionary with component names as keys and the curtailed energy + (in MWh) as values. + + Notes + ----- + This function should only be called if curtailment is enabled in the config + (check `config['curtailment_14a_enwg']['enable_curtailment']`). + + """ + # Check if curtailment is enabled + config = edisgo_obj.config + if not config["curtailment_14a_enwg"]["enable_curtailment"]: + logger.debug("§14a curtailment is disabled in config.") + return {} + + # Read config values if not provided + if max_power_kw is None: + max_power_kw = float(config["curtailment_14a_enwg"]["max_power_kw"]) + + if components_type is None: + components_type = config["curtailment_14a_enwg"]["components_type"] + if components_type == "both": + components_type = None + + if curtailment_priority is None: + curtailment_priority = config["curtailment_14a_enwg"]["curtailment_priority"] + + # Identify components to curtail + components = identify_components_for_curtailment( + edisgo_obj, + critical_components=None, # Consider all components + curtailment_priority=curtailment_priority, + ) + + if not components: + logger.debug("No components available for §14a curtailment.") + return {} + + # Apply curtailment + logger.debug( + f"Applying §14a EnWG curtailment during grid reinforcement to " + f"{len(components)} components." + ) + + curtailed_energy = apply_curtailment_14a( + edisgo_obj, + components=components, + max_power_kw=max_power_kw, + components_type=components_type, + ) + + return curtailed_energy diff --git a/edisgo/flex_opt/reinforce_grid.py b/edisgo/flex_opt/reinforce_grid.py index 6da2449f5..c52fb3192 100644 --- a/edisgo/flex_opt/reinforce_grid.py +++ b/edisgo/flex_opt/reinforce_grid.py @@ -11,6 +11,7 @@ from edisgo.flex_opt import check_tech_constraints as checks from edisgo.flex_opt import exceptions, reinforce_measures from edisgo.flex_opt.costs import grid_expansion_costs +from edisgo.flex_opt.curtailment_14a import apply_curtailment_during_reinforcement from edisgo.flex_opt.reinforce_measures import separate_lv_grid from edisgo.tools import tools from edisgo.tools.temporal_complexity_reduction import get_most_critical_time_steps @@ -88,6 +89,10 @@ def reinforce_grid( This is used in case worst-case grid reinforcement is conducted in order to reinforce MV/LV stations for LV worst-cases. Default: False. + skip_curtailment_14a : bool + If True, §14a EnWG curtailment is not applied even if enabled in the config. + This can be used to compare scenarios with and without curtailment. + Default: False. num_steps_loading : int In case `reduced_analysis` is set to True, this parameter can be used to specify the number of most critical overloading events to consider. @@ -203,6 +208,36 @@ def reinforce_grid( scale_timeseries=scale_timeseries, ) + # CHECK IF §14a EnWG CURTAILMENT SHOULD BE APPLIED + # This allows curtailing heat pumps and charging points to 4.2 kW + # instead of performing grid expansion + # curtailment_applied = False + if edisgo.config["curtailment_14a_enwg"]["enable_curtailment"] and not kwargs.get( + "skip_curtailment_14a", False + ): + logger.info("==> Checking if §14a EnWG curtailment can be applied.") + try: + curtailed_energy = apply_curtailment_during_reinforcement(edisgo) + if curtailed_energy: + logger.info( + f"==> §14a curtailment applied to {len(curtailed_energy)} " + f"components. Total curtailed energy: " + f"{sum(curtailed_energy.values()):.2f} MWh" + ) + # Re-run power flow analysis with curtailed time series + edisgo.analyze( + mode=analyze_mode, + timesteps=timesteps_pfa, + lv_grid_id=lv_grid_id, + scale_timeseries=scale_timeseries, + ) + else: + logger.info("==> No components available for §14a curtailment.") + except Exception as e: + logger.warning( + f"==> §14a curtailment failed: {e}. Proceeding with grid expansion." + ) + # REINFORCE OVERLOADED TRANSFORMERS AND LINES logger.debug("==> Check station load.") overloaded_mv_station = ( diff --git a/edisgo/io/electromobility_import.py b/edisgo/io/electromobility_import.py index 7fcc0a2e0..d248530a6 100644 --- a/edisgo/io/electromobility_import.py +++ b/edisgo/io/electromobility_import.py @@ -846,7 +846,8 @@ def distribute_private_charging_demand(edisgo_obj): designated_charging_point_capacity_df = pd.DataFrame( index=user_centric_weights_df.index, columns=["designated_charging_point_capacity"], - data=0, + data=0.0, + dtype=float, ) for destination in private_charging_df.destination.sort_values().unique(): @@ -986,7 +987,8 @@ def distribute_public_charging_demand(edisgo_obj, **kwargs): designated_charging_point_capacity_df = pd.DataFrame( index=grid_and_user_centric_weights_df.index, columns=["designated_charging_point_capacity"], - data=0, + data=0.0, + dtype=float, ) columns = [ diff --git a/edisgo/io/generators_import.py b/edisgo/io/generators_import.py index ed72c4f86..8434af840 100755 --- a/edisgo/io/generators_import.py +++ b/edisgo/io/generators_import.py @@ -341,12 +341,14 @@ def _validate_sample_geno_location(): generators_conv_mv = _import_conv_generators(session) generators_res_mv, generators_res_lv = _import_res_generators(session) - generators_mv = pd.concat( - [ - generators_conv_mv, - generators_res_mv, - ] - ) + # Filter out empty DataFrames before concatenation to avoid FutureWarning + dfs_to_concat = [ + df for df in [generators_conv_mv, generators_res_mv] if not df.empty + ] + if dfs_to_concat: + generators_mv = pd.concat(dfs_to_concat) + else: + generators_mv = pd.DataFrame() # validate that imported generators are located inside the grid district _validate_sample_geno_location() @@ -465,14 +467,20 @@ def _check_mv_generator_geom(generator_data): cap_diff_threshold = 10**-4 # get all imported generators - imported_gens = pd.concat( - [imported_generators_lv, imported_generators_mv], sort=True - ) + # Filter out empty DataFrames before concatenation to avoid FutureWarning + dfs_to_concat = [ + df for df in [imported_generators_lv, imported_generators_mv] if not df.empty + ] + if dfs_to_concat: + imported_gens = pd.concat(dfs_to_concat, sort=True) + else: + imported_gens = pd.DataFrame() logger.debug(f"{len(imported_gens)} generators imported.") - # get existing generators and append ID column - existing_gens = edisgo_object.topology.generators_df + # get existing generators and append ID column (make copy to avoid + # SettingWithCopyWarning) + existing_gens = edisgo_object.topology.generators_df.copy() existing_gens["id"] = list( map(lambda _: int(_.split("_")[-1]), existing_gens.index) ) @@ -1052,15 +1060,15 @@ def _integrate_pv_rooftop(edisgo_object, pv_rooftop_df): suffixes=("_old", ""), ).set_index("gen_name") # add building id - edisgo_object.topology.generators_df.loc[ - gens_existing.index, "building_id" - ] = gens_existing.building_id + edisgo_object.topology.generators_df.loc[gens_existing.index, "building_id"] = ( + gens_existing.building_id + ) # update plants where capacity decreased gens_decreased_cap = gens_existing.query("p_nom < p_nom_old") if len(gens_decreased_cap) > 0: - edisgo_object.topology.generators_df.loc[ - gens_decreased_cap.index, "p_nom" - ] = gens_decreased_cap.p_nom + edisgo_object.topology.generators_df.loc[gens_decreased_cap.index, "p_nom"] = ( + gens_decreased_cap.p_nom + ) # update plants where capacity increased gens_increased_cap = gens_existing.query("p_nom > p_nom_old") for gen in gens_increased_cap.index: @@ -1073,9 +1081,9 @@ def _integrate_pv_rooftop(edisgo_object, pv_rooftop_df): if voltage_level_new >= voltage_level_old: # simply update p_nom if plant doesn't need to be connected to higher # voltage level - edisgo_object.topology.generators_df.at[ - gen, "p_nom" - ] = gens_increased_cap.at[gen, "p_nom"] + edisgo_object.topology.generators_df.at[gen, "p_nom"] = ( + gens_increased_cap.at[gen, "p_nom"] + ) else: # if plant needs to be connected to higher voltage level, remove existing # plant and integrate new component based on geolocation @@ -1183,10 +1191,10 @@ def _integrate_new_pv_rooftop_to_buildings(edisgo_object, pv_rooftop_df): # add voltage level for gen in pv_rooftop_df.index: - pv_rooftop_df.at[ - gen, "voltage_level" - ] = determine_grid_integration_voltage_level( - edisgo_object, pv_rooftop_df.at[gen, "p_nom"] + pv_rooftop_df.at[gen, "voltage_level"] = ( + determine_grid_integration_voltage_level( + edisgo_object, pv_rooftop_df.at[gen, "p_nom"] + ) ) # check for duplicated generator names and choose random name for duplicates @@ -1404,9 +1412,9 @@ def _integrate_new_power_plant(edisgo_object, comp_data): if voltage_level_new >= voltage_level_old: # simply update p_nom if plant doesn't need to be connected to higher # voltage level - edisgo_object.topology.generators_df.at[ - gen, "p_nom" - ] = gens_increased_cap.at[gen, "p_nom"] + edisgo_object.topology.generators_df.at[gen, "p_nom"] = ( + gens_increased_cap.at[gen, "p_nom"] + ) else: # if plant needs to be connected to higher voltage level, remove # existing plant and integrate new component based on geolocation diff --git a/edisgo/io/powermodels_io.py b/edisgo/io/powermodels_io.py index 541e01ffc..dea05faac 100644 --- a/edisgo/io/powermodels_io.py +++ b/edisgo/io/powermodels_io.py @@ -31,6 +31,7 @@ def to_powermodels( flexible_loads=None, flexible_storage_units=None, opf_version=1, + curtailment_14a=None, ): """ Convert eDisGo representation of the network topology and timeseries to @@ -58,6 +59,12 @@ def to_powermodels( Version of optimization models to choose from. Must be one of [1, 2, 3, 4]. For more information see :func:`edisgo.opf.powermodels_opf.pm_optimize`. Default: 1. + curtailment_14a : dict or None + Dictionary with §14a EnWG curtailment settings. If provided, curtailment + constraints will be added to the OPF. Dictionary should contain: + - 'max_power_mw': Maximum power in MW (e.g., 0.0042 for 4.2 kW) + - 'components': List of component names to curtail (CPs and/or HPs) + Default: None (no curtailment). Returns ------- @@ -209,6 +216,17 @@ def to_powermodels( opf_version = 2 pm["opf_version"] = opf_version + + # Add §14a EnWG curtailment settings if provided + if curtailment_14a is not None: + logger.info("Adding §14a EnWG curtailment constraints to PowerModels dict.") + pm["curtailment_14a"] = { + "max_power_mw": curtailment_14a.get("max_power_mw", 0.0042), + "components": curtailment_14a.get("components", []), + } + else: + pm["curtailment_14a"] = None + logger.info( "Transforming components timeseries into PowerModels dictionary format." ) @@ -250,10 +268,10 @@ def from_powermodels( Base value of apparent power for per unit system. Default: 1 MVA. """ - if type(pm_results) == str: + if isinstance(pm_results, str): with open(pm_results) as f: pm = json.loads(json.load(f)) - elif type(pm_results) == dict: + elif isinstance(pm_results, dict): pm = pm_results else: raise ValueError( @@ -363,9 +381,11 @@ def from_powermodels( for flex in df2.columns: abs_error = abs(df2[flex].values - hv_flex_dict[flex]) rel_error = [ - abs_error[i] / hv_flex_dict[flex][i] - if ((abs_error > 0.01)[i] & (hv_flex_dict[flex][i] != 0)) - else 0 + ( + abs_error[i] / hv_flex_dict[flex][i] + if ((abs_error > 0.01)[i] & (hv_flex_dict[flex][i] != 0)) + else 0 + ) for i in range(len(abs_error)) ] df2[flex] = rel_error @@ -478,6 +498,70 @@ def from_powermodels( edisgo_object.opf_results.grid_slacks_t.hp_load_shedding = df elif var == "phps2": edisgo_object.opf_results.grid_slacks_t.hp_operation_slack = df + + # Extract and log §14a EnWG curtailment results + try: + cp_curtailment_14a = _result_df( + pm, + "electromobility", + "pcp_14a_curt", + timesteps, + edisgo_object.timeseries.timeindex, + s_base, + ) + # Save §14a curtailment to opf_results + edisgo_object.opf_results.curtailment_14a_cp = cp_curtailment_14a + + # Log §14a curtailment statistics for charging points + if not cp_curtailment_14a.empty: + total_cp_curt = cp_curtailment_14a.sum().sum() + if total_cp_curt > 0: + logger.info( + f"§14a Charging Point Curtailment: {total_cp_curt:.4f} MWh total across " + f"{len(cp_curtailment_14a.columns)} components" + ) + # Find components with curtailment + cp_with_curt = cp_curtailment_14a.sum()[cp_curtailment_14a.sum() > 0] + if not cp_with_curt.empty: + logger.info( + f" Components curtailed: {', '.join(cp_with_curt.index.tolist())}" + ) + for comp in cp_with_curt.index: + logger.info( + f" {comp}: {cp_with_curt[comp]:.4f} MWh " + f"(max: {cp_curtailment_14a[comp].max():.4f} MW)" + ) + else: + logger.info("§14a Charging Point Curtailment: NOT used (0 MWh)") + except (KeyError, IndexError) as e: + logger.debug(f"No §14a charging point curtailment data found: {e}") + edisgo_object.opf_results.curtailment_14a_cp = pd.DataFrame() + + try: + hp_curtailment_14a = _result_df( + pm, + "heatpumps", + "php_14a_curt", + timesteps, + edisgo_object.timeseries.timeindex, + s_base, + ) + # Save §14a curtailment to opf_results + edisgo_object.opf_results.curtailment_14a_hp = hp_curtailment_14a + + # Log §14a curtailment statistics for heat pumps + if not hp_curtailment_14a.empty: + total_hp_curt = hp_curtailment_14a.sum().sum() + if total_hp_curt > 0: + logger.info( + f"§14a Heat Pump Curtailment: {total_hp_curt:.4f} MWh total across " + f"{len(hp_curtailment_14a.columns)} components" + ) + else: + logger.info("§14a Heat Pump Curtailment: NOT used (0 MWh)") + except (KeyError, IndexError) as e: + logger.debug(f"No §14a heat pump curtailment data found: {e}") + edisgo_object.opf_results.curtailment_14a_hp = pd.DataFrame() # save line flows and currents to edisgo object for variable in ["pf", "qf", "ccm"]: @@ -664,19 +748,19 @@ def _build_gen(edisgo_obj, psa_net, pm, flexible_storage_units, s_base): idx_bus = _mapping( psa_net, edisgo_obj, - gen.bus[gen_i], + gen.bus.iloc[gen_i], flexible_storage_units=flexible_storage_units, ) pf, sign = _get_pf(edisgo_obj, pm, idx_bus, "generator") q = [ - sign * np.tan(np.arccos(pf)) * gen.p_nom[gen_i], - sign * np.tan(np.arccos(pf)) * gen.p_nom_min[gen_i], + sign * np.tan(np.arccos(pf)) * gen.p_nom.iloc[gen_i], + sign * np.tan(np.arccos(pf)) * gen.p_nom_min.iloc[gen_i], ] pm[text][str(gen_i + 1)] = { - "pg": psa_net.generators_t.p_set[gen.index[gen_i]][0] / s_base, - "qg": psa_net.generators_t.q_set[gen.index[gen_i]][0] / s_base, - "pmax": gen.p_nom[gen_i].round(20) / s_base, - "pmin": gen.p_nom_min[gen_i].round(20) / s_base, + "pg": psa_net.generators_t.p_set[gen.index[gen_i]].iloc[0] / s_base, + "qg": psa_net.generators_t.q_set[gen.index[gen_i]].iloc[0] / s_base, + "pmax": gen.p_nom.iloc[gen_i].round(20) / s_base, + "pmin": gen.p_nom_min.iloc[gen_i].round(20) / s_base, "qmax": max(q).round(20) / s_base, "qmin": min(q).round(20) / s_base, "P": 0, @@ -684,7 +768,7 @@ def _build_gen(edisgo_obj, psa_net, pm, flexible_storage_units, s_base): "vg": 1, "pf": pf, "sign": sign, - "mbase": gen.p_nom[gen_i] / s_base, + "mbase": gen.p_nom.iloc[gen_i] / s_base, "gen_bus": idx_bus, "gen_status": 1, "name": gen.index[gen_i], @@ -707,13 +791,13 @@ def _build_gen(edisgo_obj, psa_net, pm, flexible_storage_units, s_base): pf, sign = _get_pf(edisgo_obj, pm, idx_bus, "storage_unit") p_g = max( [ - psa_net.storage_units_t.p_set[inflexible_storage_units[stor_i]][0], + psa_net.storage_units_t.p_set[inflexible_storage_units[stor_i]].iloc[0], 0.0, ] ) q_g = min( [ - psa_net.storage_units_t.q_set[inflexible_storage_units[stor_i]][0], + psa_net.storage_units_t.q_set[inflexible_storage_units[stor_i]].iloc[0], 0.0, ] ) @@ -793,38 +877,38 @@ def _build_branch(edisgo_obj, psa_net, pm, flexible_storage_units, s_base): idx_f_bus = _mapping( psa_net, edisgo_obj, - branches.bus0[branch_i], + branches.bus0.iloc[branch_i], flexible_storage_units=flexible_storage_units, ) idx_t_bus = _mapping( psa_net, edisgo_obj, - branches.bus1[branch_i], + branches.bus1.iloc[branch_i], flexible_storage_units=flexible_storage_units, ) pm["branch"][str(branch_i + 1)] = { "name": branches.index[branch_i], - "br_r": branches.r_pu[branch_i] * s_base, - "r": branches.r[branch_i], - "br_x": branches.x_pu[branch_i] * s_base, + "br_r": branches.r_pu.iloc[branch_i] * s_base, + "r": branches.r.iloc[branch_i], + "br_x": branches.x_pu.iloc[branch_i] * s_base, "f_bus": idx_f_bus, "t_bus": idx_t_bus, - "g_to": branches.g_pu[branch_i] / 2 * s_base, - "g_fr": branches.g_pu[branch_i] / 2 * s_base, - "b_to": branches.b_pu[branch_i] / 2 * s_base, - "b_fr": branches.b_pu[branch_i] / 2 * s_base, - "shift": shift[branch_i], + "g_to": branches.g_pu.iloc[branch_i] / 2 * s_base, + "g_fr": branches.g_pu.iloc[branch_i] / 2 * s_base, + "b_to": branches.b_pu.iloc[branch_i] / 2 * s_base, + "b_fr": branches.b_pu.iloc[branch_i] / 2 * s_base, + "shift": shift.iloc[branch_i], "br_status": 1.0, - "rate_a": branches.s_nom[branch_i].real / s_base, + "rate_a": branches.s_nom.iloc[branch_i].real / s_base, "rate_b": 250 / s_base, "rate_c": 250 / s_base, "angmin": -np.pi / 6, "angmax": np.pi / 6, - "transformer": bool(transformer[branch_i]), + "transformer": bool(transformer.iloc[branch_i]), "storage": False, - "tap": tap[branch_i], - "length": branches.length.fillna(1)[branch_i].round(20), - "cost": branches.capital_cost[branch_i].round(20), + "tap": tap.iloc[branch_i], + "length": branches.length.fillna(1).iloc[branch_i].round(20), + "cost": branches.capital_cost.iloc[branch_i].round(20), "storage_pf": 0, "index": branch_i + 1, } @@ -912,7 +996,7 @@ def _build_load( idx_bus = _mapping( psa_net, edisgo_obj, - loads_df.bus[load_i], + loads_df.bus.iloc[load_i], flexible_storage_units=flexible_storage_units, ) if ( @@ -938,8 +1022,8 @@ def _build_load( p_d = psa_net.loads_t.p_set[loads_df.index[load_i]] q_d = psa_net.loads_t.q_set[loads_df.index[load_i]] pm["load"][str(load_i + 1)] = { - "pd": p_d[0].round(20) / s_base, - "qd": q_d[0].round(20) / s_base, + "pd": round(p_d.iloc[0], 20) / s_base, + "qd": round(q_d.iloc[0], 20) / s_base, "load_bus": idx_bus, "status": True, "pf": pf, @@ -958,19 +1042,19 @@ def _build_load( pf, sign = _get_pf(edisgo_obj, pm, idx_bus, "storage_unit") p_d = -min( [ - psa_net.storage_units_t.p_set[inflexible_storage_units[stor_i]][0], + psa_net.storage_units_t.p_set[inflexible_storage_units[stor_i]].iloc[0], np.float64(0.0), ] ) q_d = -max( [ - psa_net.storage_units_t.q_set[inflexible_storage_units[stor_i]][0], + psa_net.storage_units_t.q_set[inflexible_storage_units[stor_i]].iloc[0], np.float64(0.0), ] ) pm["load"][str(stor_i + len(loads_df.index) + 1)] = { - "pd": p_d.round(20) / s_base, - "qd": q_d.round(20) / s_base, + "pd": round(p_d, 20) / s_base, + "qd": round(q_d, 20) / s_base, "load_bus": idx_bus, "status": True, "pf": pf, @@ -1216,20 +1300,20 @@ def _build_heatpump(psa_net, pm, edisgo_obj, s_base, flexible_hps): ) ) for hp_i in np.arange(len(heat_df.index)): - idx_bus = _mapping(psa_net, edisgo_obj, heat_df.bus[hp_i]) + idx_bus = _mapping(psa_net, edisgo_obj, heat_df.bus.iloc[hp_i]) # retrieve power factor and sign from config pf, sign = _get_pf(edisgo_obj, pm, idx_bus, "heat_pump") - q = sign * np.tan(np.arccos(pf)) * heat_df.p_set[hp_i] + q = sign * np.tan(np.arccos(pf)) * heat_df.p_set.iloc[hp_i] p_d = heat_df2[heat_df.index[hp_i]] pm["heatpumps"][str(hp_i + 1)] = { - "pd": p_d[0].round(20) / s_base, # heat demand + "pd": p_d.iloc[0].round(20) / s_base, # heat demand "pf": pf, "sign": sign, "p_min": 0, - "p_max": heat_df.p_set[hp_i].round(20) / s_base, + "p_max": heat_df.p_set.iloc[hp_i].round(20) / s_base, "q_min": min(q, 0).round(20) / s_base, "q_max": max(q, 0).round(20) / s_base, - "cop": hp_cop[heat_df.index[hp_i]][0].round(20), + "cop": hp_cop[heat_df.index[hp_i]].iloc[0].round(20), "hp_bus": idx_bus, "name": heat_df.index[hp_i], "index": hp_i + 1, @@ -1325,7 +1409,7 @@ def _build_heat_storage(psa_net, pm, edisgo_obj, s_base, flexible_hps, opf_versi heat_storage_df = heat_storage_df.loc[flexible_hps] for stor_i in np.arange(len(flexible_hps)): idx_bus = _mapping( - psa_net, edisgo_obj, psa_net.loads.loc[flexible_hps].bus[stor_i] + psa_net, edisgo_obj, psa_net.loads.loc[flexible_hps].bus.iloc[stor_i] ) if ( edisgo_obj.topology.loads_df.loc[heat_storage_df.index[stor_i]].sector @@ -1338,9 +1422,9 @@ def _build_heat_storage(psa_net, pm, edisgo_obj, s_base, flexible_hps, opf_versi "ps": 0, "p_loss": p_loss, # 4% of SOC per day "energy": 0, - "capacity": heat_storage_df.capacity[stor_i].round(20) / s_base, - "charge_efficiency": heat_storage_df.efficiency[stor_i].round(20), - "discharge_efficiency": heat_storage_df.efficiency[stor_i].round(20), + "capacity": heat_storage_df.capacity.iloc[stor_i].round(20) / s_base, + "charge_efficiency": heat_storage_df.efficiency.iloc[stor_i].round(20), + "discharge_efficiency": heat_storage_df.efficiency.iloc[stor_i].round(20), "storage_bus": idx_bus, "name": heat_storage_df.index[stor_i], "soc_initial": ( diff --git a/edisgo/io/pypsa_io.py b/edisgo/io/pypsa_io.py index 023e5fd69..bf2614c72 100755 --- a/edisgo/io/pypsa_io.py +++ b/edisgo/io/pypsa_io.py @@ -355,6 +355,18 @@ def _set_slack(grid): if kwargs.get("use_seed", False) and pypsa_network.mode == "mv": set_seed(edisgo_object, pypsa_network) + # Ensure all time-varying data has float64 dtype to prevent scipy sparse errors + for component_type in ['loads_t', 'generators_t', 'storage_units_t']: + component_t = getattr(pypsa_network, component_type, None) + if component_t is not None: + for attr in ['p', 'q', 'p_set', 'q_set']: + if hasattr(component_t, attr): + df = getattr(component_t, attr) + if df is not None and isinstance(df, pd.DataFrame) and not df.empty: + # Convert all columns to float64 + for col in df.columns: + df[col] = df[col].astype('float64', errors='ignore') + return pypsa_network diff --git a/edisgo/network/timeseries.py b/edisgo/network/timeseries.py index 47c537e43..e5230d7ce 100644 --- a/edisgo/network/timeseries.py +++ b/edisgo/network/timeseries.py @@ -1804,19 +1804,21 @@ def _get_q_sign_and_power_factor_per_component( ], inplace=True, ) - self.time_series_raw.q_control = pd.concat( - [ - self.time_series_raw.q_control, - pd.DataFrame( - index=components_names, - data={ - "type": "fixed_cosphi", - "q_sign": q_sign, - "power_factor": power_factor, - }, - ), - ] + new_data = pd.DataFrame( + index=components_names, + data={ + "type": "fixed_cosphi", + "q_sign": q_sign, + "power_factor": power_factor, + }, ) + frames_to_concat = [ + df for df in [self.time_series_raw.q_control, new_data] if not df.empty + ] + if frames_to_concat: + self.time_series_raw.q_control = pd.concat(frames_to_concat) + elif not new_data.empty: + self.time_series_raw.q_control = new_data return q_sign, power_factor # set reactive power for generators diff --git a/edisgo/network/topology.py b/edisgo/network/topology.py index 1eb0a9272..a75d0d609 100755 --- a/edisgo/network/topology.py +++ b/edisgo/network/topology.py @@ -171,10 +171,10 @@ def _load_equipment_data(config=None): config = {} for voltage_level, eq_list in equipment.items(): for i in eq_list: - config[ - "equipment_{}_parameters_{}".format(voltage_level, i) - ] = "equipment-parameters_{}_{}.csv".format( - voltage_level.upper(), i + config["equipment_{}_parameters_{}".format(voltage_level, i)] = ( + "equipment-parameters_{}_{}.csv".format( + voltage_level.upper(), i + ) ) else: equipment_dir = config["system_dirs"]["equipment_dir"] @@ -1447,9 +1447,7 @@ def _get_line_data(): line_name = "Line_{}_{}".format(bus0, bus1) while line_name in self.lines_df.index: random.seed(a=line_name) - line_name = "Line_{}_{}_{}".format( - bus0, bus1, random.randint(10**8, 10**9) - ) + line_name = "Line_{}_{}_{}".format(bus0, bus1, random.randint(10**8, 10**9)) # check if all necessary data is now available if b is None: @@ -1560,7 +1558,7 @@ def remove_load(self, name): """ if name in self.loads_df.index: bus = self.loads_df.at[name, "bus"] - self._loads_df.drop(name, inplace=True) + self._loads_df = self._loads_df.drop(name) # if no other elements are connected, remove line and bus as well if self._check_bus_for_removal(bus): @@ -1583,7 +1581,7 @@ def remove_generator(self, name): """ if name in self.generators_df.index: bus = self.generators_df.at[name, "bus"] - self._generators_df.drop(name, inplace=True) + self._generators_df = self._generators_df.drop(name) # if no other elements are connected to same bus, remove line # and bus @@ -1610,7 +1608,7 @@ def remove_storage_unit(self, name): # remove storage unit and time series if name in self.storage_units_df.index: bus = self.storage_units_df.at[name, "bus"] - self._storage_units_df.drop(name, inplace=True) + self._storage_units_df = self._storage_units_df.drop(name) # if no other elements are connected, remove line and bus as well if self._check_bus_for_removal(bus): @@ -1723,9 +1721,9 @@ def update_number_of_parallel_lines(self, lines_num_parallel): ) # update number parallel lines - self._lines_df.loc[ - lines_num_parallel.index, "num_parallel" - ] = lines_num_parallel + self._lines_df.loc[lines_num_parallel.index, "num_parallel"] = ( + lines_num_parallel + ) def change_line_type(self, lines, new_line_type): """ @@ -1878,7 +1876,7 @@ def connect_to_mv(self, edisgo_object, comp_data, comp_type="generator"): power = comp_data.pop("p") # create new bus for new component - if type(comp_data["geom"]) != Point: + if not isinstance(comp_data["geom"], Point): geom = wkt_loads(comp_data["geom"]) else: geom = comp_data["geom"] @@ -3197,9 +3195,9 @@ def aggregate_lv_grid_at_station(self, lv_grid_id: int | str) -> None: self.buses_df = self.buses_df[~self.buses_df.index.isin(buses_to_drop)] self.lines_df = self.lines_df[~self.lines_df.index.isin(lines_to_drop)] self.loads_df.loc[self.loads_df.bus.isin(buses_to_drop), "bus"] = station_bus - self.generators_df.loc[ - self.generators_df.bus.isin(buses_to_drop), "bus" - ] = station_bus + self.generators_df.loc[self.generators_df.bus.isin(buses_to_drop), "bus"] = ( + station_bus + ) self.storage_units_df.loc[ self.storage_units_df.bus.isin(buses_to_drop), "bus" ] = station_bus diff --git a/edisgo/opf/eDisGo_OPF.jl/src/core/constraint_curtailment_14a.jl b/edisgo/opf/eDisGo_OPF.jl/src/core/constraint_curtailment_14a.jl new file mode 100644 index 000000000..c26fa3ebb --- /dev/null +++ b/edisgo/opf/eDisGo_OPF.jl/src/core/constraint_curtailment_14a.jl @@ -0,0 +1,97 @@ +# §14a EnWG Curtailment Constraints +# Diese Constraints implementieren §14a EnWG Abregelung als Alternative zu Netzausbau +# §14a erlaubt Netzbetreibern die Leistung auf MINIMAL 4.2 kW zu begrenzen + +""" +Add §14a EnWG curtailment constraints for charging points and heat pumps. + +§14a EnWG allows network operators to curtail devices to a MINIMUM of 4.2 kW. +This is implemented as a HARD constraint - devices cannot go below 4.2 kW! + +Correct interpretation: +- Devices normally run at their requested power (can be > 4.2 kW) +- §14a allows curtailment DOWN TO a minimum of 4.2 kW +- Below 4.2 kW is NOT allowed (hard constraint, no slack!) + +Formulation: + pcp >= curtailment_limit (HARD constraint) + pcp_14a_curt >= p_max - pcp (tracking variable for cost calculation) + +Where: +- Device power is bounded: [4.2 kW, p_max] +- If grid issues require < 4.2 kW: other constraints will be violated instead +- This ensures §14a minimum is ALWAYS respected + +Cost structure (hierarchy from cheap to expensive): +1. Normal flexibility (storage, redispatch): Factor 0.6 +2. §14a curtailment: Factor 100 (expensive! last resort before grid violations) +3. Grid violations (voltage/current): Factor 10,000 + +Examples: +- No grid issues: Device runs at p_max (e.g., 11 kW), pcp_14a_curt = 0, no §14a costs +- Grid congestion: Device curtailed to 6 kW, pcp_14a_curt = 5 kW, costs = 100 × 5 = 500 +- Critical: Device at 4.2 kW minimum, pcp_14a_curt = 6.8 kW, costs = 100 × 6.8 = 680 +- Impossible: If even 4.2 kW too much → grid violations occur (cost: 10,000) + +This ensures §14a is used only when normal flexibility is insufficient, +but prevents grid violations when possible. +""" +function constraint_curtailment_14a!(pm::AbstractPowerModel, nw::Int) + # Check if curtailment_14a is defined in the data + curtailment_data = get(PowerModels.ref(pm, nw), :curtailment_14a, nothing) + + if curtailment_data === nothing || curtailment_data == "nothing" + # No curtailment constraints + return + end + + curtailment_limit_mw = get(curtailment_data, "max_power_mw", 0.0042) # 4.2 kW MINIMUM + components = get(curtailment_data, "components", []) + + # Apply HARD minimum power constraint to charging points + for (i, cp) in PowerModels.ref(pm, nw, :electromobility) + if isempty(components) || get(cp, "name", "") in components + pcp = PowerModels.var(pm, nw, :pcp, i) + pcp_14a_curt = PowerModels.var(pm, nw, :pcp_14a_curt, i) + + # Only apply constraint if p_max > curtailment_limit + # (doesn't make sense to enforce 4.2 kW minimum on a 3 kW device) + if cp["p_max"] >= curtailment_limit_mw + # HARD CONSTRAINT: pcp >= 4.2 kW + # No slack variable - this is absolute! + # + # Examples for 11 kW device: + # pcp can be: [4.2, 11.0] kW + # pcp cannot be: < 4.2 kW (infeasible!) + # + # If grid cannot handle 4.2 kW: grid violations will occur + JuMP.@constraint(pm.model, pcp >= curtailment_limit_mw) + + # Track curtailment amount for cost calculation + # pcp_14a_curt measures how much power is curtailed due to §14a + # Example: p_max=11 kW, pcp=6 kW → pcp_14a_curt=5 kW (curtailment amount) + # MUST be equality (==) so negative costs work correctly! + JuMP.@constraint(pm.model, pcp_14a_curt == cp["p_max"] - pcp) + end + end + end + + # Apply HARD minimum power constraint to heat pumps + for (i, hp) in PowerModels.ref(pm, nw, :heatpumps) + if isempty(components) || get(hp, "name", "") in components + php = PowerModels.var(pm, nw, :php, i) + php_14a_curt = PowerModels.var(pm, nw, :php_14a_curt, i) + + # Only apply constraint if p_max > curtailment_limit + if hp["p_max"] >= curtailment_limit_mw + # HARD CONSTRAINT: php >= 4.2 kW (no slack!) + JuMP.@constraint(pm.model, php >= curtailment_limit_mw) + + # Track curtailment amount for cost calculation + # php_14a_curt measures how much power is curtailed due to §14a + # MUST be equality (==) so negative costs work correctly! + JuMP.@constraint(pm.model, php_14a_curt == hp["p_max"] - php) + end + end + end +end diff --git a/edisgo/opf/eDisGo_OPF.jl/src/core/objective.jl b/edisgo/opf/eDisGo_OPF.jl/src/core/objective.jl index 5a8470750..c9dc88734 100644 --- a/edisgo/opf/eDisGo_OPF.jl/src/core/objective.jl +++ b/edisgo/opf/eDisGo_OPF.jl/src/core/objective.jl @@ -28,7 +28,17 @@ function objective_min_losses_slacks(pm::AbstractBFModelEdisgo) phps = Dict(n => PowerModels.var(pm, n, :phps) for n in nws) phps2 = Dict(n => PowerModels.var(pm, n, :phps2) for n in nws) phss = Dict(n => PowerModels.var(pm, n, :phss) for n in nws) + pcp_14a_curt = Dict(n => PowerModels.var(pm, n, :pcp_14a_curt) for n in nws) + php_14a_curt = Dict(n => PowerModels.var(pm, n, :php_14a_curt) for n in nws) + + # Cost factors: + # - Normal flexibility (storage, redispatch): 0.6 + # - §14a curtailment: 100 (expensive! use as last resort before grid violations) + # - Grid violations (voltage/current): 10,000 (most expensive) + factor_slacks = 0.6 + factor_14a = 100 + return JuMP.@objective(pm.model, Min, (1-factor_slacks) * sum(sum(ccm[n][b] * r[n][b] for (b,i,j) in PowerModels.ref(pm, n, :arcs_from) ) for n in nws) # minimize line losses incl. storage losses + factor_slacks * sum(sum(pgc[n][i] for i in keys(PowerModels.ref(pm,1 , :gen_nd))) for n in nws) # minimize non-dispatchable curtailment @@ -36,6 +46,8 @@ function objective_min_losses_slacks(pm::AbstractBFModelEdisgo) + factor_slacks * sum(sum(pds[n][i] for i in keys(PowerModels.ref(pm,1 , :load))) for n in nws) # minimize load shedding + factor_slacks * sum(sum(pcps[n][i] for i in keys(PowerModels.ref(pm,1 , :electromobility))) for n in nws) # minimize cp load sheddin + factor_slacks * sum(sum(phps[n][i] for i in keys(PowerModels.ref(pm,1 , :heatpumps))) for n in nws) # minimize hp load shedding + + factor_14a * sum(sum(pcp_14a_curt[n][i] for i in keys(PowerModels.ref(pm, 1, :electromobility))) for n in nws) # minimize §14a CP curtailment (expensive!) + + factor_14a * sum(sum(php_14a_curt[n][i] for i in keys(PowerModels.ref(pm, 1, :heatpumps))) for n in nws) # minimize §14a HP curtailment (expensive!) + 1e4 * sum(sum(phss[n][i] + phps2[n][i] for i in keys(PowerModels.ref(pm, 1 , :heatpumps))) for n in nws) ) end @@ -69,12 +81,18 @@ function objective_min_losses_slacks_OG(pm::AbstractBFModelEdisgo) phps2 = Dict(n => PowerModels.var(pm, n, :phps2) for n in nws) phss = Dict(n => PowerModels.var(pm, n, :phss) for n in nws) phvs = Dict(n => PowerModels.var(pm, n, :phvs) for n in nws) + pcp_14a_curt = Dict(n => PowerModels.var(pm, n, :pcp_14a_curt) for n in nws) + php_14a_curt = Dict(n => PowerModels.var(pm, n, :php_14a_curt) for n in nws) + parameters = [r[1][i] for i in keys(r[1])] parameters = parameters[parameters .>0] #factor_hv_slacks = length(nws) * exp10(floor(log10(maximum(parameters)))+2) factor_hv_slacks = exp10(floor(log10(maximum(parameters)))+1) #println(factor_hv_slacks) factor_slacks = 0.6 + factor_14a = 100 + + return JuMP.@objective(pm.model, Min, (1-factor_slacks) * sum(sum(ccm[n][b]*r[n][b] for (b,i,j) in PowerModels.ref(pm, n, :arcs_from)) for n in nws) # minimize line losses + factor_slacks * sum(sum(pgc[n][i] for i in keys(PowerModels.ref(pm,1 , :gen_nd))) for n in nws) # minimize non-dispatchable curtailment @@ -82,6 +100,8 @@ function objective_min_losses_slacks_OG(pm::AbstractBFModelEdisgo) + factor_slacks * sum(sum(pds[n][i] for i in keys(PowerModels.ref(pm,1 , :load))) for n in nws) # minimize load shedding + factor_slacks * sum(sum(pcps[n][i] for i in keys(PowerModels.ref(pm,1 , :electromobility))) for n in nws) # minimize cp load shedding + factor_slacks * sum(sum(phps[n][i] for i in keys(PowerModels.ref(pm, 1 , :heatpumps))) for n in nws) # minimize hp load shedding + + factor_14a * sum(sum(pcp_14a_curt[n][i] for i in keys(PowerModels.ref(pm, 1, :electromobility))) for n in nws) # minimize §14a CP curtailment (expensive!) + + factor_14a * sum(sum(php_14a_curt[n][i] for i in keys(PowerModels.ref(pm, 1, :heatpumps))) for n in nws) # minimize §14a HP curtailment (expensive!) + factor_hv_slacks * sum(sum(phvs[n][i]^2 * flex["count"] for (i, flex) in PowerModels.ref(pm, n, :HV_requirements) if flex["name"]!= "dsm") for n in nws) # + factor_hv_slacks * 1e-1 * sum(sum(phvs[n][i]^2 * flex["count"] for (i, flex) in PowerModels.ref(pm, n, :HV_requirements) if flex["name"]== "dsm") for n in nws) # + 1e4 * sum(sum(phss[n][i] + phps2[n][i] for i in keys(PowerModels.ref(pm, 1 , :heatpumps))) for n in nws) diff --git a/edisgo/opf/eDisGo_OPF.jl/src/core/variables.jl b/edisgo/opf/eDisGo_OPF.jl/src/core/variables.jl index 5aec7f71c..f87bfef14 100644 --- a/edisgo/opf/eDisGo_OPF.jl/src/core/variables.jl +++ b/edisgo/opf/eDisGo_OPF.jl/src/core/variables.jl @@ -485,6 +485,32 @@ function variable_ev_slack(pm::AbstractBFModelEdisgo; nw::Int=nw_id_default, bou report && PowerModels.sol_component_value(pm, nw, :electromobility, :pcps, PowerModels.ids(pm, nw, :electromobility), pcps) end +"Charging point §14a curtailment tracking variable" +function variable_cp_14a_curtailment(pm::AbstractBFModelEdisgo; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true) + # Upper bound is p_max (maximum possible curtailment when pcp = 0) + # This prevents unbounded maximization with negative cost factors + pcp_14a_curt = PowerModels.var(pm, nw)[:pcp_14a_curt] = JuMP.@variable(pm.model, + [i in PowerModels.ids(pm, nw, :electromobility)], base_name="$(nw)_pcp_14a_curt", + lower_bound = 0.0, + upper_bound = PowerModels.ref(pm, nw, :electromobility, i)["p_max"] + ) + + report && PowerModels.sol_component_value(pm, nw, :electromobility, :pcp_14a_curt, PowerModels.ids(pm, nw, :electromobility), pcp_14a_curt) +end + +"Heat pump §14a curtailment tracking variable" +function variable_hp_14a_curtailment(pm::AbstractBFModelEdisgo; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true) + # Upper bound is p_max (maximum possible curtailment when php = 0) + # This prevents unbounded maximization with negative cost factors + php_14a_curt = PowerModels.var(pm, nw)[:php_14a_curt] = JuMP.@variable(pm.model, + [i in PowerModels.ids(pm, nw, :heatpumps)], base_name="$(nw)_php_14a_curt", + lower_bound = 0.0, + upper_bound = PowerModels.ref(pm, nw, :heatpumps, i)["p_max"] + ) + + report && PowerModels.sol_component_value(pm, nw, :heatpumps, :php_14a_curt, PowerModels.ids(pm, nw, :heatpumps), php_14a_curt) +end + "slack generator variables" function variable_slack_gen(pm::AbstractBFModelEdisgo; kwargs...) eDisGo_OPF.variable_slack_gen_real(pm; kwargs...) @@ -534,4 +560,9 @@ function variable_slack_HV_requirements_imaginary(pm::AbstractPowerModel; nw::In end +# Note: §14a EnWG slack variables have been removed! +# §14a is now implemented as HARD constraint - no violations allowed. +# If grid cannot handle minimum 4.2 kW, grid violations (voltage/current) occur instead. +# This ensures the 4.2 kW minimum is ALWAYS respected when §14a is active. + "" diff --git a/edisgo/opf/eDisGo_OPF.jl/src/eDisGo_OPF.jl b/edisgo/opf/eDisGo_OPF.jl/src/eDisGo_OPF.jl index 33df80742..bb0bd923b 100644 --- a/edisgo/opf/eDisGo_OPF.jl/src/eDisGo_OPF.jl +++ b/edisgo/opf/eDisGo_OPF.jl/src/eDisGo_OPF.jl @@ -18,6 +18,7 @@ include("core/types.jl") include("core/base.jl") include("core/constraint.jl") include("core/constraint_template.jl") +include("core/constraint_curtailment_14a.jl") include("core/data.jl") include("core/objective.jl") include("core/solution.jl") diff --git a/edisgo/opf/eDisGo_OPF.jl/src/prob/opf_bf.jl b/edisgo/opf/eDisGo_OPF.jl/src/prob/opf_bf.jl index aeaea8d65..5d98ec2aa 100644 --- a/edisgo/opf/eDisGo_OPF.jl/src/prob/opf_bf.jl +++ b/edisgo/opf/eDisGo_OPF.jl/src/prob/opf_bf.jl @@ -31,6 +31,10 @@ function build_mn_opf_bf_flex(pm::AbstractBFModelEdisgo) eDisGo_OPF.variable_dsm_storage_power(pm, nw=n) # Eq. (3.34), (3.35) eDisGo_OPF.variable_slack_gen(pm, nw=n) # keine Bounds für Slack Generator + # §14a tracking variables to measure curtailment amount (for cost calculation) + eDisGo_OPF.variable_cp_14a_curtailment(pm, nw=n) + eDisGo_OPF.variable_hp_14a_curtailment(pm, nw=n) + if PowerModels.ref(pm, 1, :opf_version) in(3, 4) # Nicht Teil der MA eDisGo_OPF.variable_slack_HV_requirements(pm, nw=n) if PowerModels.ref(pm, 1, :opf_version) in(3) @@ -58,6 +62,9 @@ function build_mn_opf_bf_flex(pm::AbstractBFModelEdisgo) eDisGo_OPF.constraint_hp_operation(pm, i, n) # Eq. (3.19) end + # §14a EnWG curtailment constraint + eDisGo_OPF.constraint_curtailment_14a!(pm, n) + end # CONSTRAINTS diff --git a/edisgo/opf/powermodels_opf.py b/edisgo/opf/powermodels_opf.py index 85da160a8..b03cf56ef 100644 --- a/edisgo/opf/powermodels_opf.py +++ b/edisgo/opf/powermodels_opf.py @@ -26,6 +26,7 @@ def pm_optimize( method: str = "soc", warm_start: bool = False, silence_moi: bool = False, + curtailment_14a: Optional[dict] = None, ) -> None: """ Run OPF for edisgo object in julia subprocess and write results of OPF to edisgo @@ -93,6 +94,13 @@ def pm_optimize( hence there will be no logging coming from julia subprocess in python process. Default: False. + curtailment_14a : dict or None + Dictionary with §14a EnWG curtailment settings. If provided, curtailment + constraints will be added to the OPF for charging points and heat pumps. + Dictionary should contain: + - 'max_power_mw': Maximum power in MW (default: 0.0042 for 4.2 kW) + - 'components': List of component names to apply curtailment to + Default: None (no curtailment constraints). save_heat_storage : bool Indicates whether to save results of heat storage variables from the optimization to eDisGo object. @@ -118,6 +126,7 @@ def pm_optimize( flexible_loads=flexible_loads, flexible_storage_units=flexible_storage_units, opf_version=opf_version, + curtailment_14a=curtailment_14a, ) def _convert(o): @@ -163,6 +172,14 @@ def _convert(o): hv_flex_dict=hv_flex_dict, s_base=s_base, ) + # Fix dtypes after writing OPF results to prevent scipy sparse matrix errors + for ts_df_name in ['loads_active_power', 'loads_reactive_power', + 'generators_active_power', 'generators_reactive_power', + 'storage_units_active_power', 'storage_units_reactive_power']: + ts_df = getattr(edisgo_obj.timeseries, ts_df_name, None) + if ts_df is not None and not ts_df.empty: + for col in ts_df.columns: + ts_df[col] = ts_df[col].astype('float64') elif out.rstrip().startswith("Set parameter") or out.rstrip().startswith( "Academic" ): diff --git a/edisgo/tools/tools.py b/edisgo/tools/tools.py index cf94fbc7d..053d337d0 100644 --- a/edisgo/tools/tools.py +++ b/edisgo/tools/tools.py @@ -581,6 +581,7 @@ def assign_voltage_level_to_component(df, buses_df): (either 'mv' or 'lv'). """ + df = df.copy() df["voltage_level"] = df.apply( lambda _: "lv" if buses_df.at[_.bus, "v_nom"] < 1 else "mv", axis=1,