Merge remote-tracking branch 'origin/perf/two-stage-solve' into rollup/two-stage-plus-130

Author: andig

src/optimizer/optimizer.py

@@ -87,6 +87,11 @@ def objective_scale(objective) -> float:
 # peak. A hundredth of a cent, three orders below the gap the cost stage may already leave.
 COST_BOUND_TOLERANCE = 1e-4
 
+# how far off 0 or 1 a binary may land and still count as integral. CBC's own integer tolerance
+# defaults to 1e-6 and it rounds within that before reporting, so anything above this came back
+# from a relaxation rather than from a rounding difference.
+INTEGRALITY_TOLERANCE = 1e-5
+
 # share of OPTIMIZER_TIME_LIMIT the probe gets before the solve falls back to the split. Kept well
 # under half: the probe is pure loss on a request that ends up splitting anyway, so it should be
 # long enough to catch the ordinary ones and no longer. Measured over the captured slow requests,
@@ -827,7 +832,17 @@ def _solve_preferences(self, tmpdir, deadline) -> None:
         # rather than read off the status, so a solver that reports the wrong one cannot spend
         # money.
         self.preference_stage = pulp.LpStatus[self.problem.status]
-        improved = (pulp.value(self.preference_objective) > undecided
+        # a stage that ran out of clock before it found an integer solution leaves the relaxation
+        # in the variables, and pulp reads that back like any other result. It looks like a large
+        # improvement precisely because it is one the model forbids: the binaries land between 0
+        # and 1, and every rule they gate stops holding, c_min among them. Checked here beside the
+        # other two conditions, for the same reason they are checked here rather than read off the
+        # status: a solver that reports the wrong one must not be able to spend money, and it must
+        # not be able to hand back a schedule the model does not allow either.
+        integral = self.problem.sol_status in (pulp.LpSolutionOptimal,
+                                               pulp.LpSolutionIntegerFeasible)
+        improved = (integral
+                    and pulp.value(self.preference_objective) > undecided
                     and pulp.value(self.cost_objective) >= cost - budget - COST_BOUND_TOLERANCE)
         if not improved:
             self.preference_stage += ', kept the first stage'
@@ -891,6 +906,16 @@ def _probe_then_split(self, tmpdir, deadline) -> None:
                 var.varValue = value
             self.problem.status = pulp.LpStatusOptimal
 
+    def _is_integral(self) -> bool:
+        """Whether every integer variable of the current solution came back on a whole number.
+
+        pulp stores a binary as an integer bounded to 0 and 1, so LpBinary never appears on a
+        variable and every gate in this model is covered by the integer category alone.
+        """
+        return all(abs(var.varValue - round(var.varValue)) <= INTEGRALITY_TOLERANCE
+                   for var in self.problem.variables()
+                   if var.cat == pulp.LpInteger and var.varValue is not None)
+
     def solve(self) -> Dict:
         """
         Creates the MILP model if none exists and solves the optimization problem.
@@ -921,6 +946,14 @@ def solve(self) -> Dict:
         if status == 'Optimal' and self.problem.sol_status != pulp.LpSolutionOptimal:
             status = 'Feasible'
 
+        # last line of defence. Every stage above decides for itself whether to keep what came
+        # back, so nothing should reach this point off the integers, but a schedule that breaks
+        # the model is worse than no schedule: it looks like an answer, and the caller charges a
+        # battery by it. One pass over the binaries against a solve measured in seconds.
+        if status in ('Optimal', 'Feasible') and not self._is_integral():
+            print("solver returned a fractional solution, reporting no schedule")
+            status = 'Not Solved'
+
         # grid import and export if no demand rate is active
         # if a limit is set and exceeded, this is the part that is actually imported / exported.
         # the exceeding portion is captured in 'e_imp_lim_exc' and / or 'e_exp_lim_exc'

tests/test_fractional_solution.py

@@ -0,0 +1,83 @@
+import pulp
+import pytest
+from test_objective_split import build
+
+from optimizer.optimizer import INTEGRALITY_TOLERANCE
+
+# small enough to solve in well under a second, and it has binaries to go fractional: the c_min
+# gate of the early charging case is exactly the rule a relaxed binary stops enforcing
+CASE = '012-early-charging-not-perfect'
+
+
+def binaries(optimizer):
+    # pulp stores a binary as an integer bounded to 0 and 1, LpBinary never survives on a variable
+    return [var for var in optimizer.problem.variables() if var.cat == pulp.LpInteger]
+
+
+def relax(optimizer, value=0.3):
+    """Scribble a relaxation over the binaries, the way a stage that found no integer solution
+    leaves them behind."""
+    for var in binaries(optimizer):
+        var.varValue = value
+
+
+def test_a_preference_stage_without_an_integer_solution_is_not_kept(monkeypatch):
+    # the second stage decides whether to keep its result by reading the variables, and a solver
+    # that ran out of clock before it found an integer solution leaves the relaxation in them. That
+    # point scores better on the preferences than any real schedule, because it is one the model
+    # forbids, so it has to be refused on the status rather than on its score.
+    optimizer = build(CASE)
+    optimizer.settings.probe_seconds = 0
+    optimizer.create_model()
+
+    real_solve = optimizer.problem.solve
+    calls = []
+
+    def solve(*args, **kwargs):
+        calls.append(1)
+        if len(calls) == 1:          # the cost stage, left alone
+            return real_solve(*args, **kwargs)
+        relax(optimizer)             # the preference stage, out of time and empty handed
+        optimizer.problem.status = pulp.LpStatusNotSolved
+        optimizer.problem.sol_status = pulp.LpSolutionNoSolutionFound
+        return optimizer.problem.status
+
+    monkeypatch.setattr(optimizer.problem, 'solve', solve)
+    optimizer.solve()
+
+    assert len(calls) == 2, f'the preference stage did not run, {len(calls)} solves'
+    assert optimizer.preference_stage.endswith('kept the first stage'), \
+        f'preference stage ended as {optimizer.preference_stage}'
+    for var in binaries(optimizer):
+        assert min(abs(var.varValue), abs(var.varValue - 1)) <= INTEGRALITY_TOLERANCE, \
+            f'{var.name} came back at {var.varValue}'
+
+
+def test_a_fractional_solution_is_reported_as_no_schedule(monkeypatch):
+    # the last line of defence, standing in for every stage above deciding correctly. A schedule
+    # that breaks the model is worse than no schedule: it looks like an answer and the caller
+    # charges a battery by it. Here the solve is faked wholesale, a solver claiming a solution it
+    # does not have.
+    optimizer = build(CASE)
+    optimizer.create_model()
+
+    def probe_then_split(tmpdir, deadline):
+        relax(optimizer)
+        optimizer.problem.status = pulp.LpStatusOptimal
+        optimizer.problem.sol_status = pulp.LpSolutionIntegerFeasible
+
+    monkeypatch.setattr(optimizer, '_probe_then_split', probe_then_split)
+    result = optimizer.solve()
+
+    assert result['status'] == 'Not Solved', f"reported {result['status']}"
+    assert result['objective_value'] is None
+    assert result['batteries'] == []
+
+
+@pytest.mark.parametrize('value', [0.0, 1.0, INTEGRALITY_TOLERANCE / 2, 1 - INTEGRALITY_TOLERANCE / 2])
+def test_a_solution_on_the_integers_passes(value):
+    # the guard must not fire on CBC's own rounding, which it reports within its integer tolerance
+    optimizer = build(CASE)
+    optimizer.create_model()
+    relax(optimizer, value)
+    assert optimizer._is_integral()