From a891a1b5adf9d90eea09d38d0da723affd020789 Mon Sep 17 00:00:00 2001
From: Bryn Pickering <17178478+brynpickering@users.noreply.github.com>
Date: Thu, 18 Jan 2024 14:16:38 +0100
Subject: [PATCH] Add shadow prices

---
 src/calliope/backend/backend_model.py       | 11 +++++++++++
 src/calliope/backend/pyomo_backend_model.py | 19 +++++++++++++++++++
 2 files changed, 30 insertions(+)

diff --git a/src/calliope/backend/backend_model.py b/src/calliope/backend/backend_model.py
index 08b90bbb..10528256 100644
--- a/src/calliope/backend/backend_model.py
+++ b/src/calliope/backend/backend_model.py
@@ -647,6 +647,17 @@ def get_global_expression(
             xr.DataArray: global expression array.
         """
 
+    @abstractmethod
+    def get_shadow_prices(self, name: str) -> xr.DataArray:
+        """Extract shadow prices (a.k.a. duals) from a constraint.
+
+        Args:
+            name (str): Name of constraint for which you're seeking duals.
+
+        Returns:
+            xr.DataArray: duals array.
+        """
+
     @abstractmethod
     def update_parameter(
         self, name: str, new_values: Union[xr.DataArray, SupportsFloat]
diff --git a/src/calliope/backend/pyomo_backend_model.py b/src/calliope/backend/pyomo_backend_model.py
index 217eb7b8..d7ab0789 100644
--- a/src/calliope/backend/pyomo_backend_model.py
+++ b/src/calliope/backend/pyomo_backend_model.py
@@ -59,6 +59,8 @@ def __init__(self, inputs: xr.Dataset, **kwargs):
         self._instance.constraints = pmo.constraint_dict()
         self._instance.objectives = pmo.objective_dict()
 
+        self._instance.dual = pmo.suffix(direction=pmo.suffix.IMPORT)
+
         self._add_all_inputs_as_parameters()
 
     def add_parameter(
@@ -247,6 +249,14 @@ def get_global_expression(
         else:
             return global_expression
 
+    def get_shadow_prices(self, name: str) -> xr.DataArray:
+        constraint = self.get_constraint(name, as_backend_objs=True)
+        return self._apply_func(
+            self._duals_from_pyomo_constraint,
+            constraint,
+            dual_getter=self._instance.dual,
+        )
+
     def solve(
         self,
         solver: str,
@@ -752,6 +762,15 @@ def _from_pyomo_expr(val: pmo.expression, *, eval_body: bool = False) -> Any:
             else:
                 return val.to_string()
 
+    @staticmethod
+    def _duals_from_pyomo_constraint(
+        val: pmo.constraint, *, dual_getter: pmo.suffix
+    ) -> float:
+        if pd.isnull(val):
+            return np.nan
+        else:
+            return dual_getter.get(val)
+
     @contextmanager
     def _datetime_as_string(self, data: Union[xr.DataArray, xr.Dataset]) -> Iterator:
         """Context manager to temporarily convert np.dtype("datetime64[ns]") coordinates (e.g. timesteps) to strings with a resolution of minutes.