diff --git a/doc/about.rst b/doc/about.rst index ad11ac9a..9ca4d10e 100644 --- a/doc/about.rst +++ b/doc/about.rst @@ -19,8 +19,9 @@ Research projects ==================== This software project was initially developed in the research project `open_eGo `_. -It is constantly further developed in different reserach projects, -e.g. `eGon `_ and `PoWerD `_. +It is constantly further developed in different research projects, +e.g. `eGon `_, `PoWerD `_, and +`hyBit `_. The OpenEnergy Platform diff --git a/etrago/analyze/plot.py b/etrago/analyze/plot.py index d84462b4..a0594098 100644 --- a/etrago/analyze/plot.py +++ b/etrago/analyze/plot.py @@ -46,12 +46,12 @@ if "READTHEDOCS" not in os.environ: from geoalchemy2.shape import to_shape # noqa: F401 from pyproj import Proj, transform - from shapely.geometry import LineString, Point + from shapely.geometry import LineString, Point, Polygon import geopandas as gpd import tilemapbase from etrago.execute import import_gen_from_links - from etrago.tools.utilities import find_buses_area + from etrago.tools.utilities import find_buses_area, select_buses_area __copyright__ = ( "Flensburg University of Applied Sciences, " @@ -873,6 +873,7 @@ def calc_storage_expansion_per_bus( "rural_heat_store", "central_heat_store", ], + buses=None, ): """Function that calculates storage expansion per bus and technology @@ -887,10 +888,13 @@ def calc_storage_expansion_per_bus( storage expansion per bus and technology """ - index = [(idx, "battery") for idx in network.buses.index] + if buses is None: + buses = network.buses.index + + index = [(idx, "battery") for idx in network.buses.loc[buses].index] for c in carriers: if c != "battery": - index.extend([(idx, c) for idx in network.buses.index]) + index.extend([(idx, c) for idx in network.buses.loc[buses].index]) # index.extend([(idx, 'hydrogen_storage') for idx in network.buses.index]) dist = pd.Series( @@ -900,7 +904,8 @@ def calc_storage_expansion_per_bus( if "battery" in carriers: batteries = network.storage_units[ - network.storage_units.carrier == "battery" + (network.storage_units.carrier == "battery") & + (network.storage_units.bus.isin(buses)) ] battery_distribution = ( ( @@ -921,7 +926,8 @@ def calc_storage_expansion_per_bus( ) if "H2_overground" in carriers: h2_overground = network.stores[ - network.stores.carrier == "H2_overground" + (network.stores.carrier == "H2_overground") & + (network.stores.bus.isin(buses)) ] h2_over_distribution = ( network.stores.e_nom_opt[h2_overground.index] @@ -940,7 +946,8 @@ def calc_storage_expansion_per_bus( if "H2_overground" in carriers: h2_underground = network.stores[ - network.stores.carrier == "H2_underground" + (network.stores.carrier == "H2_underground") & + (network.stores.bus.isin(buses)) ] h2_under_distribution = ( network.stores.e_nom_opt[h2_underground.index] @@ -959,7 +966,8 @@ def calc_storage_expansion_per_bus( if "rural_heat_store" in carriers: rural_heat = network.stores[ - network.stores.carrier == "rural_heat_store" + (network.stores.carrier == "rural_heat_store") & + (network.stores.bus.isin(buses)) ] rural_heat_distribution = ( network.stores.e_nom_opt[rural_heat.index] @@ -980,7 +988,8 @@ def calc_storage_expansion_per_bus( ] = rural_heat_distribution if "central_heat_store" in carriers: central_heat = network.stores[ - network.stores.carrier == "central_heat_store" + (network.stores.carrier == "central_heat_store") & + (network.stores.bus.isin(buses)) ] central_heat_distribution = ( network.stores.e_nom_opt[central_heat.index] @@ -1586,7 +1595,7 @@ def calc_ac_loading(network, timesteps): Returns ------- pandas.Series - AC line loading in MVA + AC line loading in MVA/MVA """ @@ -1609,6 +1618,36 @@ def calc_ac_loading(network, timesteps): return loading_lines / network.lines.s_nom_opt +def calc_ac_max_loading(network, timesteps): + """Calculates loading of AC-lines + + Parameters + ---------- + network : :class:`pypsa.Network + Overall container of PyPSA + timesteps : range + Defines which timesteps are considered. + + Returns + ------- + pandas.Series + AC line loading in MVA/MVA + + """ + + loading_lines = ( + network.lines_t.p0.loc[network.snapshots[timesteps]] + ).abs() + + if not network.lines_t.q0.empty: + loading_lines = np.sqrt( + loading_lines**2 + + network.lines_t.q0.loc[network.snapshots[timesteps]]** 2 + ) + + return loading_lines.max() / network.lines.s_nom_opt + + def calc_dc_loading(network, timesteps): """Calculates loading of DC-lines @@ -1624,7 +1663,7 @@ def calc_dc_loading(network, timesteps): Returns ------- pandas.Series - DC line loading in MW + DC line loading in MW/MW """ dc_links = network.links.loc[network.links.carrier == "DC", :] @@ -1649,6 +1688,41 @@ def calc_dc_loading(network, timesteps): return dc_load +def calc_dc_max_loading(network, timesteps): + """Calculates loading of DC-lines + + + Parameters + ---------- + network : :class:`pypsa.Network + Overall container of PyPSA + timesteps : range + Defines which timesteps are considered. If more than one, an + average line loading is calculated. + + Returns + ------- + pandas.Series + DC line loading in MW/MW + + """ + dc_links = network.links.loc[network.links.carrier == "DC", :] + ts = network.snapshots[timesteps] + + link_load = network.links_t.p0.loc[ts, dc_links.index].abs() + + dc_load = pd.Series(index=network.links.index, data=0.0) + dc_load.loc[dc_links.index] = ( + ( + link_load.max() / dc_links.p_nom_opt + ) + .fillna(0) + .values + ) + + return dc_load + + def plotting_colors(network): """Add color values to network.carriers @@ -2469,6 +2543,7 @@ def plot_grid( Current options: * 'line_loading': mean line loading in p.u. in selected timesteps + * 'max_line_loading': maximum line loading in p.u. in selected timesteps * 'v_nom': nominal voltage of lines * 'expansion_abs': absolute network expansion in MVA * 'expansion_rel': network expansion in p.u. of existing capacity @@ -2530,7 +2605,7 @@ def plot_grid( False, it could be assinged like this: {"H2": 50, "heat": 0.1, "battery": 10} geographical_boundaries : list, optional - Set georaphical boundaries for the plots. This parameter is overwritten + Set geographical boundaries for the plots. This parameter is overwritten when osm is used. The default is [-2.5, 16, 46.8, 58] Returns @@ -2630,6 +2705,23 @@ def plot_grid( flow[flow < 0] = -1 flow[flow > 0] = 1 + elif line_colors == "max_line_loading": + title = "Maximum line loading" + line_colors = calc_ac_max_loading(network, timesteps) + link_colors = calc_dc_max_loading(network, timesteps) + if ext_width is not False: + link_widths = link_colors.apply( + lambda x: 1 + (x / ext_width) if x != 0 else 0 + ) + line_widths = 1 + (line_colors / ext_width) + else: + link_widths = link_colors.apply(lambda x: 1 if x != 0 else 0) + line_widths = 1 + label = "max. line loading in p.u." + plot_background_grid(network, ax, geographical_boundaries, osm) + # Only active flow direction is displayed! + flow = None + elif line_colors == "v_nom": title = "Voltage levels" label = "v_nom in kV" @@ -2731,6 +2823,17 @@ def plot_grid( # Set bus colors bus_legend = False + # specify area that is plotted to select buses within this area and avoid + # very large legend entries + if osm: + x = osm["x"] + y = osm["y"] + else: + x = [geographical_boundaries[0], geographical_boundaries[1]] + y = [geographical_boundaries[2], geographical_boundaries[3]] + coords = ((x[0], y[0]), (x[0], y[1]), (x[1], y[1]), (x[0], y[1]), (x[0], y[0])) + area = Polygon(coords) + buses_area = select_buses_area(self, area).index if bus_colors == "nodal_production_balance": bus_scaling = bus_sizes bus_sizes, bus_colors = nodal_production_balance( @@ -2748,7 +2851,11 @@ def plot_grid( "battery": 10}""" ) bus_scaling = bus_sizes - bus_sizes = bus_scaling * calc_storage_expansion_per_bus(network) + bus_sizes = bus_scaling * calc_storage_expansion_per_bus( + network, + carriers=list(scaling_store_expansion.keys()), + buses=buses_area + ) for store_carrier in scaling_store_expansion.keys(): bus_sizes[ bus_sizes.index.get_level_values("carrier").str.contains( diff --git a/etrago/appl.py b/etrago/appl.py index 224e25ed..63795ddb 100644 --- a/etrago/appl.py +++ b/etrago/appl.py @@ -378,6 +378,17 @@ def run_etrago(args, json_path): * 'capacity_factor_per_gen_cntr': dict of dict of arrays Limit overall energy production country-wise for each generator by carrier. Set upper/lower limit in p.u. + * 'add_piecewise_constraint': dict(dict) + Add piecewise constraints for certain links. The dictionary should look + something like this: + { + "link_name_1": { + "breakpoints_x_pu": [0.0, 0.5, 1.0], # list(float) + "breakpoints_y_pu": [0.0, 0.4, 0.7], # list(float) + }, + } + See :func:`etrago.tools.constraints._add_piecewise_constraint` for more + information. delete_dispensable_ac_buses: bool Choose if electrical buses that are only connecting two lines should be diff --git a/etrago/execute/__init__.py b/etrago/execute/__init__.py index 5f70b80a..4fb88dae 100644 --- a/etrago/execute/__init__.py +++ b/etrago/execute/__init__.py @@ -108,6 +108,28 @@ def update_electrical_parameters(network, l_snom_pre, t_snom_pre): return l_snom_pre, t_snom_pre +def update_piecewise_link_p1(etrago): + """ + Overwrites results for links_t.p1 in case piecewise constraints were applied. + + PyPSA calculates p1 of links through p0*efficiency. In case of piecewise + constraints, the link efficiency in the links dataframe is not taken into account. + The efficiency is instead specified through the breakpoints provided in + etrago.args["extra_functionality"]["add_piecewise_constraint"]. p1 therefore needs + to be recalculated and overwritten, which is done in this function. + + Parameters + ---------- + etrago : etrago object + + """ + # overwrite p1 for piecewise links + if "add_piecewise_constraint" in etrago.args["extra_functionality"].keys(): + p1_df = pd.Series(etrago.network.model.link_out_var.get_values(), + dtype=float).unstack(0) + etrago.network.links_t.p1.loc[etrago.network.snapshots, p1_df.columns] = -p1_df + + def run_lopf(etrago, extra_functionality, method): """ Function that performs lopf with or without pyomo @@ -212,6 +234,7 @@ def run_lopf(etrago, extra_functionality, method): if etrago.network.results["Solver"][0]["Status"] != "ok": raise Exception("LOPF not solved.") + update_piecewise_link_p1(etrago) elif method["formulation"] == "linopy": status, condition = etrago.network.optimize( diff --git a/etrago/tools/constraints.py b/etrago/tools/constraints.py index 06dff351..214c6a3c 100755 --- a/etrago/tools/constraints.py +++ b/etrago/tools/constraints.py @@ -30,6 +30,8 @@ get_switchable_as_dense as get_as_dense, ) from pypsa.optimization.compat import get_var, define_constraints, linexpr +from pypsa.opf import define_nodal_balances +from pypsa.opt import LConstraint, LExpression, l_constraint import numpy as np import pandas as pd import pyomo.environ as po @@ -3657,3 +3659,166 @@ def add_chp_constraints_linopy(network, snapshots): "Link", "top_iso_fuel_line_" + i + "_" + str(snapshot), ) + + +def _add_piecewise_constraint(self, network, snapshots): + """ + Adds piecewise constraints for specified links. + + Instead of using the link's efficiency, a variable efficiency is set through this + function, that is defined through breakpoints of a piecewise linear function. For + example, if the link represents an electrolyser, the specified x-breakpoints would + represent the electrolyser's electricity demand and the y-breakpoints the resulting + hydrogen production at that respective electricity demand. + + It is assumed that the link's bus0 is the input bus and bus1 the output bus. If the + link has multiple outputs, their constraints won't be overwritten by this function. + + The piecewise constraints are parametrised through a dictionary whose keys give + the names of the links to apply a piecewise constraint to and whose values are + again dictionaries giving the piecewise function for the respective link: + + .. code-block:: python + + { + "link_name": { + "breakpoints_x_pu": [0.0, 0.5, 1.0], + "breakpoints_y_pu": [0.0, 0.4, 0.7], + }, + } + + Where: + + link_name : str + Name of link as in index of network.links of which to change constraints. + breakpoints_x_pu : list(float) + x-values of piecewise linear function, normed by p_nom of the link, so that + values are between 0 and 1. + breakpoints_y_pu : list(float) + y-values of piecewise linear function, normed by p_nom of the link, so that + values are between 0 and 1. + + The piecewise representation used here is SOS2 (special ordered set 2). + + Parameters + ---------- + self : etrago object + network : :class:`pypsa.Network` + Overall container of PyPSA + snapshots : pandas.DatetimeIndex + List of timesteps considered in the optimization + + """ + # define parameters + param_dict = self.args["extra_functionality"]["add_piecewise_constraint"] + link_names = list(param_dict.keys()) + link_p_nom = network.links.loc[link_names, "p_nom"] + out_bus = network.links.loc[link_names, "bus1"] + link_efficiency = network.links.loc[link_names, "efficiency"] + for k in link_names: + param_dict[k]["breakpoints_x"] = list( + np.array(param_dict[k]["breakpoints_x_pu"]) * link_p_nom.at[k]) + param_dict[k]["breakpoints_y"] = list( + np.array(param_dict[k]["breakpoints_y_pu"]) * link_p_nom.at[k]) + + # make new variable for link "input" + def link_p_nom_bounds(model, k, sn): + return ( + 0, link_p_nom.at[k], + ) + + network.model.link_in_var = po.Var( + link_names, network.snapshots, + bounds=link_p_nom_bounds, domain=po.Reals + ) + + # set new variable for link input equal to link_p variable + def link_input_rule(model, k, sn): + lhs = network.model.link_in_var[(k, sn)] + rhs = network.model.link_p[(k, sn)] + return lhs == rhs + + network.model.link_input_equal_link_flow = po.Constraint( + link_names, network.snapshots, rule=link_input_rule + ) + + # make new variable for link "output" + network.model.link_out_var = po.Var( + link_names, network.snapshots, + bounds=link_p_nom_bounds, domain=po.Reals + ) + + # deactivate power balance constraint for output bus + for k in param_dict.keys(): + for sn in network.snapshots: + network.model.power_balance[(out_bus.at[k], sn)].deactivate() + + # get power balances from pypsa to set up new power balance for output bus + define_nodal_balances(network, network.snapshots) + + def add_passive_branches(): + # from pypsa.opf.define_nodal_balance_constraints() + passive_branches = network.passive_branches() + for branch in passive_branches.index: + bus0 = passive_branches.at[branch, "bus0"] + bus1 = passive_branches.at[branch, "bus1"] + bt = branch[0] + bn = branch[1] + for sn in snapshots: + network._p_balance[bus0, sn].variables.append( + (-1, network.model.passive_branch_p[bt, bn, sn]) + ) + network._p_balance[bus1, sn].variables.append( + (1, network.model.passive_branch_p[bt, bn, sn]) + ) + + # add passive branches to nodal balances (is not yet done in define_nodal_balances() + # but in other pypsa function pypsa.opf.define_nodal_balance_constraints()) + add_passive_branches() + + # adapt power balance from pypsa + p_balance = {} + for k in link_names: + for sn in network.snapshots: + p_balance[(k, sn)] = network._p_balance[out_bus.at[k], sn] + # remove initial term + p_balance[(k, sn)].variables.remove( + (link_efficiency.at[k], network.model.link_p[k, sn]) + ) + # add new variable + p_balance[(k, sn)].variables.append( + (1, network.model.link_out_var[(k, sn)]) + ) + + # add new power balance to model + power_balance = { + k: LConstraint(v, "==", LExpression()) for k, v in + p_balance.items() + } + l_constraint( + network.model, + "power_balance_new", + power_balance, + link_names, + network.snapshots, + ) + # tidy up auxiliary expression + del network._p_balance + + # add piecewise constraint + def _conversion_function(model, k, sn, el_in): + ind = param_dict[k]["breakpoints_x"].index(el_in) + return param_dict[k]["breakpoints_y"][ind] + + breakpoints = {} + for k in link_names: + for sn in network.snapshots: + breakpoints[(k, sn)] = param_dict[k]["breakpoints_x"] + + network.model.piecewise_balance = po.Piecewise( + link_names, network.snapshots, + network.model.link_out_var, network.model.link_in_var, + pw_pts=breakpoints, pw_constr_type='EQ', f_rule=_conversion_function, + pw_repn="SOS2" + ) + diff --git a/etrago/tools/utilities.py b/etrago/tools/utilities.py index ccaabca1..f96ea03c 100755 --- a/etrago/tools/utilities.py +++ b/etrago/tools/utilities.py @@ -3217,20 +3217,7 @@ def find_buses_area(etrago, carrier): "not supported format supplied to 'interest_area' argument" ) - try: - buses_area = gpd.GeoDataFrame( - etrago.network.buses, geometry="geom", crs=4326 - ) - except: - buses_area = etrago.network.buses[["x", "y", "carrier"]] - buses_area["geom"] = buses_area.apply( - lambda x: Point(x["x"], x["y"]), axis=1 - ) - buses_area = gpd.GeoDataFrame( - buses_area, geometry="geom", crs=4326 - ) - - buses_area = gpd.clip(buses_area, de_areas) + buses_area = select_buses_area(etrago, de_areas) buses_area = buses_area[buses_area.carrier == carrier] else: @@ -3239,6 +3226,39 @@ def find_buses_area(etrago, carrier): return buses_area.index +def select_buses_area(etrago, area): + """ + Select buses inside given area. + + Parameters + ---------- + etrago : ETraGo object + area : shapely.Polygon + Polygon specifying the area for which to select all buses. + + Returns + ------- + geopandas.GeoDataFrame + GeoDataFrame with all buses inside given area. The dataframe is a subset of + network.buses with columns x, y, carrier and geom. + + """ + try: + buses_area = gpd.GeoDataFrame( + etrago.network.buses, geometry="geom", crs=4326 + ) + except: + buses_area = etrago.network.buses[["x", "y", "carrier"]] + buses_area["geom"] = buses_area.apply( + lambda x: Point(x["x"], x["y"]), axis=1 + ) + buses_area = gpd.GeoDataFrame( + buses_area, geometry="geom", crs=4326 + ) + + return gpd.clip(buses_area, area) + + def adjust_before_optimization(self): def check_e_initial(etrago):