Backport PR #520: Fix scheduling with efficiencies (#520)

Fix scheduling with imperfect efficiencies, which resulted in exceeding the device's lower SoC limit.

* Apply efficiencies to conversion from flow to stock change and vice versa

Signed-off-by: F.N. Claessen <felix@seita.nl>

* Reorder print statements for more convenient uncommenting

Signed-off-by: F.N. Claessen <felix@seita.nl>

* Fix test

Signed-off-by: F.N. Claessen <felix@seita.nl>

* Refactor: rename variables

Signed-off-by: F.N. Claessen <felix@seita.nl>

* changelog entry

Signed-off-by: F.N. Claessen <felix@seita.nl>

Signed-off-by: F.N. Claessen <felix@seita.nl>

documentation/changelog.rst

@@ -7,6 +7,7 @@ v0.11.3 | November XX, 2022
 
 Bugfixes
 -----------
+* Fix scheduling with imperfect efficiencies, which resulted in exceeding the device's lower SoC limit. [see `PR #520 <http://www.github.com/FlexMeasures/flexmeasures/pull/520>`_]
 * Fix scheduler for Charge Points when taking into account inflexible devices [see `PR #517 <http://www.github.com/FlexMeasures/flexmeasures/pull/517>`_]
 
 

flexmeasures/api/v1_3/tests/test_api_v1_3.py

@@ -112,7 +112,9 @@ def test_post_udi_event_and_get_device_message(
     # check targets, if applicable
     if "targets" in message:
         start_soc = message["value"] / 1000  # in MWh
-        soc_schedule = integrate_time_series(consumption_schedule, start_soc, 6)
+        soc_schedule = integrate_time_series(
+            consumption_schedule, start_soc, decimal_precision=6
+        )
         print(consumption_schedule)
         print(soc_schedule)
         for target in message["targets"]:

flexmeasures/api/v3_0/tests/test_sensor_schedules.py

@@ -90,7 +90,9 @@ def test_trigger_and_get_schedule(
     # check targets, if applicable
     if "targets" in message:
         start_soc = message["soc-at-start"] / 1000  # in MWh
-        soc_schedule = integrate_time_series(consumption_schedule, start_soc, 6)
+        soc_schedule = integrate_time_series(
+            consumption_schedule, start_soc, decimal_precision=6
+        )
         print(consumption_schedule)
         print(soc_schedule)
         for target in message["targets"]:

flexmeasures/data/models/planning/solver.py

@@ -46,8 +46,8 @@ def device_scheduler(  # noqa C901
         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)
-        derivative down efficiency: ratio of downwards flows (flow into EMS : flow out of device)
-        derivative up efficiency: ratio of upwards flows (flow into device : flow out of EMS)
+        derivative down efficiency: conversion efficiency of flow out of a device (flow out : stock decrease)
+        derivative up efficiency: conversion efficiency of flow into a device (stock increase : flow in)
     EMS constraints are on an EMS level. Handled constraints (listed by column name):
         derivative max: maximum flow
         derivative min: minimum flow
@@ -228,10 +228,12 @@ 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(
-                m.device_power_down[d, k] + m.device_power_up[d, k]
+                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)
             ),
             m.device_max[d, j],
@@ -275,12 +277,10 @@ def ems_flow_commitment_equalities(m, j):
         )
 
     def device_derivative_equalities(m, d, j):
-        """Couple device flows to EMS flows per device, applying efficiencies."""
+        """Couple device flows to EMS flows per device."""
         return (
             0,
-            m.device_power_up[d, j] / m.device_derivative_up_efficiency[d, j]
-            + m.device_power_down[d, j] * m.device_derivative_down_efficiency[d, j]
-            - m.ems_power[d, j],
+            m.device_power_up[d, j] + m.device_power_down[d, j] - m.ems_power[d, j],
             0,
         )
 
@@ -321,10 +321,7 @@ def cost_function(m):
     planned_costs = value(model.costs)
     planned_power_per_device = []
     for d in model.d:
-        planned_device_power = [
-            model.device_power_down[d, j].value + model.device_power_up[d, j].value
-            for j in model.j
-        ]
+        planned_device_power = [model.ems_power[d, j].value for j in model.j]
         planned_power_per_device.append(
             pd.Series(
                 index=pd.date_range(
@@ -335,7 +332,7 @@ def cost_function(m):
         )
 
     # model.pprint()
+    # model.display()
     # print(results.solver.termination_condition)
     # print(planned_costs)
-    # model.display()
     return planned_power_per_device, planned_costs, results

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

@@ -90,7 +90,13 @@ def test_battery_solver_day_2(add_battery_assets, roundtrip_efficiency: float):
         soc_max=soc_max,
         roundtrip_efficiency=roundtrip_efficiency,
     )
-    soc_schedule = integrate_time_series(schedule, soc_at_start, decimal_precision=6)
+    soc_schedule = integrate_time_series(
+        schedule,
+        soc_at_start,
+        up_efficiency=roundtrip_efficiency**0.5,
+        down_efficiency=roundtrip_efficiency**0.5,
+        decimal_precision=6,
+    )
 
     with pd.option_context("display.max_rows", None, "display.max_columns", 3):
         print(soc_schedule)

flexmeasures/utils/calculations.py

@@ -1,6 +1,7 @@
 """ Calculations """
+from __future__ import annotations
+
 from datetime import timedelta
-from typing import Optional
 
 import numpy as np
 import pandas as pd
@@ -37,11 +38,15 @@ def drop_nan_rows(a, b):
 
 
 def integrate_time_series(
-    s: pd.Series, s0: float, decimal_precision: Optional[int] = None
+    series: pd.Series,
+    initial_stock: float,
+    up_efficiency: float | pd.Series = 1,
+    down_efficiency: float | pd.Series = 1,
+    decimal_precision: int | None = None,
 ) -> pd.Series:
     """Integrate time series of length n and closed="left" (representing a flow)
     to a time series of length n+1 and closed="both" (representing a stock),
-    given a starting stock s0.
+    given an initial stock (i.e. the constant of integration).
     The unit of time is hours: i.e. the stock unit is flow unit times hours (e.g. a flow in kW becomes a stock in kWh).
     Optionally, set a decimal precision to round off the results (useful for tests failing over machine precision).
 
@@ -63,12 +68,24 @@ def integrate_time_series(
         2001-01-01 07:00:00    15.0
         dtype: float64
     """
-    resolution = pd.to_timedelta(s.index.freq)
+    resolution = pd.to_timedelta(series.index.freq)
+    stock_change = pd.Series(data=np.NaN, index=series.index)
+    stock_change.loc[series > 0] = series[series > 0] * (
+        up_efficiency[series > 0]
+        if isinstance(up_efficiency, pd.Series)
+        else up_efficiency
+    )
+    stock_change.loc[series <= 0] = series[series <= 0] / (
+        down_efficiency[series <= 0]
+        if isinstance(down_efficiency, pd.Series)
+        else down_efficiency
+    )
     int_s = pd.concat(
         [
-            pd.Series(s0, index=pd.date_range(s.index[0], periods=1)),
-            s.shift(1, freq=resolution).cumsum() * (resolution / timedelta(hours=1))
-            + s0,
+            pd.Series(initial_stock, index=pd.date_range(series.index[0], periods=1)),
+            stock_change.shift(1, freq=resolution).cumsum()
+            * (resolution / timedelta(hours=1))
+            + initial_stock,
         ]
     )
     if decimal_precision is not None: