goodwe-addon/goodwe_controller/tools/slew_sim.py

#!/usr/bin/env python3
"""SAFETY-39's evidence, on admin/goodwe's OWN plant model.

    python3 tools/slew_sim.py --repro     # ⭐ the 2026-09-05 excursion, replayed
    python3 tools/slew_sim.py --slew      # the derivation: slew against the rail
    python3 tools/slew_sim.py --gain      # what the slew does NOT own
    python3 tools/slew_sim.py --roster    # N=1 -> N=2, and the rail that moved
    python3 tools/slew_sim.py --all

WHY THIS IS A SCRIPT AND NOT A `test_*.py`
==========================================
Byte for byte tools/alloc_sim.py's reason, and it has not changed: this needs
admin/goodwe checked out beside this repository for `plantsim.Plant`, and the
add-on's CI job has no credentials for that repository - it is private and the
job's token is minted for this one. A `test_*.py` here is globbed and RUN by
.gitea/workflows/checks.yml, so a file that needs the sibling repo would either
fail the gate or - worse - pass SILENTLY having checked nothing.

So the rig evidence lives here and is quoted in the pull request and in
DECISIONS/SAFETY-39-slew-against-the-rail.md, and test_control.py's SAFETY-39
section grades the same properties in CI against a plant it models itself and
says it models. The two plants are INDEPENDENT - a first-order lag written in
test_control.py versus `plantsim.Plant`'s ramp-plus-lag with an efficiency and
a SoC integrator - and they agree on every sign and rank below, which is worth
more than either alone.

WHAT IS REAL IN HERE
====================
  * the control law - `compute()` and `Tuning`, imported out of THIS checkout,
    never transcribed.
  * the wire - every command goes out through `frame_check.setpoint_atomic()`,
    the firmware's own encoder, and back in through `plantsim.parse`/`decode`,
    because `Inverter.command()` is what drives the modelled unit.
  * the plant - `plantsim.Plant`, the model the rig uses.

MODELLED, and therefore not evidence about the thing modelled:
  * the meter. `grid = load - pv - sum(battery)` with no poll delay and no
    noise, sampled once per tick. The real P1 is a poll with its own age.
  * the tick. 2.0 s, which is what the live install was MEASURED running on
    2026-09-05; it is a parameter here (`--tick`) and not a constant, because
    the whole of this card is that the derivation depends on the regime.

⚠️ WHAT THIS CANNOT SHOW, SAID BEFORE THE RESULTS RATHER THAN AFTER THEM
=======================================================================
THE -5154 W BATTERY EXCURSION IS NOT REPRODUCIBLE HERE, AND NOT BECAUSE THE
TUNING IS WRONG. `plantsim.Plant.limit()` clamps the plant at `p_max`, so the
modelled battery cannot exceed its rating by one watt at any tuning, on any
step, ever. Neither can test_control.py's model. An acceptance criterion of
"battery power no longer exceeds max_w in simulation" would therefore be
UNFALSIFIABLE - it passes on the broken configuration too - so no criterion in
this card is written on it, and the ticket's own AC to that effect is the one
that had to be rewritten rather than met.

What IS reproducible, and is what everything below measures:
  * the setpoint stepping a full `slew_w` per tick, over and over (the live
    symptom that escalated 2 -> 1 -> 13 -> 43 large steps per 30 min);
  * the loop OVER-ABSORBING - taking more than the site needed - and driving
    the meter from EXPORT into IMPORT, i.e. manufacturing the error it chases.

⚠️ INVERTER B DOES NOT EXIST. `plantsim.INVERTER_B` is inverter A's numbers
COPIED, `measured=False`. Every N=2 row below is therefore a result about TWO
COPIES OF A on a 10000 W rail, NOT about the ~9600 W site that B's datasheet
and serial imply (goodwe/NOTES.md §1.1b; B is a GW5048D-ES rated 4600 W, and
that figure may itself be DC). Each such row prints the caveat with it.
"""
import argparse
import math
import os
import pathlib
import sys

HERE = pathlib.Path(__file__).resolve().parent
ADDON = HERE.parent.parent            # the goodwe-addon checkout root


def find_sim():
    """admin/goodwe's `sim/` directory, or None - and the caller says which.

    $GOODWE_SIM first, the same override test_freeze.py and alloc_sim.py take,
    then the sibling checkouts that actually occur. A missing sibling is not a
    finding about the control law and must never be reported as one.
    """
    cands = []
    if os.environ.get("GOODWE_SIM"):
        cands.append(pathlib.Path(os.environ["GOODWE_SIM"]))
    for up in (ADDON.parent, ADDON.parent.parent):
        cands += [up / "goodwe" / "sim", up / "goodwe-session" / "goodwe" / "sim"]
    for c in cands:
        if (c / "plantsim.py").exists():
            return c
    return None


def load(sim):
    sys.path[:0] = [str(sim), str(sim.parent / "master")]
    os.environ.setdefault("GOODWE_ADDON_REPO", str(ADDON))
    sys.path[:0] = [str(ADDON / "goodwe_controller")]
    import plantsim                                          # noqa: E402
    from app import control                                  # noqa: E402
    return plantsim, control


# ---------------------------------------------------------------------------
# one experiment


def step_run(plantsim, control, *, gain, slew_w, rail_w, n=1, tick_s=2.0,
             export_w=3385.0, load_w=350.0, phase=5, cycles=300, dt=0.05,
             lag_s=None):
    """One export step onto a converged loop, through the real wire.

    `rail_w` is the SITE rail - what `main.py` hands `compute()` as `max_w`
    (`Controller.site_max_w`, the roster sum). Each of the `n` modelled units
    is rated `rail_w / n` and is commanded its even share, which is the
    arrangement ARCH-01's split-invariance says any valid split is equivalent
    to at the meter.
    """
    over = {} if lag_s is None else dict(lag_s=lag_s)
    spec = dict(plantsim.INVERTER_A, rating_w=rail_w / n, **over)
    invs = [plantsim.Inverter.from_spec(name, spec) for name in "AB"[:n]]
    site = plantsim.Site(invs, load=load_w, pv=0.0)
    tuning = control.Tuning(gain=gain, max_w=rail_w, slew_w=slew_w,
                            target_grid_w=-10.0)
    prev, i_w, sat = 0.0, None, 0
    worst_import, peak, biggest, last = -1e9, 0.0, 0.0, None
    for c in range(cycles):
        if c == phase:
            site.pv = export_w + load_w
            site.grid = site.load - site.pv - site.battery
        # ⚠️ ONE UNIT'S POWER, NOT THE SITE'S SUM, at n > 1. The shipped add-on
        # binds a single `inverter_power` sensor, so a two-inverter install has
        # a loop that sees unit A only. sitesim.run_site() makes the same choice
        # and says at length why feeding the sum would be a KINDER rig than the
        # deployment.
        actual = invs[0].plant.actual
        d = control.compute(prev, site.grid, actual, tuning, sat, i_w)
        if c >= phase and last is not None:
            biggest = max(biggest, abs(d.target_w - last))
        last = d.target_w
        prev, i_w, sat = d.target_w, d.i_w, d.sat_count
        for inv in invs:
            inv.command(d.target_w / n)
        site.settle(tick_s, dt=dt)
        if c >= phase:
            worst_import = max(worst_import, site.grid)
            peak = max(peak, abs(site.battery))
    return dict(imp=worst_import, over=peak - export_w, step=biggest)


def family(plantsim, control, *, exports=None, phases=(5, 6, 7), **kw):
    """The worst case over a family of step sizes and phases.

    ⚠️ A FAMILY, NOT ONE STEP, AND THAT IS ITSELF A FINDING. A single run of
    this loop is not evidence: the 15 W deadband, the 10 W quantisation and the
    phase of the step against the tick move the answer by hundreds of watts,
    NON-MONOTONICALLY, so two adjacent slews rank either way depending on which
    step you happen to run. The first version of this derivation was a single
    step and reported gain 0.3 as WORSE than gain 0.6, which the family shows is
    backwards. Everything printed below is a worst case over the family, which
    is how admin/goodwe grades the same kind of quantity (config.yaml's "1890
    closed-loop runs", "924 runs").
    """
    exports = exports or range(1000, 4300, 200)
    worst = dict(imp=-1e9, over=-1e9, step=0.0, big=0, runs=0)
    for export_w in exports:
        for phase in phases:
            r = step_run(plantsim, control, export_w=float(export_w),
                         phase=phase, **kw)
            worst["runs"] += 1
            worst["imp"] = max(worst["imp"], r["imp"])
            worst["over"] = max(worst["over"], r["over"])
            worst["step"] = max(worst["step"], r["step"])
            worst["big"] += r["step"] >= 900.0
    return worst


def row(label, w):
    print(f"  {label:<44} step {w['step']:6.0f} W   "
          f"{w['big']:3}/{w['runs']} runs >= 900 W   "
          f"over-absorbed {w['over']:6.0f} W   meter to {w['imp']:+7.0f} W")


# ---------------------------------------------------------------------------


def repro(ps, ct, tick_s):
    print("⭐ THE 2026-09-05 EXCURSION, REPLAYED - gain 0.6, max_w 5000, "
          f"slew_w 1000, {tick_s:.0f} s tick")
    print("   Live, 11:43:55-11:45:55: grid -3385 W, setpoint 380 -> 5000 in five")
    print("   ticks, battery -5154 W, grid +2831 W. One 3385 W export step here:\n")
    print(f"   {'t':>5} {'grid in':>9} {'setpoint':>9} {'battery':>9} {'grid out':>9}  reason")
    spec = dict(ps.INVERTER_A)
    inv = ps.Inverter.from_spec("A", spec)
    site = ps.Site([inv], load=350.0, pv=0.0)
    t = ct.Tuning(gain=0.6, max_w=5000.0, slew_w=1000.0, target_grid_w=-10.0)
    prev, i_w, sat = 0.0, None, 0
    for c in range(34):
        if c == 5:
            site.pv = 3385.0 + 350.0
            site.grid = site.load - site.pv - site.battery
        g = site.grid
        d = ct.compute(prev, g, inv.plant.actual, t, sat, i_w)
        prev, i_w, sat = d.target_w, d.i_w, d.sat_count
        inv.command(d.target_w)
        site.settle(tick_s, dt=0.05)
        print(f"   {c * tick_s:5.0f} {g:9.0f} {d.target_w:9.0f} "
              f"{site.battery:9.0f} {site.grid:9.0f}  {d.reason}")
    print("\n   ⚠️ THE BATTERY NEVER PASSES 5000 W HERE AND CANNOT: plantsim.Plant")
    print("      clamps at its rating. The -5154 W figure is NOT reproducible in")
    print("      this model at any tuning, so nothing in this card is graded on it.")
    print("      What IS reproduced is the setpoint stepping 1000 W a tick to the")
    print("      rail, and the meter crossing from export into import.")


def slew_table(ps, ct, tick_s, rail):
    frac = ct.SLEW_RAIL_FRAC
    print(f"THE DERIVATION - gain 0.6, site rail {rail:.0f} W, {tick_s:.0f} s tick, "
          f"51 export steps 1.0-4.2 kW\n")
    print(f"  the constraint: the loop must not cross its whole output range "
          f"faster than the")
    print(f"  plant settles. {rail:.0f} W / ({ct.PLANT_SETTLE_S:.0f} s / "
          f"{tick_s:.0f} s per tick) = {frac * rail:.0f} W per tick.\n")
    overs = []
    for slew in (1000, 800, 700, 600, int(frac * rail), 500, 400, 300, 250):
        w = family(ps, ct, gain=0.6, slew_w=float(slew), rail_w=rail,
                   tick_s=tick_s)
        tag = "   <- the derived value" if slew == int(frac * rail) else (
              "   <- the 2026-09-05 live pair" if slew == 1000 else (
              "   <- the live interim" if slew == 400 else ""))
        row(f"slew_w {slew}{tag}", w)
        overs.append(w["over"])
    print("\n  READ IT THIS WAY. The >= 900 W step - the live symptom, which")
    print("  escalated 2 -> 1 -> 13 -> 43 per 30 min that morning - is gone from")
    print("  every run below 900 W of slew, and that is what the slew owns.")

    # ⚠️ THE NARRATION BELOW IS DERIVED FROM THE COLUMN JUST PRINTED, NOT
    # WRITTEN ABOUT IT. SAFETY-39 wrote "the over-absorption does not move, and
    # at the slowest slews it is slightly WORSE" - true of release/1.0's law,
    # which is what it was measured against, and FALSE under SAFETY-48's lead
    # cap, where the same script on the same plant prints a monotonically
    # FALLING column. A transcript whose prose contradicts its own table is a
    # failing check (MUTATION.md habit 8), and a sentence that is right under
    # one law and wrong under another has to be computed rather than typed. So:
    # which law is imported is read off `control.py` itself, and whether the
    # column falls is read off the column.
    cap = getattr(ct, "SLEW_LEAD_STEPS", None)
    falling = all(b < a for a, b in zip(overs, overs[1:]))
    law = ("release/1.0's law - the integrator's lead over the command is "
           "UNBOUNDED" if cap is None else
           f"SAFETY-48's law - the integrator may lead the command by "
           f"{cap} slew step(s)")
    print(f"\n  THE LAW THIS RAN AGAINST: {law}.")
    if cap is None and not falling:
        print("  The OVER-ABSORPTION does not move with the slew, and at the")
        print("  slowest slews it is slightly WORSE - `slew_w` bounds the")
        print("  OUTPUT while the integrator goes on winding, so a slower slew")
        print("  leaves it winding for longer. SAFETY-48 bounds that lead and")
        print("  INVERTS this column; see DECISIONS/SAFETY-48 §6.")
    elif cap is not None and falling:
        print("  The OVER-ABSORPTION FALLS MONOTONICALLY with the slew, which")
        print("  is the OPPOSITE of what SAFETY-39 measured on release/1.0")
        print("  (958/1107/1133/1058/987/1029/1150/1143/1081 W, flat and")
        print("  non-monotone). Bounding the integrator's lead is what changed")
        print("  it: the slew now owns a share of the over-absorption it did")
        print("  not own before. DECISIONS/SAFETY-48 §6 has both tables and an")
        print("  independent confirmation on test_control.py's own plant.")
    else:
        print("  ⚠️ NEITHER NARRATION FITS. The law above and the shape of the")
        print(f"  column ({'falling' if falling else 'not falling'}) are the")
        print("  combination this script has no sentence for - which means the")
        print("  measurement has moved and DECISIONS/SAFETY-39 §3 and")
        print("  DECISIONS/SAFETY-48 §6 both need re-reading before this table")
        print("  is quoted anywhere.")
    print("\n  EITHER WAY THE RIG DOES NOT PICK THE NUMBER, and it is not asked")
    print("  to: the CONSTRAINT picks it, and this table is what says what the")
    print("  constraint does and does not buy.")


def gain_table(ps, ct, tick_s, rail):
    frac = ct.SLEW_RAIL_FRAC
    print(f"WHAT THE SLEW DOES NOT OWN - slew_w at the derived "
          f"{frac * rail:.0f} W, gain swept\n")
    for gain in (0.6, 0.5, 0.4, 0.3, 0.2, 0.1):
        w = family(ps, ct, gain=gain, slew_w=frac * rail, rail_w=rail,
                   tick_s=tick_s)
        row(f"gain {gain}" + ("   <- shipped" if gain == 0.6 else ""), w)
    print("\n  `slew_w` bounds the OUTPUT. The integrator accumulates")
    print("  `gain * error` per cycle and is bounded only by `integrator_max_w`,")
    print("  so at gain 0.6 into a 3.4 kW error it moves ~2000 W per tick and")
    print(f"  reaches the {rail:.0f} W rail in under three ticks WHATEVER the slew is.")
    print("  The over-absorption is therefore a GAIN term. `gain: 0.6` is the")
    print("  next constant to re-derive; SAFETY-39 measured it and changed")
    print("  nothing, because two control-law constants moving in one release on")
    print("  a running house makes the next measurement uninterpretable.")


def roster_table(ps, ct, tick_s):
    frac = ct.SLEW_RAIL_FRAC
    print("N=1 -> N=2: THE RAIL IS WHAT MOVED\n")
    print("  ⚠️ INVERTER B DOES NOT EXIST (SPIKE-04, ~2026-09-12) AND")
    print("     plantsim.INVERTER_B IS INVERTER A COPIED, measured=False. The")
    print("     N=2 rows are therefore about TWO COPIES OF A on a 10000 W rail.")
    print("     The site B actually implies is ~9600 W (5000 + 4600), and 4600 is")
    print("     B's DATASHEET and SERIAL, not a measurement, and may be DC -")
    print("     goodwe/NOTES.md §1.1b. Read the rows as architecture, not as a")
    print("     prediction about this house.\n")
    for n, rail in ((1, 5000.0), (2, 10000.0)):
        for slew in (1000.0, frac * rail):
            w = family(ps, ct, gain=0.6, slew_w=slew, rail_w=rail, n=n,
                       tick_s=tick_s)
            lab = (f"n={n} rail {rail:.0f} slew {slew:.0f}"
                   + ("  [B = A COPIED]" if n == 2 else ""))
            row(lab, w)
    print(f"\n  The derived value at the REAL 9600 W site is "
          f"{frac * 9600:.0f} W per tick, and")
    print("  validate.py reads it off `sum(unit_caps(opts))` - the same number")
    print("  `Controller.site_max_w` builds - so adding inverter B moves the")
    print("  bound with no operator action and no second derivation.")
    print("\n  ⚠️ AND THE DIRECTION IS NOT THE ONE THE TICKET ASSUMED. A rail that")
    print("  GROWS loosens this bound: slew_w 1000 crosses a 9600 W rail in")
    print("  19.2 s and PASSES, while at 5000 W the same 1000 is the")
    print("  configuration that was measured failing. That is what the")
    print("  constraint says and it is correct - the bound is on how long the")
    print("  crossing TAKES, and a bigger rail at the same step size takes")
    print("  longer to cross.")
    print("\n  ⚠️ AND THE >= 900 W STEP COUNT STAYS AT 48/51 AT THE DERIVED VALUE")
    print("  ON THE N=2 ROW, WHICH IS NOT A FAILURE - IT IS THE THRESHOLD GOING")
    print("  STALE, on the same day and in the same way as the constant. 900 W")
    print("  was a meaningful step against a 5000 W rail; against 10000 W the")
    print("  derived step IS 1111 W, and a check written on the absolute 900")
    print("  would report the correct tuning as broken. Anything grading this")
    print("  symptom at N>=2 has to grade it as a FRACTION of the rail too.")
    print("\n  WHAT DOES NOT GET WORSE HERE, STATED BECAUSE IT WAS EXPECTED TO.")
    print("  Over the SAME family of 1.0-4.2 kW export steps, doubling the rail")
    print("  does not worsen the over-absorption (958 -> 810 W): the extra rail")
    print("  is headroom these disturbances never ask for. That is a statement")
    print("  about this family, not about the site - a bigger site sees bigger")
    print("  steps, and the over-absorption is a `gain` term either way, which")
    print("  the enforced ratio does not bound and this card did not change.")


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--repro", action="store_true")
    ap.add_argument("--slew", action="store_true")
    ap.add_argument("--gain", action="store_true")
    ap.add_argument("--roster", action="store_true")
    ap.add_argument("--all", action="store_true")
    ap.add_argument("--tick", type=float, default=None,
                    help="seconds per control cycle (default: control.LOOP_TICK_S)")
    ap.add_argument("--rail", type=float, default=5000.0)
    a = ap.parse_args()
    if not any((a.repro, a.slew, a.gain, a.roster, a.all)):
        a.all = True
    sim = find_sim()
    if sim is None:
        sys.exit("slew_sim: admin/goodwe's sim/ is not beside this checkout, so "
                 "plantsim.Plant is unavailable. Clone it as a sibling or set "
                 "$GOODWE_SIM. This is a missing repository, NOT a result about "
                 "the control law - do not report it as one.")
    ps, ct = load(sim)
    tick = a.tick if a.tick is not None else ct.LOOP_TICK_S
    print(f"slew_sim: plant from {sim}")
    print(f"          control law from {ADDON / 'goodwe_controller' / 'app'}")
    print(f"          {ps.Inverter.from_spec('A', ps.INVERTER_A).describe()}")
    print(f"          tick {tick:.1f} s, SLEW_RAIL_FRAC {ct.SLEW_RAIL_FRAC:.4f} "
          f"= {ct.LOOP_TICK_S:.0f} s / {ct.PLANT_SETTLE_S:.0f} s\n")
    if a.repro or a.all:
        repro(ps, ct, tick); print()
    if a.slew or a.all:
        slew_table(ps, ct, tick, a.rail); print()
    if a.gain or a.all:
        gain_table(ps, ct, tick, a.rail); print()
    if a.roster or a.all:
        roster_table(ps, ct, tick)


if __name__ == "__main__":
    main()