Shadow price configuration in YAML (#587)

* Shadow price configuration in YAML, allowing shadow prices to be tracked from the CLI.

---------

Co-authored-by: Bryn Pickering <[email protected]>

CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ### User-facing changes
 
+|new| Allow extracting shadow prices into results by listing constraints in `config.solve.shadow_prices`, e.g. `config.solve.shadow_prices: ["system_balance"]`  Shadow prices will be added as variables to the model results as `shadow_price_{constraintname}`, e.g. `shadow_price_system_balance`.
+
 |new| Model stores key timestamps as attributes:
 * `timestamp_model_creation`: at the start of `Model.__init__()`
 * `timestamp_build_started`: at the start of `Model.build()`

docs/advanced/shadow_prices.md

@@ -3,7 +3,38 @@
 In a linear problem, you can obtain the [shadow prices](https://en.wikipedia.org/wiki/Shadow_price) (dual variables) for each constraint from the Pyomo backend.
 This can prove particularly useful if you are linking Calliope with other models or want to dig deeper into the economic impacts of your designed energy system.
 
-You can access shadow prices for any constraint by running once you have an optimal solution by running `model.backend.shadow_prices.get("constraint_name")`, which returns an [xarray.DataArray][].
+You can access shadow prices by specifying a list of constraints for which shadow prices should be returned in the results, in your `solve` configuration:
+
+```yaml
+config:
+    solve:
+        shadow_prices: ["system_balance", ...]
+```
+
+A list of the available constraint names can be found under "Subject to" in [our base math documentation page][base-math].
+If you [define any of your own math constraints](../user_defined_math/components.md#constraints), you can also reference those by name in the list.
+
+!!! note
+
+    * Not all solvers provide access to shadow prices.
+    For instance, we know that it is possible with Gurobi and GLPK, but not with CBC.
+    Since we cannot test all Pyomo-compatible solvers, you may run into issues depending on the solver you use.
+    * You cannot access shadow prices if you have integer/binary variables in your model.
+    If you try to do so, you will receive a warning, and shadow price tracking will remain disabled.
+    * You can check the status of shadow price tracking with `model.backend.shadow_prices.is_active`.
+
+## Shadow prices when using the command-line tool
+
+By specifying constraints in the YAML configuration (see above), shadow price tracking will be activated and the shadow prices of those constraints you have listed will be available in the results dataset, prefixed with `shadow_price_`.
+For instance, listing `system_balance` in the configuration will lead to `shadow_price_system_balance` being available in the optimisation results that are saved to file on calling [`calliope run ...`](../running.md#running-with-the-command-line-tool).
+
+## Shadow prices when running in Python
+
+When running in Python, you can additionally turn on shadow price tracking by running `model.backend.shadow_prices.activate()` after `model.build()`.
+By doing that, or by having added at least one valid constraint in `config.solve.shadow_prices` (see above), shadow price tracking will be enabled.
+
+Then, you can access shadow prices for any constraint once you have an optimal solution by running `model.backend.shadow_prices.get("constraint_name")`, which returns an [xarray.DataArray][].
+As with [running in the command-line tool](#shadow-prices-when-using-the-command-line-tool), having a list of shadow prices to track in the solve configuration will lead to shadow prices being available automatically in the model [results dataset][calliope.Model.results].
 
 !!! example
 
@@ -17,13 +48,14 @@ You can access shadow prices for any constraint by running once you have an opti
     ```
 
     1. With the Pyomo backend interface, tracking shadow prices can be memory-intensive.
-    Therefore, you must actively activate tracking before solving your model.
+    Therefore, you must manually activate tracking before solving your model.
 
-!!! note
+    Or:
 
-    * Not all solvers provide an API to access duals.
-    For instance, we know it is not possible with the CBC solver.
-    * You cannot access shadow prices if you have integer/binary variables in your model.
-    If you try to do so, you will receive a None/NaN-filled array.
-    * You can check the status of shadow price tracking with `model.backend.shadow_prices.is_active`.
-    * Currently, the only way to access shadow prices is by running Calliope in Python, and not through the command-line interface.
+    ```python
+    model = calliope.examples.national_scale()
+    model.build()
+    model.solve(shadow_prices=["system_balance"])
+
+    balance_price = model.results.shadow_price_system_balance.to_series()
+    ```

src/calliope/backend/backend_model.py

@@ -19,6 +19,7 @@
     Any,
     Callable,
     Generic,
+    Iterable,
     Literal,
     Optional,
     SupportsFloat,
@@ -32,6 +33,7 @@
 from calliope import exceptions
 from calliope.attrdict import AttrDict
 from calliope.backend import helper_functions, parsing
+from calliope.exceptions import warn as model_warn
 from calliope.io import load_config
 from calliope.util.schema import (
     MATH_SCHEMA,
@@ -758,6 +760,18 @@ def to_lp(self, path: Union[str, Path]) -> None:
             path (Union[str, Path]): Path to which the LP file will be written.
         """
 
+    @property
+    @abstractmethod
+    def has_integer_or_binary_variables(self) -> bool:
+        """Confirm whether or not the built model has binary or integer decision variables.
+
+        This can be used to understand how long the optimisation may take (MILP problems are harder to solve than LP ones),
+        and to verify whether shadow prices can be tracked (they cannot be tracked in MILP problems).
+
+        Returns:
+            bool: If True, the built model has binary or integer decision variables. If False, all decision variables are continuous.
+        """
+
     @abstractmethod
     def _solve(
         self,
@@ -798,8 +812,8 @@ def _solve(
 
     def load_results(self) -> xr.Dataset:
         """
-        Evaluate backend decision variables, global expressions, and parameters (if not in inputs)
-        after a successful model run.
+        Evaluate backend decision variables, global expressions, parameters (if not in
+        inputs), and shadow_prices (if tracked), after a successful model run.
 
         Returns:
             xr.Dataset: Dataset of optimal solution results (all numeric data).
@@ -824,7 +838,14 @@ def _drop_attrs(da):
             if expr.notnull().any()
         }
 
-        results = xr.Dataset({**all_variables, **all_global_expressions}).astype(float)
+        all_shadow_prices = {
+            f"shadow_price_{constraint}": self.shadow_prices.get(constraint)
+            for constraint in self.shadow_prices.tracked
+        }
+
+        results = xr.Dataset(
+            {**all_variables, **all_global_expressions, **all_shadow_prices}
+        ).astype(float)
 
         return results
 
@@ -875,6 +896,8 @@ class ShadowPrices:
     To keep memory overhead low. Shadow price tracking is deactivated by default.
     """
 
+    _tracked: set = set()
+
     @abstractmethod
     def get(self, name) -> xr.DataArray:
         """Extract shadow prices (a.k.a. duals) from a constraint.
@@ -898,3 +921,36 @@ def deactivate(self):
     @abstractmethod
     def is_active(self) -> bool:
         "Check whether shadow price tracking is active or not"
+
+    @property
+    @abstractmethod
+    def available_constraints(self) -> Iterable:
+        "Iterable of constraints that are available to provide shadow prices on"
+
+    @property
+    def tracked(self) -> set:
+        "Constraints being tracked for automatic addition to the results dataset"
+        return self._tracked
+
+    def track_constraints(self, constraints_to_track: list):
+        """Track constraints if they are available in the built backend model.
+
+        If there is at least one available constraint to track,
+        `self.tracked` will be updated and shadow price tracking will be activated.
+
+        Args:
+            constraints_to_track (list): Constraint names to track
+        """
+        shadow_prices = set(constraints_to_track)
+        invalid_constraints = shadow_prices.difference(self.available_constraints)
+        valid_constraints = shadow_prices.intersection(self.available_constraints)
+        if invalid_constraints:
+            model_warn(
+                f"Invalid constraints {invalid_constraints} in `config.solve.shadow_prices`. "
+                "Their shadow prices will not be tracked."
+            )
+        # Only actually activate shadow price tracking if at least one valid
+        # constraint remains in the list after filtering out invalid ones
+        if valid_constraints:
+            self.activate()
+        self._tracked = valid_constraints

src/calliope/backend/pyomo_backend_model.py

@@ -27,6 +27,7 @@
 import xarray as xr
 from pyomo.common.tempfiles import TempfileManager  # type: ignore
 from pyomo.opt import SolverFactory  # type: ignore
+from pyomo.util.model_size import build_model_size_report  # type: ignore
 
 from calliope.backend import backend_model, parsing
 from calliope.exceptions import BackendError, BackendWarning
@@ -487,6 +488,14 @@ def unfix_variable(self, name: str, where: Optional[xr.DataArray] = None) -> Non
             variable_da = variable_da.where(where.fillna(0))
         self._apply_func(self._unfix_pyomo_variable, variable_da)
 
+    @property
+    def has_integer_or_binary_variables(self) -> bool:
+        model_report = build_model_size_report(self._instance)
+        binaries = model_report["activated"]["binary_variables"]
+        integers = model_report["activated"]["integer_variables"]
+        number_of_binary_and_integer_vars = binaries + integers
+        return number_of_binary_and_integer_vars > 0
+
     def _get_capacity_bound(
         self, bound: Any, name: str, references: set
     ) -> xr.DataArray:
@@ -917,7 +926,11 @@ def get(self, name: str) -> xr.DataArray:
         )
 
     def activate(self):
-        self._dual_obj.activate()
+        if self._backend_obj.has_integer_or_binary_variables:
+            warning_text = "Shadow price tracking on a model with binary or integer variables is not possible. Proceeding without activating shadow price tracking."
+            model_warn(warning_text, _class=BackendWarning)
+        else:
+            self._dual_obj.activate()
 
     def deactivate(self):
         self._dual_obj.deactivate()
@@ -926,6 +939,10 @@ def deactivate(self):
     def is_active(self) -> bool:
         return self._dual_obj.active
 
+    @property
+    def available_constraints(self) -> Iterable:
+        return self._backend_obj.constraints.data_vars
+
     @staticmethod
     def _duals_from_pyomo_constraint(
         val: pmo.constraint, *, dual_getter: pmo.suffix

src/calliope/config/config_schema.yaml

@@ -161,6 +161,14 @@ properties:
             type: number
             default: 1e-10
             description: On postprocessing the optimisation results, values smaller than this threshold will be considered as optimisation artefacts and will be set to zero.
+          shadow_prices:
+            type: array
+            uniqueItems: true
+            items:
+              type: string
+              description: Names of model constraints.
+            default: []
+            description: List of constraints for which to extract shadow prices. Shadow prices will be added as variables to the model results as `shadow_price_{constraintname}`.
 
   parameters:
     type: [object, "null"]

src/calliope/model.py

@@ -448,6 +448,9 @@ def solve(self, force: bool = False, warmstart: bool = False, **kwargs) -> None:
 
         solver_config = update_then_validate_config("solve", self.config, **kwargs)
 
+        shadow_prices = solver_config.get("shadow_prices", [])
+        self.backend.shadow_prices.track_constraints(shadow_prices)
+
         if run_mode == "operate":
             if not self._model_data.attrs["allow_operate_mode"]:
                 raise exceptions.ModelError(

src/calliope/postprocess/postprocess.py

@@ -47,6 +47,7 @@ def postprocess_model_results(
     results["total_levelised_cost"] = systemwide_levelised_cost(
         results, model_data, total=True
     )
+
     results = clean_results(results, zero_threshold)
 
     for var_data in results.data_vars.values():

tests/common/test_model/scenarios.yaml

@@ -444,3 +444,9 @@ overrides:
       demand_elec:
         add_dimensions:
           parameters: sink_use_max
+
+  shadow_prices:
+    config.solve.shadow_prices: ["system_balance", "balance_demand"]
+
+  shadow_prices_invalid_constraint:
+    config.solve.shadow_prices: ["flow_cap_max_foobar"]

tests/test_backend_pyomo.py

@@ -2477,6 +2477,12 @@ def test_unfix_variable_where(self, simple_supply):
         assert fixed.sel(techs="test_demand_elec").all()
         assert not fixed.where(where).all()
 
+    def test_has_integer_or_binary_variables_lp(self, simple_supply):
+        assert not simple_supply.backend.has_integer_or_binary_variables
+
+    def test_has_integer_or_binary_variables_milp(self, supply_milp):
+        assert supply_milp.backend.has_integer_or_binary_variables
+
 
 class TestShadowPrices:
     @pytest.fixture(scope="function")
@@ -2491,13 +2497,50 @@ def supply_milp(self):
         m.build()
         return m
 
+    @pytest.fixture(scope="function")
+    def simple_supply_with_yaml_shadow_prices(self):
+        m = build_model({}, "simple_supply,two_hours,investment_costs,shadow_prices")
+        m.build()
+        return m
+
+    @pytest.fixture(scope="function")
+    def simple_supply_yaml(self):
+        m = build_model({}, "simple_supply,two_hours,investment_costs,shadow_prices")
+        m.build()
+        return m
+
+    @pytest.fixture(scope="function")
+    def simple_supply_yaml_invalid(self):
+        m = build_model(
+            {},
+            "simple_supply,two_hours,investment_costs,shadow_prices_invalid_constraint",
+        )
+        m.build()
+        return m
+
+    @pytest.fixture(scope="function")
+    def supply_milp_yaml(self):
+        m = build_model({}, "supply_milp,two_hours,investment_costs,shadow_prices")
+        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):
+    def test_available_constraints(self, simple_supply):
+        assert set(simple_supply.backend.shadow_prices.available_constraints) == set(
+            simple_supply.backend.constraints.data_vars
+        )
+
+    def test_activate_continuous_model(self, simple_supply):
         simple_supply.backend.shadow_prices.activate()
         assert simple_supply.backend.shadow_prices.is_active
 
+    def test_activate_milp_model(self, supply_milp):
+        with pytest.warns(exceptions.BackendWarning):
+            supply_milp.backend.shadow_prices.activate()
+        assert not supply_milp.backend.shadow_prices.is_active
+
     def test_deactivate(self, simple_supply):
         simple_supply.backend.shadow_prices.activate()
         simple_supply.backend.shadow_prices.deactivate()
@@ -2516,12 +2559,6 @@ def test_get_shadow_price_some_nan(self, simple_supply):
         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:
@@ -2531,3 +2568,37 @@ def test_shadow_prices_deactivated_with_cbc(self, simple_supply):
         assert not simple_supply.backend.shadow_prices.is_active
         shadow_prices = simple_supply.backend.shadow_prices.get("system_balance")
         assert shadow_prices.isnull().all()
+
+    def test_yaml_continuous_model_tracked(self, simple_supply_yaml):
+        # before solve, there are no constraints to track
+        assert not simple_supply_yaml.backend.shadow_prices.tracked
+
+        simple_supply_yaml.solve(solver="glpk")
+
+        assert simple_supply_yaml.backend.shadow_prices.tracked == {
+            "system_balance",
+            "balance_demand",
+        }
+
+    def test_yaml_continuous_model_result(self, simple_supply_yaml):
+        m = simple_supply_yaml
+        m.solve(solver="glpk")
+        assert m.results["shadow_price_system_balance"].sum().item() == pytest.approx(
+            0.0005030505
+        )
+        assert m.results["shadow_price_balance_demand"].sum().item() == pytest.approx(
+            0.0010061011
+        )
+
+    def test_yaml_milp_model(self, supply_milp_yaml):
+        assert not supply_milp_yaml.backend.shadow_prices.is_active
+
+    def test_yaml_with_invalid_constraint(self, simple_supply_yaml_invalid):
+        m = simple_supply_yaml_invalid
+        with pytest.warns(exceptions.ModelWarning) as warning:
+            m.solve()
+        assert check_error_or_warning(
+            warning, "Invalid constraints {'flow_cap_max_foobar'}"
+        )
+        # Since we listed only one (invalid) constraint, tracking should not be active
+        assert not m.backend.shadow_prices.is_active