Feat: full SoC preference

documentation/changelog.rst

@@ -13,6 +13,7 @@ New features
 * Support saving state-of-charge schedules to sensors with ``"%"`` unit, using the ``soc-max`` flex-model field as the capacity for unit conversion [see `PR #1996 <https://www.github.com/FlexMeasures/flexmeasures/pull/1996>`_]
 * Version headers (for server and API) in API responses [see `PR #2021 <https://www.github.com/FlexMeasures/flexmeasures/pull/2021>`_]
 * Show sensor attributes on sensor page, if not empty [see `PR #2015 <https://www.github.com/FlexMeasures/flexmeasures/pull/2015>`_]
+* Separate the ``StorageScheduler``'s tie-breaking preference for a full :abbr:`SoC (state of charge)` from its reported energy costs [see `PR #2023 <https://www.github.com/FlexMeasures/flexmeasures/pull/2023>`_]
 
 Infrastructure / Support
 ----------------------

documentation/dev/setup-and-guidelines.rst

@@ -63,7 +63,7 @@ On Linux and Windows, everything will be installed using Python packages.
 On MacOS, this will install all test dependencies, and locally install the HiGHS solver.
 For this to work, make sure you have `Homebrew <https://brew.sh/>`_ installed.
 
-Besides the HiGHS solver (as the current default), the CBC solver is required for tests as well. See `The install instructions <https://github.com/coin-or/Cbc?tab=readme-ov-file#binaries`_ for more information.
+Besides the HiGHS solver (as the current default), the CBC solver is required for tests as well. See `The install instructions <https://github.com/coin-or/Cbc?tab=readme-ov-file#binaries>`_ for more information.
 
 Configuration
 ^^^^^^^^^^^^^

flexmeasures/data/models/planning/linear_optimization.py

@@ -626,11 +626,22 @@ def cost_function(m):
         if cbc_path is not None:
             solver.set_executable(cbc_path)
 
-    # disable logs for the HiGHS solver in case that LOGGING_LEVEL is INFO
-    if current_app.config["LOGGING_LEVEL"] == "INFO" and (
-        "highs" in solver_name.lower()
-    ):
-        solver.options["output_flag"] = "false"
+    # Set tight tolerance for HiGHS solver
+    profile = {}
+    if "highs" in solver_name.lower():
+        profile = {
+            "mip_rel_gap": "0",
+            "mip_abs_gap": "0",
+            "primal_feasibility_tolerance": "1e-9",
+            "dual_feasibility_tolerance": "1e-9",
+            "mip_feasibility_tolerance": "1e-9",
+        }
+        # disable logs for the HiGHS solver in case that LOGGING_LEVEL is INFO
+        if current_app.config["LOGGING_LEVEL"] == "INFO":
+            profile["output_flag"] = "false"
+
+    for option_name, option_value in profile.items():
+        solver.options[option_name] = option_value
 
     # load_solutions=False to avoid a RuntimeError exception in appsi solvers when solving an infeasible problem.
     results = solver.solve(model, load_solutions=False)

flexmeasures/data/models/planning/storage.py

@@ -198,17 +198,6 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
         start = pd.Timestamp(start).tz_convert("UTC")
         end = pd.Timestamp(end).tz_convert("UTC")
 
-        # Add tiny price slope to prefer charging now rather than later, and discharging later rather than now.
-        # We penalise future consumption and reward future production with at most 1 per thousand times the energy price spread.
-        # todo: move to flow or stock commitment per device
-        if any(prefer_charging_sooner):
-            up_deviation_prices = add_tiny_price_slope(
-                up_deviation_prices, "event_value"
-            )
-            down_deviation_prices = add_tiny_price_slope(
-                down_deviation_prices, "event_value"
-            )
-
         # Create Series with EMS capacities
         ems_power_capacity_in_mw = get_continuous_series_sensor_or_quantity(
             variable_quantity=self.flex_context.get("ems_power_capacity_in_mw"),
@@ -441,23 +430,32 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
             # Take the contracted capacity as a hard constraint
             ems_constraints["derivative min"] = ems_production_capacity
 
-        # Flow commitments per device
+        # Commitments per device
 
-        # Add tiny price slope to prefer curtailing later rather than now.
-        # The price slope is half of the slope to prefer charging sooner
-        for d, prefer_curtailing_later_d in enumerate(prefer_curtailing_later):
-            if prefer_curtailing_later_d:
+        # StockCommitment per device to prefer a full storage by penalizing not being full
+        # This corresponds to a preference for charging now rather than later, and discharging later rather than now.
+        for d, (prefer_charging_sooner_d, prefer_curtailing_later_d) in enumerate(
+            zip(prefer_charging_sooner, prefer_curtailing_later)
+        ):
+            if prefer_charging_sooner_d:
                 tiny_price_slope = (
-                    add_tiny_price_slope(up_deviation_prices, "event_value")
+                    add_tiny_price_slope(
+                        up_deviation_prices, "event_value", order="desc"
+                    )
                     - up_deviation_prices
                 )
-                tiny_price_slope *= 0.5
-                commitment = FlowCommitment(
-                    name=f"prefer curtailing device {d} later",
-                    # Prefer curtailing consumption later by penalizing later consumption
-                    upwards_deviation_price=tiny_price_slope,
-                    # Prefer curtailing production later by penalizing later production
-                    downwards_deviation_price=-tiny_price_slope,
+                if prefer_curtailing_later_d:
+                    # Use a tiny price slope to prefer a fuller SoC sooner rather than later, by lowering penalties later
+                    penalty = tiny_price_slope
+                else:
+                    # Constant penalty
+                    penalty = tiny_price_slope.iloc[0][0]
+                commitment = StockCommitment(
+                    name=f"prefer a full storage {d} sooner",
+                    quantity=(soc_max[d] - soc_at_start[d])
+                    * (timedelta(hours=1) / resolution),
+                    upwards_deviation_price=0,
+                    downwards_deviation_price=-penalty,
                     index=index,
                     device=d,
                 )
@@ -940,7 +938,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
     def convert_to_commitments(
         self,
         **timing_kwargs,
-    ) -> list[FlowCommitment]:
+    ) -> list[FlowCommitment | StockCommitment]:
         """Convert list of commitment specifications (dicts) to a list of FlowCommitments."""
         commitment_specs = self.flex_context.get("commitments", [])
         if len(commitment_specs) == 0:
@@ -1293,7 +1291,7 @@ def compute(self, skip_validation: bool = False) -> SchedulerOutputType:
 
 
 class StorageScheduler(MetaStorageScheduler):
-    __version__ = "7"
+    __version__ = "8"
     __author__ = "Seita"
 
     fallback_scheduler_class: Type[Scheduler] = StorageFallbackScheduler

flexmeasures/data/models/planning/tests/test_solver.py

@@ -270,6 +270,7 @@ def run_test_charge_discharge_sign(
             for soc_at_start_d in soc_at_start
         ],
     )
+    assert results.solver.termination_condition == "optimal"
 
     device_power_sign = pd.Series(model.device_power_sign.extract_values())[0]
     device_power_up = pd.Series(model.device_power_up.extract_values())[0]
@@ -816,8 +817,8 @@ def compute_schedule(flex_model):
 
     # soc maxima and soc minima
     soc_maxima = [
-        {"datetime": "2015-01-02T15:00:00+01:00", "value": 1.0},
-        {"datetime": "2015-01-02T16:00:00+01:00", "value": 1.0},
+        {"datetime": "2015-01-02T12:00:00+01:00", "value": 1.0},
+        {"datetime": "2015-01-02T13:00:00+01:00", "value": 1.0},
     ]
 
     soc_minima = [{"datetime": "2015-01-02T08:00:00+01:00", "value": 3.5}]
@@ -853,7 +854,7 @@ def compute_schedule(flex_model):
 
     # test for soc_maxima
     # check that the local maximum constraint is respected
-    assert soc_schedule_2.loc["2015-01-02T15:00:00+01:00"] <= 1.0
+    assert soc_schedule_2.loc["2015-01-02T13:00:00+01:00"] <= 1.0
 
     # test for soc_targets
     # check that the SOC target (at 19 pm, local time) is met
@@ -1787,10 +1788,10 @@ def test_battery_stock_delta_sensor(
         - Battery of size 2 MWh.
         - Consumption capacity of the battery is 2 MW.
         - The battery cannot discharge.
-    With these settings, the battery needs to charge at a power or greater than the usage forecast
+    With these settings, the battery needs to charge at a power equal or greater than the usage forecast
     to keep the SOC within bounds ([0, 2 MWh]).
     """
-    _, battery = get_sensors_from_db(db, add_battery_assets)
+    epex_da, battery = get_sensors_from_db(db, add_battery_assets)
     tz = pytz.timezone("Europe/Amsterdam")
     start = tz.localize(datetime(2015, 1, 1))
     end = tz.localize(datetime(2015, 1, 2))
@@ -1835,9 +1836,20 @@ def test_battery_stock_delta_sensor(
         with pytest.raises(InfeasibleProblemException):
             scheduler.compute()
     elif stock_delta_sensor is None:
-        # No usage -> the battery does not charge
+        # No usage -> the battery only charges when energy is free
+        free_hour = "2015-01-01 17:00:00+00:00"
+        prices = epex_da.search_beliefs(start, end, resolution=resolution)
+        zero_prices = prices[prices.event_value == 0]
+        assert all(
+            zero_prices.event_starts.hour == pd.Timestamp(free_hour).hour
+        ), "this test assumes a single hour of free energy from 5 to 6 PM UTC"
         schedule = scheduler.compute()
-        assert all(schedule == 0)
+        assert all(
+            schedule[~schedule.index.isin(zero_prices.event_starts)] == 0
+        ), "no charging expected when energy is not free, given no soc-usage"
+        assert all(
+            schedule[schedule.index.isin(zero_prices.event_starts)] == capacity
+        ), "max charging expected when energy is free, because of preference to have a full SoC"
     else:
         # Some usage -> the battery needs to charge
         schedule = scheduler.compute()
@@ -2235,11 +2247,14 @@ def test_battery_storage_different_units(
         battery_name="Test battery",
         power_sensor_name=power_sensor_name,
     )
-    tz = pytz.timezone("Europe/Amsterdam")
+    tz = pytz.timezone(epex_da.timezone)
 
     # transition from cheap to expensive (90 -> 100)
     start = tz.localize(datetime(2015, 1, 2, 14, 0, 0))
     end = tz.localize(datetime(2015, 1, 2, 16, 0, 0))
+    assert len(epex_da.search_beliefs(start, end)) == 2
+    assert epex_da.search_beliefs(start, end).values[0][0] == 90
+    assert epex_da.search_beliefs(start, end).values[1][0] == 100
     resolution = timedelta(minutes=15)
 
     flex_model = {

flexmeasures/data/models/planning/utils.py

@@ -2,6 +2,7 @@
 
 from packaging import version
 from datetime import date, datetime, timedelta
+from typing import Literal
 
 from flask import current_app
 import pandas as pd
@@ -67,21 +68,32 @@ def initialize_index(
 
 
 def add_tiny_price_slope(
-    orig_prices: pd.DataFrame, col_name: str = "event_value", d: float = 10**-4
+    orig_prices: pd.DataFrame,
+    col_name: str = "event_value",
+    d: float = 10**-4,
+    order: Literal["asc", "desc"] = "asc",
 ) -> pd.DataFrame:
     """Add tiny price slope to col_name to represent e.g. inflation as a simple linear price increase.
     This is meant to break ties, when multiple time slots have equal prices, in favour of acting sooner.
-    We penalise the future with at most d times the price spread (1 per thousand by default).
+    We penalise the future with at most d times the price spread (1 per thousand by default),
+    divided over the number of planning steps.
     """
     prices = orig_prices.copy()
     price_spread = prices[col_name].max() - prices[col_name].min()
     if price_spread > 0:
-        max_penalty = price_spread * d
+        max_penalty = price_spread * d / len(prices)
     else:
-        max_penalty = d
-    prices[col_name] = prices[col_name] + np.linspace(
-        0, max_penalty, prices[col_name].size
-    )
+        max_penalty = d / len(prices)
+    if order == "asc":
+        prices[col_name] = prices[col_name] + np.linspace(
+            0, max_penalty, prices[col_name].size
+        )
+    elif order == "desc":
+        prices[col_name] = prices[col_name] + np.linspace(
+            max_penalty, 0, prices[col_name].size
+        )
+    else:
+        raise ValueError(f"order must be 'asc' or 'desc', got '{order}'")
     return prices