From cc162a9e028fb7a2bac5289e27b90ab57e46f10b Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Mon, 31 Jul 2023 17:09:59 +0200
Subject: [PATCH] option for losses on bidirectional links via link splitting

---
 config/config.default.yaml        |  5 ++++
 scripts/prepare_sector_network.py | 40 +++++++++++++++++++++++++++++++
 scripts/solve_network.py          | 22 +++++++++++++++++
 3 files changed, 67 insertions(+)

diff --git a/config/config.default.yaml b/config/config.default.yaml
index b162b75d..4413b8f5 100644
--- a/config/config.default.yaml
+++ b/config/config.default.yaml
@@ -478,6 +478,11 @@ sector:
   electricity_distribution_grid: true
   electricity_distribution_grid_cost_factor: 1.0
   electricity_grid_connection: true
+  transmission_losses:
+    # per 1000 km
+    DC: 0
+    H2 pipeline: 0
+    gas pipeline: 0
   H2_network: true
   gas_network: false
   H2_retrofit: false
diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index 11406bff..8719c281 100644
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -3280,6 +3280,34 @@ def set_temporal_aggregation(n, opts, solver_name):
     return n
 
 
+def lossy_bidirectional_links(n, carrier, losses_per_thousand_km=0.0):
+    "Split bidirectional links into two unidirectional links to include transmission losses."
+
+    carrier_i = n.links.query("carrier == @carrier").index
+
+    if not losses_per_thousand_km or carrier_i.empty:
+        return
+
+    logger.info(
+        f"Specified losses for {carrier} transmission. Splitting bidirectional links."
+    )
+
+    carrier_i = n.links.query("carrier == @carrier").index
+    n.links.loc[carrier_i, "p_min_pu"] = 0
+    n.links["reversed"] = False
+    n.links.loc[carrier_i, "efficiency"] = (
+        1 - n.links.loc[carrier_i, "length"] * losses_per_thousand_km / 1e3
+    )
+    rev_links = (
+        n.links.loc[carrier_i].copy().rename({"bus0": "bus1", "bus1": "bus0"}, axis=1)
+    )
+    rev_links.capital_cost = 0
+    rev_links.reversed = True
+    rev_links.index = rev_links.index.map(lambda x: x + "-reversed")
+
+    n.links = pd.concat([n.links, rev_links], sort=False)
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
@@ -3446,6 +3474,18 @@ if __name__ == "__main__":
     if options["electricity_grid_connection"]:
         add_electricity_grid_connection(n, costs)
 
+    for k, v in options["transmission_losses"].items():
+        lossy_bidirectional_links(n, k, v)
+
+    # Workaround: Remove lines with conflicting (and unrealistic) properties
+    # cf. https://github.com/PyPSA/pypsa-eur/issues/444
+    if snakemake.config["solving"]["options"]["transmission_losses"]:
+        idx = n.lines.query("num_parallel == 0").index
+        logger.info(
+            f"Removing {len(idx)} line(s) with properties conflicting with transmission losses functionality."
+        )
+        n.mremove("Line", idx)
+
     first_year_myopic = (snakemake.params.foresight == "myopic") and (
         snakemake.params.planning_horizons[0] == investment_year
     )
diff --git a/scripts/solve_network.py b/scripts/solve_network.py
index 836544b4..a68ca074 100644
--- a/scripts/solve_network.py
+++ b/scripts/solve_network.py
@@ -494,6 +494,27 @@ def add_battery_constraints(n):
     n.model.add_constraints(lhs == 0, name="Link-charger_ratio")
 
 
+def add_lossy_bidirectional_link_constraints(n):
+    if not n.links.p_nom_extendable.any() or not "reversed" in n.links.columns:
+        return
+
+    carriers = n.links.loc[n.links.reversed, "carrier"].unique()
+
+    backward_i = n.links.query(
+        "carrier in @carriers and reversed and p_nom_extendable"
+    ).index
+    forward_i = n.links.query(
+        "carrier in @carriers and ~reversed and p_nom_extendable"
+    ).index
+
+    assert len(forward_i) == len(backward_i)
+
+    lhs = n.model["Link-p_nom"].loc[backward_i]
+    rhs = n.model["Link-p_nom"].loc[forward_i]
+
+    n.model.add_constraints(lhs == rhs, name="Link-bidirectional_sync")
+
+
 def add_chp_constraints(n):
     electric = (
         n.links.index.str.contains("urban central")
@@ -593,6 +614,7 @@ def extra_functionality(n, snapshots):
         if "EQ" in o:
             add_EQ_constraints(n, o)
     add_battery_constraints(n)
+    add_lossy_bidirectional_link_constraints(n)
     add_pipe_retrofit_constraint(n)