feat: add Binary constraint to prevent energy losses (at negative pri…

…ces) (#770) Add binary constraint to avoid energy leakages during periods with negative prices. * add highs to requirements Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * docs: add changelog entry Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fix: get results with infeasible termination status instead of RuntimeError Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * feat: add Binary constraint to prevent energy losses, and start new test against negative prices Signed-off-by: F.N. Claessen <felix@seita.nl> * fx: avoid double solving Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * style: fix HiGHS capitalization Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * remove HiGHS from requirements Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * remove dependency Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * add dependency back Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * docs: document how to install HiGHS Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * add HIghs to Dockerfile Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * remove extra lines Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fix typos Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * load solution when termination_condition!=infeasible Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * refactor: run data preparation step in a different method Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * refactor: create run_device_scheduler function to return model and results objects and using it in the device_scheduler function Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * add asserts and docstring in test_battery_solver_day_3 Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * add documentation for the constraints device_up_derivative_sign and device_down_derivative_sign Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * install HiGHS in the CI testing env Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * address some textual changes Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fx CBC capitalization Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fix grammar Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * check if there are results in a more robustly Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * little fixes Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * improve docstring of _prepare Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fix definition of M Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * improve test Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * reorder installation instructions for CBC and HiGHS Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * add production price fixture Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fix loading results twice Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * remove TODO Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * increment StorageScheduler version Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * add changelog entry Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fix conftest Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * update test_hashing with a new hash due to the change in the version of the StorageScheduler Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * feat: in play mode, allow showing any sensor on the asset page (#740) This PR uses the FlexMeasures play mode to allow showing sensors across organisations. When showing simulation results, we want to be able to show results from different scenarios (stored under different accounts) on the same asset page. For example, when comparing a scenario against some other benchmark scenario. * feat: in play mode, allow showing any sensor on the asset page Signed-off-by: F.N. Claessen <felix@seita.nl> * fix checking for config value Signed-off-by: Nicolas Höning <nicolas@seita.nl> * more consistent naming of account variable Signed-off-by: Nicolas Höning <nicolas@seita.nl> * chore: rename variable of accessible sensors for readability Signed-off-by: Nicolas Höning <nicolas@seita.nl> * docs: mention this new possibility in function docstring Signed-off-by: Nicolas Höning <nicolas@seita.nl> * fix: do not fail if sensors are not accessible, log a warning for them as well Signed-off-by: Nicolas Höning <nicolas@seita.nl> * docs: changelog entry Signed-off-by: F.N. Claessen <felix@seita.nl> * docs: changelog warning Signed-off-by: F.N. Claessen <felix@seita.nl> --------- Signed-off-by: F.N. Claessen <felix@seita.nl> Signed-off-by: Nicolas Höning <nicolas@seita.nl> Co-authored-by: Nicolas Höning <nicolas@seita.nl> * fix fixture Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fix: changelog Signed-off-by: F.N. Claessen <felix@seita.nl> * revert (the refactoring part of aa32839): prepend model to the return tuple of the device_scheduler; this prevents duplicating a large docstring at the cost of introducing a breaking change in the function signature of the device_scheduler (but no external code uses that function to the best of my knowledge) (#781) Signed-off-by: F.N. Claessen <felix@seita.nl> * clean code and latitude, longitude and account_id to get_or_create_generic_asset Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * add highspy to test.txt and test.in Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fix: typos Signed-off-by: F.N. Claessen <felix@seita.nl> * fix: test that needs to get the battery object Signed-off-by: F.N. Claessen <felix@seita.nl> * fix: typo Signed-off-by: F.N. Claessen <felix@seita.nl> * docs: clarify choice of SoC values in test Signed-off-by: F.N. Claessen <felix@seita.nl> * fix: test case 2 says it tests for no oscillation, which is only guaranteed when prefer-charging-sooner is true; otherwise, the solver will be completely indifferent to oscilations during the period with negative prices Signed-off-by: F.N. Claessen <felix@seita.nl> * style: streamline variable names in test Signed-off-by: F.N. Claessen <felix@seita.nl> * revert 48fd692: unnecessary after #695 Signed-off-by: F.N. Claessen <felix@seita.nl> * style: black Signed-off-by: F.N. Claessen <felix@seita.nl> * remove installation of highspy as it's already being installed through the test requirements. Signed-off-by: Victor Garcia Reolid <victor@seita.nl> * fix (docs): comment out cross reference to masked documentation page Signed-off-by: F.N. Claessen <felix@seita.nl> --------- Signed-off-by: Victor Garcia Reolid <victor@seita.nl> Signed-off-by: F.N. Claessen <felix@seita.nl> Signed-off-by: Victor <victor@seita.nl> Signed-off-by: Nicolas Höning <nicolas@seita.nl> Co-authored-by: Victor Garcia Reolid <victor@seita.nl> Co-authored-by: Nicolas Höning <nicolas@seita.nl>
FlexMeasures · Aug 3, 2023 · b1a85e8 · b1a85e8
1 parent 3bdff2d
commit b1a85e8
Show file tree

Hide file tree

Showing 15 changed files with 389 additions and 70 deletions.
diff --git a/.github/workflows/lint-and-test.yml b/.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'

diff --git a/documentation/changelog.rst b/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
 ----------------------
 

diff --git a/documentation/concepts/benefits.rst b/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`
diff --git a/documentation/concepts/device_scheduler.rst b/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
 ^^^^^^^^^^^^^^^^^
 

diff --git a/documentation/dev/introduction.rst b/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
 

diff --git a/flexmeasures/conftest.py b/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")

diff --git a/flexmeasures/data/models/generic_assets.py b/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
 
 

diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/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
diff --git a/flexmeasures/data/models/planning/storage.py b/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,