feat(scheduling): Storage scheduler supports losses over time

flexmeasures/api/v3_0/sensors.py

@@ -271,6 +271,7 @@ def trigger_schedule(  # noqa: C901
         at which the state of charge (soc) is 12.1 kWh, with a target state of charge of 25 kWh at 4.00pm.
         The minimum and maximum soc are set to 10 and 25 kWh, respectively.
         Roundtrip efficiency for use in scheduling is set to 98%.
+        Storage efficiency is set to 99.99%, denoting the state of charge left after each time step equal to the sensor's resolution.
         Aggregate consumption (of all devices within this EMS) should be priced by sensor 9,
         and aggregate production should be priced by sensor 10,
         where the aggregate power flow in the EMS is described by the sum over sensors 13, 14 and 15
@@ -294,6 +295,7 @@ def trigger_schedule(  # noqa: C901
                     "soc-min": 10,
                     "soc-max": 25,
                     "roundtrip-efficiency": 0.98,
+                    "storage-efficiency": 0.9999,
                 },
                 "flex-context": {
                     "consumption-price-sensor": 9,

flexmeasures/api/v3_0/tests/test_sensor_schedules.py

@@ -221,13 +221,19 @@ def test_trigger_and_get_schedule(
         roundtrip_efficiency = (
             float(message["roundtrip-efficiency"].replace("%", "")) / 100.0
         )
+        storage_efficiency = (
+            float(message["storage-efficiency"].replace("%", "")) / 100.0
+        )
         soc_targets = message.get("soc-targets")
     else:
         start_soc = message["flex-model"]["soc-at-start"] / 1000  # in MWh
         roundtrip_efficiency = (
             float(message["flex-model"]["roundtrip-efficiency"].replace("%", ""))
             / 100.0
         )
+        storage_efficiency = (
+            float(message["flex-model"]["storage-efficiency"].replace("%", "")) / 100.0
+        )
         soc_targets = message["flex-model"].get("soc-targets")
     resolution = sensor.event_resolution
     if soc_targets:
@@ -271,6 +277,7 @@ def test_trigger_and_get_schedule(
             start_soc,
             up_efficiency=roundtrip_efficiency**0.5,
             down_efficiency=roundtrip_efficiency**0.5,
+            storage_efficiency=storage_efficiency,
             decimal_precision=6,
         )
         print(consumption_schedule)

flexmeasures/api/v3_0/tests/utils.py

@@ -61,6 +61,7 @@ def message_for_trigger_schedule(
         "soc-max": 40,  # in kWh, according to soc-unit
         "soc-unit": "kWh",
         "roundtrip-efficiency": "98%",
+        "storage-efficiency": "99.99%",
     }
     if with_targets:
         if realistic_targets:

flexmeasures/cli/data_add.py

@@ -45,6 +45,7 @@
     LongitudeField,
     SensorIdField,
 )
+from flexmeasures.data.schemas.scheduling.storage import EfficiencyField
 from flexmeasures.data.schemas.sensors import SensorSchema
 from flexmeasures.data.schemas.units import QuantityField
 from flexmeasures.data.schemas.generic_assets import (
@@ -1005,11 +1006,22 @@ def create_schedule(ctx):
 @click.option(
     "--roundtrip-efficiency",
     "roundtrip_efficiency",
-    type=QuantityField("%", validate=validate.Range(min=0, max=1)),
+    type=EfficiencyField(),
     required=False,
     default=1,
     help="Round-trip efficiency (e.g. 85% or 0.85) to use for the schedule. Defaults to 100% (no losses).",
 )
+@click.option(
+    "--storage-efficiency",
+    "storage_efficiency",
+    type=EfficiencyField(),
+    required=False,
+    default=1,
+    help="Storage efficiency (e.g. 95% or 0.95) to use for the schedule,"
+    " applied over each time step equal to the sensor resolution."
+    " For example, a storage efficiency of 99 percent per (absolute) day, for scheduling a 1-hour resolution sensor, should be passed as a storage efficiency of 0.99**(1/24)."
+    " Defaults to 100% (no losses).",
+)
 @click.option(
     "--as-job",
     is_flag=True,
@@ -1029,6 +1041,7 @@ def add_schedule_for_storage(
     soc_min: ur.Quantity | None = None,
     soc_max: ur.Quantity | None = None,
     roundtrip_efficiency: ur.Quantity | None = None,
+    storage_efficiency: ur.Quantity | None = None,
     as_job: bool = False,
 ):
     """Create a new schedule for a storage asset.
@@ -1084,6 +1097,8 @@ def add_schedule_for_storage(
         soc_max = convert_units(soc_max.magnitude, str(soc_max.units), "MWh", capacity=capacity_str)  # type: ignore
     if roundtrip_efficiency is not None:
         roundtrip_efficiency = roundtrip_efficiency.magnitude / 100.0
+    if storage_efficiency is not None:
+        storage_efficiency = storage_efficiency.magnitude / 100.0
 
     scheduling_kwargs = dict(
         start=start,
@@ -1097,6 +1112,7 @@ def add_schedule_for_storage(
             "soc-max": soc_max,
             "soc-unit": "MWh",
             "roundtrip-efficiency": roundtrip_efficiency,
+            "storage-efficiency": storage_efficiency,
         },
         flex_context={
             "consumption-price-sensor": consumption_price_sensor.id,

flexmeasures/data/models/planning/linear_optimization.py

@@ -21,6 +21,7 @@
 from pyomo.opt import SolverFactory, SolverResults
 
 from flexmeasures.data.models.planning.utils import initialize_series
+from flexmeasures.utils.calculations import apply_stock_changes_and_losses
 
 infinity = float("inf")
 
@@ -31,6 +32,7 @@ def device_scheduler(  # noqa C901
     commitment_quantities: List[pd.Series],
     commitment_downwards_deviation_price: Union[List[pd.Series], List[float]],
     commitment_upwards_deviation_price: Union[List[pd.Series], List[float]],
+    initial_stock: float = 0,
 ) -> Tuple[List[pd.Series], float, SolverResults]:
     """This generic device scheduler is able to handle an EMS with multiple devices,
     with various types of constraints on the EMS level and on the device level,
@@ -43,6 +45,7 @@ def device_scheduler(  # noqa C901
         max: maximum stock assuming an initial stock of zero (e.g. in MWh or boxes)
         min: minimum stock assuming an initial stock of zero
         equal: exact amount of stock (we do this by clamping min and max)
+        efficiency: amount of stock left at the next datetime (the rest is lost)
         derivative max: maximum flow (e.g. in MW or boxes/h)
         derivative min: minimum flow
         derivative equals: exact amount of flow (we do this by clamping derivative min and derivative max)
@@ -171,6 +174,16 @@ def ems_derivative_min_select(m, j):
         else:
             return v
 
+    def device_efficiency(m, d, j):
+        """Assume perfect efficiency if no efficiency information is available."""
+        try:
+            eff = device_constraints[d]["efficiency"].iloc[j]
+        except KeyError:
+            return 1
+        if np.isnan(eff):
+            return 1
+        return eff
+
     def device_derivative_down_efficiency(m, d, j):
         """Assume perfect efficiency if no efficiency information is available."""
         try:
@@ -206,6 +219,7 @@ def device_derivative_up_efficiency(m, d, j):
     )
     model.ems_derivative_max = Param(model.j, initialize=ems_derivative_max_select)
     model.ems_derivative_min = Param(model.j, initialize=ems_derivative_min_select)
+    model.device_efficiency = Param(model.d, model.j, initialize=device_efficiency)
     model.device_derivative_down_efficiency = Param(
         model.d, model.j, initialize=device_derivative_down_efficiency
     )
@@ -228,14 +242,24 @@ def device_derivative_up_efficiency(m, d, j):
 
     # Add constraints as a tuple of (lower bound, value, upper bound)
     def device_bounds(m, d, j):
-        """Apply efficiencies to conversion from flow to stock change and vice versa."""
-        return (
-            m.device_min[d, j],
-            sum(
+        """Apply conversion efficiencies to conversion from flow to stock change and vice versa,
+        and apply storage efficiencies to stock levels from one datetime to the next."""
+        stock_changes = [
+            (
                 m.device_power_down[d, k] / m.device_derivative_down_efficiency[d, k]
                 + m.device_power_up[d, k] * m.device_derivative_up_efficiency[d, k]
-                for k in range(0, j + 1)
-            ),
+            )
+            for k in range(0, j + 1)
+        ]
+        efficiencies = [m.device_efficiency[d, k] for k in range(0, j + 1)]
+        return (
+            m.device_min[d, j],
+            [
+                stock - initial_stock
+                for stock in apply_stock_changes_and_losses(
+                    initial_stock, stock_changes, efficiencies
+                )
+            ][-1],
             m.device_max[d, j],
         )
 

flexmeasures/data/models/planning/storage.py

@@ -56,6 +56,7 @@ def compute(
         soc_min = self.flex_model.get("soc_min")
         soc_max = self.flex_model.get("soc_max")
         roundtrip_efficiency = self.flex_model.get("roundtrip_efficiency")
+        storage_efficiency = self.flex_model.get("storage_efficiency")
         prefer_charging_sooner = self.flex_model.get("prefer_charging_sooner", True)
 
         consumption_price_sensor = self.flex_context.get("consumption_price_sensor")
@@ -114,6 +115,7 @@ def compute(
             "equals",
             "max",
             "min",
+            "efficiency",
             "derivative equals",
             "derivative max",
             "derivative min",
@@ -167,6 +169,9 @@ def compute(
         )
         device_constraints[0]["derivative up efficiency"] = roundtrip_efficiency**0.5
 
+        # Apply storage efficiency (accounts for losses over time)
+        device_constraints[0]["efficiency"] = storage_efficiency
+
         # Set up EMS constraints
         columns = ["derivative max", "derivative min"]
         ems_constraints = initialize_df(columns, start, end, resolution)
@@ -181,6 +186,7 @@ def compute(
             commitment_quantities,
             commitment_downwards_deviation_price,
             commitment_upwards_deviation_price,
+            initial_stock=soc_at_start * (timedelta(hours=1) / resolution),
         )
         if scheduler_results.solver.termination_condition == "infeasible":
             # Fallback policy if the problem was unsolvable
@@ -250,7 +256,19 @@ def deserialize_flex_config(self):
             elif self.sensor.unit in ("MW", "kW"):
                 self.flex_model["soc-unit"] = self.sensor.unit + "h"
 
+        # Check for storage efficiency
+        # todo: simplify to: `if self.flex_model.get("storage-efficiency") is None:`
+        if (
+            "storage-efficiency" not in self.flex_model
+            or self.flex_model["storage-efficiency"] is None
+        ):
+            # Get default from sensor, or use 100% otherwise
+            self.flex_model["storage-efficiency"] = self.sensor.get_attribute(
+                "storage_efficiency", 1
+            )
+
         # Check for round-trip efficiency
+        # todo: simplify to: `if self.flex_model.get("roundtrip-efficiency") is None:`
         if (
             "roundtrip-efficiency" not in self.flex_model
             or self.flex_model["roundtrip-efficiency"] is None

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

@@ -11,12 +11,46 @@
 from flexmeasures.data.models.planning.utils import (
     initialize_series,
 )
-from flexmeasures.utils.calculations import integrate_time_series
+from flexmeasures.utils.calculations import (
+    apply_stock_changes_and_losses,
+    integrate_time_series,
+)
 
 
 TOLERANCE = 0.00001
 
 
+@pytest.mark.parametrize(
+    "initial_stock, stock_deltas, expected_stocks, storage_efficiency",
+    [
+        (
+            1000,
+            [100, -100, -100, 100],
+            [1000, 1089, 979.11, 870.3189, 960.615711],
+            0.99,
+        ),
+        (
+            2.5,
+            [-0.5, -0.5, -0.5, -0.5],
+            [2.5, 1.8, 1.17, 0.603, 0.0927],
+            0.9,
+        ),
+    ],
+)
+def test_storage_loss_function(
+    initial_stock, stock_deltas, expected_stocks, storage_efficiency
+):
+    stocks = apply_stock_changes_and_losses(
+        initial_stock,
+        stock_deltas,
+        storage_efficiency=storage_efficiency,
+        how="left",
+        decimal_precision=6,
+    )
+    print(stocks)
+    assert all(a == b for a, b in zip(stocks, expected_stocks))
+
+
 @pytest.mark.parametrize("use_inflexible_device", [False, True])
 def test_battery_solver_day_1(
     add_battery_assets, add_inflexible_device_forecasts, use_inflexible_device
@@ -60,14 +94,18 @@ def test_battery_solver_day_1(
 
 
 @pytest.mark.parametrize(
-    "roundtrip_efficiency",
+    "roundtrip_efficiency, storage_efficiency",
     [
-        1,
-        0.99,
-        0.01,
+        (1, 1),
+        (1, 0.999),
+        (1, 0.5),
+        (0.99, 1),
+        (0.01, 1),
     ],
 )
-def test_battery_solver_day_2(add_battery_assets, roundtrip_efficiency: float):
+def test_battery_solver_day_2(
+    add_battery_assets, roundtrip_efficiency: float, storage_efficiency: float
+):
     """Check battery scheduling results for day 2, which is set up with
     8 expensive, then 8 cheap, then again 8 expensive hours.
     If efficiency losses aren't too bad, we expect the scheduler to:
@@ -98,6 +136,7 @@ def test_battery_solver_day_2(add_battery_assets, roundtrip_efficiency: float):
             "soc-min": soc_min,
             "soc-max": soc_max,
             "roundtrip-efficiency": roundtrip_efficiency,
+            "storage-efficiency": storage_efficiency,
         },
     )
     schedule = scheduler.compute()
@@ -106,6 +145,7 @@ def test_battery_solver_day_2(add_battery_assets, roundtrip_efficiency: float):
         soc_at_start,
         up_efficiency=roundtrip_efficiency**0.5,
         down_efficiency=roundtrip_efficiency**0.5,
+        storage_efficiency=storage_efficiency,
         decimal_precision=6,
     )
 
@@ -124,22 +164,30 @@ def test_battery_solver_day_2(add_battery_assets, roundtrip_efficiency: float):
         soc_min, battery.get_attribute("min_soc_in_mwh")
     )  # Battery sold out at the end of its planning horizon
 
-    # As long as the roundtrip efficiency isn't too bad (I haven't computed the actual switch point)
-    if roundtrip_efficiency > 0.9:
+    # As long as the efficiencies aren't too bad (I haven't computed the actual switch points)
+    if roundtrip_efficiency > 0.9 and storage_efficiency > 0.9:
         assert soc_schedule.loc[start + timedelta(hours=8)] == max(
             soc_min, battery.get_attribute("min_soc_in_mwh")
         )  # Sell what you begin with
         assert soc_schedule.loc[start + timedelta(hours=16)] == min(
             soc_max, battery.get_attribute("max_soc_in_mwh")
         )  # Buy what you can to sell later
-    else:
-        # If the roundtrip efficiency is poor, best to stand idle
+    elif storage_efficiency > 0.9:
+        # If only the roundtrip efficiency is poor, best to stand idle (keep a high SoC as long as possible)
         assert soc_schedule.loc[start + timedelta(hours=8)] == battery.get_attribute(
             "soc_in_mwh"
         )
         assert soc_schedule.loc[start + timedelta(hours=16)] == battery.get_attribute(
             "soc_in_mwh"
         )
+    else:
+        # If the storage efficiency is poor, regardless of whether the roundtrip efficiency is poor, best to sell asap
+        assert soc_schedule.loc[start + timedelta(hours=8)] == max(
+            soc_min, battery.get_attribute("min_soc_in_mwh")
+        )
+        assert soc_schedule.loc[start + timedelta(hours=16)] == max(
+            soc_min, battery.get_attribute("min_soc_in_mwh")
+        )
 
 
 @pytest.mark.parametrize(
@@ -186,6 +234,9 @@ def test_charging_station_solver_day_2(target_soc, charging_station_name):
             "roundtrip_efficiency": charging_station.get_attribute(
                 "roundtrip_efficiency", 1
             ),
+            "storage_efficiency": charging_station.get_attribute(
+                "storage_efficiency", 1
+            ),
             "soc_targets": soc_targets,
         },
     )
@@ -259,6 +310,9 @@ def test_fallback_to_unsolvable_problem(target_soc, charging_station_name):
             "roundtrip_efficiency": charging_station.get_attribute(
                 "roundtrip_efficiency", 1
             ),
+            "storage_efficiency": charging_station.get_attribute(
+                "storage_efficiency", 1
+            ),
             "soc_targets": soc_targets,
         },
     )
@@ -349,6 +403,7 @@ def test_building_solver_day_2(
             "soc_min": soc_min,
             "soc_max": soc_max,
             "roundtrip_efficiency": battery.get_attribute("roundtrip_efficiency", 1),
+            "storage_efficiency": battery.get_attribute("storage_efficiency", 1),
         },
         flex_context={
             "inflexible_device_sensors": inflexible_devices.values(),