diff --git a/documentation/api/change_log.rst b/documentation/api/change_log.rst index 406d898cc..b19fdf993 100644 --- a/documentation/api/change_log.rst +++ b/documentation/api/change_log.rst @@ -5,6 +5,11 @@ API change log .. note:: The FlexMeasures API follows its own versioning scheme. This is also reflected in the URL, allowing developers to upgrade at their own pace. +v3.0-10 | 2023-06-12 +"""""""""""""""""""" + +- Introduced the ``storage-efficiency`` field to the ``flex-model``field for `/sensors//schedules/trigger` (POST). + v3.0-9 | 2023-04-26 """"""""""""""""""" diff --git a/documentation/api/notation.rst b/documentation/api/notation.rst index 9b3de79d8..801063093 100644 --- a/documentation/api/notation.rst +++ b/documentation/api/notation.rst @@ -185,27 +185,36 @@ This means that API and CLI users don't have to send the whole flex model every Here are the three types of flexibility models you can expect to be built-in: -1) For storage devices (e.g. batteries, charge points, electric vehicle batteries connected to charge points), the schedule deals with the state of charge (SOC). +1) For **storage devices** (e.g. batteries, and :abbr:`EV (electric vehicle)` batteries connected to charge points), the schedule deals with the state of charge (SOC). - The possible flexibility parameters are: + The possible flexibility parameters are: - - ``soc-at-start`` (defaults to 0) - - ``soc-unit`` (kWh or MWh) - - ``soc-min`` (defaults to 0) - - ``soc-max`` (defaults to max soc target) - - ``soc-targets`` (defaults to NaN values) - - ``roundtrip-efficiency`` (defaults to 100%) - - ``prefer-charging-sooner`` (defaults to True, also signals a preference to discharge later) + - ``soc-at-start`` (defaults to 0) + - ``soc-unit`` (kWh or MWh) + - ``soc-min`` (defaults to 0) + - ``soc-max`` (defaults to max soc target) + - ``soc-minima`` (defaults to NaN values) + - ``soc-maxima`` (defaults to NaN values) + - ``soc-targets`` (defaults to NaN values) + - ``roundtrip-efficiency`` (defaults to 100%) + - ``storage-efficiency`` (defaults to 100%) [#]_ + - ``prefer-charging-sooner`` (defaults to True, also signals a preference to discharge later) - For some examples, see the `[POST] /sensors/(id)/schedules/trigger <../api/v3_0.html#post--api-v3_0-sensors-(id)-schedules-trigger>`_ endpoint docs. + .. [#] The storage efficiency (e.g. 95% or 0.95) to use for the schedule is applied over each time step equal to the sensor resolution. For example, a storage efficiency of 95 percent per (absolute) day, for scheduling a 1-hour resolution sensor, should be passed as a storage efficiency of :math:`0.95^{1/24} = 0.997865`. -2) Shiftable process + For some examples, see the `[POST] /sensors/(id)/schedules/trigger <../api/v3_0.html#post--api-v3_0-sensors-(id)-schedules-trigger>`_ endpoint docs. + +2) For **shiftable processes** .. todo:: A simple algorithm exists, needs integration into FlexMeasures and asset type clarified. -3) Heat pumps - - .. todo:: Also work in progress, needs model for heat loss compensation. +3) For **buffer devices** (e.g. thermal energy storage systems connected to heat pumps), use the same flexibility parameters described above for storage devices. Here are some tips to model a buffer with these parameters: + + - Describe the thermal energy content in kWh or MWh. + - Set ``soc-minima`` to the accumulative usage forecast. + - Set ``roundtrip-efficiency`` to the square of the conversion efficiency. [#]_ + + .. [#] Setting a roundtrip efficiency of higher than 1 is not supported. We plan to implement a separate field for :abbr:`COP (coefficient of performance)` values. In addition, folks who write their own custom scheduler (see :ref:`plugin_customization`) might also require their custom flexibility model. That's no problem, FlexMeasures will let the scheduler decide which flexibility model is relevant and how it should be validated. diff --git a/documentation/changelog.rst b/documentation/changelog.rst index e406f77e6..8f7441b15 100644 --- a/documentation/changelog.rst +++ b/documentation/changelog.rst @@ -9,10 +9,10 @@ v0.14.0 | June XX, 2023 New features ------------- -* Add multiple maxima and minima constraints into `StorageScheduler` [see `PR #680 `_] -* Add CLI command ``flexmeasures add report`` [see `PR #659 `_] -* Add CLI command ``flexmeasures show reporters`` [see `PR #686 `_] -* Add CLI command ``flexmeasures show schedulers`` [see `PR #708 `_] +* Allow setting a storage efficiency using the new ``storage-efficiency`` field when calling `/sensors//schedules/trigger` (POST) through the API (within the ``flex-model`` field), or when calling ``flexmeasures add schedule for-storage`` through the CLI [see `PR #679 `_] +* Allow setting multiple :abbr:`SoC (state of charge)` maxima and minima constraints for the `StorageScheduler`, using the new ``soc-minima`` and ``soc-maxima`` fields when calling `/sensors//schedules/trigger` (POST) through the API (within the ``flex-model`` field) [see `PR #680 `_] +* New CLI command ``flexmeasures add report`` to calculate a custom report from sensor data and save the results to the database, with the option to export them to a CSV or Excel file [see `PR #659 `_] +* New CLI commands ``flexmeasures show reporters`` and ``flexmeasures show schedulers`` to list available reporters and schedulers, respectively, including any defined in registered plugins [see `PR #686 `_ and `PR #708 `_] Bugfixes ----------- @@ -146,7 +146,7 @@ Bugfixes * The CLI command ``flexmeasures show beliefs`` now supports plotting time series data that includes NaN values, and provides better support for plotting multiple sensors that do not share the same unit [see `PR #516 `_ and `PR #539 `_] * Fixed JSON wrapping of return message for `/sensors/data` (GET) [see `PR #543 `_] * Consistent CLI/UI support for asset lat/lng positions up to 7 decimal places (previously the UI rounded to 4 decimal places, whereas the CLI allowed more than 4) [see `PR #522 `_] -* Stop trimming the planning window in response to price availability, which is a problem when SoC targets occur outside of the available price window, by making a simplistic assumption about future prices [see `PR #538 `_] +* Stop trimming the planning window in response to price availability, which is a problem when :abbr:`SoC (state of charge)` targets occur outside of the available price window, by making a simplistic assumption about future prices [see `PR #538 `_] * Faster loading of initial charts and calendar date selection [see `PR #533 `_] Infrastructure / Support @@ -177,7 +177,7 @@ v0.11.3 | November 2, 2022 Bugfixes ----------- -* Fix scheduling with imperfect efficiencies, which resulted in exceeding the device's lower SoC limit. [see `PR #520 `_] +* Fix scheduling with imperfect efficiencies, which resulted in exceeding the device's lower :abbr:`SoC (state of charge)` limit. [see `PR #520 `_] * Fix scheduler for Charge Points when taking into account inflexible devices [see `PR #517 `_] * Prevent rounding asset lat/long positions to 4 decimal places when editing an asset in the UI [see `PR #522 `_] diff --git a/documentation/conf.py b/documentation/conf.py index ca55d2c86..ef945899a 100644 --- a/documentation/conf.py +++ b/documentation/conf.py @@ -46,7 +46,7 @@ "sphinx_rtd_theme", "sphinx.ext.intersphinx", "sphinx.ext.coverage", - "sphinx.ext.imgmath", + "sphinx.ext.mathjax", "sphinx.ext.ifconfig", "sphinx.ext.todo", "sphinx_copybutton", diff --git a/flexmeasures/api/v3_0/sensors.py b/flexmeasures/api/v3_0/sensors.py index 405bcdec5..9f61a67ed 100644 --- a/flexmeasures/api/v3_0/sensors.py +++ b/flexmeasures/api/v3_0/sensors.py @@ -273,6 +273,7 @@ def trigger_schedule( # noqa: C901 To guarantee a minimum SOC in the period prior to 4.00pm, local minima constraints are imposed (via soc-minima) at 2.00pm and 3.00pm, for 15kWh and 20kWh, 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 @@ -306,6 +307,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, diff --git a/flexmeasures/api/v3_0/tests/test_sensor_schedules.py b/flexmeasures/api/v3_0/tests/test_sensor_schedules.py index 53891c677..b2ab198f0 100644 --- a/flexmeasures/api/v3_0/tests/test_sensor_schedules.py +++ b/flexmeasures/api/v3_0/tests/test_sensor_schedules.py @@ -221,6 +221,9 @@ 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 @@ -228,6 +231,9 @@ def test_trigger_and_get_schedule( 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) diff --git a/flexmeasures/api/v3_0/tests/utils.py b/flexmeasures/api/v3_0/tests/utils.py index 2ab606a4e..41d4e8e5a 100644 --- a/flexmeasures/api/v3_0/tests/utils.py +++ b/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: diff --git a/flexmeasures/cli/data_add.py b/flexmeasures/cli/data_add.py index e66467bad..ad425e36b 100755 --- a/flexmeasures/cli/data_add.py +++ b/flexmeasures/cli/data_add.py @@ -49,6 +49,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 ( @@ -1012,11 +1013,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, @@ -1036,6 +1048,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. @@ -1091,6 +1104,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, @@ -1104,6 +1119,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, diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index 32d8b7456..48247ac79 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/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], ) diff --git a/flexmeasures/data/models/planning/storage.py b/flexmeasures/data/models/planning/storage.py index 74191180a..15ae79723 100644 --- a/flexmeasures/data/models/planning/storage.py +++ b/flexmeasures/data/models/planning/storage.py @@ -34,6 +34,7 @@ class StorageScheduler(Scheduler): "equals", "max", "min", + "efficiency", "derivative equals", "derivative max", "derivative min", @@ -73,6 +74,7 @@ def compute(self, skip_validation: bool = False) -> pd.Series | None: soc_minima = self.flex_model.get("soc_minima") soc_maxima = self.flex_model.get("soc_maxima") 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") @@ -126,7 +128,7 @@ def compute(self, skip_validation: bool = False) -> pd.Series | None: down_deviation_prices.loc[start : end - resolution]["event_value"] ] - # Set up device _constraints: only one scheduled flexible device for this EMS (at index 0), plus the forecasted inflexible devices (at indices 1 to n). + # Set up device constraints: only one scheduled flexible device for this EMS (at index 0), plus the forecasted inflexible devices (at indices 1 to n). device_constraints = [ initialize_df(StorageScheduler.COLUMNS, start, end, resolution) for i in range(1 + len(inflexible_device_sensors)) @@ -170,6 +172,9 @@ def compute(self, skip_validation: bool = False) -> pd.Series | None: ) device_constraints[0]["derivative up efficiency"] = roundtrip_efficiency**0.5 + # Apply storage efficiency (accounts for losses over time) + device_constraints[0]["efficiency"] = storage_efficiency + # check that storage constraints are fulfilled if not skip_validation: constraint_violations = validate_storage_constraints( @@ -199,6 +204,7 @@ def compute(self, skip_validation: bool = False) -> pd.Series | None: 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 @@ -268,7 +274,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 diff --git a/flexmeasures/data/models/planning/tests/test_solver.py b/flexmeasures/data/models/planning/tests/test_solver.py index 9baf56cb5..904aafd97 100644 --- a/flexmeasures/data/models/planning/tests/test_solver.py +++ b/flexmeasures/data/models/planning/tests/test_solver.py @@ -13,12 +13,46 @@ validate_storage_constraints, ) from flexmeasures.data.models.planning.utils import initialize_series, initialize_df -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 @@ -62,14 +96,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: @@ -100,6 +138,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() @@ -108,6 +147,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, ) @@ -126,22 +166,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( @@ -188,6 +236,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, }, ) @@ -261,6 +312,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, }, ) @@ -351,6 +405,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(), diff --git a/flexmeasures/data/schemas/scheduling/storage.py b/flexmeasures/data/schemas/scheduling/storage.py index e15843459..25f840900 100644 --- a/flexmeasures/data/schemas/scheduling/storage.py +++ b/flexmeasures/data/schemas/scheduling/storage.py @@ -19,6 +19,33 @@ from flexmeasures.utils.unit_utils import ur +class EfficiencyField(QuantityField): + """Field that deserializes to a Quantity with % units. Must be greater than 0% and less than or equal to 100%. + + Examples: + + >>> ef = EfficiencyField() + >>> ef.deserialize(0.9) + + >>> ef.deserialize("90%") + + >>> ef.deserialize("0%") + Traceback (most recent call last): + ... + marshmallow.exceptions.ValidationError: ['Must be greater than 0 and less than or equal to 1.'] + """ + + def __init__(self, *args, **kwargs): + super().__init__( + "%", + validate=validate.Range( + min=0, max=1, min_inclusive=False, max_inclusive=True + ), + *args, + **kwargs, + ) + + class SOCValueSchema(Schema): """ A point in time with a target value. @@ -66,11 +93,8 @@ class StorageFlexModelSchema(Schema): data_key="soc-unit", ) # todo: allow unit to be set per field, using QuantityField("%", validate=validate.Range(min=0, max=1)) soc_targets = fields.List(fields.Nested(SOCValueSchema()), data_key="soc-targets") - roundtrip_efficiency = QuantityField( - "%", - validate=validate.Range(min=0, max=1, min_inclusive=False, max_inclusive=True), - data_key="roundtrip-efficiency", - ) + roundtrip_efficiency = EfficiencyField(data_key="roundtrip-efficiency") + storage_efficiency = EfficiencyField(data_key="storage-efficiency") prefer_charging_sooner = fields.Bool(data_key="prefer-charging-sooner") def __init__(self, start: datetime, sensor: Sensor, *args, **kwargs): @@ -110,10 +134,10 @@ def post_load_sequence(self, data: dict, **kwargs) -> dict: target["value"] /= 1000.0 data["soc_unit"] = "MWh" - # Convert round-trip efficiency to dimensionless (to the (0,1] range) - if data.get("roundtrip_efficiency") is not None: - data["roundtrip_efficiency"] = ( - data["roundtrip_efficiency"].to(ur.Quantity("dimensionless")).magnitude - ) + # Convert efficiencies to dimensionless (to the (0,1] range) + efficiency_fields = ("storage_efficiency", "roundtrip_efficiency") + for field in efficiency_fields: + if data.get(field) is not None: + data[field] = data[field].to(ur.Quantity("dimensionless")).magnitude return data diff --git a/flexmeasures/data/tests/test_scheduling_jobs.py b/flexmeasures/data/tests/test_scheduling_jobs.py index a488fb98c..cdd2f1a4a 100644 --- a/flexmeasures/data/tests/test_scheduling_jobs.py +++ b/flexmeasures/data/tests/test_scheduling_jobs.py @@ -37,6 +37,10 @@ def test_scheduling_a_battery(db, app, add_battery_assets, setup_test_data): end=end, belief_time=start, resolution=resolution, + flex_model={ + "roundtrip-efficiency": "98%", + "storage-efficiency": 0.999, + }, ) print("Job: %s" % job.id) @@ -57,6 +61,9 @@ def test_scheduling_a_battery(db, app, add_battery_assets, setup_test_data): ) print([v.event_value for v in power_values]) assert len(power_values) == 96 + assert ( + sum(v.event_value for v in power_values) < -0.5 + ), "some cycling should have occurred to make a profit, resulting in overall consumption due to losses" scheduler_specs = { diff --git a/flexmeasures/data/tests/test_scheduling_jobs_fresh_db.py b/flexmeasures/data/tests/test_scheduling_jobs_fresh_db.py index ca125614f..02b966a62 100644 --- a/flexmeasures/data/tests/test_scheduling_jobs_fresh_db.py +++ b/flexmeasures/data/tests/test_scheduling_jobs_fresh_db.py @@ -42,7 +42,12 @@ def test_scheduling_a_charging_station( end=end, belief_time=start, resolution=resolution, - flex_model={"soc-at-start": soc_at_start, "soc-targets": soc_targets}, + flex_model={ + "soc-at-start": soc_at_start, + "soc-targets": soc_targets, + "roundtrip-efficiency": "100%", + "storage-efficiency": 1, + }, ) print("Job: %s" % job.id) diff --git a/flexmeasures/utils/calculations.py b/flexmeasures/utils/calculations.py index f85813fa7..e1dee05ad 100644 --- a/flexmeasures/utils/calculations.py +++ b/flexmeasures/utils/calculations.py @@ -2,6 +2,7 @@ from __future__ import annotations from datetime import timedelta +import math import numpy as np import pandas as pd @@ -37,11 +38,76 @@ def drop_nan_rows(a, b): return d[:, 0], d[:, 1] +def apply_stock_changes_and_losses( + initial: float, + changes: list[float], + storage_efficiency: float | list[float], + how: str = "linear", + decimal_precision: int | None = None, +) -> list[float]: + r"""Assign stock changes and determine losses from storage efficiency. + + The initial stock is exponentially decayed, as with each consecutive (constant-resolution) time step, + some constant percentage of the previous stock remains. For example: + + .. math:: + + 100 \rightarrow 90 \rightarrow 81 \rightarrow 72.9 \rightarrow ... + + For computing the decay of the changes, we make an assumption on how a delta :math:`d` is distributed within a given time step. + In case it happens at a constant rate, this leads to a linear stock change from one time step to the next. + + An :math:`e` is introduced when we apply exponential decay to that. + To see that, imagine we cut one time step in :math:`n` pieces (each with a stock change :math:`\frac{d}{n}` ), + apply the efficiency to each piece :math:`k` (for the corresponding fraction of the time step :math:`k/n`), + and then take the limit :math:`n \rightarrow \infty`: + + .. math:: + + \lim_{n \rightarrow \infty} \sum_{k=0}^{n}{\frac{d}{n} \eta^{k/n}} + + `which is `_: + + .. math:: + + d \cdot \frac{\eta - 1}{e^{\eta}} + + :param initial: initial stock + :param changes: stock change for each step + :param storage_efficiency: ratio of stock left after a step (constant ratio or one per step) + :param how: left, right or linear; how stock changes should be applied, which affects how losses are applied + :param decimal_precision: Optional decimal precision to round off results (useful for tests failing over machine precision) + """ + stocks = [initial] + if not isinstance(storage_efficiency, list): + storage_efficiency = [storage_efficiency] * len(changes) + for d, e in zip(changes, storage_efficiency): + s = stocks[-1] + if e == 1: + next_stock = s + d + elif how == "left": + # First apply the stock change, then apply the losses (i.e. the stock changes on the left side of the time interval in which the losses apply) + next_stock = (s + d) * e + elif how == "right": + # First apply the losses, then apply the stock change (i.e. the stock changes on the right side of the time interval in which the losses apply) + next_stock = s * e + d + elif how == "linear": + # Assume the change happens at a constant rate, leading to a linear stock change, and exponential decay, within the current interval + next_stock = s * e + d * (e - 1) / math.log(e) + else: + raise NotImplementedError(f"Missing implementation for how='{how}'.") + stocks.append(next_stock) + if decimal_precision is not None: + stocks = [round(s, decimal_precision) for s in stocks] + return stocks + + def integrate_time_series( series: pd.Series, initial_stock: float, up_efficiency: float | pd.Series = 1, down_efficiency: float | pd.Series = 1, + storage_efficiency: float | pd.Series = 1, decimal_precision: int | None = None, ) -> pd.Series: """Integrate time series of length n and inclusive="left" (representing a flow) @@ -69,25 +135,42 @@ def integrate_time_series( dtype: float64 """ resolution = pd.to_timedelta(series.index.freq) + storage_efficiency = ( + storage_efficiency + if isinstance(storage_efficiency, pd.Series) + else pd.Series(storage_efficiency, index=series.index) + ) + + # Convert from flow to stock change, applying conversion efficiencies 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] + * ( + up_efficiency[series > 0] + if isinstance(up_efficiency, pd.Series) + else up_efficiency + ) + * (resolution / timedelta(hours=1)) ) - stock_change.loc[series <= 0] = series[series <= 0] / ( - down_efficiency[series <= 0] - if isinstance(down_efficiency, pd.Series) - else down_efficiency + stock_change.loc[series <= 0] = ( + series[series <= 0] + / ( + down_efficiency[series <= 0] + if isinstance(down_efficiency, pd.Series) + else down_efficiency + ) + * (resolution / timedelta(hours=1)) + ) + + stocks = apply_stock_changes_and_losses( + initial_stock, stock_change.tolist(), storage_efficiency.tolist() ) - int_s = pd.concat( + stocks = pd.concat( [ 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, + pd.Series(stocks[1:], index=series.index).shift(1, freq=resolution), ] ) if decimal_precision is not None: - int_s = int_s.round(decimal_precision) - return int_s + stocks = stocks.round(decimal_precision) + return stocks