Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Features/#665 interface pre market market optimization #722

Merged
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions etrago/appl.py
Original file line number Diff line number Diff line change
@@ -54,7 +54,12 @@
"type": "market_grid", # type of optimization, 'lopf' or 'market_grid'
"n_iter": 1, # abort criterion of iterative optimization, 'n_iter' or 'threshold'
"pyomo": True, # set if pyomo is used for model building
"formulation": "pyomo",
"market_zones": "status_quo", # only used if type='market_grid'
"rolling_horizon": { # Define parameter of market optimization
"planning_horizon": 72, # number of snapshots in each optimization
"overlap": 24, # number of overlapping hours
},
},
"pf_post_lopf": {
"active": False, # choose if perform a pf after lopf
45 changes: 37 additions & 8 deletions etrago/cluster/electrical.py
Original file line number Diff line number Diff line change
@@ -95,7 +95,9 @@ def leader(x):
return leader


def adjust_no_electric_network(etrago, busmap, cluster_met):
def adjust_no_electric_network(
etrago, busmap, cluster_met, apply_on="grid_model"
):
"""
Adjusts the non-electric network based on the electrical network
(esp. eHV network), adds the gas buses to the busmap, and creates the
@@ -139,7 +141,16 @@ def find_de_closest(network, bus_ne):

return new_ehv_bus

network = etrago.network.copy()
if apply_on == "grid_model":
network = etrago.network.copy()
elif apply_on == "market_model":
network = etrago.network_tsa.copy()
else:
logger.warning(
"""Parameter apply_on must be either 'grid_model' or 'market_model'
"""
)

# network2 is supposed to contain all the not electrical or gas buses
# and links
network2 = network.copy(with_time=False)
@@ -457,7 +468,7 @@ def ehv_clustering(self):
logger.info("Network clustered to EHV-grid")


def select_elec_network(etrago):
def select_elec_network(etrago, apply_on="grid_model"):
"""
Selects the electric network based on the clustering settings specified
in the Etrago object.
@@ -475,7 +486,15 @@ def select_elec_network(etrago):
n_clusters : int
number of clusters used in the clustering process.
"""
elec_network = etrago.network.copy()
if apply_on == "grid_model":
elec_network = etrago.network.copy()
elif apply_on == "market_model":
elec_network = etrago.network_tsa.copy()
else:
logger.warning(
"""Parameter apply_on must be either 'grid_model' or 'market_model'
"""
)
settings = etrago.args["network_clustering"]
if settings["cluster_foreign_AC"]:
elec_network.buses = elec_network.buses[
@@ -626,7 +645,7 @@ def unify_foreign_buses(etrago):
return busmap_foreign


def preprocessing(etrago):
def preprocessing(etrago, apply_on="grid_model"):
"""
Preprocesses an Etrago object to prepare it for network clustering.

@@ -647,7 +666,16 @@ def preprocessing(etrago):
The Series object with the foreign bus mapping data.
"""

network = etrago.network
if apply_on == "grid_model":
network = etrago.network
elif apply_on == "market_model":
network = etrago.network_tsa
else:
logger.warning(
"""Parameter apply_on must be either 'grid_model' or 'market_model'
"""
)

settings = etrago.args["network_clustering"]

# problem our lines have no v_nom. this is implicitly defined by the
@@ -736,7 +764,7 @@ def preprocessing(etrago):
else:
busmap_foreign = pd.Series(name="foreign", dtype=str)

network_elec, n_clusters = select_elec_network(etrago)
network_elec, n_clusters = select_elec_network(etrago, apply_on=apply_on)

if settings["method"] == "kmedoids-dijkstra":
lines_col = network_elec.lines.columns
@@ -777,6 +805,7 @@ def postprocessing(
medoid_idx=None,
aggregate_generators_carriers=None,
aggregate_links=True,
apply_on="grid_model",
):
"""
Postprocessing function for network clustering.
@@ -846,7 +875,7 @@ def postprocessing(
)

network, busmap = adjust_no_electric_network(
etrago, busmap, cluster_met=method
etrago, busmap, cluster_met=method, apply_on=apply_on
)

# merge busmap for foreign buses with the German buses
2 changes: 1 addition & 1 deletion etrago/cluster/snapshot.py
Original file line number Diff line number Diff line change
@@ -811,7 +811,7 @@ def skip_snapshots(self):
if (
self.args["temporal_disaggregation"]["active"]
and not self.args["snapshot_clustering"]["active"]
):
) or self.args["method"]["type"] == "market_grid":
self.network_tsa = self.network.copy()

n_skip = self.args["skip_snapshots"]
26 changes: 22 additions & 4 deletions etrago/execute/grid_optimization.py
Original file line number Diff line number Diff line change
@@ -225,7 +225,11 @@ def add_redispatch_generators(self):
p_max_pu_all.loc[:, gens_redispatch].mul(
self.network.generators.loc[gens_redispatch, "p_nom"]
)
- (self.market_model.generators_t.p.loc[:, gens_redispatch])
- (
self.market_model.generators_t.p.loc[
self.network.snapshots, gens_redispatch
]
)
)
.clip(lower=0.0)
.values
@@ -251,7 +255,11 @@ def add_redispatch_generators(self):
self.network.links_t.p_max_pu.loc[:, links_redispatch + " ramp_up"] = (
(
self.network.links.loc[links_redispatch, "p_nom"]
- (self.market_model.links_t.p0.loc[:, links_redispatch])
- (
self.market_model.links_t.p0.loc[
self.network.snapshots, links_redispatch
]
)
)
.clip(lower=0.0)
.values
@@ -277,7 +285,13 @@ def add_redispatch_generators(self):
# (disaggregated) generators in the market model
self.network.generators_t.p_min_pu.loc[
:, gens_redispatch + " ramp_down"
] = (-(self.market_model.generators_t.p.loc[:, gens_redispatch])).values
] = (
-(
self.market_model.generators_t.p.loc[
self.network.snapshots, gens_redispatch
].clip(lower=0.0)
)
).values

# Add ramp down links to the network for the grid optimization
# Marginal cost are currently only the management fee of 4 EUR/MWh,
@@ -299,7 +313,11 @@ def add_redispatch_generators(self):
# Ramp down can be at maximum as high as the feed-in of the
# (disaggregated) links in the market model
self.network.links_t.p_min_pu.loc[:, links_redispatch + " ramp_down"] = (
-(self.market_model.links_t.p0.loc[:, links_redispatch])
-(
self.market_model.links_t.p0.loc[
self.network.snapshots, links_redispatch
].clip(lower=0.0)
)
).values

# Check if the network contains any problems
114 changes: 108 additions & 6 deletions etrago/execute/market_optimization.py
Original file line number Diff line number Diff line change
@@ -48,6 +48,7 @@ def market_optimization(self):
logger.info("Start building pre market model")
build_market_model(self)
logger.info("Start solving pre market model")

self.pre_market_model.lopf(
solver_name=self.args["solver"],
solver_options=self.args["solver_options"],
@@ -60,13 +61,23 @@ def market_optimization(self):
build_shortterm_market_model(self)
logger.info("Start solving short-term UC market model")

self.market_model.optimize.optimize_with_rolling_horizon(
# Set 'linopy' as formulation to make sure that constraints are added
method_args = self.args["method"]["formulation"]
self.args["method"]["formulation"] = "linopy"

optimize_with_rolling_horizon(
self.market_model,
self.pre_market_model,
snapshots=None,
horizon=168,
overlap=144,
horizon=self.args["method"]["rolling_horizon"]["planning_horizon"],
overlap=self.args["method"]["rolling_horizon"]["overlap"],
solver_name=self.args["solver"],
extra_functionality=Constraints(self.args, False).functionality,
)

# Reset formulation to previous setting of args
self.args["method"]["formulation"] = method_args

# Export results of market model
if self.args["csv_export"]:
path = self.args["csv_export"]
@@ -75,6 +86,95 @@ def market_optimization(self):
self.market_model.export_to_csv_folder(path + "/market")


def optimize_with_rolling_horizon(
n, pre_market, snapshots=None, horizon=2, overlap=0, **kwargs
):
"""
Optimizes the network in a rolling horizon fashion.

Parameters
----------
n : pypsa.Network
snapshots : list-like
Set of snapshots to consider in the optimization. The default is None.
horizon : int
Number of snapshots to consider in each iteration. Defaults to 100.
overlap : int
Number of snapshots to overlap between two iterations. Defaults to 0.
**kwargs:
Keyword argument used by `linopy.Model.solve`, such as `solver_name`,

Returns
-------
None
"""
if snapshots is None:
snapshots = n.snapshots

if horizon <= overlap:
raise ValueError("overlap must be smaller than horizon")

starting_points = range(0, len(snapshots), horizon - overlap)
for i, start in enumerate(starting_points):
end = min(len(snapshots), start + horizon)
sns = snapshots[start:end]
logger.info(
f"""Optimizing network for snapshot horizon
[{sns[0]}:{sns[-1]}] ({i+1}/{len(starting_points)})."""
)

if i:
if not n.stores.empty:
n.stores.e_initial = n.stores_t.e.loc[snapshots[start - 1]]

# Select seasonal stores
seasonal_stores = n.stores.index[
n.stores.carrier.isin(
["central_heat_store", "H2_overground", "CH4"]
)
]

# Set e_initial from pre_market model for seasonal stores
n.stores.e_initial[seasonal_stores] = (
pre_market.stores_t.e.loc[
snapshots[start - 1], seasonal_stores
]
)

# Set e at the end of the horizon
# by setting e_max_pu and e_min_pu
n.stores_t.e_max_pu.loc[
snapshots[end - 1], seasonal_stores
] = pre_market.stores_t.e.loc[
snapshots[end - 1], seasonal_stores
].div(
pre_market.stores.e_nom_opt[seasonal_stores]
)
n.stores_t.e_min_pu.loc[
snapshots[end - 1], seasonal_stores
] = pre_market.stores_t.e.loc[
snapshots[end - 1], seasonal_stores
].div(
pre_market.stores.e_nom_opt[seasonal_stores]
)
n.stores_t.e_min_pu.fillna(0.0, inplace=True)
n.stores_t.e_max_pu.fillna(1.0, inplace=True)

if not n.storage_units.empty:
n.storage_units.state_of_charge_initial = (
n.storage_units_t.state_of_charge.loc[snapshots[start - 1]]
)

status, condition = n.optimize(sns, **kwargs)

if status != "ok":
logger.warning(
f"""Optimization failed with status {status}
and condition {condition}"""
)
return n


def build_market_model(self):
"""Builds market model based on imported network from eTraGo

@@ -91,8 +191,9 @@ def build_market_model(self):
"""

# use existing preprocessing to get only the electricity system

net, weight, n_clusters, busmap_foreign = preprocessing(self)
net, weight, n_clusters, busmap_foreign = preprocessing(
self, apply_on="market_model"
)

# Define market regions based on settings.
# Currently the only option is 'status_quo' which means that the current
@@ -139,6 +240,7 @@ def build_market_model(self):
medoid_idx,
aggregate_generators_carriers=[],
aggregate_links=False,
apply_on="market_model",
)

self.update_busmap(busmap)
@@ -165,7 +267,7 @@ def build_market_model(self):
)
# net.buses.loc[net.buses.carrier == 'AC', 'carrier'] = "DC"

net.generators_t.p_max_pu = self.network.generators_t.p_max_pu
net.generators_t.p_max_pu = self.network_tsa.generators_t.p_max_pu

self.pre_market_model = net

175 changes: 170 additions & 5 deletions etrago/tools/constraints.py
Original file line number Diff line number Diff line change
@@ -29,7 +29,7 @@
expand_series,
get_switchable_as_dense as get_as_dense,
)
from pypsa.linopt import define_constraints, define_variables, get_var, linexpr
from pypsa.optimization.compat import get_var, define_constraints, linexpr
import numpy as np
import pandas as pd
import pyomo.environ as po
@@ -1327,7 +1327,84 @@ def read_max_gas_generation(self):

return arg

def add_ch4_constraints_linopy(self, network, snapshots):
"""
Add CH4 constraints for optimization with linopy
Functionality that limits the dispatch of CH4 generators. In
Germany, there is one limitation specific for biogas and one
limitation specific for natural gas (natural gas only in eGon2035).
Abroad, each generator has its own limitation contains in the
column e_nom_max.
Parameters
----------
network : :class:`pypsa.Network`
Overall container of PyPSA
snapshots : pandas.DatetimeIndex
List of timesteps considered in the optimization
Returns
-------
None.
"""
scn_name = self.args["scn_name"]
n_snapshots = self.args["end_snapshot"] - self.args["start_snapshot"] + 1

# Add constraint for Germany
arg = read_max_gas_generation(self)
gas_carrier = arg.keys()

carrier_names = {
"eGon2035": {"CH4": "CH4_NG", "biogas": "CH4_biogas"},
"eGon2035_lowflex": {"CH4": "CH4_NG", "biogas": "CH4_biogas"},
"eGon100RE": {"biogas": "CH4"},
}
for c in gas_carrier:
gens = network.generators.index[
(network.generators.carrier == carrier_names[scn_name][c])
& (
network.generators.bus.astype(str).isin(
network.buses.index[network.buses.country == "DE"]
)
)
]
if not gens.empty:
factor = arg[c]
generation = (
get_var(network, "Generator", "p")
.loc[snapshots, gens] *
network.snapshot_weightings.generators
)
define_constraints(
network, generation, "<=", factor * (n_snapshots / 8760),
"Genertor", "max_flh_DE_" + c
)

# Add contraints for neigbouring countries
gen_abroad = network.generators[
(network.generators.carrier == "CH4")
& (
network.generators.bus.astype(str).isin(
network.buses.index[network.buses.country != "DE"]
)
)
& (network.generators.e_nom_max != np.inf)
]
for g in gen_abroad.index:
factor = network.generators.e_nom_max[g]

generation = (
get_var(network, "Generator", "p")
.loc[snapshots, g] *
network.snapshot_weightings.generators
)
define_constraints(
network, generation, "<=", factor * (n_snapshots / 8760),
"Genertor", "max_flh_abroad_" + str(g).replace(" ", "_")
)


def add_ch4_constraints(self, network, snapshots):
"""
Add CH4 constraints for optimization with pyomo
@@ -2751,20 +2828,25 @@ def functionality(self, network, snapshots):
"""
if "CH4" in network.buses.carrier.values:
if self.args["method"]["pyomo"]:
if self.args["method"]["formulation"] == "pyomo":
add_chp_constraints(network, snapshots)
if (self.args["scn_name"] != "status2019") & (
len(network.snapshots) > 1500
):
add_ch4_constraints(self, network, snapshots)
elif self.args["method"]["formulation"] == "linopy":
if (self.args["scn_name"] != "status2019") & (
len(network.snapshots) > 1500
):
add_ch4_constraints_linopy(self, network, snapshots)
add_chp_constraints_linopy(network, snapshots)
else:
add_chp_constraints_nmp(network)
if self.args["scn_name"] != "status2019":
add_ch4_constraints_nmp(self, network, snapshots)

for constraint in self.args["extra_functionality"].keys():
try:
type(network.model)
if self.args["method"]["formulation"] == "pyomo":
try:
eval("_" + constraint + "(self, network, snapshots)")
logger.info(
@@ -2776,7 +2858,13 @@ def functionality(self, network, snapshots):
+ ". New constraints can be defined in"
+ " etrago/tools/constraint.py."
)
except:
elif self.args["method"]["formulation"] == "linopy":
logger.warning(
"Constraint {} not defined for linopy formulation".format(constraint)
+ ". New constraints can be defined in"
+ " etrago/tools/constraint.py."
)
else:
try:
eval("_" + constraint + "_nmp(self, network, snapshots)")
logger.info(
@@ -3047,3 +3135,80 @@ def top_iso_fuel_line(model, snapshot):
"top_iso_fuel_line_" + str(i),
Constraint(list(snapshots), rule=top_iso_fuel_line),
)

def add_chp_constraints_linopy(network, snapshots):
"""
Limits the dispatch of combined heat and power links based on
T.Brown et. al : Synergies of sector coupling and transmission
reinforcement in a cost-optimised, highly renewable European energy system,
2018
Parameters
----------
network : pypsa.Network
Network container
snapshots : pandas.DataFrame
Timesteps to optimize
Returns
-------
None.
"""

# backpressure limit
c_m = 0.75

# marginal loss for each additional generation of heat
c_v = 0.15
electric_bool = network.links.carrier == "central_gas_CHP"
heat_bool = network.links.carrier == "central_gas_CHP_heat"

electric = network.links.index[electric_bool]
heat = network.links.index[heat_bool]

network.links.loc[heat, "efficiency"] = (
network.links.loc[electric, "efficiency"] / c_v
).values.mean()

ch4_nodes_with_chp = network.buses.loc[
network.links.loc[electric, "bus0"].values
].index.unique()

for i in ch4_nodes_with_chp:
elec_chp = network.links[
(network.links.carrier == "central_gas_CHP")
& (network.links.bus0 == i)
].index

heat_chp = network.links[
(network.links.carrier == "central_gas_CHP_heat")
& (network.links.bus0 == i)
].index

dispatch_heat = (
c_m * get_var(network, "Link", "p").loc[snapshots, heat_chp] *
network.links.loc[heat_chp, "efficiency"]).sum()
dispatch_elec = (
get_var(network, "Link", "p").loc[snapshots, elec_chp] *
network.links.loc[elec_chp, "efficiency"]).sum()

define_constraints(
network,
(dispatch_heat - dispatch_elec),
"<=",
0,
"Link", "backpressure_" + i
)

define_constraints(
network,
get_var(network, "Link", "p").loc[snapshots, heat_chp].sum() +
get_var(network, "Link", "p").loc[snapshots, elec_chp].sum(),
"<=",
network.links[
(network.links.carrier == "central_gas_CHP")
& (network.links.bus0 == i)
].p_nom.sum(),
"Link", "top_iso_fuel_line_" + i
)