fix(genetic): same plan on every CPU (fixed-order dot products, stable sort of equal prices) (#1518)

* fix(genetic): sum the self-consumption dot products in a fixed order

GENETIC gave different plans for the same input and seed on different
CPUs. SelfConsumptionProbabilityInterpolator computed its two dot
products with np.dot, which calls BLAS; OpenBLAS picks a ddot kernel per
CPU at run time (and macOS uses Accelerate), and the kernels sum in
different orders, so the result differs in the last bit. The optimizer
compares fitness values exactly, so one bit is enough for a different
plan.

fixed_order_dot() adds the rounded products in index order with Python
floats, without BLAS: the same result on every CPU. The vectors have one
entry per 50 W load bin and the results are cached per (load, PV) pair,
so the run time does not change noticeably. Plans change once in the
last bits compared with the previous np.dot result of a given CPU.

* fix(genetic): sort equal prices of the terminal value window in slot order

np.argsort with the default kind (quicksort) leaves the order of equal
values to the platform: NumPy sorts with SIMD kernels on x86-64
(x86-simd-sort) and with introsort elsewhere. Slots with the same price
then fill the terminal value curve in another order, the battery
capacity cuts it elsewhere, and the same input and seed gave another
plan on x86-64 than on aarch64 (seen as a native-vs-Python compare
difference on x86-64 in
test_plan_keeps_reserve_for_a_peak_beyond_the_control_horizon).

The window is now sorted with kind="stable": slots with the same price
keep their slot order on every platform. Same goal as the dot product
change in this PR: the same input gives the same plan on every CPU.

Assisted-by: Claude Code (claude-opus-5-5)
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

src/akkudoktoreos/optimization/genetic/terminalvalue.py

@@ -276,7 +276,10 @@ def build_terminal_value_curve(
     # LCOS is charged on delivered DC energy; the curve is in AC energy.
     lcos_per_wh_ac = (lcos_euro_per_kwh / 1000.0) / max(dc_to_ac_efficiency, 1e-9)
 
-    order = np.argsort(-prices)
+    # Stable: slots with the same price keep their order. The default (quicksort) leaves the
+    # order of equal prices to the platform - NumPy sorts with SIMD kernels on x86-64 - and that
+    # order decides where the curve is cut at the battery capacity, i.e. the plan.
+    order = np.argsort(-prices, kind="stable")
     energy_points: list[float] = [0.0]
     value_points: list[float] = [0.0]
     marginals: list[float] = []

src/akkudoktoreos/prediction/interpolator.py

@@ -14,6 +14,29 @@ from akkudoktoreos.utils.gridinterpolator import (
 )
 
 
+def fixed_order_dot(a: np.ndarray, b: np.ndarray) -> float:
+    """Return the dot product of two vectors, summed from the first element to the last.
+
+    ``np.dot`` calls BLAS, whose kernel is chosen per CPU at run time (OpenBLAS
+    ``DYNAMIC_ARCH``, Accelerate). The kernels sum in different orders and some use fused
+    multiply-add, so the last bit of the result depends on the machine. The optimizer
+    compares fitness values exactly, so one bit is enough to give a different plan for the
+    same input and seed. Here every product is rounded and added in index order, with
+    Python floats: the same result on every CPU.
+
+    Args:
+        a: First vector.
+        b: Second vector, same length as ``a``.
+
+    Returns:
+        ``sum(a[i] * b[i])`` in index order.
+    """
+    total = 0.0
+    for x, y in zip(a.tolist(), b.tolist(), strict=True):
+        total += x * y
+    return total
+
+
 class SelfConsumptionProbabilityInterpolator:
     LOAD_DISTRIBUTION_CACHE_SIZE = 4096
     DIRECT_CONSUMPTION_CACHE_SIZE = 8192
@@ -166,15 +189,16 @@ class SelfConsumptionProbabilityInterpolator:
             return 0.0
 
         load_levels_w, probabilities = self._load_distribution(mean_load_power_w)
-        modeled_mean_load_w = float(np.dot(probabilities, load_levels_w))
+        # fixed_order_dot instead of np.dot: same result on every CPU (see there).
+        modeled_mean_load_w = fixed_order_dot(probabilities, load_levels_w)
         if modeled_mean_load_w <= 0.0:
             return 0.0
 
         # Preserve the requested mean load while keeping the conditional shape
         # from the probability table.
         normalized_load_levels_w = load_levels_w * (mean_load_power_w / modeled_mean_load_w)
-        expected_direct_power_w = float(
-            np.dot(probabilities, np.minimum(normalized_load_levels_w, pv_power_w))
+        expected_direct_power_w = fixed_order_dot(
+            probabilities, np.minimum(normalized_load_levels_w, pv_power_w)
         )
         # Scalar bounds: min/max instead of np.clip (same result, no array round trip).
         return min(max(expected_direct_power_w, 0.0), min(mean_load_power_w, pv_power_w))

tests/test_interpolator.py

@@ -2,7 +2,10 @@ import numpy as np
 import pytest
 
 from akkudoktoreos.core.cache import CacheEnergyManagementStore
-from akkudoktoreos.prediction.interpolator import get_eos_load_interpolator
+from akkudoktoreos.prediction.interpolator import (
+    fixed_order_dot,
+    get_eos_load_interpolator,
+)
 
 
 def test_quarter_hour_energy_is_converted_back_to_same_mean_power():
@@ -133,3 +136,37 @@ def test_expected_direct_consumption_cache_is_bounded(monkeypatch):
     assert interpolator.calculate_expected_direct_consumption(
         800.0, 500.0
     ) == interpolator._expected_direct_consumption(800.0, 500.0)
+
+
+def test_fixed_order_dot_sums_in_index_order():
+    """The sum runs from the first product to the last, whatever BLAS would do."""
+    ones = np.ones(3)
+    # 1e16 + 1 rounds back to 1e16, so the order decides whether the 1 survives.
+    assert fixed_order_dot(np.array([1e16, 1.0, -1e16]), ones) == 0.0
+    assert fixed_order_dot(np.array([1e16, -1e16, 1.0]), ones) == 1.0
+
+    rng = np.random.default_rng(7)
+    a, b = rng.random(97), rng.random(97) * 5000.0
+    expected = 0.0
+    for x, y in zip(a.tolist(), b.tolist()):
+        expected += x * y
+    assert fixed_order_dot(a, b) == expected
+    assert fixed_order_dot(a, b) == pytest.approx(float(np.dot(a, b)))
+
+
+def test_expected_direct_consumption_does_not_use_blas(monkeypatch):
+    """The expected direct consumption must not depend on the CPU's BLAS kernel."""
+    interpolator = get_eos_load_interpolator()
+    load_levels_w, probabilities = interpolator._load_distribution(800.0)
+    modeled_mean_load_w = 0.0
+    for p, level in zip(probabilities.tolist(), load_levels_w.tolist()):
+        modeled_mean_load_w += p * level
+    expected = 0.0
+    for p, level in zip(probabilities.tolist(), load_levels_w.tolist()):
+        expected += p * min(level * (800.0 / modeled_mean_load_w), 1200.0)
+
+    def no_blas(*args, **kwargs):
+        raise AssertionError("np.dot must not be used: its rounding depends on the CPU")
+
+    monkeypatch.setattr(np, "dot", no_blas)
+    assert interpolator._expected_direct_consumption(800.0, 1200.0) == min(expected, 800.0)

tests/test_terminalvalue.py

@@ -127,3 +127,27 @@ def test_trailing_window_takes_the_end_of_the_horizon():
     # A window longer than the horizon yields what there is.
     assert list(trailing_window(values, end_slot=2, window_slots=5)) == [0.0, 1.0]
     assert list(trailing_window(None, end_slot=8, window_slots=3)) == []
+
+
+def test_equal_prices_keep_their_slot_order():
+    """Slots with the same price fill the curve in slot order, on every platform.
+
+    The capacity cuts the curve inside the equal-price slots, so their order decides the
+    breakpoints; NumPy's default sort leaves it to the platform (SIMD kernels on x86-64).
+    """
+    prices = np.array([0.0003] * 40 + [0.0004] + [0.0003] * 40)
+    load = np.array([100.0 + slot for slot in range(81)])
+    curve = _curve(
+        prices_euro_per_wh=prices,
+        load_wh=load,
+        pv_wh=np.zeros(81),
+        feed_in_euro_per_wh=np.zeros(81),
+        max_energy_wh=1000.0,
+    )
+
+    # The 0.40 slot (index 40, 140 Wh) first, then the 0.30 slots from index 0 on.
+    expected = [0.0, 140.0]
+    for slot in range(81):
+        if slot != 40 and expected[-1] < 1000.0:
+            expected.append(min(expected[-1] + 100.0 + slot, 1000.0))
+    assert curve.energy_wh == pytest.approx(expected)