Add shadow prices (duals) (#540)

src/calliope/backend/backend_model.py

@@ -564,6 +564,7 @@ def __init__(self, inputs: xr.Dataset, instance: T, **kwargs) -> None:
         """
         super().__init__(inputs, **kwargs)
         self._instance = instance
+        self.shadow_prices: ShadowPrices
 
     @abstractmethod
     def get_parameter(self, name: str, as_backend_objs: bool = True) -> xr.DataArray:
@@ -843,3 +844,34 @@ def _rebuild_reference(self, reference: str) -> None:
         obj_type = self._dataset[reference].attrs["obj_type"]
         self.delete_component(reference, obj_type)
         getattr(self, "add_" + obj_type.removesuffix("s"))(name=reference)
+
+
+class ShadowPrices:
+    """Object containing methods to interact with the backend object "shadow prices" tracker, which can be used to access duals for constraints.
+
+    To keep memory overhead low. Shadow price tracking is deactivated by default.
+    """
+
+    @abstractmethod
+    def get(self, name) -> 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 activate(self):
+        "Activate shadow price tracking."
+
+    @abstractmethod
+    def deactivate(self):
+        "Deactivate shadow price tracking."
+
+    @property
+    @abstractmethod
+    def is_active(self) -> bool:
+        "Check whether shadow price tracking is active or not"

src/calliope/backend/pyomo_backend_model.py

@@ -64,6 +64,9 @@ def __init__(self, inputs: xr.Dataset, **kwargs) -> None:
         self._instance.constraints = pmo.constraint_dict()
         self._instance.objectives = pmo.objective_dict()
 
+        self._instance.dual = pmo.suffix(direction=pmo.suffix.IMPORT)
+        self.shadow_prices = PyomoShadowPrices(self._instance.dual, self)
+
         self._add_all_inputs_as_parameters()
 
     def add_parameter(
@@ -261,6 +264,11 @@ def _solve(
         warmstart: bool = False,
         **solve_config,
     ) -> xr.Dataset:
+        if solver == "cbc" and self.shadow_prices.is_active:
+            model_warn(
+                "Switching off shadow price tracker as constraint duals cannot be accessed from the CBC solver"
+            )
+            self.shadow_prices.deactivate()
         opt = SolverFactory(solver, solver_io=solver_io)
 
         if solver_options:
@@ -272,6 +280,7 @@ def _solve(
             solve_kwargs.update({"symbolic_solver_labels": True, "keepfiles": True})
             os.makedirs(save_logs, exist_ok=True)
             TempfileManager.tempdir = save_logs  # Sets log output dir
+
         if warmstart and solver in ["glpk", "cbc"]:
             model_warn(
                 "The chosen solver, {}, does not support warmstart, which may "
@@ -862,3 +871,35 @@ def __init__(self, **kwds):
 
     def getname(self, *args, **kwargs):
         return self._update_name(pmo.constraint.getname(self, *args, **kwargs))
+
+
+class PyomoShadowPrices(backend_model.ShadowPrices):
+    def __init__(self, dual_obj: pmo.suffix, backend_obj: PyomoBackendModel):
+        self._dual_obj = dual_obj
+        self._backend_obj = backend_obj
+        self.deactivate()
+
+    def get(self, name: str) -> xr.DataArray:
+        constraint = self._backend_obj.get_constraint(name, as_backend_objs=True)
+        return self._backend_obj._apply_func(
+            self._duals_from_pyomo_constraint, constraint, dual_getter=self._dual_obj
+        )
+
+    def activate(self):
+        self._dual_obj.activate()
+
+    def deactivate(self):
+        self._dual_obj.deactivate()
+
+    @property
+    def is_active(self) -> bool:
+        return self._dual_obj.active
+
+    @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)

tests/test_backend_pyomo.py

@@ -2470,3 +2470,58 @@ def test_unfix_variable_where(self, simple_supply):
         simple_supply.backend.unfix_variable("flow_cap")  # reset
         assert fixed.sel(techs="test_demand_elec").all()
         assert not fixed.where(where).all()
+
+
+class TestShadowPrices:
+    @pytest.fixture(scope="function")
+    def simple_supply(self):
+        m = build_model({}, "simple_supply,two_hours,investment_costs")
+        m.build()
+        return m
+
+    @pytest.fixture(scope="function")
+    def supply_milp(self):
+        m = build_model({}, "supply_milp,two_hours,investment_costs")
+        m.build()
+        return m
+
+    def test_default_to_deactivated(self, simple_supply):
+        assert not simple_supply.backend.shadow_prices.is_active
+
+    def test_activate(self, simple_supply):
+        simple_supply.backend.shadow_prices.activate()
+        assert simple_supply.backend.shadow_prices.is_active
+
+    def test_deactivate(self, simple_supply):
+        simple_supply.backend.shadow_prices.activate()
+        simple_supply.backend.shadow_prices.deactivate()
+        assert not simple_supply.backend.shadow_prices.is_active
+
+    def test_get_shadow_price(self, simple_supply):
+        simple_supply.backend.shadow_prices.activate()
+        simple_supply.solve(solver="glpk")
+        shadow_prices = simple_supply.backend.shadow_prices.get("system_balance")
+        assert shadow_prices.notnull().all()
+
+    def test_get_shadow_price_some_nan(self, simple_supply):
+        simple_supply.backend.shadow_prices.activate()
+        simple_supply.solve(solver="glpk")
+        shadow_prices = simple_supply.backend.shadow_prices.get("balance_demand")
+        assert shadow_prices.notnull().any()
+        assert shadow_prices.isnull().any()
+
+    def test_get_shadow_price_empty_milp(self, supply_milp):
+        supply_milp.backend.shadow_prices.activate()
+        supply_milp.solve(solver="glpk")
+        shadow_prices = supply_milp.backend.shadow_prices.get("system_balance")
+        assert shadow_prices.isnull().all()
+
+    def test_shadow_prices_deactivated_with_cbc(self, simple_supply):
+        simple_supply.backend.shadow_prices.activate()
+        with pytest.warns(exceptions.ModelWarning) as warning:
+            simple_supply.solve(solver="cbc")
+
+        assert check_error_or_warning(warning, "Switching off shadow price tracker")
+        assert not simple_supply.backend.shadow_prices.is_active
+        shadow_prices = simple_supply.backend.shadow_prices.get("system_balance")
+        assert shadow_prices.isnull().all()