diff --git a/CHANGELOG.md b/CHANGELOG.md index cdecbcd0..066ba618 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -84,6 +84,12 @@ compatibility changes and a short test procedure. These changes are not yet a ta - Lower idle memory: matplotlib (PDF reports), statsmodels (price fallback forecast), pvlib, bs4/lxml, scipy (self-consumption table), the GENETIC0 optimizer and the timezone lookup data are imported on first use instead of at server start. Results are unchanged. +- Faster GENETIC runs with unchanged plans (about 60 % less run time on a Raspberry Pi 4): + table lookups, the appliance starts and the battery state decoding are worked out once per + run instead of once per candidate, the slot simulation and the battery bookkeeping work on + plain numbers instead of single NumPy elements, the terminal value report is built for the + final solution only, offspring are created without `deepcopy`, and identical genomes of a + generation are evaluated once also when the fitness cache is off. ### Fixed diff --git a/src/akkudoktoreos/core/cache.py b/src/akkudoktoreos/core/cache.py index 1ad3fc8e..073866b1 100644 --- a/src/akkudoktoreos/core/cache.py +++ b/src/akkudoktoreos/core/cache.py @@ -84,6 +84,9 @@ class CacheEnergyManagementStore(SingletonMixin): last_value: ClassVar[Any] = None hit_count: ClassVar[int] = 0 miss_count: ClassVar[int] = 0 + # Counts clear() calls. Callers that keep their own per-run cache compare + # it to drop their entries together with this store. + generation: ClassVar[int] = 0 def __init__(self, *args: Any, **kwargs: Any) -> None: """Initializes the `CacheEnergyManagementStore` instance with default parameters. @@ -188,6 +191,7 @@ class CacheEnergyManagementStore(SingletonMixin): """ if hasattr(self.cache, "clear") and callable(getattr(self.cache, "clear")): CacheEnergyManagementStore.cache.clear() + CacheEnergyManagementStore.generation += 1 CacheEnergyManagementStore.last_event = None CacheEnergyManagementStore.last_key = None CacheEnergyManagementStore.last_value = None diff --git a/src/akkudoktoreos/devices/genetic/battery.py b/src/akkudoktoreos/devices/genetic/battery.py index 2bef2450..640c1bfc 100644 --- a/src/akkudoktoreos/devices/genetic/battery.py +++ b/src/akkudoktoreos/devices/genetic/battery.py @@ -211,11 +211,17 @@ class Battery: self.max_charge_power_w = self.capacity_wh # TODO this should not be equal capacity_wh self.discharge_array = np.full(self.prediction_hours, 0) self.charge_array = np.full(self.prediction_hours, 0) - self._discharged_raw_wh_per_slot = np.zeros(self.prediction_hours, dtype=float) - self._charged_raw_wh_per_slot = np.zeros(self.prediction_hours, dtype=float) + # Discharged energy per slot, a plain list: it is read and written one + # slot at a time for every candidate plan, where single NumPy elements + # (and the NumPy scalars they return) are slow. + self._discharged_raw_wh_per_slot: list[float] = [0.0] * self.prediction_hours + self._charged_raw_wh_per_slot: list[float] = [0.0] * self.prediction_hours # Optional per-slot cap on the raw charge energy from all sources. It is # unbounded unless an inverter restricts a slot (see limit_slot_charge). - self._charge_limit_raw_wh_per_slot = np.full(self.prediction_hours, np.inf) + # Like the charged energy it is kept in a plain list: both are read and + # written one slot at a time for every candidate plan, where single + # NumPy elements (and the NumPy scalars they return) are slow. + self._charge_limit_raw_wh_per_slot: list[float] = [float("inf")] * self.prediction_hours self.soc_wh = (self.initial_soc_percentage / 100) * self.capacity_wh self.min_soc_wh = (self.min_soc_percentage / 100) * self.capacity_wh self.max_soc_wh = (self.max_soc_percentage / 100) * self.capacity_wh @@ -238,6 +244,18 @@ class Battery: # Yield values before idx in reverse (descending) return (charge_rates_fast[j] for j in range(idx - 1, -1, -1)) + def _charge_rates_descending(self) -> tuple[float, ...]: + """Return the charge rates from highest to lowest as plain numbers. + + Built once per ``charge_rates`` array instead of searching and slicing + the array for every charge request that does not fit. + """ + rates = self.charge_rates + if rates is not getattr(self, "_charge_rates_desc_source", None): + self._charge_rates_desc = tuple(float(rate) for rate in rates[::-1]) + self._charge_rates_desc_source = rates + return self._charge_rates_desc + def to_dict(self) -> dict[str, Any]: """Converts the object to a dictionary representation.""" return { @@ -257,11 +275,11 @@ class Battery: """Resets the battery state to its initial values.""" self.soc_wh = (self.initial_soc_percentage / 100) * self.capacity_wh self.soc_wh = min(self.soc_wh, self.max_soc_wh) # Only clamp to max - self.discharge_array = np.full(self.prediction_hours, 0) - self.charge_array = np.full(self.prediction_hours, 0) - self._discharged_raw_wh_per_slot = np.zeros(self.prediction_hours, dtype=float) - self._charged_raw_wh_per_slot = np.zeros(self.prediction_hours, dtype=float) - self._charge_limit_raw_wh_per_slot = np.full(self.prediction_hours, np.inf) + self.discharge_array = np.zeros(self.prediction_hours, dtype=np.int_) + self.charge_array = np.zeros(self.prediction_hours, dtype=np.int_) + self._discharged_raw_wh_per_slot = [0.0] * self.prediction_hours + self._charged_raw_wh_per_slot = [0.0] * self.prediction_hours + self._charge_limit_raw_wh_per_slot = [float("inf")] * self.prediction_hours def limit_slot_charge(self, hour: int, raw_wh: float) -> None: """Cap the raw energy the battery may take in one slot, from all sources. @@ -295,6 +313,10 @@ class Battery: raw_soc_available_wh = max(self.soc_wh - self.min_soc_wh, 0.0) return min(raw_power_remaining_wh, raw_soc_available_wh) * self.discharging_efficiency + def discharge_released(self, hour: int) -> bool: + """Whether the plan lets the battery discharge in this slot.""" + return self.discharge_array[hour] != 0 + def discharged_energy_wh(self, hour: int) -> float: """Return DC energy delivered by the battery in one optimization slot.""" return self._discharged_raw_wh_per_slot[hour] * self.discharging_efficiency @@ -464,9 +486,13 @@ class Battery: raw_request_wh = max_charge_per_slot_wh_fast * charge_factor raw_charge_wh = max(self.max_soc_wh - soc_wh_fast, 0.0) / charging_efficiency_fast if raw_request_wh > raw_charge_wh: - # Use a lower charge factor - lower_charge_factors = self._lower_charge_rates_desc(charge_factor) - for charge_factor in lower_charge_factors: + # Use a lower charge factor: the rates below the requested + # one, highest first (as _lower_charge_rates_desc yields them). + requested_charge_factor = charge_factor + for lower_charge_factor in self._charge_rates_descending(): + if lower_charge_factor >= requested_charge_factor: + continue + charge_factor = lower_charge_factor raw_request_wh = max_charge_per_slot_wh_fast * charge_factor if raw_request_wh <= raw_charge_wh: self.charge_array[hour] = charge_factor diff --git a/src/akkudoktoreos/devices/genetic/inverter.py b/src/akkudoktoreos/devices/genetic/inverter.py index 09453bae..48db3dc2 100644 --- a/src/akkudoktoreos/devices/genetic/inverter.py +++ b/src/akkudoktoreos/devices/genetic/inverter.py @@ -1,3 +1,4 @@ +from bisect import bisect_left from typing import Optional from loguru import logger @@ -159,8 +160,31 @@ class Inverter: """ if self.dc_to_ac_efficiency_curve is None: return self.dc_to_ac_efficiency + curve: list[tuple[float, float]] = self.dc_to_ac_efficiency_curve load_fraction = ac_wh / self.max_power_wh if self.max_power_wh > 0 else 0.0 - return interpolate_efficiency_curve(self.dc_to_ac_efficiency_curve, load_fraction) + # This runs for every battery discharge of every candidate plan. Find + # the curve segment by bisection on the load fractions, split off once + # per curve, instead of walking the points. The segment and the + # interpolation are those of interpolate_efficiency_curve(). + if curve is not getattr(self, "_curve_source", None): + self._curve_fractions = [point[0] for point in curve] + self._curve_efficiencies = [point[1] for point in curve] + self._curve_source = curve + fractions = self._curve_fractions + if load_fraction != load_fraction or len(fractions) != len(curve): + # NaN, or a curve changed in place: the plain walk. + return interpolate_efficiency_curve(curve, load_fraction) + efficiencies = self._curve_efficiencies + index = bisect_left(fractions, load_fraction) + if index == 0: + return efficiencies[0] + if index == len(fractions): + return efficiencies[-1] + lower_fraction = fractions[index - 1] + lower_efficiency = efficiencies[index - 1] + return lower_efficiency + (load_fraction - lower_fraction) / ( + fractions[index] - lower_fraction + ) * (efficiencies[index] - lower_efficiency) def ac_charge_factor(self, factor: float) -> float: """Return the AC charge factor the inverter can actually execute. @@ -220,23 +244,82 @@ class Inverter: ac_energy = dc_energy / self.ac_to_dc_efficiency return ac_energy, battery_losses + (ac_energy - dc_energy) + def _delivered_ac_wh( + self, + battery_discharge_dc: float, + dc_request: float, + requested_ac_wh: float, + dc_to_ac_efficiency: float, + ) -> float: + """Return the AC energy of a discharge that delivered ``battery_discharge_dc``. + + If the battery delivered less than requested (low state of charge or its + power limit), the conversion runs at a lower AC load than requested. With + an efficiency curve the efficiency must then be the one at the AC energy + that is actually delivered, not at the request. + """ + curve = self.dc_to_ac_efficiency_curve + if curve is not None and battery_discharge_dc < dc_request * (1.0 - 1e-12): + return solve_ac_for_dc_energy( + curve, battery_discharge_dc, self.max_power_wh, requested_ac_wh + ) + return battery_discharge_dc * dc_to_ac_efficiency + def _discharge_battery_to_ac(self, requested_ac_wh: float, hour: int) -> tuple[float, float]: """Discharge battery energy and convert it to AC energy.""" - if not self.battery or requested_ac_wh <= 0.0: + battery = self.battery + if not battery or requested_ac_wh <= 0.0: return 0.0, 0.0 + if type(battery) is not Battery: + # Another battery implementation: ask it, step by step. + if not battery.discharge_released(hour): + return 0.0, 0.0 + dc_to_ac_efficiency = self.dc_to_ac_efficiency_at(requested_ac_wh) + dc_request = requested_ac_wh / dc_to_ac_efficiency + battery_discharge_dc, discharge_losses = battery.discharge_energy(dc_request, hour) + battery_discharge_ac = self._delivered_ac_wh( + battery_discharge_dc, dc_request, requested_ac_wh, dc_to_ac_efficiency + ) + inverter_discharge_losses = battery_discharge_dc - battery_discharge_ac + return battery_discharge_ac, discharge_losses + inverter_discharge_losses + # This runs for every slot of every candidate plan in which the battery + # could serve load or export. For the plain Battery the release check, + # the efficiency and Battery.discharge_energy() are computed in one go, + # without the nested calls. The arithmetic is that of + # Battery.discharge_energy(); keep both in step. + if battery.discharge_array[hour] == 0: + # Not released in this slot: the battery delivers nothing. + return 0.0, 0.0 # With an efficiency curve, the efficiency is taken at the requested AC # energy of this conversion. - dc_to_ac_efficiency = self.dc_to_ac_efficiency_at(requested_ac_wh) + if self.dc_to_ac_efficiency_curve is None: + dc_to_ac_efficiency = self.dc_to_ac_efficiency + else: + dc_to_ac_efficiency = self.dc_to_ac_efficiency_at(requested_ac_wh) dc_request = requested_ac_wh / dc_to_ac_efficiency - battery_discharge_dc, discharge_losses = self.battery.discharge_energy(dc_request, hour) + + soc_wh = battery.soc_wh + min_soc_wh = battery.min_soc_wh + discharged_raw_wh = battery._discharged_raw_wh_per_slot + discharging_efficiency = battery.discharging_efficiency + raw_available_wh = max(soc_wh - min_soc_wh, 0.0) + max_raw_wh = max( + battery.max_charge_power_w * battery.slot_duration_h - discharged_raw_wh[hour], 0.0 + ) + max_deliverable_wh = min(raw_available_wh, max_raw_wh) * discharging_efficiency + battery_discharge_dc = min(dc_request, max_deliverable_wh) + raw_used_wh = battery_discharge_dc / discharging_efficiency + soc_wh -= raw_used_wh + battery.soc_wh = max(soc_wh, min_soc_wh) + discharged_raw_wh[hour] += raw_used_wh + discharge_losses = raw_used_wh - battery_discharge_dc + + # Same rule as _delivered_ac_wh(), without the call: a limited discharge + # is converted with the efficiency at the AC energy actually delivered. if self.dc_to_ac_efficiency_curve is not None and battery_discharge_dc < dc_request * ( 1.0 - 1e-12 ): - # The battery delivered less than requested (low state of charge or - # its power limit). The conversion then runs at a lower AC load than - # requested, so the efficiency must be the one at the AC energy that - # is actually delivered, not at the request. battery_discharge_ac = solve_ac_for_dc_energy( self.dc_to_ac_efficiency_curve, battery_discharge_dc, @@ -275,75 +358,120 @@ class Inverter: battery and the inverter allow, which is the behaviour when no export rates are configured. """ + # This method runs once per slot of every candidate plan. The bounds + # below are written as conditional expressions instead of min()/max() + # calls: "0.0 if x < 0.0 else x" is exactly max(x, 0.0) and + # "b if b < a else a" is exactly min(a, b), NaN and signed zero included. losses = 0.0 grid_export = 0.0 - generation = max(float(generation), 0.0) - consumption = max(float(consumption), 0.0) + generation = float(generation) + if generation < 0.0: + generation = 0.0 + consumption = float(consumption) + if consumption < 0.0: + consumption = 0.0 + battery = self.battery + max_power_wh = self.max_power_wh # Convert interval energy [Wh] to mean power [W] for the probability # lookup, then convert its expected direct power back to slot energy. if generation > 0.0 and consumption > 0.0: + slot_duration_h = self.slot_duration_h expected_direct_power_w = ( self.self_consumption_predictor.calculate_expected_direct_consumption( - consumption / self.slot_duration_h, - generation / self.slot_duration_h, + consumption / slot_duration_h, + generation / slot_duration_h, ) ) - direct_pv_energy = expected_direct_power_w * self.slot_duration_h + direct_pv_energy = expected_direct_power_w * slot_duration_h else: direct_pv_energy = 0.0 # Direct PV is bounded by both input energies and by the AC energy the # inverter can move during this slot. - direct_pv_energy = min( - max(direct_pv_energy, 0.0), - generation, - consumption, - self.max_power_wh, - ) - remaining_load = max(consumption - direct_pv_energy, 0.0) - pv_surplus = max(generation - direct_pv_energy, 0.0) - remaining_inverter_ac_capacity = max(self.max_power_wh - direct_pv_energy, 0.0) + if direct_pv_energy < 0.0: + direct_pv_energy = 0.0 + if generation < direct_pv_energy: + direct_pv_energy = generation + if consumption < direct_pv_energy: + direct_pv_energy = consumption + if max_power_wh < direct_pv_energy: + direct_pv_energy = max_power_wh + remaining_load = consumption - direct_pv_energy + if remaining_load < 0.0: + remaining_load = 0.0 + pv_surplus = generation - direct_pv_energy + if pv_surplus < 0.0: + pv_surplus = 0.0 + remaining_inverter_ac_capacity = max_power_wh - direct_pv_energy + if remaining_inverter_ac_capacity < 0.0: + remaining_inverter_ac_capacity = 0.0 # Load gaps and PV surplus may both occur within the same coarse slot. # Cover the load gap first; this preserves the existing chronological # approximation and can create headroom for later PV charging. battery_discharge_ac = 0.0 - if remaining_load > 0.0 and self.battery and remaining_inverter_ac_capacity > 0.0: - requested_ac_wh = min(remaining_load, remaining_inverter_ac_capacity) + if remaining_load > 0.0 and battery and remaining_inverter_ac_capacity > 0.0: + requested_ac_wh = ( + remaining_inverter_ac_capacity + if remaining_inverter_ac_capacity < remaining_load + else remaining_load + ) battery_discharge_ac, battery_discharge_losses = self._discharge_battery_to_ac( requested_ac_wh, hour ) - remaining_load = max(remaining_load - battery_discharge_ac, 0.0) - remaining_inverter_ac_capacity = max( - remaining_inverter_ac_capacity - battery_discharge_ac, 0.0 - ) + remaining_load = remaining_load - battery_discharge_ac + if remaining_load < 0.0: + remaining_load = 0.0 + remaining_inverter_ac_capacity = remaining_inverter_ac_capacity - battery_discharge_ac + if remaining_inverter_ac_capacity < 0.0: + remaining_inverter_ac_capacity = 0.0 losses += battery_discharge_losses grid_import = remaining_load - # Charge from the probabilistic PV surplus on the DC path. Stored energy - # plus charge losses equals the PV energy accepted by the battery. - remaining_surplus = pv_surplus - if remaining_surplus > 0.0 and self.battery: - charged_energy, charge_losses = self.battery.charge_energy(remaining_surplus, hour) - remaining_surplus = max(remaining_surplus - charged_energy - charge_losses, 0.0) - losses += charge_losses + # Without PV surplus (night slots, or all PV consumed directly) there is + # nothing to charge, export or curtail. (!= instead of >: a NaN forecast + # keeps propagating into the result as before.) + if pv_surplus != 0.0: + # Charge from the probabilistic PV surplus on the DC path. Stored + # energy plus charge losses equals the PV energy accepted by the + # battery. + remaining_surplus = pv_surplus + if remaining_surplus > 0.0 and battery: + charged_energy, charge_losses = battery.charge_energy(remaining_surplus, hour) + remaining_surplus = remaining_surplus - charged_energy - charge_losses + if remaining_surplus < 0.0: + remaining_surplus = 0.0 + losses += charge_losses - pv_grid_export = min(remaining_surplus, remaining_inverter_ac_capacity) - grid_export += pv_grid_export - remaining_inverter_ac_capacity = max(remaining_inverter_ac_capacity - pv_grid_export, 0.0) - # PV which can neither charge the battery nor pass through the inverter - # is curtailed and reported as a loss. - losses += max(remaining_surplus - pv_grid_export, 0.0) + pv_grid_export = ( + remaining_inverter_ac_capacity + if remaining_inverter_ac_capacity < remaining_surplus + else remaining_surplus + ) + grid_export += pv_grid_export + remaining_inverter_ac_capacity = remaining_inverter_ac_capacity - pv_grid_export + if remaining_inverter_ac_capacity < 0.0: + remaining_inverter_ac_capacity = 0.0 + # PV which can neither charge the battery nor pass through the + # inverter is curtailed and reported as a loss. + curtailed = remaining_surplus - pv_grid_export + if curtailed < 0.0: + curtailed = 0.0 + losses += curtailed - if allow_battery_grid_export and self.battery and remaining_inverter_ac_capacity > 0.0: - export_factor = min(max(float(battery_grid_export_factor), 0.0), 1.0) + if allow_battery_grid_export and battery and remaining_inverter_ac_capacity > 0.0: + export_factor = float(battery_grid_export_factor) + if export_factor < 0.0: + export_factor = 0.0 + if 1.0 < export_factor: + export_factor = 1.0 # Upper bounds of the export. With an efficiency curve, the # efficiency is estimated at the DC energy of each bound; the # conversion itself in _discharge_battery_to_ac() is exact and is # limited by the energy the battery can actually deliver. - remaining_battery_dc = self.battery.remaining_discharge_energy_wh(hour) + remaining_battery_dc = battery.remaining_discharge_energy_wh(hour) remaining_battery_ac = remaining_battery_dc * self.dc_to_ac_efficiency_at( remaining_battery_dc ) @@ -351,11 +479,13 @@ class Inverter: # power, so it stays a plain power setpoint ("export at 50 %") that # does not silently grow when self-consumption used less of the slot. # At factor 1.0 this bound never binds; behaviour is unchanged. - rated_export_dc = self.battery.rated_discharge_energy_wh() * export_factor + rated_export_dc = battery.rated_discharge_energy_wh() * export_factor rated_export_ac = rated_export_dc * self.dc_to_ac_efficiency_at(rated_export_dc) - export_capacity = min( - remaining_inverter_ac_capacity, remaining_battery_ac, rated_export_ac - ) + export_capacity = remaining_inverter_ac_capacity + if remaining_battery_ac < export_capacity: + export_capacity = remaining_battery_ac + if rated_export_ac < export_capacity: + export_capacity = rated_export_ac battery_export_ac, battery_export_losses = self._discharge_battery_to_ac( export_capacity, hour ) diff --git a/src/akkudoktoreos/optimization/genetic/genetic.py b/src/akkudoktoreos/optimization/genetic/genetic.py index 6c3ad56e..b85f1870 100644 --- a/src/akkudoktoreos/optimization/genetic/genetic.py +++ b/src/akkudoktoreos/optimization/genetic/genetic.py @@ -124,6 +124,18 @@ def _pack_genes(values: list[int]) -> PackedGenes: gene fits nearly always; any value outside 0-255 switches the whole genome to 8 bytes per gene. The leading tag byte keeps both encodings apart. """ + # Nearly every genome is narrow: let bytes() do the range check in C. It + # rejects exactly the values the check below rejects (outside 0-255). + try: + return ( + b"\x00", + *[ + bytes(values[start : start + _NARROW_CHUNK_GENES]) + for start in range(0, len(values), _NARROW_CHUNK_GENES) + ], + ) + except (TypeError, ValueError): + pass narrow = all(0 <= value <= 255 for value in values) width = _NARROW_CHUNK_GENES if narrow else _WIDE_CHUNK_GENES chunks: list[bytes] = [] @@ -136,6 +148,19 @@ def _pack_genes(values: list[int]) -> PackedGenes: return (b"\x00" if narrow else b"\x01", *chunks) +def _slot_values(values: Any) -> Any: + """Return per-slot values as plain Python numbers for the simulation loop. + + Reading single elements of a NumPy array yields NumPy scalars, whose + arithmetic is several times slower than that of Python numbers. The values + are the same IEEE doubles (or integers), so the results do not change. + Other element types are left as they are. + """ + if isinstance(values, np.ndarray) and (values.dtype == np.float64 or values.dtype.kind in "iu"): + return values.tolist() + return values + + def _unpack_genes(packed: PackedGenes) -> list[int]: """Decode genes encoded by ``_pack_genes``.""" if packed[0] == b"\x00": @@ -533,21 +558,27 @@ class GeneticSimulation(PydanticBaseModel): len(load_energy_array_fast), self.prediction_hours or len(load_energy_array_fast) ) total_hours = end_hour - start_hour + if total_hours < 0: + raise ValueError("negative dimensions are not allowed") - # Pre-allocate arrays for the results, optimized for speed - loads_energy_per_hour = np.full((total_hours), np.nan) - feedin_energy_per_hour = np.full((total_hours), np.nan) - consumption_energy_per_hour = np.full((total_hours), np.nan) - costs_per_hour = np.full((total_hours), np.nan) - revenue_per_hour = np.full((total_hours), np.nan) - losses_wh_per_hour = np.full((total_hours), np.nan) - electricity_price_per_hour = np.full((total_hours), np.nan) - feed_in_tariff_per_hour = np.full((total_hours), np.nan) + # Per-slot results are collected in plain lists (element writes into + # NumPy arrays are slow) and turned into arrays after the loop. + nan = float("nan") + loads_energy_per_hour: Any = [nan] * total_hours + feedin_energy_per_hour: Any = [nan] * total_hours + consumption_energy_per_hour: Any = [nan] * total_hours + costs_per_hour: Any = [nan] * total_hours + revenue_per_hour: Any = [nan] * total_hours + losses_wh_per_hour: Any = [nan] * total_hours + electricity_price_per_hour: Any = [nan] * total_hours + feed_in_tariff_per_hour: Any = [nan] * total_hours + soc_per_hour: Any + soc_ev_per_hour: Any + home_appliance_wh_per_hour: Any # Set initial state if battery_fast: - # Pre-allocate arrays for the results, optimized for speed - soc_per_hour = np.full((total_hours), np.nan) + soc_per_hour = [nan] * total_hours soc_per_hour[0] = battery_fast.current_soc_percentage() @@ -605,8 +636,7 @@ class GeneticSimulation(PydanticBaseModel): ac_charging_possible = False if ev_fast: - # Pre-allocate arrays for the results, optimized for speed - soc_ev_per_hour = np.full((total_hours), np.nan) + soc_ev_per_hour = [nan] * total_hours soc_ev_per_hour[0] = ev_fast.current_soc_percentage() # Fill the charge array of the ev @@ -623,48 +653,115 @@ class GeneticSimulation(PydanticBaseModel): if home_appliances_fast: home_appliance_enabled = True - # Pre-allocate the aggregate appliance load array (sum over all - # devices). Each appliance already carries its own resampled load - # curve, built from the decoded start(s) before this call. - home_appliance_wh_per_hour = np.full((total_hours), np.nan) + # The aggregate appliance load (sum over all devices). Each + # appliance already carries its own resampled load curve, built + # from the decoded start(s) before this call. + home_appliance_wh_per_hour = [nan] * total_hours + appliance_loads_fast = [ + appliance.get_load_for_hour for appliance in home_appliances_fast + ] else: home_appliance_enabled = False # Default return if no home appliance is available home_appliance_wh_per_hour = np.full((total_hours), 0) + appliance_loads_fast = [] + + # The per-slot inputs as plain numbers. They are final here: the loop + # itself changes only the charge arrays of the devices, and those after + # the slot value was read. + load_energy_array_fast = _slot_values(load_energy_array_fast) + pv_prediction_wh_fast = _slot_values(pv_prediction_wh_fast) + elect_price_hourly_fast = _slot_values(elect_price_hourly_fast) + elect_revenue_per_hour_arr_fast = _slot_values(elect_revenue_per_hour_arr_fast) + bat_grid_export_hours_fast = _slot_values(bat_grid_export_hours_fast) + ac_charge_hours_fast = _slot_values(ac_charge_hours_fast) + ev_charge_hours_fast = _slot_values(ev_charge_hours_fast) + if battery_fast: + battery_lcos_per_wh_factor = battery_fast.levelized_cost_of_storage_kwh + + # Aggregate appliance load per slot, summed once per candidate in the + # order the loop used to add the appliances (0.0 + first + second ...). + # Anything unusual (another appliance class, a curve that does not + # cover the simulated slots) keeps the per-slot lookups and their + # range check. + appliance_load_fast: Optional[list[float]] = None + if home_appliance_enabled and start_hour >= 0: + appliance_load_fast = [0.0] * end_hour + for appliance in home_appliances_fast: + curve = getattr(appliance, "load_curve", None) + if ( + type(appliance).get_load_for_hour is not HomeAppliance.get_load_for_hour + or not isinstance(curve, np.ndarray) + or curve.dtype != np.float64 + or curve.ndim != 1 + or len(curve) < end_hour + or appliance.prediction_hours < end_hour + ): + appliance_load_fast = None + break + appliance_load_fast = [ + total + load for total, load in zip(appliance_load_fast, curve.tolist()) + ] + + # State of charge in percent, read directly where the device uses the + # plain Battery formula (two method calls per slot otherwise). + battery_soc_direct = bool( + battery_fast + and type(battery_fast).current_soc_percentage is Battery.current_soc_percentage + ) + ev_soc_direct = bool( + ev_fast and type(ev_fast).current_soc_percentage is Battery.current_soc_percentage + ) + battery_lcos_direct = False + if battery_fast is not None and type(battery_fast) is Battery: + battery_lcos_direct = True + battery_discharged_raw_fast = battery_fast._discharged_raw_wh_per_slot + battery_discharging_efficiency_fast = battery_fast.discharging_efficiency for hour in range(start_hour, end_hour): hour_idx = hour - start_hour # Accumulate loads and PV generation consumption = load_energy_array_fast[hour] - losses_wh_per_hour[hour_idx] = 0.0 + slot_losses = 0.0 # Home appliances (sum the per-slot load of all flexible consumers) if home_appliance_enabled: - ha_load = 0.0 - for appliance in home_appliances_fast: - ha_load += appliance.get_load_for_hour(hour) + if appliance_load_fast is not None: + ha_load = appliance_load_fast[hour] + else: + ha_load = 0.0 + for appliance_load in appliance_loads_fast: + ha_load += appliance_load(hour) consumption += ha_load home_appliance_wh_per_hour[hour_idx] = ha_load # E-Auto handling if ev_fast: - soc_ev_per_hour[hour_idx] = ev_fast.current_soc_percentage() # save begin state - if ev_charge_hours_fast[hour] > 0: + # save begin state + if ev_soc_direct: + soc_ev_per_hour[hour_idx] = (ev_fast.soc_wh / ev_fast.capacity_wh) * 100 + else: + soc_ev_per_hour[hour_idx] = ev_fast.current_soc_percentage() + ev_charge_factor = ev_charge_hours_fast[hour] + if ev_charge_factor > 0: stored_energy_ev, verluste_eauto = ev_fast.charge_energy( - wh=None, hour=hour, charge_factor=ev_charge_hours_fast[hour] + wh=None, hour=hour, charge_factor=ev_charge_factor ) # The inverter/grid must supply the EV charger's raw input, # not only the energy stored after charging losses. consumption += stored_energy_ev + verluste_eauto - losses_wh_per_hour[hour_idx] += verluste_eauto + slot_losses += verluste_eauto # Save battery SOC before inverter processing = true begin-of-interval state. # Must be recorded here (before DC charge/discharge) so the displayed SOC at # timestamp T reflects what the battery actually had at the START of interval T, # not the post-DC result. Consistent with the EV SOC convention above. if battery_fast: - soc_per_hour[hour_idx] = battery_fast.current_soc_percentage() + if battery_soc_direct: + soc_per_hour[hour_idx] = (battery_fast.soc_wh / battery_fast.capacity_wh) * 100 + else: + soc_per_hour[hour_idx] = battery_fast.current_soc_percentage() # Process inverter logic energy_feedin_grid_actual = energy_consumption_grid_actual = losses = eigenverbrauch = ( @@ -675,7 +772,8 @@ class GeneticSimulation(PydanticBaseModel): ac_charge_factor = 0.0 if battery_fast and ac_charging_possible: ac_charge_factor = ac_charge_hours_fast[hour] - if inverter_fast: + # No AC charging in this slot (the usual case): nothing to cap. + if inverter_fast and ac_charge_factor != 0: ac_charge_factor = inverter_fast.ac_charge_factor(ac_charge_factor) if inverter_fast: @@ -703,8 +801,8 @@ class GeneticSimulation(PydanticBaseModel): energy_produced, consumption, hour, - allow_battery_grid_export=battery_grid_export_allowed, - battery_grid_export_factor=battery_grid_export_factor, + battery_grid_export_allowed, + battery_grid_export_factor, ) else: hourly_feed_in_tariff = elect_revenue_per_hour_arr_fast[hour] @@ -724,7 +822,7 @@ class GeneticSimulation(PydanticBaseModel): ac_energy = stored + ac_charge_losses consumption += ac_energy energy_consumption_grid_actual += ac_energy - losses_wh_per_hour[hour_idx] += ac_charge_losses + slot_losses += ac_charge_losses # Update hourly arrays if ( @@ -732,12 +830,12 @@ class GeneticSimulation(PydanticBaseModel): and hourly_feed_in_tariff < 0.0 and energy_feedin_grid_actual > 0.0 ): - losses_wh_per_hour[hour_idx] += energy_feedin_grid_actual + slot_losses += energy_feedin_grid_actual energy_feedin_grid_actual = 0.0 feedin_energy_per_hour[hour_idx] = energy_feedin_grid_actual consumption_energy_per_hour[hour_idx] = energy_consumption_grid_actual - losses_wh_per_hour[hour_idx] += losses + losses_wh_per_hour[hour_idx] = slot_losses + losses loads_energy_per_hour[hour_idx] = consumption hourly_electricity_price = elect_price_hourly_fast[hour] electricity_price_per_hour[hour_idx] = hourly_electricity_price @@ -749,15 +847,36 @@ class GeneticSimulation(PydanticBaseModel): # not charged on input energy, internal discharge losses, or the # downstream DC-to-AC inverter loss. battery_lcos_cost = 0.0 - if battery_fast: + if battery_lcos_direct: + # Battery.discharged_energy_wh(), without the call per slot. battery_lcos_cost = ( - battery_fast.discharged_energy_wh(hour) - * battery_fast.levelized_cost_of_storage_kwh + battery_discharged_raw_fast[hour] + * battery_discharging_efficiency_fast + * battery_lcos_per_wh_factor / 1000.0 ) + elif battery_fast: + battery_lcos_cost = ( + battery_fast.discharged_energy_wh(hour) * battery_lcos_per_wh_factor / 1000.0 + ) costs_per_hour[hour_idx] = grid_cost + battery_lcos_cost revenue_per_hour[hour_idx] = energy_feedin_grid_actual * hourly_feed_in_tariff + loads_energy_per_hour = np.array(loads_energy_per_hour, dtype=float) + feedin_energy_per_hour = np.array(feedin_energy_per_hour, dtype=float) + consumption_energy_per_hour = np.array(consumption_energy_per_hour, dtype=float) + costs_per_hour = np.array(costs_per_hour, dtype=float) + revenue_per_hour = np.array(revenue_per_hour, dtype=float) + losses_wh_per_hour = np.array(losses_wh_per_hour, dtype=float) + electricity_price_per_hour = np.array(electricity_price_per_hour, dtype=float) + feed_in_tariff_per_hour = np.array(feed_in_tariff_per_hour, dtype=float) + if battery_fast: + soc_per_hour = np.array(soc_per_hour, dtype=float) + if ev_fast: + soc_ev_per_hour = np.array(soc_ev_per_hour, dtype=float) + if home_appliance_enabled: + home_appliance_wh_per_hour = np.array(home_appliance_wh_per_hour, dtype=float) + total_cost = np.nansum(costs_per_hour) total_losses = np.nansum(losses_wh_per_hour) total_revenue = np.nansum(revenue_per_hour) @@ -785,6 +904,7 @@ class GeneticSimulation(PydanticBaseModel): class GeneticOptimization(OptimizationBase): """GENETIC algorithm to solve energy optimization.""" + APPLIANCE_START_CACHE_SIZE = 4096 WARM_START_COPIES = 10 WARM_START_MUTATIONS = 50 EDUCATED_GUESS_TARGET = 100 @@ -826,8 +946,19 @@ class GeneticOptimization(OptimizationBase): # evaluation fails with an IndexError. The snapshot is taken in # optimize_ems() and released when the run ends; new values apply to the # next run. + # Energy-shift targets per source slot of the current run + # (see _energy_shift_target_slots). + _energy_shift_cache: Any = None + # Decoded appliance starts and load curves per appliance gene combination + # of the current layout (see _apply_appliance_starts). + _appliance_start_cache: Any = None + # decode_charge_discharge() per battery state of the current layout. + _decode_table: Any = None _run_genetic_cfg: Any = None _run_prediction_hours: Optional[float] = None + # control_end_slot of the frozen configuration. It is read several hundred + # thousand times per run and cannot change while the snapshot is in place. + _run_control_end_slot: Optional[int] = None @property def _genetic_cfg(self) -> Any: @@ -881,6 +1012,8 @@ class GeneticOptimization(OptimizationBase): @property def control_end_slot(self) -> int: """Exclusive control end in run-relative device arrays.""" + if self._run_control_end_slot is not None: + return self._run_control_end_slot return self._control_start_slot() + self.control_slots def _control_start_slot(self) -> int: @@ -1232,6 +1365,24 @@ class GeneticOptimization(OptimizationBase): logger.warning(self._terminal_value_reason) return curve + def _terminal_value_credit(self, parameters: GeneticOptimizationParameters) -> float: + """Credit for the energy left in the battery, without the report. + + The same credit as ``_terminal_value()`` returns. Every candidate needs + the amount only; the report object (component values, diagnostics, + curves) is built once, for the final solution. + """ + battery = self.simulation.battery + if battery is None: + return 0.0 + energy_wh = battery.current_energy_content() + if self.simulation.inverter: + energy_wh *= self.simulation.inverter.reference_dc_to_ac_efficiency + curve = getattr(self, "_terminal_value_curve", None) + if curve is not None and curve.energy_wh: + return curve.value(energy_wh) + return energy_wh * parameters.ems.price_per_wh_battery + def _terminal_value( self, parameters: GeneticOptimizationParameters, @@ -1511,9 +1662,33 @@ class GeneticOptimization(OptimizationBase): """Build every appliance's load curve from the decoded starts.""" if not self.simulation.home_appliances: return + # The decoded starts, the canonical gene values and the load curves + # depend on the appliance genes only, and a run sees few different + # combinations of them among tens of thousands of candidates. Each + # combination is decoded and built once per run; afterwards the + # appliances get the stored starts and curves (never changed in place: + # a rebuild starts from a new array). + appliances = self.simulation.home_appliances + cache = self._appliance_start_cache + if cache is None or cache[0] is not self.appliance_layout: + entries: dict[tuple[int, ...], Any] = {} + cache = (self.appliance_layout, entries) + self._appliance_start_cache = cache + key = tuple(int(value) for value in appliance_gene_values) + entry = cache[1].get(key) + if entry is not None and len(entry[1]) == len(appliances): + appliance_gene_values[:] = entry[0] + for appliance, (start_hours, load_curve) in zip(appliances, entry[1]): + appliance.start_hours = list(start_hours) + appliance.load_curve = load_curve + return starts_per_appliance = self._decode_appliance_starts(appliance_gene_values) - for appliance_index, appliance in enumerate(self.simulation.home_appliances): + built = [] + for appliance_index, appliance in enumerate(appliances): appliance.build_load_curve(starts_per_appliance.get(appliance_index, [])) + built.append((tuple(appliance.start_hours), appliance.load_curve)) + if len(cache[1]) < self.APPLIANCE_START_CACHE_SIZE: + cache[1][key] = (list(appliance_gene_values), built) def _start_solution_matches_layout(self, start_solution: list[float]) -> bool: """Check that a start solution's appliance tail fits the current layout. @@ -1778,7 +1953,41 @@ class GeneticOptimization(OptimizationBase): def decode_charge_discharge( self, discharge_hours_bin: np.ndarray ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: - """Decode the input array into charge, self-consumption discharge and export arrays.""" + """Decode the input array into charge, self-consumption discharge and export arrays. + + Every slot is decoded on its own from its state. The four results per + state are therefore worked out once per state layout + (``_decode_states``) and looked up for the slots of a candidate. + """ + states = np.asarray(discharge_hours_bin) + if states.ndim == 1 and states.dtype.kind in "iu": + layout = self._battery_state_layout() + key = ( + layout, + tuple(self.bat_possible_charge_values), + tuple(self.bat_possible_grid_export_values), + ) + table = self._decode_table + if table is None or table[0] != key: + decoded = self._decode_states(np.arange(max(layout.total_states, 1))) + table = (key, decoded) + self._decode_table = table + lookup = table[1] + # States outside the table (never produced by the optimizer) keep + # the explicit decoding. + if states.size and 0 <= states.min() and states.max() < lookup[0].size: + return ( + lookup[0][states], + lookup[1][states], + lookup[2][states], + lookup[3][states], + ) + return self._decode_states(discharge_hours_bin) + + def _decode_states( + self, discharge_hours_bin: np.ndarray + ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Decode states into charge, self-consumption discharge and export arrays.""" discharge_hours_bin_np = np.array(discharge_hours_bin) # Battery AC charge uses its own charge-level list (bat_possible_charge_values). len_bat = len(self.bat_possible_charge_values) @@ -1864,30 +2073,54 @@ class GeneticOptimization(OptimizationBase): source_slot: int, ) -> list[int]: """Return later idle slots where retained battery energy avoids costly import.""" - try: - prices = np.asarray(self.simulation.elect_price_hourly, dtype=float) - feed_in = np.asarray(self.simulation.elect_revenue_per_hour_arr, dtype=float) - pv = np.asarray(self.simulation.pv_prediction_wh, dtype=float) - load = np.asarray(self.simulation.load_energy_array, dtype=float) - except Exception: - return [] - if any(values.size < self.control_end_slot for values in (prices, feed_in, pv, load)): + # Which later slots are worth it, and in which order, depends on the + # forecasts and the source slot only - not on the individual. That + # part is worked out once per source slot and run (this is called for + # every export slot of every energy-shift mutation); the individual + # only filters the ordered slots. The cache is tied to the forecast + # arrays and the control end, and dropped at the start of a run. + simulation = self.simulation + sources = ( + simulation.elect_price_hourly, + simulation.elect_revenue_per_hour_arr, + simulation.pv_prediction_wh, + simulation.load_energy_array, + ) + control_end_slot = self.control_end_slot + cache = self._energy_shift_cache + if ( + cache is None + or cache[0] != control_end_slot + or any(known is not current for known, current in zip(cache[1], sources)) + ): + try: + arrays = tuple(np.asarray(values, dtype=float) for values in sources) + except Exception: + return [] + ordered_by_source: dict[int, list[int]] = {} + cache = (control_end_slot, sources, arrays, ordered_by_source) + self._energy_shift_cache = cache + prices, feed_in, pv, load = cache[2] + if any(values.size < control_end_slot for values in (prices, feed_in, pv, load)): return [] + ordered = cache[3].get(source_slot) + if ordered is None: + source_tariff = float(feed_in[source_slot]) + candidates = [ + slot + for slot in range(source_slot + 1, control_end_slot) + if load[slot] > pv[slot] and prices[slot] > source_tariff + ] + ordered = sorted( + candidates, + key=lambda slot: (float(prices[slot]), float(load[slot] - pv[slot])), + reverse=True, + ) + cache[3][source_slot] = ordered + len_bat = len(self.bat_possible_charge_values) - source_tariff = float(feed_in[source_slot]) - candidates = [ - slot - for slot in range(source_slot + 1, self.control_end_slot) - if 0 <= int(individual[slot]) < len_bat - and load[slot] > pv[slot] - and prices[slot] > source_tariff - ] - return sorted( - candidates, - key=lambda slot: (float(prices[slot]), float(load[slot] - pv[slot])), - reverse=True, - ) + return [slot for slot in ordered if 0 <= int(individual[slot]) < len_bat] def _mutate_energy_shift(self, individual: list[int]) -> bool: """Move battery energy from a weak export into later expensive self-consumption.""" @@ -2142,21 +2375,26 @@ class GeneticOptimization(OptimizationBase): if abs(ev_possible_charge_values[zero_charge_index]) > 1e-12: return False - _, ev_charge_indices, _ = self.split_individual(individual) - if ev_charge_indices is None: - return False - ev_soc = np.asarray(simulation_result.get("EAuto_SoC_pro_Stunde", []), dtype=float) start_slot = self._control_start_slot() result_slots = min(ev_soc.size, self.control_end_slot - start_slot) if result_slots <= 0: return False + # Nearly every candidate has no slot that begins at full SoC: decide + # that before the genome is split. + full = ev_soc[:result_slots] >= 100.0 - 1e-9 + if not full.any(): + return False + + _, ev_charge_indices, _ = self.split_individual(individual) + if ev_charge_indices is None: + return False changed = False - for offset in range(result_slots): - slot = start_slot + offset + for offset in np.flatnonzero(full): + slot = start_slot + int(offset) charge_index = int(ev_charge_indices[slot]) - if ev_soc[offset] >= 100.0 - 1e-9 and ev_possible_charge_values[charge_index] > 0.0: + if ev_possible_charge_values[charge_index] > 0.0: ev_charge_indices[slot] = zero_charge_index changed = True @@ -2590,28 +2828,73 @@ class GeneticOptimization(OptimizationBase): ) self._stop_evaluation_pool() pending = [individual for individual in invalid if not individual.fitness.valid] - fitnesses = self.toolbox.map(self.toolbox.evaluate, pending) - for individual, fitness in zip(pending, fitnesses): - individual.fitness.values = fitness + if getattr(self, "_fitness_cache_enabled", False): + # evaluate() answers repeated genomes from the cache. + fitnesses = self.toolbox.map(self.toolbox.evaluate, pending) + for individual, fitness in zip(pending, fitnesses): + individual.fitness.values = fitness + return len(invalid) + # Without the cache, identical genomes of this batch are still + # evaluated once only (see _group_identical). + groups = self._group_identical(pending) + firsts = [group[0] for group in groups] + fitnesses = self.toolbox.map(self.toolbox.evaluate, firsts) + for group, fitness in zip(groups, fitnesses): + first = group[0] + first.fitness.values = fitness + for individual in group[1:]: + self._copy_evaluation(first, individual) return len(invalid) + def _group_identical(self, individuals: list[Any]) -> list[list[Any]]: + """Group individuals with identical genomes, in order of first occurrence. + + Crossover and mutation regularly produce the same genome several times + within one generation. The evaluation is deterministic, so one + evaluation per group gives every member the result it would get on its + own. Unlike the fitness cache this needs no memory beyond the batch. + """ + groups: list[list[Any]] = [] + by_genome: dict[tuple[Any, ...], list[Any]] = {} + for individual in individuals: + genome = tuple(individual) + group = by_genome.get(genome) + if group is None: + group = by_genome[genome] = [] + groups.append(group) + else: + self._duplicate_evaluations_skipped = ( + getattr(self, "_duplicate_evaluations_skipped", 0) + 1 + ) + group.append(individual) + return groups + + @staticmethod + def _copy_evaluation(source: Any, target: Any) -> None: + """Give ``target`` the evaluated genome and results of ``source``.""" + target[:] = source + if hasattr(source, "extra_data"): + target.extra_data = source.extra_data + elif hasattr(target, "extra_data"): + del target.extra_data + target.fitness.values = source.fitness.values + def _evaluate_parallel(self, pool: Any, invalid: list[Any]) -> None: """Evaluate individuals in the worker processes, same results as evaluate(). - Cache lookups and stores stay in this process. With the cache active, - identical genomes are evaluated once, as a serial run would hit the - cache for the repeats. Results are applied only when the whole batch + Cache lookups and stores stay in this process. Identical genomes are + evaluated once: with the cache active as a serial run would hit the + cache for the repeats, without it by grouping the batch. Results are applied only when the whole batch has returned, so a failed batch leaves nothing half-assigned. """ cache_enabled = getattr(self, "_fitness_cache_enabled", False) groups: list[list[Any]] = [] keys: list[Optional[PackedGenes]] = [] by_key: dict[PackedGenes, list[Any]] = {} - for individual in invalid: - if not cache_enabled: - groups.append([individual]) - keys.append(None) - continue + if not cache_enabled: + groups = self._group_identical(invalid) + keys = [None] * len(groups) + for individual in invalid if cache_enabled else (): key = self._fitness_key(individual) cached = self._fitness_cache.get(key) if cached is not None: @@ -2844,15 +3127,21 @@ class GeneticOptimization(OptimizationBase): ) -> list[Any]: """Create offspring where crossover and mutation can both be applied.""" offspring: list[Any] = [] + # Offspring only take the genes of their parents: fitness, extra_data and + # immigrant protection are dropped by _invalidate_individual() below. A + # genome is a flat list of ints, so a new Individual from the parent's + # genes is the same child as a deepcopy that is invalidated afterwards - + # without copy.deepcopy walking every gene (8 % of a run). + clone = creator.Individual for _ in range(count): - child = self.toolbox.clone(random.choice(population)) # noqa: S311 + child = clone(random.choice(population)) # noqa: S311 crossed = False if ( len(child) > 1 and len(population) > 1 and random.random() < self.CROSSOVER_PROBABILITY # noqa: S311 ): - partner = self.toolbox.clone(random.choice(population)) # noqa: S311 + partner = clone(random.choice(population)) # noqa: S311 child, _ = self.toolbox.mate(child, partner) crossed = True @@ -3220,14 +3509,25 @@ class GeneticOptimization(OptimizationBase): def _fitness_key(self, individual: list[int]) -> PackedGenes: """Return the fitness-relevant genome, excluding elapsed control slots.""" start_slot = self._control_start_slot() - relevant = list(individual[start_slot : self.control_end_slot]) - if self.optimize_ev: - ev_start = self.control_end_slot + start_slot - relevant.extend(individual[ev_start : self.control_end_slot * 2]) + control_end_slot = self.control_end_slot n_appliance_genes = self.appliance_layout.n_genes - if n_appliance_genes > 0: - relevant.extend(individual[-n_appliance_genes:]) - return _pack_genes([int(value) for value in relevant]) + if start_slot == 0 and len(individual) == ( + control_end_slot * (2 if self.optimize_ev else 1) + n_appliance_genes + ): + # The whole genome is relevant (the usual case): no copy needed. + relevant = individual + else: + relevant = list(individual[start_slot:control_end_slot]) + if self.optimize_ev: + ev_start = control_end_slot + start_slot + relevant.extend(individual[ev_start : control_end_slot * 2]) + if n_appliance_genes > 0: + relevant.extend(individual[-n_appliance_genes:]) + try: + # Genes are plain ints; bytes() rejects anything else. + return _pack_genes(relevant) + except (TypeError, ValueError): + return _pack_genes([int(value) for value in relevant]) def _ev_soc_at_deadline(self, simulation_result: dict[str, Any], start_slot: int) -> float: """EV state of charge the target is checked against [%]. @@ -3363,8 +3663,7 @@ class GeneticOptimization(OptimizationBase): # replaces the most expensive hour after the horizon, the last one # replaces nothing. A scalar cannot express that (see terminalvalue.py). if self.simulation.battery: - restwert_akku, _ = self._terminal_value(parameters) - gesamtbilanz += -restwert_akku + gesamtbilanz += -self._terminal_value_credit(parameters) # --- AC charging break-even penalty --- # Penalise AC charging decisions that cannot be economically justified given the @@ -3616,6 +3915,8 @@ class GeneticOptimization(OptimizationBase): self._fitness_cache_hits = 0 self._fitness_cache_misses = 0 self._fitness_cache_enabled = max_entries != 0 + self._duplicate_evaluations_skipped = 0 + self._energy_shift_cache = None self._start_evaluation_pool() local_evaluations = 0 local_improvements = 0 @@ -3667,7 +3968,11 @@ class GeneticOptimization(OptimizationBase): cache_hit_rate = self._fitness_cache_hits / cache_lookups if cache_lookups > 0 else 0.0 cache_keys = len(self._fitness_cache) if max_entries == 0: - logger.info("Fitness cache: disabled (fitness_cache_max_entries = 0).") + logger.info( + "Fitness cache: disabled (fitness_cache_max_entries = 0); " + "{} identical genomes within a generation evaluated once.", + getattr(self, "_duplicate_evaluations_skipped", 0), + ) elif max_entries is None: logger.info( "Fitness cache: {} hits, {} misses, {:.1%} hit rate, {} keys.", @@ -3741,11 +4046,13 @@ class GeneticOptimization(OptimizationBase): self.config.validate_optimization_horizons() self._run_genetic_cfg = self.config.optimization.genetic.model_copy(deep=True) self._run_prediction_hours = self.config.prediction.hours + self._run_control_end_slot = self._control_start_slot() + self.control_slots try: return self._optimize_ems_run(parameters, start_hour, worst_case, ngen, individuals) finally: self._run_genetic_cfg = None self._run_prediction_hours = None + self._run_control_end_slot = None def _optimize_ems_run( self, diff --git a/src/akkudoktoreos/optimization/genetic/tailvalue.py b/src/akkudoktoreos/optimization/genetic/tailvalue.py index 316bcb32..f5b68c21 100644 --- a/src/akkudoktoreos/optimization/genetic/tailvalue.py +++ b/src/akkudoktoreos/optimization/genetic/tailvalue.py @@ -48,9 +48,9 @@ def _simulate_action( """Apply one tail action from one stored-energy state.""" dc, discharge, ac_rate, export = action bat.soc_wh = float(energy_wh) - bat._charged_raw_wh_per_slot.fill(0) - bat._discharged_raw_wh_per_slot.fill(0) - bat._charge_limit_raw_wh_per_slot.fill(np.inf) + bat._charged_raw_wh_per_slot[:] = [0.0] * len(bat._charged_raw_wh_per_slot) + bat._discharged_raw_wh_per_slot[:] = [0.0] * len(bat._discharged_raw_wh_per_slot) + bat._charge_limit_raw_wh_per_slot[:] = [float("inf")] * len(bat._charge_limit_raw_wh_per_slot) ac_enabled = inv.ac_to_dc_efficiency > 0 and ( inv.max_ac_charge_power_w is None or inv.max_ac_charge_power_w > 0 ) diff --git a/src/akkudoktoreos/prediction/interpolator.py b/src/akkudoktoreos/prediction/interpolator.py index d1a52026..92bbafad 100644 --- a/src/akkudoktoreos/prediction/interpolator.py +++ b/src/akkudoktoreos/prediction/interpolator.py @@ -5,7 +5,7 @@ from typing import Optional import numpy as np -from akkudoktoreos.core.cache import cache_energy_management +from akkudoktoreos.core.cache import CacheEnergyManagementStore, cache_energy_management from akkudoktoreos.core.coreabc import SingletonMixin from akkudoktoreos.utils.gridinterpolator import ( SwitchableGridInterpolator, @@ -15,6 +15,9 @@ from akkudoktoreos.utils.gridinterpolator import ( class SelfConsumptionProbabilityInterpolator: + LOAD_DISTRIBUTION_CACHE_SIZE = 4096 + DIRECT_CONSUMPTION_CACHE_SIZE = 8192 + def __init__(self, filepath: str | Path): self.filepath = filepath # The table is a pickled SciPy RegularGridInterpolator. It is read @@ -27,6 +30,13 @@ class SelfConsumptionProbabilityInterpolator: self.load_power_max_w = float(self.interpolator.grid[0][-1]) self.minute_load_levels_w = np.asarray(self.interpolator.grid[1], dtype=float) self.minute_load_max_w = float(self.interpolator.grid[1][-1]) + # _load_distribution() per (backend, bounded mean load). The table is + # static, so an entry never goes stale; the size is bounded. + self._load_distributions: dict[tuple[str, float], tuple[np.ndarray, np.ndarray]] = {} + # calculate_expected_direct_consumption() per (mean load, PV power), + # valid for one energy management run (see there). + self._direct_consumption: dict[tuple[float, float], float] = {} + self._direct_consumption_generation = CacheEnergyManagementStore.generation def _load_distribution(self, mean_load_power_w: float) -> tuple[np.ndarray, np.ndarray]: """Return the conditional minute-load distribution for a mean load. @@ -36,9 +46,18 @@ class SelfConsumptionProbabilityInterpolator: numerical deviations, so negative masses are removed and the result is normalized explicitly. """ - bounded_mean_load_w = float( - np.clip(mean_load_power_w, self.load_power_min_w, self.load_power_max_w) + # The distribution depends on the mean load only, not on the PV power. + # An optimization asks for the same few loads (one per forecast slot, + # plus appliance combinations) with ever-changing PV values, so the + # table lookup over all minute-load bins is done once per load. Callers + # must not modify the returned arrays. + bounded_mean_load_w = min( + max(float(mean_load_power_w), self.load_power_min_w), self.load_power_max_w ) + key = (str(configured_grid_interpolator_backend()), bounded_mean_load_w) + cached = self._load_distributions.get(key) + if cached is not None: + return cached points = np.column_stack( ( np.full(self.minute_load_levels_w.shape, bounded_mean_load_w), @@ -48,8 +67,14 @@ class SelfConsumptionProbabilityInterpolator: probabilities = np.maximum(np.asarray(self.interpolator(points), dtype=float), 0.0) probability_sum = float(probabilities.sum()) if probability_sum <= 0.0: - return self.minute_load_levels_w, probabilities - return self.minute_load_levels_w, probabilities / probability_sum + distribution = (self.minute_load_levels_w, probabilities) + else: + distribution = (self.minute_load_levels_w, probabilities / probability_sum) + if bounded_mean_load_w == bounded_mean_load_w: # not NaN: NaN never matches a key + if len(self._load_distributions) >= self.LOAD_DISTRIBUTION_CACHE_SIZE: + self._load_distributions.clear() + self._load_distributions[key] = distribution + return distribution def _generate_points( self, mean_load_power_w: float, pv_power_w: float @@ -91,7 +116,6 @@ class SelfConsumptionProbabilityInterpolator: probabilities = self.interpolator(points) return float(np.clip(probabilities.sum(), 0.0, 1.0)) - @cache_energy_management def calculate_expected_direct_consumption( self, mean_load_power_w: float, pv_power_w: float ) -> float: @@ -106,6 +130,13 @@ class SelfConsumptionProbabilityInterpolator: forecast mean exactly. This compensates for discretization and the finite upper table boundary while retaining the distribution shape. + The results are cached until the start of the next energy management + run, like ``cache_energy_management`` does. The optimizer asks several + hundred thousand times per run, so the cache is a plain dictionary of + this instance instead of the shared store: no key built from the + callable, no callback per call, and room for all slots of a run instead + of the last 100 results. + Args: mean_load_power_w: Mean load power of the forecast interval [W]. pv_power_w: Mean PV power of the forecast interval [W]. @@ -113,6 +144,22 @@ class SelfConsumptionProbabilityInterpolator: Returns: Expected direct PV-to-load power [W]. """ + cache = self._direct_consumption + if self._direct_consumption_generation != CacheEnergyManagementStore.generation: + cache.clear() + self._direct_consumption_generation = CacheEnergyManagementStore.generation + key = (mean_load_power_w, pv_power_w) + cached = cache.get(key) + if cached is not None: + return cached + result = self._expected_direct_consumption(mean_load_power_w, pv_power_w) + if len(cache) >= self.DIRECT_CONSUMPTION_CACHE_SIZE: + cache.clear() + cache[key] = result + return result + + def _expected_direct_consumption(self, mean_load_power_w: float, pv_power_w: float) -> float: + """Calculate the expected direct PV-to-load power in watts, uncached.""" mean_load_power_w = max(float(mean_load_power_w), 0.0) pv_power_w = max(float(pv_power_w), 0.0) if mean_load_power_w == 0.0 or pv_power_w == 0.0: @@ -129,7 +176,8 @@ class SelfConsumptionProbabilityInterpolator: expected_direct_power_w = float( np.dot(probabilities, np.minimum(normalized_load_levels_w, pv_power_w)) ) - return float(np.clip(expected_direct_power_w, 0.0, min(mean_load_power_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)) # def calculate_self_consumption(self, load_1h_power: float, pv_power: float) -> float: # """Calculate the PV self-consumption rate using RegularGridInterpolator. diff --git a/tests/test_genetic_parallel_evaluation.py b/tests/test_genetic_parallel_evaluation.py index bb6e19fe..c3db9176 100644 --- a/tests/test_genetic_parallel_evaluation.py +++ b/tests/test_genetic_parallel_evaluation.py @@ -1,6 +1,7 @@ """Parallel fitness evaluation gives the same result as evaluating in one process.""" import multiprocessing +from types import SimpleNamespace import pytest @@ -203,3 +204,47 @@ def test_worker_cpu_set(workers, allowed, expected): from akkudoktoreos.optimization.genetic.genetic import worker_cpu_set assert worker_cpu_set(workers, allowed) == expected + + +class _Fitness: + def __init__(self) -> None: + self.values: tuple[float, ...] = () + + @property + def valid(self): + return bool(self.values) + + +class _Individual(list): + extra_data: tuple[float, float, float] + + def __init__(self, genes): + super().__init__(genes) + self.fitness = _Fitness() + + +def test_identical_genomes_of_a_batch_are_evaluated_once_without_cache(): + """Repeats get the repaired genome and results of the one evaluation.""" + optimization = GeneticOptimization.__new__(GeneticOptimization) + evaluated: list[list[int]] = [] + + def evaluate(individual): + evaluated.append(list(individual)) + individual[0] = 9 # an evaluation may repair the genome in place + individual.extra_data = (float(sum(individual)), 0.0, 0.0) + return (float(sum(individual)),) + + optimization.toolbox = SimpleNamespace(map=map, evaluate=evaluate) + population = [_Individual(genes) for genes in ([1, 2, 3], [4, 5, 6], [1, 2, 3], [7, 8, 9])] + done = _Individual([0, 0, 0]) + done.fitness.values = (1.0,) + + assert optimization._evaluate_invalid(population + [done]) == 4 + + assert evaluated == [[1, 2, 3], [4, 5, 6], [7, 8, 9]] + assert optimization._duplicate_evaluations_skipped == 1 + assert population[2] == population[0] == [9, 2, 3] + assert population[2].fitness.values == population[0].fitness.values == (14.0,) + assert population[2].extra_data == population[0].extra_data + assert population[1].fitness.values == (20.0,) + assert done == [0, 0, 0] and done.fitness.values == (1.0,) diff --git a/tests/test_interpolator.py b/tests/test_interpolator.py index df3cbfb5..ccb9a0b9 100644 --- a/tests/test_interpolator.py +++ b/tests/test_interpolator.py @@ -1,6 +1,7 @@ import numpy as np import pytest +from akkudoktoreos.core.cache import CacheEnergyManagementStore from akkudoktoreos.prediction.interpolator import get_eos_load_interpolator @@ -103,3 +104,32 @@ def test_genetic_inverter_boundary_flows_remain_nonnegative(load, pv): inverter = Inverter(InverterParameters(device_id="boundary", max_power_wh=25000)) flows = inverter.process_energy(generation=pv, consumption=load, hour=0) assert all(np.isfinite(value) and value >= 0.0 for value in flows) + + +def test_expected_direct_consumption_cache_matches_uncached_and_follows_the_run_cache(): + """The per-run cache returns the computed value and is dropped with the EMS cache.""" + interpolator = get_eos_load_interpolator() + CacheEnergyManagementStore().clear() + + first = interpolator.calculate_expected_direct_consumption(800.0, 1200.0) + assert first == interpolator._expected_direct_consumption(800.0, 1200.0) + assert interpolator._direct_consumption[(800.0, 1200.0)] == first + assert interpolator.calculate_expected_direct_consumption(800.0, 1200.0) == first + + CacheEnergyManagementStore().clear() + assert interpolator.calculate_expected_direct_consumption(900.0, 1200.0) > 0.0 + assert (800.0, 1200.0) not in interpolator._direct_consumption + + +def test_expected_direct_consumption_cache_is_bounded(monkeypatch): + interpolator = get_eos_load_interpolator() + CacheEnergyManagementStore().clear() + monkeypatch.setattr(interpolator, "DIRECT_CONSUMPTION_CACHE_SIZE", 3) + + for pv_power_w in (500.0, 600.0, 700.0, 800.0, 900.0): + interpolator.calculate_expected_direct_consumption(800.0, pv_power_w) + + assert len(interpolator._direct_consumption) <= 3 + assert interpolator.calculate_expected_direct_consumption( + 800.0, 500.0 + ) == interpolator._expected_direct_consumption(800.0, 500.0)