Drop the free level deviation term again, leaving the peak alone

The measurements that justified it do not hold up: on the cases where the two
terms differ the search is truncated by the time limit, so the profile that
comes out is the incumbent CBC happens to hold, not what the term ordered. On
one production request the whole peak penalty is 8e-4 against a remaining gap
of 4.7e4, seven orders below anything the solver reacts to.

What is left is the ramp removal, which stands on its own.

Author: andig

README.md

@@ -37,11 +37,9 @@ The optimizer buys all 29 kWh of grid energy in the three cheapest hours of the
 
 Cheapest is not always kindest to the grid connection. The same house on a *flat* tariff — where nothing but the strategy decides when to charge — takes its 29 kWh at the full 11 kW in the last two hours before the EV deadline, and pushes the midday surplus out in two short bursts.
 
-`attenuate_grid_peaks` penalizes the highest grid power over the horizon and the distance of every single step from a common level, on the import and the feed-in side. The same energy then arrives as a flat 3.6 kW plateau, and the feed-in peak drops from 1.3 kW to 0.25 kW. Nothing costs more — the connection just sees a calmer profile. `attenuate_demand_peaks` and `attenuate_feedin_peaks` do the same for one side only.
+`attenuate_grid_peaks` penalizes the highest grid power over the horizon, on the import and the feed-in side. The same energy then arrives as a flat 3.6 kW plateau, and the feed-in peak drops from 1.3 kW to 0.25 kW. Nothing costs more — the connection just sees a calmer profile. `attenuate_demand_peaks` and `attenuate_feedin_peaks` do the same for one side only.
 
-The second term is what keeps the profile level where the peak is out of reach. A load spike the schedule cannot touch — an oven, a heat pump defrost — fixes the horizon maximum, and a peak penalty alone has nothing left to win below it: the cheapest way to charge is then flat out against the spike until the goal is met. Pricing the distance from a level instead spreads the same energy over the whole window, which is what a minimum square deviation from constant grid power asks for.
-
-The level is placed freely, which leaves one gap: on a side that mostly rests at zero the level settles at zero too, and the term then only adds up the energy through that side — a plateau, a jagged profile and a single spike of the same energy score alike, and nothing below the peak tells them apart.
+The maximum is a single value out of the horizon, which leaves one gap: a load spike the schedule cannot touch — an oven, a heat pump defrost — fixes it, and the penalty then has nothing left to win below it. Charging flat out against the spike scores the same as spreading the same energy over the window, and the solver may pick either.
 
 <picture>
   <source media="(prefers-color-scheme: dark)" srcset="docs/img/example-peak-dark.svg">

docs/comparison_objective_terms.md

@@ -14,7 +14,6 @@
 | | Final SOC value | `bat.p_a` | *(raw, now also feeds `penalty_base` directly)* | ~1.7e2 (up to whatever `p_a` is configured) | actual currency |
 | | Demand-rate charge | `grid.prc_p_exc_imp` | *(raw; `prc_p_exc_imp·3600/min(dt)` now also feeds `penalty_base`, only when `p_max_imp` makes the rate active)* | unbounded, but tier 1 now tracks it | actual currency |
 | **3 — tie-break nudges** | Peak leveling | `prc_p_peak = penalty_base·1e-3` | ×1e-3 | 0.16 – 300 | shape, not cost |
-| | Level deviation | `prc_p_dev = penalty_base·1e-4`, per step weighted by `dt[t]/Σdt` | ×1e-4·dt[t]/Σdt | 2.5e-5 – 2.5e-2 | shape, not cost |
 | | Charge-before-export | `min_import_price·2e-5·(T−t)` | ×2e-5·(T−t) | ~1e-3 – tens | shape, not cost |
 | | Discharge-before-import | `min_import_price·5e-6·(T−t)` | ×5e-6·(T−t) | ~1e-4 – tens | shape, not cost |
 | | Charging priority | `min_import_price·5e-5·(T−t)·c_priority`, on `c` and on `d` | ×5e-5·(T−t)·c_priority | ~1e-3 – tens | shape, not cost |
@@ -24,19 +23,15 @@
 - The SOC penalty is charged per time step, so a violation lasting the whole horizon is paid for
   once per step while one confined to a single step is paid for once. `penalty_base` is sized
   against the single-step case, which is the weaker of the two.
-- The two leveling terms are ordered per W, not per horizon. Peak leveling is charged once on a
-  single value, level deviation is a time average, so neither grows with the horizon and the order
-  between them holds at any sampling.
-- Level deviation carries the smallest coefficients in the model, a decade below peak leveling and
-  smaller again once the time average is spread over the steps. Level and deviation variables are
-  therefore bounded by the power the side can physically carry: unbounded they relax into the same
-  range as the big M of 1e6 and, on the largest stored case, cost the solver a factor of ten.
-- Level deviation has a blind spot. The level is free, and on a side that mostly rests at zero it
-  settles there; with `p_grid ≥ 0` the term is then `Σ p_grid[t]·dt[t]`, the energy through the
-  side, which the energy balance already fixes. A plateau, a jagged profile and a single spike of
-  the same energy score identically and only the peak separates them. A ramp leveling term at
-  `penalty_base·1e-5` used to order them; it was dropped because it priced transitions only and
-  cost around a third of the solve time of a leveling case. `test_peak_leveling` pins both halves.
+- Peak leveling is charged once on a single value, not per step, so its weight does not grow with
+  the horizon or change with the sampling.
+- Peak leveling has a blind spot. The maximum is one value out of the horizon, so once a load the
+  schedule cannot touch pins it, nothing orders the steps below it: a plateau, a jagged profile and
+  a single spike of the same energy score identically. A ramp leveling term at `penalty_base·1e-5`
+  used to order them; it was dropped because it priced transitions only, so it could not tell a
+  straight climb from a plateau either, and the second absolute value system over the same grid
+  power cost around a fifth of the solve time of a leveling case. `test_peak_leveling` pins both
+  halves.
 - The tier 3 terms keyed off `min_import_price` invert sign when market prices go negative, the
-  reason peak and level deviation leveling use `penalty_base` instead.
+  reason peak leveling uses `penalty_base` instead.
 

src/optimizer/optimizer.py

@@ -169,21 +169,12 @@ def __init__(self, strategy: OptimizationStrategy, grid: GridConfig, batteries:
         # solar power
         self.prc_e_grid_exp_pen = self.penalty_base * 10e1
 
-        # weights of the grid peak leveling strategies. Capping the horizon maximum alone leaves the
-        # profile below the cap arbitrary: once the maximum is pinned by a load the schedule cannot
-        # touch, the term has nothing left to win below it. The deviation from a free level closes
-        # that gap, it is what a mean square deviation from constant grid power would measure and it
-        # prices every single time step rather than one value out of the horizon.
-        # Per unit of W the peak outweighs the deviation, so lowering the peak still wins a direct
-        # trade against leveling the profile underneath it.
-        #
-        # A step to step ramp term sat between these two until the deviation term arrived. It priced
-        # the transitions only, so it could not tell a straight climb from a plateau, and carrying
-        # both absolute value systems over the same grid power cost around a third of the solve time
-        # of a leveling case for no measurable gain in deviation. See the blind spot noted in
-        # _add_energy_balance_constraints for what its removal gives up.
+        # weight of the grid peak leveling strategies. A step to step ramp term at penalty_base·1e-5
+        # sat next to it until this change. It priced the transitions only, so it could not tell a
+        # straight climb from a plateau, and the second absolute value system over the same grid
+        # power cost around a fifth of the solve time of a leveling case. See the blind spot noted
+        # in _add_energy_balance_constraints for what its removal gives up.
         self.prc_p_peak = self.penalty_base * 1e-3
-        self.prc_p_dev = self.penalty_base * 1e-4
 
         # grid sides leveled by the active peak attenuation strategy, empty for all other strategies
         self.peak_sides = PEAK_STRATEGY_SIDES.get(strategy.charging_strategy, ())
@@ -281,23 +272,9 @@ def _setup_variables(self):
             self.variables['p_max_imp_exc'] = pulp.LpVariable("p_max_imp_exc", lowBound=0)
 
         # highest grid power over the whole horizon (W) per side, used by the peak attenuation
-        # strategies.
-        # p_{side}_lvl is the constant the profile is leveled towards and p_{side}_dev the distance
-        # of each step from it (W). The level is a variable and not the mean of the profile, so the
-        # optimizer places it where the profile is cheapest to level, and nothing forces the schedule
-        # towards a level picked up front.
+        # strategies
         for side in self.peak_sides:
-            # the profile of a side cannot exceed what the horizon can physically carry there.
-            # bounding level and distance by it keeps the relaxation of the new rows from wandering
-            # off into big M territory, the same reason the flow direction constraints below prefer
-            # a natural cap over the global big M.
-            p_side_max = self._peak_side_bound(side)
             self.variables[f'p_{side}_peak'] = pulp.LpVariable(f"p_{side}_peak", lowBound=0)
-            self.variables[f'p_{side}_lvl'] = pulp.LpVariable(f"p_{side}_lvl", lowBound=0, upBound=p_side_max)
-            self.variables[f'p_{side}_dev'] = [
-                pulp.LpVariable(f"p_{side}_dev_{t}", lowBound=0, upBound=p_side_max)
-                for t in self.time_steps
-            ]
 
         # Binary variable: power flow direction to / from grid variables
         # these variables
@@ -326,30 +303,6 @@ def _setup_variables(self):
                 for t in self.time_steps
             ]
 
-    def _peak_side_bound(self, side):
-        """
-        Highest grid power the horizon can put on a leveled side (W), or the big M where there
-        is no natural cap.
-
-        Import is the demand of a step plus everything the batteries may pull from the grid,
-        export the production plus everything they may push back into it, both taken over the
-        step that carries the most. The flow direction constraints keep the two sides apart and
-        make those caps hold - except for the excess variable of a configured limit, which stays
-        outside them. With a limit on the opposite side that side can feed this one past its
-        natural cap, the same case where the flow direction constraints fall back to the big M.
-        """
-        if side == 'imp':
-            if self.grid.p_max_exp is not None:
-                return self.M
-            series = self.time_series.gt
-            cap = sum(bat.c_max for bat in self.batteries if bat.charge_from_grid)
-        else:
-            if self.grid.p_max_imp is not None:
-                return self.M
-            series = self.time_series.ft
-            cap = sum(bat.d_max for bat in self.batteries if bat.discharge_to_grid)
-        return max(series[t] * 3600. / self.time_series.dt[t] for t in self.time_steps) + cap
-
     def _setup_target_function(self):
         """
         Gather all target function contributions and instantiate the objective
@@ -438,21 +391,13 @@ def _setup_target_function(self):
 
         # level the grid profile to unload the public grid from peaks. attenuate_demand_peaks levels
         # grid import, attenuate_feedin_peaks levels grid export, attenuate_grid_peaks levels both.
-        # the penalty sits on the horizon maximum and on the distance of every step from a free
-        # level instead of on charge power, so the optimizer spreads charging at partial power over
-        # several time steps rather than running one step at full power, and keeps the profile below
-        # the cap leveled too.
+        # the penalty sits on the horizon maximum instead of on charge power, so the optimizer
+        # spreads charging at partial power over several time steps rather than running one step at
+        # full power.
         # penalty_base is used instead of min_import_price because negative market prices would turn
         # this penalty into a reward for peaks.
-        horizon = float(sum(self.time_series.dt))
         for side in self.peak_sides:
             objective += - self.variables[f'p_{side}_peak'] * self.prc_p_peak
-            # weighted by step length and divided by the horizon, which makes the term the time
-            # average of the distance from the level. A step average would let the same day weigh
-            # differently depending on how finely it is sampled, and would count a 15 min excursion
-            # like one lasting an hour.
-            objective += - pulp.lpSum(self.variables[f'p_{side}_dev'][t] * self.time_series.dt[t]
-                                      for t in self.time_steps) / horizon * self.prc_p_dev
 
         # prefer discharging batteries completely before importing from grid
         if self.strategy.discharging_strategy == 'discharge_before_import':
@@ -559,17 +504,14 @@ def _add_energy_balance_constraints(self):
                 self.problem += (self.variables['e_exp_lim_exc'][t]
                                  <= self.M * (1 - self.variables['z_exp_lim'][t]))
 
-        # track the horizon maximum and the distance from the level of the total grid power for every
-        # side the strategy levels. Both include the portion beyond p_max_imp / p_max_exp, so they
-        # stay correct in demand rate mode and when a limit is violated.
+        # track the horizon maximum of the total grid power for every side the strategy levels. It
+        # includes the portion beyond p_max_imp / p_max_exp, so it stays correct in demand rate mode
+        # and when a limit is violated.
         #
-        # Known blind spot: the level is free, and on a side that mostly rests at zero it settles at
-        # zero. With p_grid >= 0 the term is then sum(p_grid[t] * dt[t]), the energy through the
-        # side, which the energy balance already fixes - so it scores a plateau, a jagged profile and
-        # a single spike of the same energy identically and orders none of them. Only the peak still
-        # separates them there. The step to step ramp used to, and test_peak_leveling pins what that
-        # costs; on the stored corpus it is worth at most 8 W of deviation, and it is not worth a
-        # third of the solve time.
+        # Known blind spot: the maximum is a single value out of the horizon, so once it is pinned by
+        # a load the schedule cannot touch, nothing orders the steps below it any more - a plateau, a
+        # jagged profile and a single spike of the same energy all score alike. A step to step ramp
+        # term used to order them, and test_peak_leveling pins what dropping it costs.
         for side in self.peak_sides:
             grid_var, lim_exc_var = PEAK_SIDE_VARIABLES[side]
             # total grid power of this side per time step (W)
@@ -582,10 +524,6 @@ def _add_energy_balance_constraints(self):
 
             for t in self.time_steps:
                 self.problem += p_grid[t] <= self.variables[f'p_{side}_peak']
-            # distance from the level: p_dev[t] >= |p_grid[t] - p_lvl|
-            for t in self.time_steps:
-                self.problem += self.variables[f'p_{side}_dev'][t] >= p_grid[t] - self.variables[f'p_{side}_lvl']
-                self.problem += self.variables[f'p_{side}_dev'][t] >= self.variables[f'p_{side}_lvl'] - p_grid[t]
 
         # if demand rate is applied, the maximum grid import power value
         # of all time steps drives the demand rate charge

test_cases/027-attenuate-feedin-peaks-below-cap.json

@@ -69,10 +69,10 @@
             {
                 "charging_power": [
                     50.0,
-                    42.5,
-                    42.5,
-                    42.5,
-                    42.5
+                    50.0,
+                    50.0,
+                    50.0,
+                    20.0
                 ],
                 "discharging_power": [
                     0.0,
@@ -83,9 +83,9 @@
                 ],
                 "state_of_charge": [
                     50.0,
-                    92.5,
-                    135.0,
-                    177.5,
+                    100.0,
+                    150.0,
+                    200.0,
                     220.0
                 ]
             }
@@ -99,10 +99,10 @@
         ],
         "grid_export": [
             950.0,
-            57.5,
-            57.5,
-            57.5,
-            57.5
+            50.0,
+            50.0,
+            50.0,
+            80.0
         ],
         "flow_direction": [
             1,