From 3e7384750bcff8ae47511b4f09851d2262746a93 Mon Sep 17 00:00:00 2001 From: Felix Claessen <30658763+Flix6x@users.noreply.github.com> Date: Wed, 2 Nov 2022 21:07:37 +0100 Subject: [PATCH] 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 * Reorder print statements for more convenient uncommenting Signed-off-by: F.N. Claessen * Fix test Signed-off-by: F.N. Claessen * Refactor: rename variables Signed-off-by: F.N. Claessen * changelog entry Signed-off-by: F.N. Claessen Signed-off-by: F.N. Claessen --- documentation/changelog.rst | 1 + flexmeasures/api/v1_3/tests/test_api_v1_3.py | 4 ++- .../api/v3_0/tests/test_sensor_schedules.py | 4 ++- flexmeasures/data/models/planning/solver.py | 21 ++++++------- .../data/models/planning/tests/test_solver.py | 8 ++++- flexmeasures/utils/calculations.py | 31 ++++++++++++++----- 6 files changed, 47 insertions(+), 22 deletions(-) diff --git a/documentation/changelog.rst b/documentation/changelog.rst index 174fca2c4..ba4038efa 100644 --- a/documentation/changelog.rst +++ b/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 `_] * Fix scheduler for Charge Points when taking into account inflexible devices [see `PR #517 `_] diff --git a/flexmeasures/api/v1_3/tests/test_api_v1_3.py b/flexmeasures/api/v1_3/tests/test_api_v1_3.py index 54429774a..0ab3221e9 100644 --- a/flexmeasures/api/v1_3/tests/test_api_v1_3.py +++ b/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"]: diff --git a/flexmeasures/api/v3_0/tests/test_sensor_schedules.py b/flexmeasures/api/v3_0/tests/test_sensor_schedules.py index 9c9ac912a..6d9a4516d 100644 --- a/flexmeasures/api/v3_0/tests/test_sensor_schedules.py +++ b/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"]: diff --git a/flexmeasures/data/models/planning/solver.py b/flexmeasures/data/models/planning/solver.py index e16af9488..42c8acf29 100644 --- a/flexmeasures/data/models/planning/solver.py +++ b/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 diff --git a/flexmeasures/data/models/planning/tests/test_solver.py b/flexmeasures/data/models/planning/tests/test_solver.py index 15396e66b..0295079eb 100644 --- a/flexmeasures/data/models/planning/tests/test_solver.py +++ b/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) diff --git a/flexmeasures/utils/calculations.py b/flexmeasures/utils/calculations.py index fef972f82..84cb8dd79 100644 --- a/flexmeasures/utils/calculations.py +++ b/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: