feat: add Binary constraint to prevent energy losses (at negative prices)

.github/workflows/lint-and-test.yml

@@ -49,7 +49,6 @@ jobs:
             ${{ runner.os }}-pip-
       - run: |
           ci/setup-postgres.sh
-          sudo apt-get -y install coinor-cbc
       - name: Install FlexMeasures & exact dependencies for tests
         run: make install-for-test
         if: github.event_name == 'push' && steps.cache.outputs.cache-hit != 'true'

documentation/changelog.rst

@@ -8,6 +8,8 @@ v0.15.0 | July XX, 2023
 
 .. warning:: Upgrading to this version requires running ``flexmeasures db upgrade`` (you can create a backup first with ``flexmeasures db-ops dump``).
 
+.. warning:: Upgrading to this version requires installing the LP/MILP solver HiGHS using ``pip install highspy``.
+
 .. warning:: If your server is running in play mode (``FLEXMEASURES_MODE = "play"``), users will be able to see sensor data from any account [see `PR #740 <https://www.github.com/FlexMeasures/flexmeasures/pull/740>`_].
 
 New features
@@ -27,6 +29,8 @@ New features
 Bugfixes
 -----------
 
+* Add binary constraint to avoid energy leakages during periods with negative prices [see `PR #770 <https://www.github.com/FlexMeasures/flexmeasures/pull/770>`_]
+
 Infrastructure / Support
 ----------------------
 

documentation/concepts/benefits.rst

@@ -34,4 +34,6 @@ The platform operator (as ESCo or Aggregator) and asset owners can share the pro
 FlexMeasures plans on providing basic accounting for this.
 
 
-.. note:: Read more on flexibility opportunities and activations, as well as profit sharing on :ref:`benefits_of_flex`
+..
+    <This cross reference has been commented out while we rewrite the benefits_of_flex page to fit out new data model and UI>
+    note:: Read more on flexibility opportunities and activations, as well as profit sharing on :ref:`benefits_of_flex`

documentation/concepts/device_scheduler.rst

@@ -48,6 +48,7 @@ Symbol                              Variable in the Code
 :math:`P^{ems}_{min}(j)`              ems_derivative_min                                 Minimum flow of the EMS during time period :math:`j`.
 :math:`P^{ems}_{max}(j)`              ems_derivative_max                                 Maximum flow of the EMS during time period :math:`j`.
 :math:`Commitment(c,j)`               commitment_quantity                                Commitment c (at EMS level) over time step :math:`j`.
+:math:`M`                             M                                                  Large constant number, upper bound of :math:`Power_{up}(d,j)` and :math:`|Power_{down}(d,j)|`
 ================================ ================================================ ==============================================================================================================  
 
 
@@ -62,6 +63,7 @@ Symbol                              Variable in the Code
 :math:`P_{up}(d,j)`                   device_power_up                                    Upwards power of device :math:`d` during time period :math:`j`.
 :math:`P_{down}(d,j)`                 device_power_down                                  Downwards power of device :math:`d` during time period :math:`j`.
 :math:`P^{ems}(j)`                    ems_power                                          Aggregated power of all the devices during time period :math:`j`.
+:math:`\sigma(d,j)`                   device_power_sign                                  Upwards power activation if :math:`\sigma(d,j)=1`, downwards power activation otherwise.
 ================================ ================================================ ==============================================================================================================  
 
 Cost function
@@ -154,6 +156,22 @@ Device bounds
     0 \leq P_{up}(d,j)\leq max(P_{max}(d,j),0)
 
 
+Upwards/Downwards activation selection
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Avoid simultaneous upwards and downwards activation during the same time period.
+
+.. math:: 
+  :name: device_up_derivative_sign
+
+    P_{up}(d,j) \leq M \cdot \sigma(d,j)
+
+.. math:: 
+  :name: device_down_derivative_sign
+
+    -P_{down}(d,j) \leq M \cdot (1-\sigma(d,j))
+
+
 Grid constraints
 ^^^^^^^^^^^^^^^^^
 

documentation/dev/introduction.rst

@@ -48,7 +48,14 @@ Go into the ``flexmeasures`` folder and install all dependencies including the o
    $ cd flexmeasures
    $ make install-for-dev
 
-:ref:`Install the LP solver <install-lp-solver>`. On Unix the Cbc LP solver can be installed with:
+:ref:`Install the LP solver <install-lp-solver>`. On Linux, the HiGHS solver can be installed with:
+
+.. code-block:: bash
+
+   $ pip install highspy
+
+
+Alternatively, the CBC solver can be installed with:
 
 .. code-block:: bash
 

flexmeasures/conftest.py

@@ -249,24 +249,27 @@ def create_test_markets(db) -> dict[str, Sensor]:
         name="epex",
         generic_asset_type=day_ahead,
     )
-    epex_da = Sensor(
-        name="epex_da",
-        generic_asset=epex,
-        event_resolution=timedelta(hours=1),
-        unit="EUR/MWh",
-        knowledge_horizon=(
-            x_days_ago_at_y_oclock,
-            {"x": 1, "y": 12, "z": "Europe/Paris"},
-        ),
-        attributes=dict(
-            daily_seasonality=True,
-            weekly_seasonality=True,
-            yearly_seasonality=True,
-        ),
-    )
-    db.session.add(epex_da)
-    db.session.flush()  # assign an id, so it can be used to set a market_id attribute on a GenericAsset or Sensor
-    return {"epex_da": epex_da}
+    price_sensors = {}
+    for sensor_name in ("epex_da", "epex_da_production"):
+        price_sensor = Sensor(
+            name=sensor_name,
+            generic_asset=epex,
+            event_resolution=timedelta(hours=1),
+            unit="EUR/MWh",
+            knowledge_horizon=(
+                x_days_ago_at_y_oclock,
+                {"x": 1, "y": 12, "z": "Europe/Paris"},
+            ),
+            attributes=dict(
+                daily_seasonality=True,
+                weekly_seasonality=True,
+                yearly_seasonality=True,
+            ),
+        )
+        db.session.add(price_sensor)
+        price_sensors[sensor_name] = price_sensor
+    db.session.flush()  # assign an id, so the markets can be used to set a market_id attribute on a GenericAsset or Sensor
+    return price_sensors
 
 
 @pytest.fixture(scope="module")
@@ -526,7 +529,7 @@ def create_beliefs(db: SQLAlchemy, setup_markets, setup_sources) -> int:
 def add_market_prices(
     db: SQLAlchemy, setup_assets, setup_markets, setup_sources
 ) -> dict[str, Sensor]:
-    """Add two days of market prices for the EPEX day-ahead market."""
+    """Add three days of market prices for the EPEX day-ahead market."""
 
     # one day of test data (one complete sine curve)
     time_slots = initialize_index(
@@ -568,7 +571,46 @@ def add_market_prices(
         for dt, val in zip(time_slots, values)
     ]
     db.session.add_all(day2_beliefs)
-    return {"epex_da": setup_markets["epex_da"]}
+
+    # the third day of test data (8 hours with negative prices, followed by 16 expensive hours)
+    time_slots = initialize_index(
+        start=pd.Timestamp("2015-01-03").tz_localize("Europe/Amsterdam"),
+        end=pd.Timestamp("2015-01-04").tz_localize("Europe/Amsterdam"),
+        resolution="1H",
+    )
+
+    # consumption prices
+    values = [-10] * 8 + [100] * 16
+    day3_beliefs = [
+        TimedBelief(
+            event_start=dt,
+            belief_horizon=timedelta(hours=0),
+            event_value=val,
+            source=setup_sources["Seita"],
+            sensor=setup_markets["epex_da"],
+        )
+        for dt, val in zip(time_slots, values)
+    ]
+    db.session.add_all(day3_beliefs)
+
+    # production prices = consumption prices - 40
+    values = [-50] * 8 + [60] * 16
+    day3_beliefs_production = [
+        TimedBelief(
+            event_start=dt,
+            belief_horizon=timedelta(hours=0),
+            event_value=val,
+            source=setup_sources["Seita"],
+            sensor=setup_markets["epex_da_production"],
+        )
+        for dt, val in zip(time_slots, values)
+    ]
+    db.session.add_all(day3_beliefs_production)
+
+    yield {
+        "epex_da": setup_markets["epex_da"],
+        "epex_da_production": setup_markets["epex_da_production"],
+    }
 
 
 @pytest.fixture(scope="module")

flexmeasures/data/models/generic_assets.py

@@ -581,6 +581,7 @@ def create_generic_asset(generic_asset_type: str, **kwargs) -> GenericAsset:
     ).one_or_none()
     if generic_asset_type is None:
         raise ValueError(f"Cannot find GenericAssetType {asset_type_name} in database.")
+
     new_generic_asset = GenericAsset(
         name=kwargs["name"],
         generic_asset_type_id=generic_asset_type.id,
@@ -591,6 +592,7 @@ def create_generic_asset(generic_asset_type: str, **kwargs) -> GenericAsset:
             setattr(new_generic_asset, arg, kwargs[arg])
     db.session.add(new_generic_asset)
     db.session.flush()  # generates the pkey for new_generic_asset
+
     return new_generic_asset
 
 

flexmeasures/data/models/planning/linear_optimization.py

@@ -12,6 +12,7 @@
     Reals,
     NonNegativeReals,
     NonPositiveReals,
+    Binary,
     Constraint,
     Objective,
     minimize,
@@ -33,7 +34,7 @@ def device_scheduler(  # noqa C901
     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]:
+) -> Tuple[List[pd.Series], float, SolverResults, ConcreteModel]:
     """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,
     and with multiple market commitments on the EMS level.
@@ -92,6 +93,15 @@ def device_scheduler(  # noqa C901
                 % (resolution, resolution_c)
             )
 
+    # Compute a good value for M
+    M = 0.1
+    for device_constraint in device_constraints:
+        M = max(
+            M,
+            device_constraint["derivative max"].max(),
+            -device_constraint["derivative min"].min(),
+        )
+
     # Turn prices per commitment into prices per commitment flow
     if len(commitment_downwards_deviation_price) != 0:
         if all(
@@ -234,6 +244,7 @@ def device_derivative_up_efficiency(m, d, j):
         model.d, model.j, domain=NonPositiveReals, initialize=0
     )
     model.device_power_up = Var(model.d, model.j, domain=NonNegativeReals, initialize=0)
+    model.device_power_sign = Var(model.d, model.j, domain=Binary, initialize=0)
     model.commitment_downwards_deviation = Var(
         model.c, model.j, domain=NonPositiveReals, initialize=0
     )
@@ -287,6 +298,14 @@ def device_up_derivative_bounds(m, d, j):
             max(0, m.device_derivative_max[d, j]),
         )
 
+    def device_up_derivative_sign(m, d, j):
+        """Derivative up if sign points up, derivative not up if sign points down."""
+        return m.device_power_up[d, j] <= M * m.device_power_sign[d, j]
+
+    def device_down_derivative_sign(m, d, j):
+        """Derivative down if sign points down, derivative not down if sign points up."""
+        return -m.device_power_down[d, j] <= M * (1 - m.device_power_sign[d, j])
+
     def ems_derivative_bounds(m, j):
         return m.ems_derivative_min[j], sum(m.ems_power[:, j]), m.ems_derivative_max[j]
 
@@ -319,6 +338,12 @@ def device_derivative_equalities(m, d, j):
     model.device_power_up_bounds = Constraint(
         model.d, model.j, rule=device_up_derivative_bounds
     )
+    model.device_power_up_sign = Constraint(
+        model.d, model.j, rule=device_up_derivative_sign
+    )
+    model.device_power_down_sign = Constraint(
+        model.d, model.j, rule=device_down_derivative_sign
+    )
     model.ems_power_bounds = Constraint(model.j, rule=ems_derivative_bounds)
     model.ems_power_commitment_equalities = Constraint(
         model.j, rule=ems_flow_commitment_equalities
@@ -366,4 +391,4 @@ def cost_function(m):
     # model.display()
     # print(results.solver.termination_condition)
     # print(planned_costs)
-    return planned_power_per_device, planned_costs, results
+    return planned_power_per_device, planned_costs, results, model

flexmeasures/data/models/planning/storage.py

@@ -26,7 +26,7 @@
 
 class StorageScheduler(Scheduler):
 
-    __version__ = "1"
+    __version__ = "2"
     __author__ = "Seita"
 
     COLUMNS = [
@@ -50,13 +50,16 @@ def compute_schedule(self) -> pd.Series | None:
 
         return self.compute()
 
-    def compute(self, skip_validation: bool = False) -> pd.Series | None:
-        """Schedule a battery or Charge Point based directly on the latest beliefs regarding market prices within the specified time window.
-        For the resulting consumption schedule, consumption is defined as positive values.
+    def _prepare(self, skip_validation: bool = False) -> tuple:
+        """This function prepares the required data to compute the schedule:
+            - price data
+            - device constraint
+            - ems constraints
 
         :param skip_validation: If True, skip validation of constraints specified in the data.
-        :returns:               The computed schedule.
+        :returns:               Input data for the scheduler
         """
+
         if not self.config_deserialized:
             self.deserialize_config()
 
@@ -201,7 +204,41 @@ def compute(self, skip_validation: bool = False) -> pd.Series | None:
             ems_constraints["derivative min"] = ems_capacity * -1
             ems_constraints["derivative max"] = ems_capacity
 
-        ems_schedule, expected_costs, scheduler_results = device_scheduler(
+        return (
+            sensor,
+            start,
+            end,
+            resolution,
+            soc_at_start,
+            device_constraints,
+            ems_constraints,
+            commitment_quantities,
+            commitment_downwards_deviation_price,
+            commitment_upwards_deviation_price,
+        )
+
+    def compute(self, skip_validation: bool = False) -> pd.Series | None:
+        """Schedule a battery or Charge Point based directly on the latest beliefs regarding market prices within the specified time window.
+        For the resulting consumption schedule, consumption is defined as positive values.
+
+        :param skip_validation: If True, skip validation of constraints specified in the data.
+        :returns:               The computed schedule.
+        """
+
+        (
+            sensor,
+            start,
+            end,
+            resolution,
+            soc_at_start,
+            device_constraints,
+            ems_constraints,
+            commitment_quantities,
+            commitment_downwards_deviation_price,
+            commitment_upwards_deviation_price,
+        ) = self._prepare(skip_validation=skip_validation)
+
+        ems_schedule, expected_costs, scheduler_results, _ = device_scheduler(
             device_constraints,
             ems_constraints,
             commitment_quantities,