From 533179adee981b4b853a66a1d9f29da13b2c4797 Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Sat, 22 Aug 2026 20:06:26 +0200 Subject: [PATCH 1/3] fix(mpc): replay every plan and represent a true zero PV cap Go fallback used to skip ValidatePlan, and forward simulation snapped SoC onto the DP grid, so reported energy did not replay. Python serialized full curtailment as pv_limit_w=0, the same sentinel as no cap, and applied the PV-charge bonus in every mode. Share BatteryEnergyDeltaWh, replay aggregate Go DP trajectories, start forward sim from real SoC, and add pv_curtail_active so a zero cap is a distinct contract. Gate the PV-charge bonus to passive_arbitrage. Leftover-PV and battery-to-EV residual constraints run per scenario. Signed-off-by: Fredrik Ahlgren --- .changeset/planning-physics-replay.md | 5 + go/internal/loadpoint/site_power.go | 24 +++ go/internal/loadpoint/site_power_test.go | 78 +++++++- .../loadpoint/testdata/site_physics.json | 85 +++++++++ go/internal/mpc/diagnose.go | 85 +++++---- go/internal/mpc/external_optimizer.go | 83 +++++---- go/internal/mpc/mpc.go | 109 +++++------ go/internal/mpc/service.go | 13 ++ go/internal/mpc/validate_dp_test.go | 175 ++++++++++++++++++ optimizer/ftw_optimizer/direct_highs.py | 7 +- optimizer/ftw_optimizer/model.py | 49 ++++- optimizer/ftw_optimizer/multistage.py | 22 +-- optimizer/ftw_optimizer/progressive.py | 14 +- optimizer/ftw_optimizer/recourse.py | 8 +- optimizer/tests/test_model.py | 43 ++++- optimizer/tests/test_site_physics.py | 65 +++++++ 16 files changed, 691 insertions(+), 174 deletions(-) create mode 100644 .changeset/planning-physics-replay.md create mode 100644 go/internal/loadpoint/testdata/site_physics.json create mode 100644 go/internal/mpc/validate_dp_test.go create mode 100644 optimizer/tests/test_site_physics.py diff --git a/.changeset/planning-physics-replay.md b/.changeset/planning-physics-replay.md new file mode 100644 index 000000000..7f38c72b0 --- /dev/null +++ b/.changeset/planning-physics-replay.md @@ -0,0 +1,5 @@ +--- +"ftw": patch +--- + +Go fallback plans are replayed against the same site-power and battery-energy identities as the mathematical optimizer before they can become the live plan. A true zero PV cap is now a distinct `pv_curtail_active` flag, so full curtailment is no longer serialized as “no cap”. A trajectory that cannot be reconstructed from the request is kept off dispatch. diff --git a/go/internal/loadpoint/site_power.go b/go/internal/loadpoint/site_power.go index cb4f00913..292ff1c9d 100644 --- a/go/internal/loadpoint/site_power.go +++ b/go/internal/loadpoint/site_power.go @@ -61,3 +61,27 @@ func BatteryDischargeFeedsEV(batteryW, evW, loadW, pvW float64) bool { func PlannedSurplusForEVW(loadW, pvW, batteryW, gridW float64) float64 { return -pvW - loadW - PlannedPVSoakW(batteryW, gridW) } + +// BatteryEnergyDeltaWh is the cell-side energy change for a site-signed +// AC battery power over dtH hours. Charge (powerW > 0) lands ηc of the +// AC energy in the cells. Discharge (powerW < 0) draws AC / ηd from the +// cells, so a 1000 W discharge at 0.95 efficiency removes ~1053 Wh/h. +func BatteryEnergyDeltaWh(powerW, dtH, chargeEff, dischargeEff float64) float64 { + if powerW >= 0 { + return powerW * dtH * chargeEff + } + return powerW * dtH / dischargeEff +} + +// EffectivePVW is the site-signed PV used in grid replay after an +// optional curtailment cap. Inactive means the forecast stands. +// Active with pvLimitW = 0 is a true zero cap (no generation this slot). +func EffectivePVW(pvW, pvLimitW float64, curtailActive bool) float64 { + if !curtailActive { + return pvW + } + if pvLimitW < 0 { + pvLimitW = 0 + } + return -pvLimitW +} diff --git a/go/internal/loadpoint/site_power_test.go b/go/internal/loadpoint/site_power_test.go index c98495dd0..040299bc0 100644 --- a/go/internal/loadpoint/site_power_test.go +++ b/go/internal/loadpoint/site_power_test.go @@ -1,6 +1,11 @@ package loadpoint -import "testing" +import ( + "encoding/json" + "math" + "os" + "testing" +) func TestGridWIncludesEV(t *testing.T) { // 500 house, 8 kW PV, 10 kW battery charge, 4.14 kW EV → import. @@ -51,6 +56,18 @@ func TestBatteryDischargeFeedsEV(t *testing.T) { } } +func TestEffectivePVWActiveZeroIsTrueCap(t *testing.T) { + if got := EffectivePVW(-5000, 0, false); got != -5000 { + t.Errorf("inactive = %.0f, want -5000", got) + } + if got := EffectivePVW(-5000, 0, true); got != 0 { + t.Errorf("active zero = %.0f, want 0", got) + } + if got := EffectivePVW(-5000, 2000, true); got != -2000 { + t.Errorf("partial cap = %.0f, want -2000", got) + } +} + func TestPlannedSurplusForEVWSkipsGridFundedCharge(t *testing.T) { // leftover 7500, battery soaking 2000 of it. if got := PlannedSurplusForEVW(500, -8000, 2000, 0); got != 5500 { @@ -65,3 +82,62 @@ func TestPlannedSurplusForEVWSkipsGridFundedCharge(t *testing.T) { t.Errorf("soak+EV: got %.0f, want 3500", got) } } + +type sitePhysicsFixture struct { + Flows []struct { + Name string `json:"name"` + LoadW float64 `json:"load_w"` + PVW float64 `json:"pv_w"` + BatteryW float64 `json:"battery_w"` + EVW float64 `json:"ev_w"` + GridW float64 `json:"grid_w"` + LeftoverW float64 `json:"leftover_w"` + HouseResidualW float64 `json:"house_residual_w"` + FeedsEV bool `json:"feeds_ev"` + } `json:"flows"` + EnergySteps []struct { + Name string `json:"name"` + PowerW float64 `json:"power_w"` + DtH float64 `json:"dt_h"` + ChargeEff float64 `json:"charge_eff"` + DischargeEff float64 `json:"discharge_eff"` + DeltaWh float64 `json:"delta_wh"` + } `json:"energy_steps"` +} + +func loadSitePhysicsFixture(t *testing.T) sitePhysicsFixture { + t.Helper() + raw, err := os.ReadFile("testdata/site_physics.json") + if err != nil { + t.Fatal(err) + } + var fixture sitePhysicsFixture + if err := json.Unmarshal(raw, &fixture); err != nil { + t.Fatal(err) + } + return fixture +} + +func TestSitePhysicsTable(t *testing.T) { + fixture := loadSitePhysicsFixture(t) + for _, row := range fixture.Flows { + if got := GridW(row.LoadW, row.PVW, row.BatteryW, row.EVW); got != row.GridW { + t.Errorf("%s: GridW = %.6g, want %.6g", row.Name, got, row.GridW) + } + if got := PVLeftoverAfterHouseW(row.LoadW, row.PVW); got != row.LeftoverW { + t.Errorf("%s: leftover = %.6g, want %.6g", row.Name, got, row.LeftoverW) + } + if got := HouseResidualW(row.LoadW, row.PVW); got != row.HouseResidualW { + t.Errorf("%s: residual = %.6g, want %.6g", row.Name, got, row.HouseResidualW) + } + if got := BatteryDischargeFeedsEV(row.BatteryW, row.EVW, row.LoadW, row.PVW); got != row.FeedsEV { + t.Errorf("%s: feeds EV = %v, want %v", row.Name, got, row.FeedsEV) + } + } + for _, row := range fixture.EnergySteps { + got := BatteryEnergyDeltaWh(row.PowerW, row.DtH, row.ChargeEff, row.DischargeEff) + if math.Abs(got-row.DeltaWh) > 1e-9 { + t.Errorf("%s: delta = %.12g, want %.12g", row.Name, got, row.DeltaWh) + } + } +} diff --git a/go/internal/loadpoint/testdata/site_physics.json b/go/internal/loadpoint/testdata/site_physics.json new file mode 100644 index 000000000..8adfbdf7c --- /dev/null +++ b/go/internal/loadpoint/testdata/site_physics.json @@ -0,0 +1,85 @@ +{ + "flows": [ + { + "name": "pv leftover after house", + "load_w": 500, + "pv_w": -8000, + "battery_w": 0, + "ev_w": 0, + "grid_w": -7500, + "leftover_w": 7500, + "house_residual_w": 0, + "feeds_ev": false + }, + { + "name": "house residual after weak pv", + "load_w": 2000, + "pv_w": -500, + "battery_w": 0, + "ev_w": 0, + "grid_w": 1500, + "leftover_w": 0, + "house_residual_w": 1500, + "feeds_ev": false + }, + { + "name": "leftover pv into ev and grid into battery", + "load_w": 500, + "pv_w": -8000, + "battery_w": 10000, + "ev_w": 4140, + "grid_w": 6640, + "leftover_w": 7500, + "house_residual_w": 0, + "feeds_ev": false + }, + { + "name": "battery discharge covers house only", + "load_w": 500, + "pv_w": 0, + "battery_w": -400, + "ev_w": 4000, + "grid_w": 4100, + "leftover_w": 0, + "house_residual_w": 500, + "feeds_ev": false + }, + { + "name": "battery discharge feeds ev", + "load_w": 500, + "pv_w": 0, + "battery_w": -4000, + "ev_w": 4000, + "grid_w": 500, + "leftover_w": 0, + "house_residual_w": 500, + "feeds_ev": true + } + ], + "energy_steps": [ + { + "name": "1 kW charge for 1 h at 0.95", + "power_w": 1000, + "dt_h": 1, + "charge_eff": 0.95, + "discharge_eff": 0.95, + "delta_wh": 950 + }, + { + "name": "1 kW discharge for 1 h at 0.95", + "power_w": -1000, + "dt_h": 1, + "charge_eff": 0.95, + "discharge_eff": 0.95, + "delta_wh": -1052.6315789473683 + }, + { + "name": "idle", + "power_w": 0, + "dt_h": 0.25, + "charge_eff": 0.95, + "discharge_eff": 0.95, + "delta_wh": 0 + } + ] +} diff --git a/go/internal/mpc/diagnose.go b/go/internal/mpc/diagnose.go index a3757ae8b..7cd39331a 100644 --- a/go/internal/mpc/diagnose.go +++ b/go/internal/mpc/diagnose.go @@ -36,13 +36,14 @@ type DiagnosticSlot struct { WeatherRowAvailableAtMs int64 `json:"weather_row_available_at_ms,omitempty"` // Outputs - BatteryW float64 `json:"battery_w"` - GridW float64 `json:"grid_w"` - SoC float64 `json:"soc"` // SoC at END of slot - CostOre float64 `json:"cost_ore"` // raw (un-blended) slot cost - Reason string `json:"reason"` - EMSMode string `json:"ems_mode"` - PVLimitW float64 `json:"pv_limit_w,omitempty"` + BatteryW float64 `json:"battery_w"` + GridW float64 `json:"grid_w"` + SoC float64 `json:"soc"` // 0–1 at END of slot + CostOre float64 `json:"cost_ore"` + Reason string `json:"reason"` + EMSMode string `json:"ems_mode"` + PVLimitW float64 `json:"pv_limit_w,omitempty"` + PVCurtailActive bool `json:"pv_curtail_active,omitempty"` // EV outputs — present only when the plan included a loadpoint. // `omitempty` + the web renderer's `lpActive` gate mean plans @@ -53,12 +54,12 @@ type DiagnosticSlot struct { // against `LOAD 1.6 kW` and reasonably assumed the battery was // exporting — reality was `LOAD 1.6 + EV 4.0 = 5.6 kW covered`, // grid ≈ 0. See issue #174. - LoadpointW float64 `json:"loadpoint_w,omitempty"` + LoadpointW float64 `json:"loadpoint_w,omitempty"` LoadpointSoC float64 `json:"loadpoint_soc,omitempty"` - LoadpointPowerW map[string]float64 `json:"loadpoint_power_w,omitempty"` + LoadpointPowerW map[string]float64 `json:"loadpoint_power_w,omitempty"` LoadpointSoCByID map[string]float64 `json:"loadpoint_soc_by_id,omitempty"` - StoragePowerW map[string]float64 `json:"storage_power_w,omitempty"` - StorageEnergyWh map[string]float64 `json:"storage_energy_wh,omitempty"` + StoragePowerW map[string]float64 `json:"storage_power_w,omitempty"` + StorageEnergyWh map[string]float64 `json:"storage_energy_wh,omitempty"` } // DiagnosticParams is a JSON-friendly subset of the Params struct — @@ -66,9 +67,9 @@ type DiagnosticSlot struct { // without pulling the whole internal struct. type DiagnosticParams struct { Mode Mode `json:"mode"` - InitialSoC float64 `json:"initial_soc"` - SoCMin float64 `json:"soc_min"` - SoCMax float64 `json:"soc_max"` + InitialSoC float64 `json:"initial_soc"` + SoCMin float64 `json:"soc_min"` + SoCMax float64 `json:"soc_max"` PVChargeBonusOreKwh float64 `json:"pv_charge_bonus_ore_kwh,omitempty"` SoCLevels int `json:"soc_levels"` ActionLevels int `json:"action_levels"` @@ -170,15 +171,16 @@ func buildDiagnostic(plan *Plan, slots []Slot, p Params, zone string, WeatherRowAvailableAtMs: slot.WeatherRowAvailableAtMs, BatteryW: action.BatteryW, GridW: action.GridW, - SoC: action.SoC, + SoC: action.SoC, CostOre: action.CostOre, Reason: action.Reason, EMSMode: action.EMSMode, PVLimitW: action.PVLimitW, + PVCurtailActive: action.PVCurtailActive, LoadpointW: action.LoadpointW, - LoadpointSoC: action.LoadpointSoC, + LoadpointSoC: action.LoadpointSoC, LoadpointPowerW: action.LoadpointPowerW, - LoadpointSoCByID: action.LoadpointSoCByID, + LoadpointSoCByID: action.LoadpointSoCByID, StoragePowerW: action.StoragePowerW, StorageEnergyWh: action.StorageEnergyWh, } @@ -202,9 +204,9 @@ func buildDiagnostic(plan *Plan, slots []Slot, p Params, zone string, OptimizerInput: append(json.RawMessage(nil), plan.OptimizerInput...), Params: DiagnosticParams{ Mode: p.Mode, - InitialSoC: p.InitialSoC, - SoCMin: p.SoCMin, - SoCMax: p.SoCMax, + InitialSoC: p.InitialSoC, + SoCMin: p.SoCMin, + SoCMax: p.SoCMax, PVChargeBonusOreKwh: p.PVChargeBonusOreKwh, SoCLevels: p.SoCLevels, ActionLevels: p.ActionLevels, @@ -324,9 +326,9 @@ func planFromDiagnostic(d *Diagnostic) (*Plan, []Slot, Params, time.Time, bool) } params := Params{ Mode: d.Params.Mode, - InitialSoC: d.Params.InitialSoC, - SoCMin: d.Params.SoCMin, - SoCMax: d.Params.SoCMax, + InitialSoC: d.Params.InitialSoC, + SoCMin: d.Params.SoCMin, + SoCMax: d.Params.SoCMax, PVChargeBonusOreKwh: d.Params.PVChargeBonusOreKwh, SoCLevels: d.Params.SoCLevels, ActionLevels: d.Params.ActionLevels, @@ -377,26 +379,27 @@ func planFromDiagnostic(d *Diagnostic) (*Plan, []Slot, Params, time.Time, bool) WeatherRowAvailableAtMs: ds.WeatherRowAvailableAtMs, }) action := Action{ - SlotStartMs: ds.SlotStartMs, - SlotLenMin: lenMin, - PriceOre: ds.PriceOre, - SpotOre: ds.SpotOre, - PVW: ds.PVW, - LoadW: ds.LoadW, - BatteryW: ds.BatteryW, - GridW: ds.GridW, + SlotStartMs: ds.SlotStartMs, + SlotLenMin: lenMin, + PriceOre: ds.PriceOre, + SpotOre: ds.SpotOre, + PVW: ds.PVW, + LoadW: ds.LoadW, + BatteryW: ds.BatteryW, + GridW: ds.GridW, SoC: ds.SoC, - CostOre: ds.CostOre, - Confidence: ds.Confidence, - Reason: ds.Reason, - EMSMode: ds.EMSMode, - PVLimitW: ds.PVLimitW, - LoadpointW: ds.LoadpointW, + CostOre: ds.CostOre, + Confidence: ds.Confidence, + Reason: ds.Reason, + EMSMode: ds.EMSMode, + PVLimitW: ds.PVLimitW, + PVCurtailActive: ds.PVCurtailActive, + LoadpointW: ds.LoadpointW, LoadpointSoC: ds.LoadpointSoC, - LoadpointPowerW: ds.LoadpointPowerW, + LoadpointPowerW: ds.LoadpointPowerW, LoadpointSoCByID: ds.LoadpointSoCByID, - StoragePowerW: ds.StoragePowerW, - StorageEnergyWh: ds.StorageEnergyWh, + StoragePowerW: ds.StoragePowerW, + StorageEnergyWh: ds.StorageEnergyWh, } if identified && ds.SlotEndMs > 0 { slotEndMs, err := checkedSlotEndMs(action.SlotStartMs, action.SlotLenMin) @@ -424,7 +427,7 @@ func planFromDiagnostic(d *Diagnostic) (*Plan, []Slot, Params, time.Time, bool) Mode: params.Mode, HorizonSlots: horizon, CapacityWh: params.CapacityWh, - InitialSoC: params.InitialSoC, + InitialSoC: params.InitialSoC, TotalCostOre: d.TotalCostOre, Actions: actions, Solver: d.Solver, diff --git a/go/internal/mpc/external_optimizer.go b/go/internal/mpc/external_optimizer.go index 502aaf237..21e1d9fa4 100644 --- a/go/internal/mpc/external_optimizer.go +++ b/go/internal/mpc/external_optimizer.go @@ -302,19 +302,20 @@ type externalPlan struct { } type externalAction struct { - SlotStartMs int64 `json:"slot_start_ms"` - SlotLenMin int `json:"slot_len_min"` - BatteryW float64 `json:"battery_w"` - GridW float64 `json:"grid_w"` - SoCPct float64 `json:"soc_pct"` - CostOre float64 `json:"cost_ore"` - PVLimitW float64 `json:"pv_limit_w"` - StoragePowerW map[string]float64 `json:"storage_power_w"` - StorageEnergy map[string]float64 `json:"storage_energy_wh"` - FlexPowerW map[string]float64 `json:"flex_power_w"` - FlexEnergyWh map[string]float64 `json:"flex_energy_wh"` - ThermalPowerW map[string]float64 `json:"thermal_power_w"` - ThermalState map[string]float64 `json:"thermal_state"` + SlotStartMs int64 `json:"slot_start_ms"` + SlotLenMin int `json:"slot_len_min"` + BatteryW float64 `json:"battery_w"` + GridW float64 `json:"grid_w"` + SoCPct float64 `json:"soc_pct"` + CostOre float64 `json:"cost_ore"` + PVLimitW float64 `json:"pv_limit_w"` + PVCurtailActive bool `json:"pv_curtail_active,omitempty"` + StoragePowerW map[string]float64 `json:"storage_power_w"` + StorageEnergy map[string]float64 `json:"storage_energy_wh"` + FlexPowerW map[string]float64 `json:"flex_power_w"` + FlexEnergyWh map[string]float64 `json:"flex_energy_wh"` + ThermalPowerW map[string]float64 `json:"thermal_power_w"` + ThermalState map[string]float64 `json:"thermal_state"` } func (o *ExternalOptimizer) Optimize(ctx context.Context, slots []Slot, p Params) (Plan, error) { @@ -527,6 +528,7 @@ func (r externalResponse) toPlan(slots []Slot, p Params) Plan { BatteryW: candidate.BatteryW, GridW: candidate.GridW, SoC: candidate.SoCPct / 100, CostOre: candidate.CostOre, PVLimitW: candidate.PVLimitW, + PVCurtailActive: candidate.PVCurtailActive, StoragePowerW: candidate.StoragePowerW, StorageEnergyWh: candidate.StorageEnergy, } @@ -605,7 +607,7 @@ func ValidatePlan(slots []Slot, p Params, plan *Plan) error { return fmt.Errorf("slot %d battery_w %.3f exceeds bounds", i, a.BatteryW) } dtH := float64(slot.LenMin) / 60 - if len(p.Storages) > 0 { + if len(p.Storages) > 0 && len(a.StoragePowerW) > 0 { var totalPowerW, totalEnergyWh float64 for _, storage := range p.Storages { powerW, powerOK := a.StoragePowerW[storage.ID] @@ -619,12 +621,8 @@ func ValidatePlan(slots []Slot, p Params, plan *Plan) error { if powerW > storage.MaxChargeW+2 || powerW < -storage.MaxDischargeW-2 { return fmt.Errorf("slot %d storage %s power %.3f exceeds bounds", i, storage.ID, powerW) } - energyWh := storageEnergy[storage.ID] - if powerW >= 0 { - energyWh += powerW * dtH * storage.ChargeEfficiency - } else { - energyWh += powerW * dtH / storage.DischargeEfficiency - } + energyWh := storageEnergy[storage.ID] + loadpoint.BatteryEnergyDeltaWh( + powerW, dtH, storage.ChargeEfficiency, storage.DischargeEfficiency) energyToleranceWh := math.Max(1, storage.CapacityWh*0.0002) if energyWh < -energyToleranceWh || energyWh > storage.CapacityWh+energyToleranceWh || math.Abs(reportedEnergyWh-energyWh) > energyToleranceWh { return fmt.Errorf("slot %d storage %s energy %.3f inconsistent with replay %.3f", i, storage.ID, reportedEnergyWh, energyWh) @@ -644,10 +642,14 @@ func ValidatePlan(slots []Slot, p Params, plan *Plan) error { return fmt.Errorf("slot %d aggregate battery_w %.3f, want %.3f", i, a.BatteryW, totalPowerW) } soc = totalEnergyWh / p.CapacityWh - } else if a.BatteryW >= 0 { - soc += a.BatteryW * dtH * p.ChargeEfficiency / p.CapacityWh } else { - soc += a.BatteryW * dtH / p.DischargeEfficiency / p.CapacityWh + // Go DP publishes an aggregate trajectory. Replay that fleet + // as one battery; per-storage maps are required only when present. + if p.CapacityWh <= 0 { + return fmt.Errorf("slot %d capacity_wh must be positive to replay aggregate energy", i) + } + soc += loadpoint.BatteryEnergyDeltaWh( + a.BatteryW, dtH, p.ChargeEfficiency, p.DischargeEfficiency) / p.CapacityWh } lowerRecovery := math.Max(0, p.SoCMin-soc) upperRecovery := math.Max(0, soc-p.SoCMax) @@ -680,12 +682,33 @@ func ValidatePlan(slots []Slot, p Params, plan *Plan) error { } totalLoadpointW += powerW } + if a.PVLimitW < 0 || (a.PVLimitW > 0 && a.PVLimitW > -slot.PVW+2) { + return fmt.Errorf("slot %d pv_limit_w %.3f exceeds forecast generation %.3f", i, a.PVLimitW, -slot.PVW) + } effectivePVW := slot.PVW - if a.PVLimitW > 0 { - if a.PVLimitW > -slot.PVW+2 { - return fmt.Errorf("slot %d pv_limit_w %.3f exceeds forecast generation %.3f", i, a.PVLimitW, -slot.PVW) + if a.PVCurtailActive { + // Applied cap, including a true zero. GridW must already + // include it; this is the optimizer encoding. + effectivePVW = loadpoint.EffectivePVW(slot.PVW, a.PVLimitW, true) + wantGridW := loadpoint.GridW(slot.LoadW, effectivePVW, a.BatteryW, totalLoadpointW) + if math.Abs(a.GridW-wantGridW) > 2 { + return fmt.Errorf("slot %d grid balance %.3f, want %.3f", i, a.GridW, wantGridW) + } + } else { + // Go DP writes a positive PVLimitW as a dispatch hint and + // leaves GridW on the uncurtailed identity. An optimizer + // that applied a positive cap (legacy, no flag) matches + // the curtailed identity instead. + uncurtailedGridW := loadpoint.GridW(slot.LoadW, slot.PVW, a.BatteryW, totalLoadpointW) + if math.Abs(a.GridW-uncurtailedGridW) > 2 { + if a.PVLimitW > 0 { + effectivePVW = loadpoint.EffectivePVW(slot.PVW, a.PVLimitW, true) + } + wantGridW := loadpoint.GridW(slot.LoadW, effectivePVW, a.BatteryW, totalLoadpointW) + if math.Abs(a.GridW-wantGridW) > 2 { + return fmt.Errorf("slot %d grid balance %.3f, want %.3f", i, a.GridW, wantGridW) + } } - effectivePVW = -a.PVLimitW } for lpIdx, lp := range activeLoadpoints { powerW := a.LoadpointPowerW[lp.ID] @@ -702,11 +725,7 @@ func ValidatePlan(slots []Slot, p Params, plan *Plan) error { return fmt.Errorf("slot %d battery discharge feeds loadpoint %s", i, lp.ID) } } - wantGridW := loadpoint.GridW(slot.LoadW, effectivePVW, a.BatteryW, totalLoadpointW) - if math.Abs(a.GridW-wantGridW) > 2 { - return fmt.Errorf("slot %d grid balance %.3f, want %.3f", i, a.GridW, wantGridW) - } - baseGridW := slot.LoadW + effectivePVW + totalLoadpointW + baseGridW := loadpoint.GridW(slot.LoadW, effectivePVW, 0, totalLoadpointW) if !modeAllows(p.Mode, baseGridW, a.GridW, a.BatteryW) { return fmt.Errorf("slot %d violates mode %s: baseline_grid_w=%.9f grid_w=%.9f battery_w=%.9f", i, p.Mode, baseGridW, a.GridW, a.BatteryW) diff --git a/go/internal/mpc/mpc.go b/go/internal/mpc/mpc.go index 4117edec2..313ab5bf2 100644 --- a/go/internal/mpc/mpc.go +++ b/go/internal/mpc/mpc.go @@ -293,13 +293,15 @@ type Action struct { EMSMode string `json:"ems_mode"` // effective EMS mode for this slot (set by SlotAt post-processing) // PVLimitW is the recommended cap on PV inverter output (W, positive). - // 0 = no curtailment. Set by post-processing when exporting would - // cost money (negative export revenue after fees). Includes house - // load + battery charge + any planned EV loadpoint charge so that - // curtailment does not starve loads the plan itself scheduled. - // Consumed by the control loop only when the driver advertises - // `supports_pv_curtail`. - PVLimitW float64 `json:"pv_limit_w,omitempty"` + // When PVCurtailActive is false, 0 means no cap (a dispatch hint may + // still use a positive PVLimitW without rewriting GridW). When + // PVCurtailActive is true, 0 is a real zero cap already applied to + // GridW. Includes house load + battery charge + any planned EV + // loadpoint charge so that curtailment does not starve loads the + // plan itself scheduled. Consumed by the control loop only when + // the driver advertises `supports_pv_curtail`. + PVLimitW float64 `json:"pv_limit_w,omitempty"` + PVCurtailActive bool `json:"pv_curtail_active,omitempty"` // LoadpointW is the EV charger power (W, positive = charging) the // DP picked for this slot. Zero when no loadpoint was in Params @@ -469,6 +471,26 @@ func finite(v float64) bool { return !math.IsNaN(v) && !math.IsInf(v, 0) } +func gridIndex(value, min, step float64, n int) int { + if n <= 1 || step <= 0 { + return 0 + } + i := int(math.Round((value - min) / step)) + if i < 0 { + return 0 + } + if i >= n { + return n - 1 + } + return i +} + +func operatingBoundWorsens(from, to, min, max float64) bool { + const eps = 1e-9 + return math.Max(0, min-to) > math.Max(0, min-from)+eps || + math.Max(0, to-max) > math.Max(0, from-max)+eps +} + func sanitizeOptimizeSlots(slots []Slot) []Slot { out := make([]Slot, 0, len(slots)) for _, s := range slots { @@ -685,12 +707,7 @@ func Optimize(slots []Slot, p Params) Plan { battW := actionAt(ba) // Battery SoC transition (independent of EV). - var dBattWh float64 - if battW >= 0 { - dBattWh = +battW * dtH * p.ChargeEfficiency - } else { - dBattWh = +battW * dtH / p.DischargeEfficiency - } + dBattWh := loadpoint.BatteryEnergyDeltaWh(battW, dtH, p.ChargeEfficiency, p.DischargeEfficiency) battSoc2 := soc + dBattWh/p.CapacityWh if battSoc2 < p.SoCMin-1e-9 || battSoc2 > p.SoCMax+1e-9 { continue @@ -971,28 +988,16 @@ func Optimize(slots []Slot, p Params) Plan { InitialSoC: p.InitialSoC, Actions: make([]Action, 0, N), } - fIdx := (p.InitialSoC - p.SoCMin) / socStep - si := int(math.Round(fIdx)) - if si < 0 { - si = 0 - } - if si >= S { - si = S - 1 - } - soc := socAt(si) - // Initial EV SoC index. + // Policy is stored on the SoC grid; energy is not. Integrate from + // the actual initial SoC so reported trajectories replay. Clamp + // only the policy lookup index onto the operating grid. + soc := p.InitialSoC + si := gridIndex(soc, p.SoCMin, socStep, S) ei := 0 var evSoc float64 if evActive { - f := (lp.InitialSoC - lp.SoCMin) / evSocStep - ei = int(math.Round(f)) - if ei < 0 { - ei = 0 - } - if ei >= EL { - ei = EL - 1 - } - evSoc = evSocAt(ei) + evSoc = lp.InitialSoC + ei = gridIndex(evSoc, lp.SoCMin, evSocStep, EL) } var totalCost float64 for t := 0; t < N; t++ { @@ -1003,27 +1008,19 @@ func Optimize(slots []Slot, p Params) Plan { ea := pol % EA actW := actionAt(ba) evW := evActionW(ea) - // Battery SoC transition. - var dSoCWh float64 - if actW >= 0 { - dSoCWh = +actW * dtH * p.ChargeEfficiency - } else { - dSoCWh = +actW * dtH / p.DischargeEfficiency - } + dSoCWh := loadpoint.BatteryEnergyDeltaWh(actW, dtH, p.ChargeEfficiency, p.DischargeEfficiency) soc2 := soc + dSoCWh/p.CapacityWh - if soc2 < p.SoCMin { - soc2 = p.SoCMin - } - if soc2 > p.SoCMax { - soc2 = p.SoCMax + if operatingBoundWorsens(soc, soc2, p.SoCMin, p.SoCMax) { + actW = 0 + soc2 = soc } - // EV SoC transition (no-op when !evActive since evW = 0). var evSoc2 float64 if evActive { dEvWh := evW * dtH * evChargeEff evSoc2 = evSoc + dEvWh/lp.CapacityWh - if evSoc2 > lp.SoCMax { - evSoc2 = lp.SoCMax + if evSoc2 > lp.SoCMax+1e-9 { + evW = 0 + evSoc2 = evSoc } } gridW := loadpoint.GridW(slot.LoadW, slot.PVW, actW, evW) @@ -1053,24 +1050,10 @@ func Optimize(slots []Slot, p Params) Plan { } plan.Actions = append(plan.Actions, a) soc = soc2 - fIdx = (soc - p.SoCMin) / socStep - si = int(math.Round(fIdx)) - if si < 0 { - si = 0 - } - if si >= S { - si = S - 1 - } + si = gridIndex(soc, p.SoCMin, socStep, S) if evActive { evSoc = evSoc2 - f := (evSoc - lp.SoCMin) / evSocStep - ei = int(math.Round(f)) - if ei < 0 { - ei = 0 - } - if ei >= EL { - ei = EL - 1 - } + ei = gridIndex(evSoc, lp.SoCMin, evSocStep, EL) } } plan.TotalCostOre = totalCost diff --git a/go/internal/mpc/service.go b/go/internal/mpc/service.go index 2e3e1e2cd..820a6e865 100644 --- a/go/internal/mpc/service.go +++ b/go/internal/mpc/service.go @@ -1473,6 +1473,19 @@ func (s *Service) runReplan(request replanRequest) *Plan { "err", err) return s.Latest() } + if err := ValidatePlan(slots, p, &plan); err != nil { + engine := "go-dp" + if plan.Solver != nil && plan.Solver.Engine != "" { + engine = plan.Solver.Engine + } + slog.Error("mpc: rejected plan that failed physical replay", + "generation", request.generation, + "mode", p.Mode, + "reason", request.reason, + "engine", engine, + "err", err) + return s.Latest() + } // Tag each action with the effective EMS mode so the UI can render // a mode-band showing which strategy drives each slot. diff --git a/go/internal/mpc/validate_dp_test.go b/go/internal/mpc/validate_dp_test.go new file mode 100644 index 000000000..083e3f403 --- /dev/null +++ b/go/internal/mpc/validate_dp_test.go @@ -0,0 +1,175 @@ +package mpc + +import ( + "math" + "testing" +) + +func TestOptimizePlansPassValidatePlan(t *testing.T) { + cases := []struct { + name string + slots []Slot + p Params + }{ + { + name: "self_consumption flat load", + slots: flatLoadSlots([]float64{100, 200, 50, 300}), + p: func() Params { + p := baseParams(ModeSelfConsumption) + p.InitialSoC = 0.80 + return p + }(), + }, + { + name: "self_consumption pv surplus", + slots: []Slot{ + {StartMs: 0, LenMin: 60, PriceOre: 100, Confidence: 1, LoadW: 2000, PVW: -3500}, + }, + p: baseParams(ModeSelfConsumption), + }, + { + name: "passive_arbitrage", + slots: flatLoadSlots([]float64{40, 200, 40, 250}), + p: baseParams(ModePassiveArbitrage), + }, + { + name: "arbitrage", + slots: flatLoadSlots([]float64{20, 250, 20, 300}), + p: baseParams(ModeArbitrage), + }, + } + for _, tc := range cases { + t.Run(tc.name, func(t *testing.T) { + plan := Optimize(tc.slots, tc.p) + if len(plan.Actions) == 0 { + t.Fatal("Optimize returned no actions") + } + if err := ValidatePlan(tc.slots, tc.p, &plan); err != nil { + t.Fatalf("ValidatePlan: %v", err) + } + }) + } +} + +func TestValidatePlanAcceptsGoDPCurtailHint(t *testing.T) { + slots := []Slot{{ + StartMs: 1, LenMin: 60, PriceOre: 100, SpotOre: -50, Confidence: 1, + LoadW: 500, PVW: -5000, + }} + p := baseParams(ModeSelfConsumption) + p.InitialSoC = 0.90 + plan := Optimize(slots, p) + if len(plan.Actions) != 1 { + t.Fatalf("got %d actions", len(plan.Actions)) + } + if plan.Actions[0].PVLimitW <= 0 { + t.Fatalf("expected a curtail hint, got pv_limit_w=%f grid_w=%f", plan.Actions[0].PVLimitW, plan.Actions[0].GridW) + } + uncurtailed := plan.Actions[0].LoadW + plan.Actions[0].PVW + plan.Actions[0].BatteryW + if math.Abs(plan.Actions[0].GridW-uncurtailed) > 2 { + t.Fatalf("Go DP GridW = %f, want uncurtailed %f", plan.Actions[0].GridW, uncurtailed) + } + if err := ValidatePlan(slots, p, &plan); err != nil { + t.Fatalf("ValidatePlan rejected DP curtail hint: %v", err) + } +} + +func TestValidatePlanRejectsFuseViolatingIdle(t *testing.T) { + slots := []Slot{{ + StartMs: 1, LenMin: 60, PriceOre: 100, Confidence: 1, + LoadW: 0, PVW: -8000, + Limits: PowerLimits{MaxExportW: 100}, + }} + p := baseParams(ModeSelfConsumption) + p.MaxChargeW = 0 + p.MaxDischargeW = 0 + p.InitialSoC = 0.90 + plan := Optimize(slots, p) + if len(plan.Actions) != 1 { + t.Fatalf("got %d actions", len(plan.Actions)) + } + err := ValidatePlan(slots, p, &plan) + if err == nil { + t.Fatal("ValidatePlan accepted idle export past MaxExportW") + } +} + +func TestValidatePlanAcceptsActiveZeroPVCap(t *testing.T) { + slots := []Slot{{ + StartMs: 1, LenMin: 60, PriceOre: 100, SpotOre: -100, Confidence: 1, + LoadW: 0, PVW: -5000, + }} + p := baseParams(ModeArbitrage) + p.InitialSoC = 0.95 + plan := Plan{ + Mode: p.Mode, HorizonSlots: 1, CapacityWh: p.CapacityWh, InitialSoC: 0.95, + Actions: []Action{{ + SlotStartMs: 1, SlotLenMin: 60, + BatteryW: 0, GridW: 0, SoC: 0.95, CostOre: 0, + PVLimitW: 0, PVCurtailActive: true, + }}, + } + if err := ValidatePlan(slots, p, &plan); err != nil { + t.Fatalf("active zero cap: %v", err) + } + + plan.Actions[0].PVCurtailActive = false + if err := ValidatePlan(slots, p, &plan); err == nil { + t.Fatal("zero grid without pv_curtail_active must not replay as uncurtailed PV") + } +} + +func TestValidatePlanReplaysAggregateWhenStorageMapsEmpty(t *testing.T) { + slots := flatLoadSlots([]float64{100, 200}) + p := baseParams(ModeSelfConsumption) + p.InitialSoC = 0.80 + p.Storages = []StorageAssetSpec{{ + ID: "home", CapacityWh: p.CapacityWh, + InitialEnergyWh: p.CapacityWh * p.InitialSoC, + MinEnergyWh: p.CapacityWh * p.SoCMin, + MaxEnergyWh: p.CapacityWh * p.SoCMax, + MaxChargeW: p.MaxChargeW, MaxDischargeW: p.MaxDischargeW, + ChargeEfficiency: p.ChargeEfficiency, DischargeEfficiency: p.DischargeEfficiency, + }} + plan := Optimize(slots, p) + if len(plan.Actions[0].StoragePowerW) != 0 { + t.Fatal("Go DP should not invent per-storage maps") + } + if err := ValidatePlan(slots, p, &plan); err != nil { + t.Fatalf("aggregate replay: %v", err) + } +} + +func TestOptimizeReplaysFromActualSoCBelowMinimum(t *testing.T) { + slots := flatLoadSlots([]float64{100, 100}) + p := baseParams(ModeSelfConsumption) + p.InitialSoC = 0.08 // below SoCMin 0.10 + plan := Optimize(slots, p) + if len(plan.Actions) == 0 { + t.Fatal("Optimize returned no actions") + } + if err := ValidatePlan(slots, p, &plan); err != nil { + t.Fatalf("out-of-band start must replay: %v", err) + } + if plan.Actions[0].SoC < p.InitialSoC-1e-9 { + t.Fatalf("first SoC %.4f worsened below start %.4f", plan.Actions[0].SoC, p.InitialSoC) + } +} + +func TestOptimizeDoesNotWorsenOperatingBoundNearFloor(t *testing.T) { + slots := []Slot{{ + StartMs: 0, LenMin: 60, PriceOre: 300, Confidence: 1, LoadW: 2000, PVW: 0, + }} + p := baseParams(ModeSelfConsumption) + p.InitialSoC = p.SoCMin + 0.001 + plan := Optimize(slots, p) + if len(plan.Actions) != 1 { + t.Fatalf("got %d actions", len(plan.Actions)) + } + if err := ValidatePlan(slots, p, &plan); err != nil { + t.Fatalf("near-floor plan must replay: %v", err) + } + if plan.Actions[0].SoC < p.SoCMin-1e-9 && plan.Actions[0].SoC < p.InitialSoC-1e-9 { + t.Fatalf("SoC %.4f left the band and worsened start %.4f", plan.Actions[0].SoC, p.InitialSoC) + } +} diff --git a/optimizer/ftw_optimizer/direct_highs.py b/optimizer/ftw_optimizer/direct_highs.py index f7400cd07..90dd3ae0d 100644 --- a/optimizer/ftw_optimizer/direct_highs.py +++ b/optimizer/ftw_optimizer/direct_highs.py @@ -12,6 +12,7 @@ from .deadline import SolveCancelled, SolveDeadline, SolveDeadlineExceeded from .model import ( _arbitrage_spread_ore_kwh, + _pv_curtail_output, _solver_options, _storage_starts_above_maximum, ) @@ -724,6 +725,7 @@ def _response( raw_total_cost += raw_cost curtailed_w = max(0.0, float(solution[base_vars.curtail[t]])) pv_forecast = prepared.base_pv if shared else base.pv + pv_limit_w, pv_curtail_active = _pv_curtail_output(pv_forecast[t], curtailed_w) actions.append( { "slot_start_ms": int(slot.get("start_ms", 0)), @@ -732,9 +734,8 @@ def _response( "grid_w": grid_w, "soc_pct": stored_wh / total_capacity * 100.0, "cost_ore": raw_cost, - "pv_limit_w": max(0.0, -pv_forecast[t] - curtailed_w) - if curtailed_w > 1e-5 - else 0.0, + "pv_limit_w": pv_limit_w, + "pv_curtail_active": pv_curtail_active, "storage_power_w": storage_power, "storage_energy_wh": storage_energy, "flex_power_w": {}, diff --git a/optimizer/ftw_optimizer/model.py b/optimizer/ftw_optimizer/model.py index f566ce8b6..0250638f2 100644 --- a/optimizer/ftw_optimizer/model.py +++ b/optimizer/ftw_optimizer/model.py @@ -232,6 +232,25 @@ def _arbitrage_spread_ore_kwh(settings: dict[str, Any], mode: str) -> float: return spread +def _pv_charge_bonus_ore_kwh(settings: dict[str, Any], mode: str) -> float: + """Return the PV-charge bonus only for passive_arbitrage. + + Parse in every mode so a malformed value still fails at the contract + boundary. Go DP applies this bias only in passive_arbitrage. + """ + + bonus = max( + 0.0, + finite_number( + settings.get("pv_charge_bonus_ore_kwh", 0), + "settings.pv_charge_bonus_ore_kwh", + ), + ) + if mode != "passive_arbitrage": + return 0.0 + return bonus + + def _requires_direction_binary(formulation: str, relaxation_unsafe: bool) -> bool: """Keep mutually exclusive physical flows when a relaxation can profit.""" @@ -264,6 +283,19 @@ def _storage_relaxation_is_unsafe( ) +def _pv_curtail_output(forecast_pv_w: float, curtailed_w: float) -> tuple[float, bool]: + """Return (pv_limit_w, pv_curtail_active) for one slot. + + Active with pv_limit_w = 0 is a true zero cap. Inactive with 0 is + release. The two must not share a sentinel. + """ + + curtailed_w = max(0.0, float(curtailed_w)) + if curtailed_w <= 1e-5: + return 0.0, False + return max(0.0, -float(forecast_pv_w) - curtailed_w), True + + def _export_price(slot: dict[str, Any], settings: dict[str, Any]) -> float: flat = finite_number(settings.get("export_ore_per_kwh", 0), "settings.export_ore_per_kwh") if flat > 0: @@ -509,13 +541,7 @@ def solve( formulation = settings.get("formulation", "auto") if formulation not in {"auto", "milp", "relaxed"}: raise ProtocolError("settings.formulation must be auto, milp, or relaxed") - pv_charge_bonus_ore = max( - 0.0, - finite_number( - settings.get("pv_charge_bonus_ore_kwh", 0), - "settings.pv_charge_bonus_ore_kwh", - ), - ) + pv_charge_bonus_ore = _pv_charge_bonus_ore_kwh(settings, mode) constraints: list[cp.Constraint] = [] discrete = False @@ -927,8 +953,9 @@ def solve( # that exceeds the slot's actual leftover. constraints.append(flex.power <= house_surplus + 50.0) if bool(flex.spec.get("no_storage_to_load", False)) and storages: - house_residual = np.maximum(0.0, base_load + base_pv) - constraints.append(total_discharge <= house_residual + max_site_power * (1 - active)) + for scenario in scenarios: + house_residual = np.maximum(0.0, scenario["load"] + scenario["pv"]) + constraints.append(total_discharge <= house_residual + max_site_power * (1 - active)) if storage_discharge_active is not None: # EV charging may coexist with house-covering discharge, but not # with battery-driven site export. @@ -1065,6 +1092,7 @@ def run_problem(problem: cp.Problem, solver_name: str) -> None: raw_cost = price[t] * max(grid_kwh, 0.0) - export_price[t] * max(-grid_kwh, 0.0) raw_total_cost += raw_cost curtailed_w = max(0.0, float(curtail.value[t])) + pv_limit_w, pv_curtail_active = _pv_curtail_output(base_pv[t], curtailed_w) actions.append( { "slot_start_ms": int(slot.get("start_ms", 0)), @@ -1073,7 +1101,8 @@ def run_problem(problem: cp.Problem, solver_name: str) -> None: "grid_w": grid_w, "soc_pct": (stored_wh / total_capacity * 100.0) if total_capacity > 0 else 0.0, "cost_ore": raw_cost, - "pv_limit_w": max(0.0, -base_pv[t] - curtailed_w) if curtailed_w > 1e-5 else 0.0, + "pv_limit_w": pv_limit_w, + "pv_curtail_active": pv_curtail_active, "storage_power_w": storage_power, "storage_energy_wh": storage_energy, "flex_power_w": flex_power, diff --git a/optimizer/ftw_optimizer/multistage.py b/optimizer/ftw_optimizer/multistage.py index 496601da8..aba7b8479 100644 --- a/optimizer/ftw_optimizer/multistage.py +++ b/optimizer/ftw_optimizer/multistage.py @@ -17,6 +17,8 @@ ReplayConsistencyError, _STORAGE_INITIAL_ABOVE_MAXIMUM_KEY, _arbitrage_spread_ore_kwh, + _pv_charge_bonus_ore_kwh, + _pv_curtail_output, _canonicalize_storage_payload, _export_price, _mode, @@ -137,13 +139,7 @@ def assign(self, prepared: PreparedMultistage) -> None: self.import_coeff.value = prepared.effective_import * prepared.dt_h / 1000.0 self.export_coeff.value = prepared.effective_export * prepared.dt_h / 1000.0 self.strict_coeff.value = 2.0 * np.maximum(prepared.effective_import, 0.0) * prepared.dt_h / 1000.0 - self.pv_bonus.value = max( - 0.0, - finite_number( - prepared.settings.get("pv_charge_bonus_ore_kwh", 0), - "settings.pv_charge_bonus_ore_kwh", - ), - ) + self.pv_bonus.value = _pv_charge_bonus_ore_kwh(prepared.settings, prepared.mode) spread = _arbitrage_spread_ore_kwh(prepared.settings, prepared.mode) for i, spec in enumerate(prepared.storages): initial = finite_number(spec.get("initial_energy_wh"), f"storages[{i}].initial_energy_wh") @@ -420,13 +416,7 @@ def _prepare(payload: dict[str, Any]) -> PreparedMultistage: formulation = str(settings.get("formulation", "auto")) if formulation not in {"auto", "milp", "relaxed"}: raise ProtocolError("settings.formulation must be auto, milp, or relaxed") - pv_charge_bonus = max( - 0.0, - finite_number( - settings.get("pv_charge_bonus_ore_kwh", 0), - "settings.pv_charge_bonus_ore_kwh", - ), - ) + pv_charge_bonus = _pv_charge_bonus_ore_kwh(settings, mode) unsafe_meter_split = bool(np.any(effective_import < effective_export - 1e-9)) base_load = np.asarray( [finite_number(slot.get("load_w", 0), f"slots[{i}].load_w") for i, slot in enumerate(slots)] @@ -986,6 +976,7 @@ def _response( raw_cost = prepared.price[t] * max(grid_kwh, 0.0) - prepared.export_price[t] * max(-grid_kwh, 0.0) raw_total_cost += raw_cost curtailed_w = max(0.0, float(curtail_values[t])) + pv_limit_w, pv_curtail_active = _pv_curtail_output(base.pv[t], curtailed_w) actions.append( { "slot_start_ms": int(slot.get("start_ms", 0)), @@ -994,7 +985,8 @@ def _response( "grid_w": grid_w, "soc_pct": stored_wh / total_capacity * 100.0, "cost_ore": raw_cost, - "pv_limit_w": max(0.0, -base.pv[t] - curtailed_w) if curtailed_w > 1e-5 else 0.0, + "pv_limit_w": pv_limit_w, + "pv_curtail_active": pv_curtail_active, "storage_power_w": storage_power, "storage_energy_wh": storage_energy, "flex_power_w": {}, diff --git a/optimizer/ftw_optimizer/progressive.py b/optimizer/ftw_optimizer/progressive.py index 2233fc471..c6681a9f2 100644 --- a/optimizer/ftw_optimizer/progressive.py +++ b/optimizer/ftw_optimizer/progressive.py @@ -10,7 +10,13 @@ from . import SCHEMA_VERSION from .deadline import SolveDeadline -from .model import OPTIMAL_STATUSES, _arbitrage_spread_ore_kwh, _solver_options +from .model import ( + OPTIMAL_STATUSES, + _arbitrage_spread_ore_kwh, + _pv_charge_bonus_ore_kwh, + _pv_curtail_output, + _solver_options, +) from .protocol import ProtocolError, finite_number if TYPE_CHECKING: @@ -53,7 +59,7 @@ def ph_eligible(prepared: "PreparedMultistage") -> tuple[bool, str]: return False, "mode is not unconstrained arbitrage" if prepared.economic_cvar_weight > 0: return False, "economic CVaR couples scenario subproblems" - if finite_number(settings.get("pv_charge_bonus_ore_kwh", 0), "settings.pv_charge_bonus_ore_kwh") != 0: + if _pv_charge_bonus_ore_kwh(settings, prepared.mode) != 0: return False, "PV charge bonus can incentivize simultaneous cycling" if np.any(prepared.effective_import < -1e-9): return False, "negative import prices require a discrete cycling guard" @@ -358,6 +364,7 @@ def _response( raw_cost = prepared.price[t] * max(grid_kwh, 0.0) - prepared.export_price[t] * max(-grid_kwh, 0.0) raw_total_cost += raw_cost curtailed_w = max(0.0, float(base_problem.curtail.value[t])) + pv_limit_w, pv_curtail_active = _pv_curtail_output(base.pv[t], curtailed_w) actions.append( { "slot_start_ms": int(slot.get("start_ms", 0)), @@ -366,7 +373,8 @@ def _response( "grid_w": grid_w, "soc_pct": stored_wh / total_capacity * 100.0, "cost_ore": raw_cost, - "pv_limit_w": max(0.0, -base.pv[t] - curtailed_w) if curtailed_w > 1e-5 else 0.0, + "pv_limit_w": pv_limit_w, + "pv_curtail_active": pv_curtail_active, "storage_power_w": storage_power, "storage_energy_wh": storage_energy, "flex_power_w": {}, diff --git a/optimizer/ftw_optimizer/recourse.py b/optimizer/ftw_optimizer/recourse.py index 25c30a99b..e9a89805f 100644 --- a/optimizer/ftw_optimizer/recourse.py +++ b/optimizer/ftw_optimizer/recourse.py @@ -14,6 +14,8 @@ OPTIMAL_STATUSES, ReplayConsistencyError, _arbitrage_spread_ore_kwh, + _pv_charge_bonus_ore_kwh, + _pv_curtail_output, _canonicalize_storage_payload, _export_price, _mode, @@ -152,7 +154,7 @@ def solve_storage_recourse( expected_pv_bonus: cp.Expression = cp.Constant(0.0) strict_sc_penalty: cp.Expression = cp.Constant(0.0) worst_service_slack = cp.Variable(nonneg=True, name="worst_service_slack") - bonus_ore = max(0.0, finite_number(settings.get("pv_charge_bonus_ore_kwh", 0), "settings.pv_charge_bonus_ore_kwh")) + bonus_ore = _pv_charge_bonus_ore_kwh(settings, mode) arbitrage_spread = _arbitrage_spread_ore_kwh(settings, mode) unsafe_cycle = _storage_relaxation_is_unsafe( eff_import, @@ -383,6 +385,7 @@ def run_problem(problem: cp.Problem, solver_name: str) -> None: raw_cost = price[t] * max(grid_kwh, 0.0) - export_price[t] * max(-grid_kwh, 0.0) raw_total_cost += raw_cost curtailed_w = max(0.0, float(base_vars["curtail"].value[t])) + pv_limit_w, pv_curtail_active = _pv_curtail_output(base["pv"][t], curtailed_w) actions.append( { "slot_start_ms": int(slot.get("start_ms", 0)), @@ -391,7 +394,8 @@ def run_problem(problem: cp.Problem, solver_name: str) -> None: "grid_w": grid_w, "soc_pct": (stored_wh / total_capacity * 100.0) if total_capacity > 0 else 0.0, "cost_ore": raw_cost, - "pv_limit_w": max(0.0, -base["pv"][t] - curtailed_w) if curtailed_w > 1e-5 else 0.0, + "pv_limit_w": pv_limit_w, + "pv_curtail_active": pv_curtail_active, "storage_power_w": storage_power, "storage_energy_wh": storage_energy, "flex_power_w": {}, diff --git a/optimizer/tests/test_model.py b/optimizer/tests/test_model.py index de00889db..a8edd2d4b 100644 --- a/optimizer/tests/test_model.py +++ b/optimizer/tests/test_model.py @@ -17,6 +17,8 @@ OPTIMAL_STATUSES, _arbitrage_spread_ore_kwh, _canonicalize_storage_payload, + _pv_charge_bonus_ore_kwh, + _pv_curtail_output, _requires_direction_binary, _storage_relaxation_is_unsafe, ) @@ -30,6 +32,25 @@ from ftw_optimizer.worker import handle, handshake +def test_pv_charge_bonus_matches_go_dp_mode_gate() -> None: + settings = {"pv_charge_bonus_ore_kwh": 30} + assert _pv_charge_bonus_ore_kwh(settings, "passive_arbitrage") == 30 + assert _pv_charge_bonus_ore_kwh(settings, "arbitrage") == 0 + assert _pv_charge_bonus_ore_kwh(settings, "self_consumption") == 0 + assert _pv_charge_bonus_ore_kwh(settings, "cheap_charge") == 0 + + +def test_pv_curtail_output_distinguishes_zero_cap_from_release() -> None: + limit, active = _pv_curtail_output(-5000, 0) + assert (limit, active) == (0.0, False) + limit, active = _pv_curtail_output(-5000, 5000) + assert active is True + assert limit == 0.0 + limit, active = _pv_curtail_output(-5000, 2000) + assert active is True + assert limit == 3000.0 + + def test_cvxpy_user_limit_is_not_an_accepted_solution() -> None: assert cp.OPTIMAL in OPTIMAL_STATUSES assert cp.OPTIMAL_INACCURATE in OPTIMAL_STATUSES @@ -602,6 +623,7 @@ def test_shared_direct_highs_matches_strict_pv_surplus_and_limit() -> None: assert math.isclose(action["battery_w"], 0, abs_tol=1e-6) assert math.isclose(action["grid_w"], -100, abs_tol=1e-6) assert math.isclose(action["pv_limit_w"], 600, abs_tol=1e-6) + assert action["pv_curtail_active"] is True assert_storage_replays(request, direct) assert_storage_replays(reference_request, reference) @@ -799,7 +821,9 @@ def test_shared_auto_falls_back_at_each_direct_eligibility_boundary() -> None: cases.append(("unsafe-cycle", negative_import)) pv_charge_bonus = base_request() - pv_charge_bonus["settings"]["pv_charge_bonus_ore_kwh"] = 1 + pv_charge_bonus["settings"].update( + {"mode": "passive_arbitrage", "pv_charge_bonus_ore_kwh": 1} + ) cases.append(("pv-charge-bonus", pv_charge_bonus)) meter_split = base_request() @@ -1579,14 +1603,20 @@ def test_multistage_auto_keeps_binary_guards_for_unsafe_incentives() -> None: assert response["solver"]["formulation"] == "multistage-milp" shared_bonus = base_request() - shared_bonus["settings"]["pv_charge_bonus_ore_kwh"] = 1 + shared_bonus["settings"].update( + {"mode": "passive_arbitrage", "pv_charge_bonus_ore_kwh": 1} + ) response = handle(shared_bonus) assert response["ok"], response assert response["solver"]["formulation"] == "milp" recourse_bonus = base_request() recourse_bonus["settings"].update( - {"scenario_policy": "recourse", "pv_charge_bonus_ore_kwh": 1} + { + "mode": "passive_arbitrage", + "scenario_policy": "recourse", + "pv_charge_bonus_ore_kwh": 1, + } ) response = handle(recourse_bonus) assert response["ok"], response @@ -1594,7 +1624,11 @@ def test_multistage_auto_keeps_binary_guards_for_unsafe_incentives() -> None: pv_bonus = base_request() pv_bonus["settings"].update( - {"scenario_policy": "multistage", "pv_charge_bonus_ore_kwh": 1} + { + "mode": "passive_arbitrage", + "scenario_policy": "multistage", + "pv_charge_bonus_ore_kwh": 1, + } ) response = handle(pv_bonus) assert response["ok"], response @@ -1748,6 +1782,7 @@ def test_relaxed_formulation_guards_pv_bonus_storage_cycles( ) -> None: request = _relaxed_flow_guard_request(scenario_policy) request["slots"][0].update({"price_ore": 0, "spot_ore": 0, "pv_w": -5000}) + request["settings"]["mode"] = "passive_arbitrage" request["settings"]["pv_charge_bonus_ore_kwh"] = 100 response = handle(request) diff --git a/optimizer/tests/test_site_physics.py b/optimizer/tests/test_site_physics.py new file mode 100644 index 000000000..a324a019d --- /dev/null +++ b/optimizer/tests/test_site_physics.py @@ -0,0 +1,65 @@ +from __future__ import annotations + +import json +from pathlib import Path + + +FIXTURE = ( + Path(__file__).resolve().parents[2] + / "go" + / "internal" + / "loadpoint" + / "testdata" + / "site_physics.json" +) + + +def grid_w(load_w: float, pv_w: float, battery_w: float, ev_w: float) -> float: + return load_w + pv_w + battery_w + ev_w + + +def leftover_w(load_w: float, pv_w: float) -> float: + return max(0.0, -(load_w + pv_w)) + + +def house_residual_w(load_w: float, pv_w: float) -> float: + return max(0.0, load_w + pv_w) + + +def battery_discharge_feeds_ev( + battery_w: float, ev_w: float, load_w: float, pv_w: float +) -> bool: + if ev_w <= 0 or battery_w >= 0: + return False + return -battery_w > house_residual_w(load_w, pv_w) + 50 + + +def battery_energy_delta_wh( + power_w: float, dt_h: float, charge_eff: float, discharge_eff: float +) -> float: + if power_w >= 0: + return power_w * dt_h * charge_eff + return power_w * dt_h / discharge_eff + + +def test_site_physics_table_matches_go_kernel() -> None: + fixture = json.loads(FIXTURE.read_text()) + for row in fixture["flows"]: + assert grid_w(row["load_w"], row["pv_w"], row["battery_w"], row["ev_w"]) == row[ + "grid_w" + ], row["name"] + assert leftover_w(row["load_w"], row["pv_w"]) == row["leftover_w"], row["name"] + assert house_residual_w(row["load_w"], row["pv_w"]) == row["house_residual_w"], row[ + "name" + ] + assert ( + battery_discharge_feeds_ev( + row["battery_w"], row["ev_w"], row["load_w"], row["pv_w"] + ) + is row["feeds_ev"] + ), row["name"] + for row in fixture["energy_steps"]: + got = battery_energy_delta_wh( + row["power_w"], row["dt_h"], row["charge_eff"], row["discharge_eff"] + ) + assert abs(got - row["delta_wh"]) < 1e-9, row["name"] From f4cc89fc8f336cca2b307f21a4e03c2bbb3c99c9 Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Sun, 23 Aug 2026 16:13:34 +0200 Subject: [PATCH 2/3] fix(mpc): ignore sub-watt grid-limit residue --- .changeset/optimizer-grid-limit-residue.md | 5 +++ go/internal/mpc/external_optimizer.go | 8 ++++- go/internal/mpc/external_optimizer_test.go | 42 ++++++++++++++++++++++ 3 files changed, 54 insertions(+), 1 deletion(-) create mode 100644 .changeset/optimizer-grid-limit-residue.md diff --git a/.changeset/optimizer-grid-limit-residue.md b/.changeset/optimizer-grid-limit-residue.md new file mode 100644 index 000000000..b009af7f0 --- /dev/null +++ b/.changeset/optimizer-grid-limit-residue.md @@ -0,0 +1,5 @@ +--- +"ftw": patch +--- + +Accept sub-watt solver residue at a slot's grid limit so an optimizer plan at the configured fuse ceiling does not trigger Go planner fallback. Larger import and export violations still fail validation. diff --git a/go/internal/mpc/external_optimizer.go b/go/internal/mpc/external_optimizer.go index 21e1d9fa4..59c256641 100644 --- a/go/internal/mpc/external_optimizer.go +++ b/go/internal/mpc/external_optimizer.go @@ -564,6 +564,11 @@ func (o *ExternalOptimizer) Health(ctx context.Context) (OptimizerRuntimeInfo, e return o.transport.Health(ctx) } +// solverGridLimitToleranceW admits only sub-watt feasibility residue from the +// mathematical optimizer. The Go planner still observes the exact slot limits, +// and dispatch keeps its separate fuse guard. +const solverGridLimitToleranceW = 0.1 + // ValidatePlan independently replays a candidate plan against the canonical // site sign convention and current constraints. Solver output is untrusted at // this boundary: NaN, stale slot alignment, energy drift, illegal EV steps, or @@ -730,7 +735,8 @@ func ValidatePlan(slots []Slot, p Params, plan *Plan) error { return fmt.Errorf("slot %d violates mode %s: baseline_grid_w=%.9f grid_w=%.9f battery_w=%.9f", i, p.Mode, baseGridW, a.GridW, a.BatteryW) } - if !slot.Limits.allowsImport(a.GridW) || !slot.Limits.allowsExport(a.GridW) { + if (slot.Limits.MaxImportW > 0 && a.GridW > slot.Limits.MaxImportW+solverGridLimitToleranceW) || + (slot.Limits.MaxExportW > 0 && a.GridW < -slot.Limits.MaxExportW-solverGridLimitToleranceW) { return fmt.Errorf("slot %d grid_w %.3f violates grid limits", i, a.GridW) } gridKWh := a.GridW * dtH / 1000 diff --git a/go/internal/mpc/external_optimizer_test.go b/go/internal/mpc/external_optimizer_test.go index 2267c5832..7b009f878 100644 --- a/go/internal/mpc/external_optimizer_test.go +++ b/go/internal/mpc/external_optimizer_test.go @@ -151,6 +151,48 @@ func TestValidatePlanAcceptsSubWattSolverResidueInPassiveMode(t *testing.T) { } } +func TestValidatePlanGridLimitAllowsOnlySubWattSolverResidue(t *testing.T) { + const limitW = 11040.0 + p := Params{ + Mode: ModeArbitrage, CapacityWh: 10000, + SoCMinPct: 10, SoCMaxPct: 95, InitialSoCPct: 50, + MaxChargeW: 5000, MaxDischargeW: 5000, + ChargeEfficiency: 1, DischargeEfficiency: 1, + } + tests := []struct { + name string + gridW float64 + wantErr bool + }{ + {name: "import solver residue", gridW: limitW + 0.000001}, + {name: "export solver residue", gridW: -limitW - 0.000001}, + {name: "import real violation", gridW: limitW + 1, wantErr: true}, + {name: "export real violation", gridW: -limitW - 1, wantErr: true}, + } + for _, tc := range tests { + t.Run(tc.name, func(t *testing.T) { + slot := Slot{ + StartMs: 1, LenMin: 15, PriceOre: 100, SpotOre: 50, Confidence: 1, + Limits: PowerLimits{MaxImportW: limitW, MaxExportW: limitW}, + } + if tc.gridW > 0 { + slot.LoadW = tc.gridW + } else { + slot.PVW = tc.gridW + } + costOre := SlotGridCostOre(slot, tc.gridW*0.25/1000, p) + plan := Plan{TotalCostOre: costOre, Actions: []Action{{ + SlotStartMs: 1, SlotLenMin: 15, GridW: tc.gridW, + SoCPct: 50, CostOre: costOre, + }}} + err := ValidatePlan([]Slot{slot}, p, &plan) + if (err != nil) != tc.wantErr { + t.Fatalf("ValidatePlan() error = %v, wantErr %v", err, tc.wantErr) + } + }) + } +} + func TestValidatePlanModeErrorIncludesPowerValues(t *testing.T) { slots := []Slot{{StartMs: 1, LenMin: 15, PriceOre: 100, Confidence: 1, LoadW: 0}} p := Params{ From 1f00ebe2e5af4ea1548fa2b66df4ad185ef6eede Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Sun, 23 Aug 2026 18:26:19 +0200 Subject: [PATCH 3/3] fix(mpc): replay Go DP from live SoC and clip at the operating band Policy lookup stays on the SoC grid. Forward simulation now integrates from the actual initial SoC, including starts outside the configured band, so ValidatePlan and the published trajectory match. An action that would worsen recovery is clipped to the remaining headroom instead of being replaced with idle or snapped after the fact. Signed-off-by: Fredrik Ahlgren --- go/internal/mpc/external_optimizer_test.go | 4 +-- go/internal/mpc/mpc.go | 32 ++++++++++++++++++++-- 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/go/internal/mpc/external_optimizer_test.go b/go/internal/mpc/external_optimizer_test.go index 7b009f878..c0340929a 100644 --- a/go/internal/mpc/external_optimizer_test.go +++ b/go/internal/mpc/external_optimizer_test.go @@ -155,7 +155,7 @@ func TestValidatePlanGridLimitAllowsOnlySubWattSolverResidue(t *testing.T) { const limitW = 11040.0 p := Params{ Mode: ModeArbitrage, CapacityWh: 10000, - SoCMinPct: 10, SoCMaxPct: 95, InitialSoCPct: 50, + SoCMin: 0.1, SoCMax: 0.95, InitialSoC: 0.5, MaxChargeW: 5000, MaxDischargeW: 5000, ChargeEfficiency: 1, DischargeEfficiency: 1, } @@ -183,7 +183,7 @@ func TestValidatePlanGridLimitAllowsOnlySubWattSolverResidue(t *testing.T) { costOre := SlotGridCostOre(slot, tc.gridW*0.25/1000, p) plan := Plan{TotalCostOre: costOre, Actions: []Action{{ SlotStartMs: 1, SlotLenMin: 15, GridW: tc.gridW, - SoCPct: 50, CostOre: costOre, + SoC: 0.5, CostOre: costOre, }}} err := ValidatePlan([]Slot{slot}, p, &plan) if (err != nil) != tc.wantErr { diff --git a/go/internal/mpc/mpc.go b/go/internal/mpc/mpc.go index 313ab5bf2..18e096718 100644 --- a/go/internal/mpc/mpc.go +++ b/go/internal/mpc/mpc.go @@ -491,6 +491,32 @@ func operatingBoundWorsens(from, to, min, max float64) bool { math.Max(0, to-max) > math.Max(0, from-max)+eps } +// clipBatteryPowerToBand reduces a DP action so continuous SoC does not +// worsen operating-bound recovery. Policy is looked up on the grid; energy +// is not, so a charge that lands on max from the nearest grid point can +// overshoot from the real SoC. +func clipBatteryPowerToBand(soc, powerW, dtH, capacityWh, etaC, etaD, min, max float64) float64 { + if capacityWh <= 0 || dtH <= 0 { + return 0 + } + delta := loadpoint.BatteryEnergyDeltaWh(powerW, dtH, etaC, etaD) / capacityWh + if !operatingBoundWorsens(soc, soc+delta, min, max) { + return powerW + } + if powerW > 0 { + headroom := max - soc + if headroom <= 0 || etaC <= 0 { + return 0 + } + return headroom * capacityWh / (dtH * etaC) + } + headroom := soc - min + if headroom <= 0 || etaD <= 0 { + return 0 + } + return -headroom * capacityWh * etaD / dtH +} + func sanitizeOptimizeSlots(slots []Slot) []Slot { out := make([]Slot, 0, len(slots)) for _, s := range slots { @@ -1006,10 +1032,10 @@ func Optimize(slots []Slot, p Params) Plan { pol := Policy[t][si][ei] ba := pol / EA ea := pol % EA - actW := actionAt(ba) + actW := clipBatteryPowerToBand(soc, actionAt(ba), dtH, p.CapacityWh, + p.ChargeEfficiency, p.DischargeEfficiency, p.SoCMin, p.SoCMax) evW := evActionW(ea) - dSoCWh := loadpoint.BatteryEnergyDeltaWh(actW, dtH, p.ChargeEfficiency, p.DischargeEfficiency) - soc2 := soc + dSoCWh/p.CapacityWh + soc2 := soc + loadpoint.BatteryEnergyDeltaWh(actW, dtH, p.ChargeEfficiency, p.DischargeEfficiency)/p.CapacityWh if operatingBoundWorsens(soc, soc2, p.SoCMin, p.SoCMax) { actW = 0 soc2 = soc