Fall back scheduler heuristics (#267)

Introduce a fallback policy for charging schedules of batteries and Charge Points, in cases where the solver is presented with an infeasible problem.


* Implement and test fallback policy for infeasible scheduler results

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

* Reform old Asset properties into Sensor properties

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

* Simplify

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

flexmeasures/api/v1/implementations.py

@@ -274,15 +274,15 @@ def create_connection_and_value_groups(  # noqa: C901
                 return unrecognized_connection_group()
 
             # Validate the sign of the values (following USEF specs with positive consumption and negative production)
-            if sensor.get_attribute("is_pure_consumer") and any(
+            if sensor.get_attribute("is_strictly_non_positive") and any(
                 v < 0 for v in value_group
             ):
                 extra_info = (
                     "Connection %s is registered as a pure consumer and can only receive non-negative values."
                     % sensor.entity_address
                 )
                 return power_value_too_small(extra_info)
-            elif sensor.get_attribute("is_pure_producer") and any(
+            elif sensor.get_attribute("is_strictly_non_negative") and any(
                 v > 0 for v in value_group
             ):
                 extra_info = (

flexmeasures/api/v2_0/implementations/sensors.py

@@ -328,15 +328,15 @@ def post_power_data(
                 return unrecognized_connection_group()
 
             # Validate the sign of the values (following USEF specs with positive consumption and negative production)
-            if sensor.get_attribute("is_pure_consumer") and any(
+            if sensor.get_attribute("is_strictly_non_positive") and any(
                 v < 0 for v in event_values
             ):
                 extra_info = (
                     "Connection %s is registered as a pure consumer and can only receive non-negative values."
                     % sensor.entity_address
                 )
                 return power_value_too_small(extra_info)
-            elif sensor.get_attribute("is_pure_producer") and any(
+            elif sensor.get_attribute("is_strictly_non_negative") and any(
                 v > 0 for v in event_values
             ):
                 extra_info = (

flexmeasures/data/models/planning/battery.py

@@ -10,6 +10,7 @@
     initialize_series,
     add_tiny_price_slope,
     get_prices,
+    fallback_charging_policy,
 )
 
 
@@ -93,13 +94,19 @@ def schedule_battery(
     columns = ["derivative max", "derivative min"]
     ems_constraints = initialize_df(columns, start, end, resolution)
 
-    ems_schedule, expected_costs = device_scheduler(
+    ems_schedule, expected_costs, scheduler_results = device_scheduler(
         device_constraints,
         ems_constraints,
         commitment_quantities,
         commitment_downwards_deviation_price,
         commitment_upwards_deviation_price,
     )
-    battery_schedule = ems_schedule[0]
+    if scheduler_results.solver.termination_condition == "infeasible":
+        # Fallback policy if the problem was unsolvable
+        battery_schedule = fallback_charging_policy(
+            sensor, device_constraints[0], start, end, resolution
+        )
+    else:
+        battery_schedule = ems_schedule[0]
 
     return battery_schedule

flexmeasures/data/models/planning/charging_station.py

@@ -10,6 +10,7 @@
     initialize_series,
     add_tiny_price_slope,
     get_prices,
+    fallback_charging_policy,
 )
 
 
@@ -82,13 +83,14 @@ def schedule_charging_station(
     ) - soc_at_start * (
         timedelta(hours=1) / resolution
     )  # Lacking information about the battery's nominal capacity, we use the highest target value as the maximum state of charge
-    if sensor.get_attribute("is_pure_consumer", False):
+
+    if sensor.get_attribute("is_strictly_non_positive"):
         device_constraints[0]["derivative min"] = 0
     else:
         device_constraints[0]["derivative min"] = (
             sensor.get_attribute("capacity_in_mw") * -1
         )
-    if sensor.get_attribute("is_pure_producer", False):
+    if sensor.get_attribute("is_strictly_non_negative"):
         device_constraints[0]["derivative max"] = 0
     else:
         device_constraints[0]["derivative max"] = sensor.get_attribute("capacity_in_mw")
@@ -97,13 +99,19 @@ def schedule_charging_station(
     columns = ["derivative max", "derivative min"]
     ems_constraints = initialize_df(columns, start, end, resolution)
 
-    ems_schedule, expected_costs = device_scheduler(
+    ems_schedule, expected_costs, scheduler_results = device_scheduler(
         device_constraints,
         ems_constraints,
         commitment_quantities,
         commitment_downwards_deviation_price,
         commitment_upwards_deviation_price,
     )
-    charging_station_schedule = ems_schedule[0]
+    if scheduler_results.solver.termination_condition == "infeasible":
+        # Fallback policy if the problem was unsolvable
+        charging_station_schedule = fallback_charging_policy(
+            sensor, device_constraints[0], start, end, resolution
+        )
+    else:
+        charging_station_schedule = ems_schedule[0]
 
     return charging_station_schedule

flexmeasures/data/models/planning/solver.py

@@ -16,7 +16,7 @@
 )
 from pyomo.environ import UnknownSolver  # noqa F401
 from pyomo.environ import value
-from pyomo.opt import SolverFactory
+from pyomo.opt import SolverFactory, SolverResults
 
 from flexmeasures.data.models.planning.utils import initialize_series
 
@@ -29,7 +29,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]],
-) -> Tuple[List[pd.Series], float]:
+) -> Tuple[List[pd.Series], float, SolverResults]:
     """Schedule devices given constraints on a device and EMS level, and given a list of commitments by the EMS.
     The commitments are assumed to be with regards to the flow of energy to the device (positive for consumption,
     negative for production). The solver minimises the costs of deviating from the commitments.
@@ -223,7 +223,9 @@ def cost_function(m):
     model.costs = Objective(rule=cost_function, sense=minimize)
 
     # Solve
-    SolverFactory(current_app.config.get("FLEXMEASURES_LP_SOLVER")).solve(model)
+    results = SolverFactory(current_app.config.get("FLEXMEASURES_LP_SOLVER")).solve(
+        model
+    )
 
     planned_costs = value(model.costs)
     planned_power_per_device = []
@@ -239,6 +241,7 @@ def cost_function(m):
         )
 
     # model.pprint()
+    # print(results.solver.termination_condition)
     # print(planned_costs)
     # input()
-    return planned_power_per_device, planned_costs
+    return planned_power_per_device, planned_costs, results

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

@@ -78,14 +78,16 @@ def test_battery_solver_day_2(add_battery_assets):
 def test_charging_station_solver_day_2(target_soc, charging_station_name):
     """Starting with a state of charge 1 kWh, within 2 hours we should be able to reach
     any state of charge in the range [1, 5] kWh for a unidirectional station,
-    or [0, 5] for a bidirectional station."""
+    or [0, 5] for a bidirectional station, given a charging capacity of 2 kW.
+    """
     soc_at_start = 1
     duration_until_target = timedelta(hours=2)
 
     epex_da = Sensor.query.filter(Sensor.name == "epex_da").one_or_none()
     charging_station = Sensor.query.filter(
         Sensor.name == charging_station_name
     ).one_or_none()
+    assert charging_station.get_attribute("capacity_in_mw") == 2
     assert Sensor.query.get(charging_station.get_attribute("market_id")) == epex_da
     start = as_server_time(datetime(2015, 1, 2))
     end = as_server_time(datetime(2015, 1, 3))
@@ -113,3 +115,61 @@ def test_charging_station_solver_day_2(target_soc, charging_station_name):
     print(consumption_schedule.head(12))
     print(soc_schedule.head(12))
     assert abs(soc_schedule.loc[target_soc_datetime] - target_soc) < 0.00001
+
+
+@pytest.mark.parametrize(
+    "target_soc, charging_station_name",
+    [
+        (9, "Test charging station"),
+        (15, "Test charging station"),
+        (5, "Test charging station (bidirectional)"),
+        (15, "Test charging station (bidirectional)"),
+    ],
+)
+def test_fallback_to_unsolvable_problem(target_soc, charging_station_name):
+    """Starting with a state of charge 10 kWh, within 2 hours we should be able to reach
+    any state of charge in the range [10, 14] kWh for a unidirectional station,
+    or [6, 14] for a bidirectional station, given a charging capacity of 2 kW.
+    Here we test target states of charge outside that range, ones that we should be able
+    to get as close to as 1 kWh difference.
+    We want our scheduler to handle unsolvable problems like these with a sensible fallback policy.
+    """
+    soc_at_start = 10
+    duration_until_target = timedelta(hours=2)
+    expected_gap = 1
+
+    epex_da = Sensor.query.filter(Sensor.name == "epex_da").one_or_none()
+    charging_station = Sensor.query.filter(
+        Sensor.name == charging_station_name
+    ).one_or_none()
+    assert charging_station.get_attribute("capacity_in_mw") == 2
+    assert Sensor.query.get(charging_station.get_attribute("market_id")) == epex_da
+    start = as_server_time(datetime(2015, 1, 2))
+    end = as_server_time(datetime(2015, 1, 3))
+    resolution = timedelta(minutes=15)
+    target_soc_datetime = start + duration_until_target
+    soc_targets = pd.Series(
+        np.nan, index=pd.date_range(start, end, freq=resolution, closed="right")
+    )
+    soc_targets.loc[target_soc_datetime] = target_soc
+    consumption_schedule = schedule_charging_station(
+        charging_station, start, end, resolution, soc_at_start, soc_targets
+    )
+    soc_schedule = integrate_time_series(
+        consumption_schedule, soc_at_start, decimal_precision=6
+    )
+
+    # Check if constraints were met
+    assert (
+        min(consumption_schedule.values)
+        >= charging_station.get_attribute("capacity_in_mw") * -1
+    )
+    assert max(consumption_schedule.values) <= charging_station.get_attribute(
+        "capacity_in_mw"
+    )
+    print(consumption_schedule.head(12))
+    print(soc_schedule.head(12))
+    assert (
+        abs(abs(soc_schedule.loc[target_soc_datetime] - target_soc) - expected_gap)
+        < 0.00001
+    )

flexmeasures/data/models/planning/utils.py

@@ -111,3 +111,64 @@ def get_prices(
                 "Prices partially unknown for planning window."
             )
     return price_df, query_window
+
+
+def fallback_charging_policy(
+    sensor: Sensor,
+    device_constraints: pd.DataFrame,
+    start: datetime,
+    end: datetime,
+    resolution: timedelta,
+) -> pd.Series:
+    """This fallback charging policy is to just start charging or discharging, or do neither,
+    depending on the first target state of charge and the capabilities of the Charge Point.
+    Note that this ignores any cause of the infeasibility and,
+    while probably a decent policy for Charge Points,
+    should not be considered a robust policy for other asset types.
+    """
+    charge_power = (
+        sensor.get_attribute("capacity_in_mw")
+        if sensor.get_attribute("is_consumer")
+        else 0
+    )
+    discharge_power = (
+        -sensor.get_attribute("capacity_in_mw")
+        if sensor.get_attribute("is_producer")
+        else 0
+    )
+
+    charge_schedule = initialize_series(charge_power, start, end, resolution)
+    discharge_schedule = initialize_series(discharge_power, start, end, resolution)
+    idle_schedule = initialize_series(0, start, end, resolution)
+    if (
+        device_constraints["max"].first_valid_index() is not None
+        and device_constraints["max"][device_constraints["max"].first_valid_index()] < 0
+    ):
+        # start discharging to try and bring back the soc below the next max constraint
+        return discharge_schedule
+    if (
+        device_constraints["min"].first_valid_index() is not None
+        and device_constraints["min"][device_constraints["min"].first_valid_index()] > 0
+    ):
+        # start charging to try and bring back the soc above the next min constraint
+        return charge_schedule
+    if (
+        device_constraints["equals"].first_valid_index() is not None
+        and device_constraints["equals"][
+            device_constraints["equals"].first_valid_index()
+        ]
+        > 0
+    ):
+        # start charging to get as close as possible to the next target
+        return charge_schedule
+    if (
+        device_constraints["equals"].first_valid_index() is not None
+        and device_constraints["equals"][
+            device_constraints["equals"].first_valid_index()
+        ]
+        < 0
+    ):
+        # start discharging to get as close as possible to the next target
+        return discharge_schedule
+    # stand idle
+    return idle_schedule

flexmeasures/data/models/time_series.py

@@ -85,14 +85,30 @@ def location(self) -> Optional[Tuple[float, float]]:
             return self.latitude, self.longitude
         return None
 
+    @property
+    def is_strictly_non_positive(self) -> bool:
+        """Return True if this sensor strictly records non-positive values."""
+        return self.get_attribute("is_consumer", False) and not self.get_attribute(
+            "is_producer", True
+        )
+
+    @property
+    def is_strictly_non_negative(self) -> bool:
+        """Return True if this sensor strictly records non-negative values."""
+        return self.get_attribute("is_producer", False) and not self.get_attribute(
+            "is_consumer", True
+        )
+
     def get_attribute(self, attribute: str, default: Any = None) -> Any:
         """Looks for the attribute on the Sensor.
         If not found, looks for the attribute on the Sensor's GenericAsset.
         If not found, returns the default.
         """
+        if hasattr(self, attribute):
+            return getattr(self, attribute)
         if attribute in self.attributes:
             return self.attributes[attribute]
-        elif attribute in self.generic_asset.attributes:
+        if attribute in self.generic_asset.attributes:
             return self.generic_asset.attributes[attribute]
         return default