From 3754643e8192c57452c84d7a8a3639abd17a49a3 Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Fri, 6 Aug 2021 12:46:03 +0200
Subject: [PATCH] add co2 sequestration potential global constraint

---
 scripts/prepare_sector_network.py |  9 +--------
 scripts/solve_network.py          | 19 +++++++++++++++++++
 2 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index ca956e44..659fbf37 100644
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -45,7 +45,7 @@ def define_spatial(nodes):
         spatial.co2.vents = nodes + " co2 vent"
     else:
         spatial.co2.nodes = ["co2 stored"]
-        spatial.co2.locations = "EU"
+        spatial.co2.locations = ["EU"]
         spatial.co2.vents = ["co2 vent"]
 
     spatial.co2.df = pd.DataFrame(vars(spatial.co2), index=nodes)
@@ -379,13 +379,6 @@ def add_co2_tracking(n, options):
         bus=spatial.co2.nodes
     )
 
-    # TODO if nodally resolved total allowed sequestration needs to
-    # be an extra_functionality constraint
-    # (best to implement it this way for either case)
-    # let e_nom_max represent geological potential
-    # don't forget to log duals of extra functionality!
-    # options['co2_sequestration_potential'] * 1e6
-
     if options['co2_vent']:
 
         n.madd("Link",
diff --git a/scripts/solve_network.py b/scripts/solve_network.py
index 8c0313f1..3796bc4d 100644
--- a/scripts/solve_network.py
+++ b/scripts/solve_network.py
@@ -3,6 +3,7 @@
 import pypsa
 
 import numpy as np
+import pandas as pd
 
 from pypsa.linopt import get_var, linexpr, define_constraints
 
@@ -150,9 +151,27 @@ def add_chp_constraints(n):
         define_constraints(n, lhs, "<=", 0, 'chplink', 'backpressure')
 
 
+def add_co2_sequestration_limit(n, sns):
+    
+    co2_stores = n.stores.loc[n.stores.carrier=='co2 stored'].index
+
+    if co2_stores.empty or ('Store', 'e') not in n.variables.index:
+        return
+    
+    vars_final_co2_stored = get_var(n, 'Store', 'e').loc[sns[-1], co2_stores]
+    
+    lhs = linexpr((1, vars_final_co2_stored)).sum()
+    rhs = n.config["sector"].get("co2_sequestration_potential", 200) * 1e6
+    
+    name = 'co2_sequestration_limit'
+    define_constraints(n, lhs, "<=", rhs, 'GlobalConstraint',
+                       'mu', axes=pd.Index([name]), spec=name)
+
+
 def extra_functionality(n, snapshots):
     add_chp_constraints(n)
     add_battery_constraints(n)
+    add_co2_sequestration_limit(n, snapshots)
 
 
 def solve_network(n, config, opts='', **kwargs):