diff --git a/src/akkudoktoreos/optimization/genetic/terminalvalue.py b/src/akkudoktoreos/optimization/genetic/terminalvalue.py index f8c45990..d0a33985 100644 --- a/src/akkudoktoreos/optimization/genetic/terminalvalue.py +++ b/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] = [] diff --git a/src/akkudoktoreos/prediction/interpolator.py b/src/akkudoktoreos/prediction/interpolator.py index 92bbafad..1e2c1a0d 100644 --- a/src/akkudoktoreos/prediction/interpolator.py +++ b/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)) diff --git a/tests/test_interpolator.py b/tests/test_interpolator.py index ccb9a0b9..581310d7 100644 --- a/tests/test_interpolator.py +++ b/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) diff --git a/tests/test_terminalvalue.py b/tests/test_terminalvalue.py index acc2f002..c91b6701 100644 --- a/tests/test_terminalvalue.py +++ b/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)