fix: relax decimal resolution in constraint validation (#731)

* fix: consider decimal precision when validating equality constraints

Signed-off-by: Victor Garcia Reolid <victor@seita.nl>

* test: add test cases to check for the equality validation considering decimal resolution

Signed-off-by: Victor Garcia Reolid <victor@seita.nl>

* fix: add self.round_to_decimals to create_constraint_violations_message

Signed-off-by: Victor Garcia Reolid <victor@seita.nl>

* fix: fix and modernize type annotations

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

* fix(docs): update docstring with split and new parameters

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

* docs: changelog entry

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>
Co-authored-by: F.N. Claessen <felix@seita.nl>

documentation/changelog.rst

@@ -16,6 +16,15 @@ Infrastructure / Support
 ----------------------
 
 
+v0.14.1 | June XX, 2023
+============================
+
+Bugfixes
+-----------
+
+* Relax constraint validation of `StorageScheduler` to accommodate violations caused by floating point precision [see `PR #731 <https://www.github.com/FlexMeasures/flexmeasures/pull/731>`_]
+
+
 v0.14.0 | June 15, 2023
 ============================
 

flexmeasures/data/models/planning/storage.py

@@ -3,7 +3,6 @@
 import re
 import copy
 from datetime import datetime, timedelta
-from typing import List, Dict
 
 import pandas as pd
 import numpy as np
@@ -187,7 +186,13 @@ def compute(self, skip_validation: bool = False) -> pd.Series | None:
 
             if len(constraint_violations) > 0:
                 # TODO: include hints from constraint_violations into the error message
-                raise ValueError("The input data yields an infeasible problem.")
+                message = create_constraint_violations_message(
+                    constraint_violations, self.round_to_decimals
+                )
+                raise ValueError(
+                    "The input data yields an infeasible problem. Constraint validation has found the following issues:\n"
+                    + message
+                )
 
         # Set up EMS constraints
         ems_constraints = initialize_df(
@@ -383,8 +388,25 @@ def ensure_soc_min_max(self):
                     )
 
 
+def create_constraint_violations_message(constraint_violations: list) -> str:
+    """Create a human-readable message with the constraint_violations.
+
+    :param constraint_violations: list with the constraint violations
+    :return: human-readable message
+    """
+    message = ""
+
+    for c in constraint_violations:
+        message += f"t={c['dt']} | {c['violation']}\n"
+
+    if len(message) > 1:
+        message = message[:-1]
+
+    return message
+
+
 def build_device_soc_values(
-    soc_values: List[Dict[str, datetime | float]] | pd.Series,
+    soc_values: list[dict[str, datetime | float]] | pd.Series,
     soc_at_start: float,
     start_of_schedule: datetime,
     end_of_schedule: datetime,
@@ -453,9 +475,9 @@ def add_storage_constraints(
     end: datetime,
     resolution: timedelta,
     soc_at_start: float,
-    soc_targets: List[Dict[str, datetime | float]] | pd.Series | None,
-    soc_maxima: List[Dict[str, datetime | float]] | pd.Series | None,
-    soc_minima: List[Dict[str, datetime | float]] | pd.Series | None,
+    soc_targets: list[dict[str, datetime | float]] | pd.Series | None,
+    soc_maxima: list[dict[str, datetime | float]] | pd.Series | None,
+    soc_minima: list[dict[str, datetime | float]] | pd.Series | None,
     soc_max: float,
     soc_min: float,
 ) -> pd.DataFrame:
@@ -578,25 +600,33 @@ def validate_storage_constraints(
     # 1) min >= soc_min
     soc_min = (soc_min - soc_at_start) * timedelta(hours=1) / resolution
     _constraints["soc_min(t)"] = soc_min
-    constraint_violations += validate_constraint(_constraints, "soc_min(t) <= min(t)")
+    constraint_violations += validate_constraint(
+        _constraints, "soc_min(t)", "<=", "min(t)"
+    )
 
     # 2) max <= soc_max
     soc_max = (soc_max - soc_at_start) * timedelta(hours=1) / resolution
     _constraints["soc_max(t)"] = soc_max
-    constraint_violations += validate_constraint(_constraints, "max(t) <= soc_max(t)")
+    constraint_violations += validate_constraint(
+        _constraints, "max(t)", "<=", "soc_max(t)"
+    )
 
     ########################################
     # B. Validation in the same time frame #
     ########################################
 
     # 1) min <= max
-    constraint_violations += validate_constraint(_constraints, "min(t) <= max(t)")
+    constraint_violations += validate_constraint(_constraints, "min(t)", "<=", "max(t)")
 
     # 2) min <= equals
-    constraint_violations += validate_constraint(_constraints, "min(t) <= equals(t)")
+    constraint_violations += validate_constraint(
+        _constraints, "min(t)", "<=", "equals(t)"
+    )
 
     # 3) equals <= max
-    constraint_violations += validate_constraint(_constraints, "equals(t) <= max(t)")
+    constraint_violations += validate_constraint(
+        _constraints, "equals(t)", "<=", "max(t)"
+    )
 
     ##########################################
     # C. Validation in different time frames #
@@ -609,32 +639,38 @@ def validate_storage_constraints(
 
     # 1) equals(t) - equals(t-1) <= derivative_max(t)
     constraint_violations += validate_constraint(
-        _constraints, "equals(t) - equals(t-1) <= derivative_max(t) * factor_w_wh(t)"
+        _constraints,
+        "equals(t) - equals(t-1)",
+        "<=",
+        "derivative_max(t) * factor_w_wh(t)",
     )
 
     # 2) derivative_min(t) <= equals(t) - equals(t-1)
     constraint_violations += validate_constraint(
-        _constraints, "derivative_min(t) * factor_w_wh(t) <= equals(t) - equals(t-1)"
+        _constraints,
+        "derivative_min(t) * factor_w_wh(t)",
+        "<=",
+        "equals(t) - equals(t-1)",
     )
 
     # 3) min(t) - max(t-1) <= derivative_max(t)
     constraint_violations += validate_constraint(
-        _constraints, "min(t) - max(t-1) <= derivative_max(t) * factor_w_wh(t)"
+        _constraints, "min(t) - max(t-1)", "<=", "derivative_max(t) * factor_w_wh(t)"
     )
 
     # 4) max(t) - min(t-1) >= derivative_min(t)
     constraint_violations += validate_constraint(
-        _constraints, "derivative_min(t) * factor_w_wh(t) <= max(t) - min(t-1)"
+        _constraints, "derivative_min(t) * factor_w_wh(t)", "<=", "max(t) - min(t-1)"
     )
 
     # 5) equals(t) - max(t-1) <= derivative_max(t)
     constraint_violations += validate_constraint(
-        _constraints, "equals(t) - max(t-1) <= derivative_max(t) * factor_w_wh(t)"
+        _constraints, "equals(t) - max(t-1)", "<=", "derivative_max(t) * factor_w_wh(t)"
     )
 
     # 6) derivative_min(t) <= equals(t) - min(t-1)
     constraint_violations += validate_constraint(
-        _constraints, "derivative_min(t) * factor_w_wh(t) <= equals(t) - min(t-1)"
+        _constraints, "derivative_min(t) * factor_w_wh(t)", "<=", "equals(t) - min(t-1)"
     )
 
     return constraint_violations
@@ -658,33 +694,85 @@ def get_pattern_match_word(word: str) -> str:
     return regex + re.escape(word) + regex
 
 
+def sanitize_expression(expression: str, columns: list) -> tuple[str, list]:
+    """Wrap column in commas to accept arbitrary column names (e.g. with spaces).
+
+    :param expression: expression to sanitize
+    :param columns: list with the name of the columns of the input data for the expression.
+    :return: sanitized expression and columns (variables) used in the expression
+    """
+
+    _expression = copy.copy(expression)
+    columns_involved = []
+
+    for column in columns:
+
+        if re.search(get_pattern_match_word(column), _expression):
+            columns_involved.append(column)
+
+        _expression = re.sub(get_pattern_match_word(column), f"`{column}`", _expression)
+
+    return _expression, columns_involved
+
+
 def validate_constraint(
-    constraints_df: pd.DataFrame, constraint_expression: str
+    constraints_df: pd.DataFrame,
+    lhs_expression: str,
+    inequality: str,
+    rhs_expression: str,
+    round_to_decimals: int | None = 6,
 ) -> list[dict]:
     """Validate the feasibility of a given set of constraints.
 
-    :param constraints_df: DataFrame with the constraints
-    :param constraint_expression: inequality expression following pd.eval format.
-                                  No need to use the syntax `column` to reference
-                                  column, just use the column name.
-    :return: List of constraint violations, specifying their time, constraint and violation.
+    :param constraints_df:      DataFrame with the constraints
+    :param lhs_expression:      left-hand side of the inequality expression following pd.eval format.
+                                No need to use the syntax `column` to reference
+                                column, just use the column name.
+    :param inequality:          inequality operator, one of ('<=', '<', '>=', '>', '==', '!=').
+    :param rhs_expression:      right-hand side of the inequality expression following pd.eval format.
+                                No need to use the syntax `column` to reference
+                                column, just use the column name.
+    :param round_to_decimals:   Number of decimals to round off to before validating constraints.
+    :return:                    List of constraint violations, specifying their time, constraint and violation.
     """
 
-    columns_involved = []
+    constraint_expression = f"{lhs_expression} {inequality} {rhs_expression}"
 
-    eval_expression = copy.copy(constraint_expression)
+    constraints_df_columns = list(constraints_df.columns)
 
-    for column in constraints_df.columns:
-        if re.search(get_pattern_match_word(column), eval_expression):
-            columns_involved.append(column)
+    lhs_expression, columns_lhs = sanitize_expression(
+        lhs_expression, constraints_df_columns
+    )
+    rhs_expression, columns_rhs = sanitize_expression(
+        rhs_expression, constraints_df_columns
+    )
 
-        eval_expression = re.sub(
-            get_pattern_match_word(column), f"`{column}`", eval_expression
-        )
+    columns_involved = columns_lhs + columns_rhs
+
+    lhs = constraints_df.fillna(0).eval(lhs_expression).round(round_to_decimals)
+    rhs = constraints_df.fillna(0).eval(rhs_expression).round(round_to_decimals)
+
+    condition = None
+
+    inequality = inequality.strip()
+
+    if inequality == "<=":
+        condition = lhs <= rhs
+    elif inequality == "<":
+        condition = lhs < rhs
+    elif inequality == ">=":
+        condition = lhs >= rhs
+    elif inequality == ">":
+        condition = lhs > rhs
+    elif inequality == "==":
+        condition = lhs == rhs
+    elif inequality == "!=":
+        condition = lhs != rhs
+    else:
+        raise ValueError(f"Inequality `{inequality} not supported.")
 
     time_condition_fails = constraints_df.index[
-        ~constraints_df.fillna(0).eval(eval_expression)
-        & ~constraints_df[columns_involved].isna().any(axis=1)
+        ~condition & ~constraints_df[columns_involved].isna().any(axis=1)
     ]
 
     constraint_violations = []
@@ -695,7 +783,7 @@ def validate_constraint(
         for column in constraints_df.columns:
             value_replaced = re.sub(
                 get_pattern_match_word(column),
-                f"{column} [{constraints_df.loc[dt, column]}]",
+                f"{column} [{constraints_df.loc[dt, column]}] ",
                 value_replaced,
             )
 
@@ -730,7 +818,7 @@ def prepend_serie(serie: pd.Series, value) -> pd.Series:
 ####################
 @deprecated(build_device_soc_values, "0.14")
 def build_device_soc_targets(
-    targets: List[Dict[str, datetime | float]] | pd.Series,
+    targets: list[dict[str, datetime | float]] | pd.Series,
     soc_at_start: float,
     start_of_schedule: datetime,
     end_of_schedule: datetime,

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

@@ -640,6 +640,26 @@ def test_add_storage_constraints(
 @pytest.mark.parametrize(
     "value_min1, value_equals1, value_max1, value_min2, value_equals2, value_max2, expected_constraint_type_violations",
     [
+        (1, np.nan, 9, 1, np.nan, 9, []),  # base case
+        (1, np.nan, 10, 1, np.nan, 10, []),  # exact equality
+        (
+            1,
+            np.nan,
+            10 + 0.5e-6,
+            1,
+            np.nan,
+            10,
+            [],
+        ),  # equality considering the precision (6 decimal figures)
+        (
+            1,
+            np.nan,
+            10 + 1e-5,
+            1,
+            np.nan,
+            10,
+            ["max(t) <= soc_max(t)"],
+        ),  # difference of 0.5e-5 > 1e-6
         (1, np.nan, 9, 2, np.nan, 20, ["max(t) <= soc_max(t)"]),
         (-1, np.nan, 9, 1, np.nan, 9, ["soc_min(t) <= min(t)"]),
         (1, 10, 9, 1, np.nan, 9, ["equals(t) <= max(t)"]),