Powerflow performance pass, and a cache that cannot outlive a failure by BDonnot · Pull Request #189 · Grid2op/lightsim2grid

BDonnot · 2026-09-04T07:31:41Z

A profiling-driven pass over the C++ powerflow path, then the correctness work that came out of it. Every performance item was measured before it was kept: a git worktree baseline at the previous commit, matched build configs, callgrind instruction counts (deterministic) plus a dedicated wall-clock instrument, and full-precision result dumps compared bit-for-bit. Several candidates measured worse and were reverted rather than committed — they are not in this branch.

Rebased on top of #188, which this branch merges (2dfaa18).

Performance

Measured on case9241pegase unless stated, NR + KLU, callgrind Ir.

change	effect
Feed `Contrib` to `setFromTriplets` directly (`4aa6a85`)	drops a whole redundant `cntrb` → `triplets` copy
Build J's CSC arrays directly in `build_J_sparsity` (`4c4398d`)	no triplets, no temporary transposed matrix; counting sort with duplicate collapse folded into the transpose
Same treatment for the DC slack-free matrix (`08e8b60`)	built in place, two passes over `ref_mat`
`fillYbus` via the same assembly	reuses the `setFromTriplets` alternative rather than copy-pasting it
Shunt and branch results on real/imag parts (`93088bc`, `77f09e0`)	avoids the NaN-recovery branch `std::complex::operator*` carries after every product; bit-identical for finite inputs
Fuse `v_kv_from_vpu` + `v_deg_from_va` (`780baf5`)	one walk instead of two
One-pass `_get_amps` with no temporaries (`e8e2e1e`)	`compute_results` −6.7% Ir, −14% wall
dS assembly gets both derivatives from one product (`4260235`)
Stop zeroing J, stop rebinding a `Ref` per coefficient (`4fab621`)
FDPF stops rediscovering `(Vm, Va)` with a `hypot` and an `atan2` (`bf2eaf2`)
Bounds-check ids we generated only in debug builds (`3ae9427`, `1bd626d`, `1d68a29`)	see below

Internal vs user-facing bounds checks

_check_in_range (always active) is now distinguished from _check_in_range_internal (#ifndef NDEBUG). Only ids the library generated itself are compiled out; no check reachable from the python API is. Verified from the binaries rather than by reading: message-string counts in Release vs Debug libraries, plus a release-built probe confirming 6/6 user-facing rejections still throw.

5a76ac5 also removed the last throw sites from inside the bus-counting bracket, which recovered the full 4.9M instructions per solve that #188's try/catch had cost — not through virtualization (disassembly showed 1196 vs 1205 instructions, 4 indirect calls each) but through the unwind edge, which changes GCC's scheduling of unrelated inlined code in the same TU. Bisected to that single construct; reverting only the try/catch restored the old figure to the instruction.

Correctness

A powerflow that throws can no longer leave a cache claiming to be in sync

A solve rebuilds its family's cache in place — bus labelling, Ybus/Sbus, the PV/PQ split, slack weights, a factorization. A throw part-way through left a mixture of the old grid and the new one behind, with the change flags still saying whatever they said before the call.

ac_pf / dc_pf now run the whole solve against a copy of the change tracking and hold the grid itself at "everything changed, both families" for the duration; the copy becomes the grid's change tracking only at a publication statement after compute_results, which a throw never reaches. Deliberately not a try/catch — an unwind edge through this code is not free, as above — while the copy is 24 bools: +21 instructions per solve, against 206M.

Input validation at the top of ac_pf / dc_pf (Vinit size, max_iter/tol, a droop grid handed to a solver that cannot do droop) stays outside that bracket: it runs before anything is touched, so a rejected call still costs the caller's cache nothing.

build_solver_input / build_dc_solver_input got the same treatment from the other side: they retire this grid's cache before copying the caller's labelling into it rather than after.

The per-bus element counts are rebuilt when the cache cannot be reused

Everything a solve builds is derived from the elements, here and now — except the counts, which are carried forward +1/-1 from every mutation the grid has ever seen. init_bus_status() only established them when they had never been counted; an armed-but-drifted count was inherited for the life of the grid.

That matters more than a stale cache would, because since #188 connectivity is the counts: a phantom bus enlarges the solved system and shifts every bus id after it, and nothing downstream can notice — an off-by-one count reads exactly like a real one.

They are now rebuilt from the elements whenever the answer to "may I reuse what was built for this grid?" is no. Hooked to that, not to redo_all: the latter also fires on has_dimension_changed(), which an ordinary topology change raises on every grid2op step, and an O(all elements) recount per step is precisely the cost these counts exist to avoid.

path	before	after
cache reuse (grid2op steady state)	619,932,638	619,932,635 (−3 over three solves)
full rebuild	863,811,268	879,669,208 (+5.29M/solve, +1.8%)

All of the delta is recompute_bus_element_counts (15,857,772 of 15,857,940).

Together with the point above, this is what closes the loop on a throw out of a mutator: the backend calls tell_solver_need_reset() on any exception, so the next solve recounts from the elements and the corruption dies there.

Tests

New C++ cases, each checked against the parent commit's library to confirm it actually fails there:

A plugin solver that converges on its first call and throws on every later one, run twice on the same grid — nothing is swapped between the two calls, so the state after the second is the throw's doing. Before: the grid still answered nothing_changed() on both families. Also: a solver throwing from compute_pf (AC and DC), a throw raised after the algorithm ran, and the success path still publishing exactly what it consumed.
A drift staged by hand — no public API can produce one any more — and checked that it dies at the next invalidation: a phantom bus (lost decrement, changes connectivity), a count one too low (lost increment, latent until the bus empties one element early), the DC family, and the negative case that an ordinary reusing solve does not pay for the recount.
New test_fdpf_algorithm.cpp, plus additions to test_powerflow_algorithm.cpp.

229 test cases / 590,856 assertions pass in C++17 Release, C++14 Release, and Debug with assertions active. 16/16 grid × algorithm result dumps are bit-identical across the performance work.

Python: 507 passed / 97 skipped / 461 subtests. The remaining failures and errors in this container are all missing optional dependencies (grid2op, pypowsybl), not assertion failures.

Notes for review

LSGrid's member layout changes, so anything casting an LSGrid across a module boundary (gpusim2grid) must be rebuilt against these headers — which the plugin ABI policy in docs/solver_plugin.rst already requires. Nothing changes for python.
process_results gained a parameter (protected; ac_pf/dc_pf are its only callers) and the _mark_cache_valid(bool) shorthand is gone — it marked the member control, which is now exactly the bug the working-copy protocol removes.
Still open, and deliberately not in this branch: fillYbus position-caching for the need_recompute_ybus-without-pattern-change path (14.1M of setFromTriplets), cheapening get_free_vm_slack_solver_buses (5.61M, called twice per solve) and fill_voltage_control_solver_data (4.10M), and the bus power mismatch computed twice per solve (~2.1M).

🤖 Generated with Claude Code

https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy

Generated by Claude Code

The powerflow inner loops asked Eigen for a few things it can only do by allocating, or by redoing work, and paid for it once per iteration -- or, in Gauss-Seidel, once per bus per sweep. Gauss-Seidel (~13x, and ~2x for the synchronous variant, on a 121-bus mesh with identical sweep counts and bit-identical voltages): `Ybus.row(k) * V_` reads a ROW of a matrix that arrives column-major, and a row of a column-major sparse matrix is not stored anywhere -- getting it means scanning every column. Once per pq bus and twice per pv bus, that made a sweep cost O(nb_bus * nnz) instead of O(nnz). `Ybus.coeff(k, k)` was a binary search inside column k, every time. Both now come from caches built once per compute_pf: a row-major copy of Ybus, and its diagonal. `BaseAlgo::_evaluate_Fx` computed its sparse * dense product three times (twice in the multi-slack overload). The mismatch was kept as a lazy expression, and a lazy expression holding a sparse * dense product re-evaluates it at every evaluation -- once per result segment assigned. The Newton-Raphson iteration (~7-14% off the algorithm-side cost, i.e. everything but the linear solver) stops allocating and freeing a full vector per call in its residual, trial-voltage and dS assembly: those buffers now live as long as the system. Every `Ybus * V` lands in a named buffer with .noalias() rather than inside a coefficient-wise expression Eigen has to allocate a temporary to evaluate, and NRAlgo scales its step in place instead of handing apply_step a `coeff * F` expression its Eigen::Ref parameter had to materialise. `dS_dVm` / `dS_dVa` are plain value arrays now, not two full SparseMatrix copies of Ybus: nothing ever indexed them by (row, column), so each topology change is two structure copies lighter -- the "TODO speed: copy only the sparsity pattern" that sat in init_topology. All of it is value-preserving: the C++ suite passes unchanged (C++14, C++23 and a Debug build with assertions), and every solver family returns voltages identical to the last bit on the benchmark grids. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

has_converged rebuilds V from (Vm, Va) and then has to put that pair back in canonical form -- magnitude >= 0, angle wrapped -- because the Q iteration (`Vm_(pq) -= q_`) can overshoot a magnitude past zero and the P iteration accumulates into Va_ without ever wrapping it. It did that with `Vm_ = V_.array().abs(); Va_ = V_.array().arg();`: a hypot and an atan2 per bus, twice per iteration, to recompute two numbers the solver already holds. V is built as Vm_ * exp(i.Va_), so its modulus is |Vm_| and its argument is Va_, plus half a turn where the magnitude went negative. Worth ~1.25-1.35x on the whole solve, both flavours, 49 to 1600 buses. The complex voltage the solve consumes is built exactly as before and is unchanged bit for bit; Vm_ and Va_ move in their last bits, which shifts a converged solution by ~1e-13 -- four orders inside the 1e-9 tolerance, with identical iteration counts on all 80 grid / loading / flavour combinations checked and no change in the distance to a Newton-Raphson reference. The canonicalisation is a named static because its two interesting branches turn out to be unreachable from any converging solve: a trajectory that overshoots a magnitude ends up diverging, and a diverged solve clears its state. They needed a test that does not go through a powerflow. Which exposed that the Fast-Decoupled family had no answer-level coverage at all -- it reached the suite only through test_cache_reuse.cpp (caching, not numbers) and test_plugin_registration.cpp (names). test_fdpf_algorithm .cpp now pins that both flavours converge to the Newton-Raphson solution of the same grid, that the reported Vm / Va stay consistent with the reported V, and that canonicalise_vm_va matches the abs()/arg() pair over Vm in [-3, 3] x Va in [-10, 10] rad. The one deliberate difference is Vm == 0 exactly, where the phase of a zero phasor is undefined and the solve cannot converge anyway. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

The synthetic meshes the earlier perf work used are not a good proxy for the grids lightsim2grid actually runs: they say the algorithm side (everything but the linear solver) is ~20% of a solve, where case118 and the PEGASE cases put it at 28-55%. Two scripts, split so the expensive part stays out of the measurement. benchmark_matpower_to_binary.py reads a MATPOWER case -- which goes through pandas and dwarfs the powerflow -- and writes a ".lsb". benchmark_binary_powerflow.py reloads that binary and runs N powerflows; python only holds the initial-voltage array, so timing or profiling it measures C++ and essentially nothing else. It reports the solver's own timers (get_timers_jacobian), which split a solve into its C++ phases without needing a profiler, and its docstring carries the callgrind recipe for function-level attribution on top -- --toggle-collect on ac_pf so start-up and load_binary stay out, and a note that a normal (stripped) install still resolves the exported template instantiations, which is enough to name the hot functions. Both modes matter and stress different code: --rebuild re-derives the bus labelling, Ybus/Sbus, the Jacobian sparsity and the symbolic factorization on every solve (a grid2op step that changed topology), --reuse keeps them (a TimeSeries workload). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Two findings from profiling case118 and the PEGASE cases under callgrind. fill_J zeroed the whole value array before filling because every write accumulated -- several contributions being able to land on one coefficient. Only the feature entries actually need that. The four dS families live at (p_row, theta_col), (p_row, vm_col), (q_row, theta_col) and (q_row, vm_col); a row is a P equation or a Q equation and a column a theta unknown or a vm unknown, never both, so the families are pairwise disjoint, and within a family the ledger hands every bus its own row and column, so no two dS entries ever claim the same coefficient. They can assign. Feature entries still accumulate, because a component may legitimately add to a dS coefficient -- an hvdc droop slope adds to the dP/dtheta of its end buses -- so their positions, a handful against every nonzero of J, are the only ones zeroed. A megabyte no longer gets rewritten at every factorisation, and the dS passes no longer read each coefficient before writing it. ~1.2x on the fill for case9241pegase. find_J_pos bound an Eigen::Ref to J on every call, and build_J_sparsity calls it four times the Ybus nonzeros -- 170k times per solve on a 9241-bus grid, the most-called function in a powerflow. The Ref does not copy (it aliases a compressed, same-Options matrix) and costs about nothing to build; paying for one per coefficient to read two pointers that never change across the loop is what added up. It takes the index arrays now, read once outside the loop. ~1.1x on the pre-processing of a rebuild. Values are unchanged to the bit: case118 and the three PEGASE cases, across NR_KLU / NRSing_KLU / NR_SparseLU / FDPF_XB_KLU, rebuild and cache-reuse, byte-identical dumps at full precision. The two properties the first change rests on are asserted in debug builds where the layout is decided, and pinned behaviourally by two tests: a second fill with nothing changed must not move J, and filling at one voltage state then another must match a system that only ever saw the second -- both with and without a droop hvdc line, the case whose feature entries land on dS coefficients. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

fill_internal_variables is the most expensive thing lightsim2grid itself does in a Newton solve -- 13.7% of one on case9241pegase, against 63% for KLU -- and callgrind put two thirds of it inside <complex>. Writing Y for Ybus(i, j), the formulas derived from pandapower are dS_dVm(i, j) = conj(Y . Vnorm_j) . V_i dS_dVa(i, j) = -conj(Y . V_j) . i . V_i and Vnorm_j is just V_j / |V_j|, so the first is the second's product divided by a real. With B = conj(Y . V_j) . V_i -- two complex products -- dS_dVm is B / |V_j|, a real scaling, and dS_dVa is -i . B, a swap and a sign flip. The diagonal's two extra terms share a product the same way: conj(Y . V_i - Ibus_i) . V_i is B - conj(Ibus_i) . V_i, and its dS_dVm term conj(Ibus_i) . Vnorm_i is that same value over |V_i|. The arithmetic is written on real and imaginary parts rather than left to std::complex, whose operator* carries a NaN-recovery branch after every product; that half was worth 1.19-1.31x on its own and was bit-identical. Two of the four nb_bus scratch vectors go away with the rewrite: the unit phasors and the conj(Ibus) products are each used once per bus, on the diagonal, so they are built there rather than carried through the solve. 1.26-1.39x on the phase, on every grid from 118 to 9241 buses, rebuild or cache-reuse; -36% of its instruction count. Values move by less than one ulp relative (1.9e-16 on case9241pegase): the rounding that differs is dividing by |V_j| after the product rather than before. This is the Jacobian, not the residual -- it steers the Newton step, it does not define the answer -- and across case118 and the three PEGASE cases on NR / NRSing / SparseLU, iteration counts are unchanged everywhere and converged voltages agree to 7.3e-13, five orders inside the 1e-8 tolerance. FDPF builds no Jacobian and stays bit-identical. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

…ser's Profiling a powerflow on case9241pegase found GenericContainer:: _check_in_range and the _get_bus it guards reached 128k times per solve -- four passes over both ends of every branch, in fillYbus, compute_results and reconnect_connected_buses. Every one of those call sites is a loop shaped for(el_id = 0; el_id < nb(); ++el_id), where the bound is a property of the loop: the check can only fire if the bug is here rather than in the caller. Both were also out-of-line calls rather than inlined compares, because the error paths build an ostringstream. So there are two checks now. _check_in_range is unchanged and still always throws -- it is for ids that crossed the python boundary, which are never trusted, and where the alternative is an out-of-bounds read on a bit-packed std::vector<bool>. _check_in_range_internal is its debug-only twin, for ids this library produced itself; it feeds the new _get_bus_internal and the get_bus_side_{1,2}_internal accessors that the branch containers' own loops now call. The assertion and sanitizer CI builds keep both, since USE_DEBUG_ASSERTS clears NDEBUG. Nothing reachable from python lost a check: get_bus_load / get_bus_gen / get_bus1_powerline and friends, deactivate / reactivate / change_bus, change_ratio / change_shift, set_regulated_bus, set_status_droop and update_slack_weights_by_id all still go through the always-on form. test_LSGrid_out_of_bounds.py -- the suite written for exactly this contract -- passes unchanged (10 tests, 82 subtests), and a release build still raises out_of_range for a bad id on each user-facing accessor. Worth ~2.9% of everything lightsim2grid itself does in a rebuild solve (14.7M instructions of 919M). Results bit-identical on case118 and the three PEGASE cases. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

The generic dS pass of NRSystem::build_J_sparsity records one Contrib per Jacobian coefficient a dS matrix feeds, and that vector has to outlive the matrix build because the dS value maps are resolved from it afterwards. The code then allocated a second, equally large vector of Eigen::Triplet and refilled it with the same (row, col) pairs plus a literal 0., only because setFromTriplets wants row()/col()/value() and Contrib spelled them jrow()/jcol() with no value at all. On case9241pegase that is 2.07 MB duplicating the 2.07 MB already held. Contrib now also answers to Eigen's triplet protocol (value() is the zero: this pass builds the sparsity pattern only, fill_J writes the numbers). The feature entries the components declare are appended to the same vector so a single range describes the whole pattern, and the map-resolution loops run over the leading dS entries. Worth 12.9M instructions of 908M on a case9241pegase rebuild solve (-1.4% of everything the library does, C++ side): the 8.5M spent constructing triplets is gone, and set_from_triplets itself drops 12.5M -> 6.6M because the zero is a constant the compiler folds into the store rather than a value loaded per entry. Eigen receives the same pattern in the same order, so results are bit-identical: verified on case118 and the three PEGASE cases across NR, NRSing, NRRefactorRetry, FDPF XB and FDPF BX with SparseLU and KLU. The unit tests pass in the C++17 and C++14 builds, and the Debug build exercises the build_J_sparsity writer-multiplicity assertions. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

setFromTriplets runs the same two counting sorts this function needs -- bucket the entries by J row, then transpose into columns, the transpose being what sorts each column -- but it carries a double per entry through a temporary SparseMatrix, collapses duplicates in a pass of its own, and materialises the transpose as a second matrix. And because it hands back only a matrix, locating each contribution in the result took 4 * nnz(Ybus) binary searches afterwards: ~170k per rebuild on a 9241-bus grid, 12.1M instructions, with Eigen's own machinery accounting for 23.5M more. Running the two sorts here carries a 4-byte entry id instead of the value, folds the duplicate collapse into the transpose, and knows each coefficient's position at the moment it writes its row index -- so the dS maps and the feature positions are filled straight from that walk and the binary searches are gone. Duplicate entries stay supported, which is what rules out the cheaper sorted-triplet path: a feature entry may legitimately land on a dS coefficient (an hvdc droop slope adds to the dP/dtheta of its end buses). 1.5-1.6x on the phase that contains it: pre_proc, which also covers update_state and init_topology, goes 5.06-5.21 -> 3.36-3.39 ms per solve on case9241pegase, 1.39-1.46 -> 0.88-0.90 on case2869pegase and 0.60-0.63 -> 0.38-0.39 on case1354pegase. The function itself goes 58.7M -> 40.4M instructions for three rebuilds (1.45x), and a rebuild solve drops 2% of its total instruction count. The invariants Eigen used to provide are asserted in the debug build (J compressed, every column filled exactly to its end, row indices strictly increasing). The build was also cross-checked entry by entry against setFromTriplets + find_J_pos over 409 sparsity builds of the C++ test suite -- which is what exercises multi-slack, voltage control and hvdc, and where the 856 duplicate entries turn up -- and on the four benchmark grids: same nnz, same outer array, same inner array, same position for every contribution. Results are bit-identical on case118 and the three PEGASE cases across NR, NRSing, NRRefactorRetry, FDPF XB, FDPF BX and DC with both SparseLU and KLU. Tests pass in the C++17, C++14 and Debug builds. Two variants measured worse and were dropped: holding the scratch buffers as members to skip the per-rebuild allocation (the extra indirection in the two hot walks cost more than the allocations saved), and leaving J's value array uninitialised. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Nothing in remove_slack_buses ever needed sorting or merging. mat_bus_id_ is a monotone compaction -- fill_mat_bus_id hands out consecutive ids in bus order, skipping the slack buses -- and the inner iterator walks each column of the source in increasing row order, so the coefficients that survive the row/column deletion come out already in the order a compressed column-major matrix stores them, one per coefficient. setFromTriplets cannot know that, so it re-derived that exact order the hard way: bucket by row, transpose back into columns, both through a temporary SparseMatrix carrying a double per entry. That cost 5% of a DC solve. The matrix is now written in place in two passes over the source, one to size each column (straight into the outer array, which a prefix sum turns into the column starts) and one to fill it. The triplet vector -- half a megabyte on case9241pegase -- goes with it, and res_mat is resized rather than reassigned, so a rebuild at constant size re-uses what the last one allocated. 3.09x on the function (11.4M -> 3.7M instructions for three rebuilds of case9241pegase) and -4.9% on a whole DC solve (dc_pf 158.5M -> 150.8M). Wall clock per rebuild solve: 12.97 -> 12.53 ms on case9241pegase, 3.39 -> 3.31 on case2869pegase, 1.53 -> 1.45 on case1354pegase. The sortedness this rests on was measured rather than assumed: an assertion over every consecutive pair found the old triplet list strictly ordered column-major with no duplicates in all 136 builds of the C++ test suite and on the four benchmark grids under both DC solvers. The new build was then cross-checked against the old one coefficient by coefficient -- same nnz, same outer array, same inner array, same values bit for bit -- over the same 136 builds and the same grids. What stays in the code is the cheaper permanent form of that check, in debug builds: compressed, each column filled exactly to the size the first pass computed, row indices strictly increasing (which is also what would catch mat_bus_id_ losing its monotonicity). Results are bit-identical -- voltages, line flows and generator P/Q to 17 digits -- on case118 and the three PEGASE cases under DC_KLU and DC_SparseLU. The AC algorithms are untouched. Tests pass in the C++17, C++14 and Debug builds. LSGrid::fillBdc, the other setFromTriplets in the DC path (6.6% of a DC solve), is deliberately left alone: its entries are unsorted and heavily duplicated, which is the shape Eigen's assembler is already good at. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

s = E * conj(y * E) with y = -(p + i.q) / sn_mva is two complex-times-complex products, each followed by std::complex's NaN-recovery branch, on each of case9241pegase's 7327 shunts on every solve. Written out on real and imaginary parts the branches are gone. -264k instructions of the 24.4M LSGrid::compute_results spends on a solve (-1.1% of the post-processing), measured under callgrind on three cache-reuse solves of case9241pegase. Bit-identical: the recovery path only fires on a non-finite product, and (-1 * x) / s and -(x / s) agree exactly in IEEE 754. Checked against the previous build on every element result -- P/Q/V/theta/amps of lines, trafos, loads, gens, shunts, sgens, storages and hvdc, at 17 digits -- over case118 and the three PEGASE cases under NR, NRSing and FDPF with SparseLU and KLU. Tests pass in the C++17, C++14 and Debug builds. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

OneSideContainer::compute_results called them back to back on the same elements -- for loads, static gens, storages, shunts, generators and SVCs -- and everything before the last line of each was the same work: read the element's bus, map it through id_grid_to_solver, check that neither is the deactivated bus. Only the final assignment differed: Vm times the bus' nominal kV against Va times 180/pi. The walk, the two gathers and the two checks now happen once, and both results are written from them. -862k instructions of the 24.1M LSGrid::compute_results spends on a case9241pegase solve (-3.6%), and -2.9% on the function's wall time (min of 11 batches of 2000 calls, four runs, the two ranges not overlapping). Bit-identical, checked on every element result -- P/Q/V/theta/amps of lines, trafos, loads, gens, shunts, sgens, storages and hvdc, at 17 digits -- over case118 and the three PEGASE cases under NR, NRSing and FDPF with SparseLU and KLU. Tests pass in the C++17, C++14 and Debug builds. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

_get_amps read a = sqrt(p^2 + q^2) / (sqrt(3) . v) as four separate expressions: the sum of squares into a vector, a square root over that vector, a full copy of v, then a scan of the copy replacing the zeros. That is two full-length heap allocations and three extra passes over memory, on every one of the four calls a solve makes -- both ends of the powerlines and of the transformers. The guard that stops a disconnected element's zero voltage dividing by zero is now a ternary inside the single loop. -14% on LSGrid::compute_results' wall time (min of 11 batches of 2000 calls: 1164 -> 1001 microseconds, four runs) and -1.55M instructions of 23.3M (-6.7%). The gain is much larger than the instruction count suggests, because allocations cost time rather than instructions. Same arithmetic in the same order, so results are bit-identical: verified on every element result -- P/Q/V/theta/amps of lines, trafos, loads, gens, shunts, sgens, storages and hvdc, at 17 digits -- over case118 and the three PEGASE cases under NR, NRSing and FDPF with SparseLU and KLU. Tests pass in the C++17, C++14 and Debug builds. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Each of case9241pegase's 16049 branches costs six complex-times-complex products -- y11.Ehv, y12.Elv, y22.Elv, y21.Ehv, then Ehv.conj(I_hvlv) and Elv.conj(I_lvhv) -- and std::complex follows every one of them with a branch that re-derives the result if it came out NaN. That was a third (4.24M) of everything the branch flow loop did. -963k instructions of the 21.7M LSGrid::compute_results spends on a solve (-4.4%) and -2.9% on the function's wall time (min of 11 batches of 2000 calls, four runs, the two ranges not overlapping). Bit-identical: the products are grouped exactly as std::complex groups them, conj is an exact sign flip, and the recovery path only fires on a non-finite product. Verified on every element result -- P/Q/V/theta/amps of lines, trafos, loads, gens, shunts, sgens, storages and hvdc, at 17 digits -- over case118 and the three PEGASE cases under NR, NRSing and FDPF with SparseLU and KLU. Tests pass in the C++17, C++14 and Debug builds. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

fillYbus, fillSbus, fillBdc and hack_Sbus_for_dc_phase_shifter are reached only from LSGrid::_build_into_cache. The bus id they test comes out of the container's own bus_id_, for an element the loop has already established is connected; the solver id comes out of the id_grid_to_solver LSGrid built two steps earlier. A deactivated bus there is not a caller error, it is an inconsistency between an element's status and its bus -- which check_grid() validates and which no public method can produce on its own. Same reasoning that put _check_in_range behind NDEBUG, applied to the same kind of check. Each was a std::ostringstream built inline in an assembly loop, so they cost the surrounding code its registers as well as the compare. 18 blocks across 7 containers. -240,759 instructions on the branch fillYbus (-3.2%), -394,671 on LSGrid::fillYbus, -359,922 on a whole case9241pegase AC rebuild solve and -169,974 on a DC one. Too small to separate from run-to-run noise in wall clock; the instruction counts are exact and reproducible. Nothing a user can reach lost a check, and the assertion builds keep every one of them -- verified from the binaries: the Release and C++14 libraries carry none of the message strings, the Debug library carries them all. Results are bit-identical on case118 and the three PEGASE cases across 16 AC configurations (NR / NRSing / FDPF with SparseLU and KLU) and 8 DC ones, on every element result at 17 digits. Tests pass in the C++17, C++14 and Debug builds. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

The same treatment the assembly path just got, and worth considerably more here: nine blocks, all inside loops that run over every element. compute_results_tsc_rxha_no_amps has four, in a loop over all 16049 branches of case9241pegase; v_kv_theta_from_vpu two, over every load, static gen, storage, shunt, generator and SVC; ShuntContainer::_compute_results two; and _get_results_back_to_orig_nodes one, per bus. Each was an std::ostringstream built inline in the loop. Same argument as before: these run only from LSGrid::process_results, on a bus the container just read out of its own bus_id_ for an element already established as connected, and on the id_me_to_solver this very solve built. LSGrid::compute_results drops -1,591,221 instructions (-7.65%) and -7.4% of its wall time (min of 11 batches of 2000 calls: 954 -> 884 microseconds, four runs, the two ranges not overlapping). Left alone deliberately: the guard in HvdcLineContainer::compute_results that rejects a half-open droop line. That is a state invariant, not an id check, and it is the documented alternative to indexing id_grid_to_solver with the open side's -1. Verified from the binaries: the Release and C++14 libraries carry none of the nine messages, the Debug library carries them all. Results bit-identical on case118 and the three PEGASE cases across 16 AC configurations and 8 DC ones, every element result at 17 digits. Tests pass in the C++17, C++14 and Debug builds. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

One conflict, in TwoSidesContainer_rxh_A.hpp: 3ae9427 had switched reconnect_connected_buses to get_bus_side_*_internal, and #188 deletes that function outright. Resolved in #188's favour -- the function is gone, so the _internal edit to it is moot. The other 26 _internal call sites, the branch flow arithmetic on parts, and the NDEBUG gates in fillYbus and compute_results_tsc_rxha_no_amps all survive the merge untouched. 227 test cases / 590,769 assertions pass in the C++17, C++14 and Debug builds. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

_apply_and_track_buses brackets a mutation between "take this element's contribution away" and "put it back". The bus-id validation lived inside that bracket, in _generic_change_bus, so a call the grid was going to refuse was rejected with the contribution already removed -- which is why it needed a try / catch(...) putting it back on the way out. That catch cost far more than the path it protected. An unwind edge through GenericContainer.hpp made GCC keep every std::vector<bool> access in fillYbus live across it, in a function that never calls any of this. Bisecting #188's nine commits against that one number found it exactly: eight at 7,463,013 and the ninth at 8,859,264, and deleting only the catch while keeping everything else in that commit restored 7,463,013 to the instruction. Checking first is both simpler and cheaper. GenericContainer::_check_new_bus_id is called by the four change_bus entry points before they enter the bracket, so a refused call never touches the counts at all rather than touching them and undoing it with a restore that has to reason about a half-applied mutation. It is always active: the id comes from the caller. -3,899,955 instructions on a case9241pegase rebuild solve (-0.45% of everything), of which the whole of the branch fillYbus' share: 12,085,107 -> 8,185,194, -32.3%. pre_process_solver 33,503,851 -> 29,603,932. An exception from deeper inside a mutation can still leave the counts short. That is deliberate: such a grid must be rebuilt and its caches dropped, not carried on with. #188's own coverage caught the first attempt at this, where the check had landed inside the lambda instead of before the bracket on the HVDC and two-sided paths -- 36 failing assertions, all of them the refused-change_bus cases that commit added. All 227 test cases / 590,769 assertions pass in the C++17, C++14 and Debug builds; results bit-identical on case118 and the three PEGASE cases across 16 AC and 8 DC configurations. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Validating the bus id before the bracket removed the biggest unwind edge through GenericContainer.hpp, but not all of them. deactivate_no_bus_tracking, reactivate_no_bus_tracking and change_bus_no_bus_tracking each re-checked the element id, and _generic_deactivate / _generic_reactivate / _generic_change_bus checked it again underneath them -- three layers of the same check, only the outermost of which a user can reach. The inner two are now _check_in_range_internal, the debug-only form. Every one of the 18 _apply_and_track_buses call sites is reached through a public entry point that has already raised for a bad id, and a throw from inside the bracket is precisely what this layer exists to avoid. That takes the branch fillYbus back to exactly what it cost before #188 landed -- 7,222,254 instructions, the same figure to the digit. Over both commits: 12,085,107 -> 7,222,254 on that function (-40.2%), pre_process_solver 33,503,851 -> 28,640,988, and -4,862,863 on a whole case9241pegase rebuild solve (-0.56%). No user-facing check was lost, verified from a release build rather than by reading: it still raises for change_bus_load(0, -1), change_bus_load(0, nb+1000), the generator equivalents, change_bus_load(999999, 0) and deactivate_load(999999) -- six out of six. The Debug library carries all ten internal messages, the Release library the four user-facing ones. Two stale comments went with it: GeneratorContainer::_change_bus and SvcContainer::_change_bus both claimed their IndexError came from _generic_change_bus "which the caller runs *after* this function". change_bus_no_bus_tracking raises first, and has for some time. 227 test cases / 590,769 assertions pass in the C++17, C++14 and Debug builds -- including #188's refused-change_bus sweep, which is what makes this safe to do. Results bit-identical across 16 AC and 8 DC configurations on case118 and the three PEGASE cases. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

A solve rebuilds its family's cache in place -- the bus labelling, Ybus / Sbus, the PV / PQ split, the slack weights, the algorithm's factorization -- so a throw part way through `ac_pf` / `dc_pf` left a mixture of the old grid and the new one behind, with the change flags still saying whatever they said before the call. Nothing that reads those flags could tell a half-built cache from a whole one. So do not try to undo it. `ac_pf` / `dc_pf` now run the whole solve against a COPY of the change tracking and hold the grid itself at "everything changed, both families" for the duration; the copy is what `pre_process_solver` reads, what the algorithm is told, and what `process_results` marks as in sync. It becomes the grid's change tracking again only at a publication statement placed after `compute_results`, which a throw never reaches. An interrupted powerflow therefore leaves `need_reset_solver()` true on both families and the next one rebuilds from scratch -- by construction, whatever a future step throws and wherever, with no unwind path, no catch block and no restore code. Deliberately not a try / catch: an unwind edge through this code is not free (the one that used to guard the bus-counting bracket cost 4.9M instructions per solve without ever running). `DualAlgoControl` is 24 bools -- three register-sized stores each way. The input validation at the top of `ac_pf` / `dc_pf` (Vinit size, max_iter / tol, a droop grid handed to a solver that cannot do droop) stays OUTSIDE that bracket: it runs before anything is touched, so a rejected call still costs the caller's cache nothing. Two supporting changes: - `_build_into_cache` no longer recomputes `redo_all` after `init_bus_status()`. That existed because `init_bus_status()` could raise `has_dimension_changed()` half way through AND `solver_control` was bound to the grid's own member, so it saw it. Neither half is true since the bus connectivity became the per-bus element counts: `init_bus_status()` is `_ensure_bus_counts()`, logically const, deciding nothing about change flags -- and the powerflow now hands over a copy anyway. Dead code, and a stale comment that contradicted the new protocol. - `_build_foreign_cache` retires this grid's cache BEFORE copying the caller's labelling into it rather than after. Same end state (the retirement writes `built_for_nb_bus` and the flags, the publication writes seven containers, neither reads the other), but a throw part way through those allocating vector assignments can no longer leave a half-published mixture behind a control that still says the cache is up to date. `built_for_nb_bus` is deliberately NOT zeroed when the bracket opens: that would make the next solve's "nothing changed" claim be checked against a cache it has just been told is inconsistent, and the debug assertion in `_pre_process_own_cache` would fire on the ordinary success path. Left alone it is simply never consulted -- `need_reset_solver()` makes `_build_into_cache` redo every step, resets the algorithm, and short-circuits both that assertion and `unset_changes()`. Measured on case9241pegase / NR_KLU (callgrind, 3 solves): cache-reuse 619,932,575 -> 619,932,638 (+63 Ir total, +21 per solve) full rebuild 863,811,316 -> 863,811,268 (-48 Ir) Wall clock indistinguishable (min of 5: 29.911 vs 29.834 ms/solve). 16 grid x algorithm result dumps bit-identical. Covered by four new sections in test_cache_reuse.cpp: a solver that throws from compute_pf (AC and DC), a throw raised after the algorithm ran (the wrong-sized voltage rejection in process_results), and the other half of the protocol -- a solve that does not throw still publishes exactly what it consumed and leaves the other family alone. 228 test cases / 590,796 assertions pass in C++17 Release, C++14 Release and Debug (assertions active). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

The sections added with the working-copy protocol each swapped the solver between a good powerflow and a throwing one, so "the grid asks for a rebuild" was true partly because `change_algorithm` had said so. This one does not swap anything: a plugin that converges on its first call and throws on every later one, run twice on the same grid. Call 1 succeeds, so both families end it with a live cache and `nothing_changed()`. Call 2 throws. Whatever the grid says afterwards was said by the throw -- and it says `need_reset_solver()` on both families, the same state a refused mutation has to leave behind. Checked against the previous commit's library with the same probe source: there the grid still answered `nothing_changed()` on both families after the throw, and would have solved from the half-rebuilt cache. 228 test cases / 590,809 assertions pass in C++17 Release, C++14 Release and Debug (assertions active). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Everything a solve builds is derived from the elements, here and now -- except the per-bus element counts, which are carried forward +1 / -1 from every mutation the grid has ever seen. `init_bus_status()` only *established* them when they had never been counted; an armed-but-drifted count was inherited forever. That matters more than a stale cache would, because connectivity IS the counts: a phantom bus enlarges the solved system and shifts every bus id after it, and nothing downstream can notice -- an off-by-one count reads exactly like a real one. So the counts are now rebuilt from the elements whenever the answer to "may I reuse what was built for this grid?" is no: `force_full_rebuild` (the caller's reuse policy, which the batch entry points always pass) or `need_reset_solver()` -- what prevent_cache_reuse(), a copy, a set_state, a throw out of a powerflow and a throw out of a mutator the caller reported all leave behind. That is the one moment a drift can be repaired rather than inherited, and it is on the path already paying for a full rebuild of everything else. Hooked to "cannot reuse", NOT to `redo_all`: that also fires on has_dimension_changed(), which an ordinary topology change raises on every grid2op step, and an O(all elements) recount per step is exactly the cost these counts exist to avoid. case9241pegase, callgrind: cache-reuse path 619,932,638 -> 619,932,635 (-3 over three solves) full-rebuild path 863,811,268 -> 879,669,208 (+5.29M/solve, +1.8%) all of the delta attributed to recompute_bus_element_counts (15,857,772 of 15,857,940). 16/16 result dumps bit-identical. Four new cases stage a drift by hand -- no public API can produce one since the throw sites left the counting bracket -- and check it dies: a phantom bus (a lost decrement, which changes connectivity), a count one too low (a lost increment, which does not, and is latent until the bus empties one element early), the same for the DC family, and the negative case that an ordinary reusing solve does not pay for the recount. Checked against the previous commit's library: there the phantom survives the invalidated solve. 229 test cases / 590,856 assertions pass in C++17 Release, C++14 Release and Debug (assertions active). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

The two branches look like two operations and are one at two strengths: init_bus_status() is _ensure_bus_counts() -- recount only if the counts were never armed -- plus a debug assertion. The direct recount does strictly more, and the assertion is all it skips, which right afterwards cannot fail. Comment only. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

BDonnot

Some more refactoring.

Do not hesitate to A/B test and not just assume everything here needs to be implemented.

There are some things that change the architecture (like having both mismatch, the real one eg used as rhs for fdpf and nr, and the complex one that is used on every DC Algo at the algorithm level so that each can reuse it (and that it can be reused afterwards).
I would like ro have the cost (or benefit) from this change particularly (and if it does not cost too much to have it implemented)

BDonnot · 2026-09-04T10:54:16Z

+            // Everything a caller can get wrong -- the element id, the bus id -- is
+            // checked by the mutators BEFORE they call in here (_check_in_range and
+            // _check_new_bus_id), so a call the grid refuses never reaches this
+            // bracket and the counts it would have left short are never touched.
+            //
+            // This used to be a try / catch that put the contribution back on the way
+            // out of an exception. It was not free: an unwind edge through this header
+            // made GCC keep every std::vector<bool> access in fillYbus live across it,
+            // for 4.9M instructions per rebuild solve of case9241pegase, in a function
+            // that never calls any of this. An exception from deeper inside a mutation
+            // can still leave the counts short; that is a grid which must be rebuilt
+            // and its caches dropped, not one to carry on solving with.


The new defense mechanism should be explained here (something that throws put the grid in a state where everything needs to be computed)

BDonnot · 2026-09-04T10:57:22Z

+    // fail: recompute_bus_element_counts() ends with recount_connected_buses(), and
+    // `connected_bus_count_is_exact()` is that same "count the non-empty buses" loop
+    // compared against what it just wrote.
+    if (force_full_rebuild || solver_control.need_reset_solver()){


I'm not a huge fan.

need_reset_solver is not "need_reset_cache". Most of the time the cache is not poisoned so it forces its rebuild for nothing. Maybe the addition of another flag in solver_control above, like "cache_maybe_poisoned" or something would be clearer and avoid the unnecessary rebuild of a consistent cache when solver needs to be reset.

BDonnot · 2026-09-04T13:23:05Z

+    // Hence the order: zero the feature positions, let the dS pass overwrite
+    // the ones it shares, then add the feature values on top. Values are
+    // unchanged: `0 + x` is `x`.
+    for (int pos : feature_pos_)


There could even be two types of feature (type constexpr, known at compile time) the one that needs 0ing and the other. For example distributed slack does not need to be zeroed. Once set values does not change I think.

BDonnot · 2026-09-04T13:28:31Z

+template <typename... Rest>
+inline RealVect NRSystem<Base, Rest...>::mismatch() const
+{
+    RealVect res(static_cast<Eigen::Index>(total_state_variables()));


Mismatch could be cached in BaseAlgorithm. Fdpf reuse the same strategy.
The size of the mismatch is not the same here and in fdpf however which is fine. Derived class could be responsible for the mismatch size.
This Would avoid other system calls here I believe.

BDonnot · 2026-09-04T13:28:38Z

+template <typename... Rest>
+inline RealVect NRSystem<Base, Rest...>::mismatch() const
+{
+    RealVect res(static_cast<Eigen::Index>(total_state_variables()));


Mismatch could be cached in BaseAlgorithm. Fdpf reuse the same strategy.
The size of the mismatch is not the same here and in fdpf however which is fine. Derived class could be responsible for the mismatch size.
This Would avoid other system calls here I believe.

BDonnot · 2026-09-04T13:30:02Z

+inline void NRSystem<Base, Rest...>::_reconstruct_V_into(
+    CplxVect& V_out,
+    const Eigen::Ref<const RealVect>& Va,
+    const Eigen::Ref<const RealVect>& Vm)
+{
+    // V = Vm * (cos(Va) + i.sin(Va)), straight into the caller's vector: the
+    // same expression as before, but assigned instead of returned, so the
+    // per-call nb_bus complex temporary is gone (Eigen resizes V_out only when
+    // the dimension actually changed, so nothing is allocated after the first
+    // call of a topology).
    const cplx_type m_i = BaseConstants::my_i;
-    return Vm.array() * (Va.array().cos().template cast<cplx_type>()
-                         + m_i * Va.array().sin().template cast<cplx_type>());
+    V_out = Vm.array() * (Va.array().cos().template cast<cplx_type>()
+                          + m_i * Va.array().sin().template cast<cplx_type>());
 }


The fdpf formula might be reused here I believe. And it has been proven faster. Needs to be A/B tested though

BDonnot · 2026-09-04T13:32:31Z

+    // trial of a line search).
+    assert(res.size() == static_cast<Eigen::Index>(total_state_variables()));
+    ybus_v_cache_.noalias() = Ybus_ref_ * V_t;
+    mis_cache_ = V_t.array() * ybus_v_cache_.array().conjugate()


The more I read this the more I think both real and complex mismatch should be defined in BaseAlgorithm (as members) and filled by the derived classes.
So that it's even usable later on the code (eg for post processing)

…sets Review of #189: need_reset_solver is not need_reset_cache. It says the solver-side data is stale -- data this grid re-derives from the elements every time it rebuilds -- which says nothing about the per-bus element counts. Gating the recount on it charged an O(all elements) walk to every caller who merely wanted a fresh solve: `allow_cache_reuse(false)` on every single solve, a `tell_solver_need_reset()` after an ordinary change, every batch sweep. AlgoControl gains `cache_maybe_poisoned()`, raised by construction, by `tell_all_changed()` (a reset, a copy, a set_state, a powerflow that threw part way through) and by the new public `LSGrid::tell_bus_counts_maybe_poisoned()` -- what a caller who changed bus membership behind LSGrid's back, or who caught an exception out of a mutator, uses to say so. `prevent_cache_reuse()` / `tell_solver_need_reset()` deliberately do not raise it. The implication runs one way only: a poisoning claim raises need_reset_solver with it. The counts decide which buses exist, so everything solver-side is built on them; recounting while re-stamping the rest would repair the counts and then solve the old bus set anyway -- the same wrong grid with tidier bookkeeping. The converse is exactly what this commit stops assuming. case9241pegase, callgrind, against the previous commit: cache-reuse path 619,932,635 -> 619,932,650 (+15 over three solves) full-rebuild path 879,669,208 -> 863,811,432 (-15,857,776) which returns the full recount cost added two commits ago to every path that was paying it without needing it, while keeping the repair. 16/16 result dumps bit-identical. The four staged-drift cases now claim poisoning rather than `prevent_cache_reuse()`, and one of them pins the new distinction directly: an ordinary solve, then a `prevent_cache_reuse()` solve, keep the drift; only the poisoning claim repairs it. They also check the repaired count reaches the bus labelling, not just the counts -- without the implication above it did not, and the DC case segfaulted. CHANGELOG gains a TODO for promoting `canonicalise_vm_va` to BaseAlgo and for the "wrap the angle once at the end" question, both out of scope here. 229 test cases / 590,862 assertions pass in C++17 Release, C++14 Release and Debug (assertions active). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Review of #189: the Newton-Raphson and the FDPF do the same thing here and had grown two implementations that ended up doing different things. Both can drive a magnitude past zero mid-solve -- the FDPF Q iteration (`Vm_(pq) -= q_`), the NR step (`Vm_(bus) += dx(col)`) -- and both accumulate into Va_ without wrapping. The FDPF repaired it with a sign flip and a half turn, unconditionally, every call. NRSystem::apply_step repaired the same overshoot with `Vm_ = V_.abs(); Va_ = V_.arg();` -- the hypot and atan2 pair the cheap form exists to replace -- but only behind `if (Vm_.minCoeff() < 0)`, and wrapped the angle as an accidental side effect of atan2's range. So: one implementation in BaseConstants (where the constants it uses already live, and which NRSystem already includes -- no new coupling, and BaseAlgo inherits it so every caller keeps resolving), split into the two halves, which have different justifications: - fix_negative_vm MUST run per FDPF iteration, because `mis_ /= Vm_` follows it and a negative magnitude flips that bus's P and Q. It now runs behind the NR's guard, which is where the win is. - wrap_va need NOT: nothing inside either solve reads Va_ except its own accumulation and the cos / sin that rebuild V, both indifferent to a multiple of 2.pi. Once per solve, and asked for before it is done. Benchmarked on case9241pegase against five alternatives -- shared with no guard, guarded per call, guarded plus a wrap guard, and the wrap moved to the end with and without its own guard (Ir per three solves, cache reused): FDPF_XB FDPF_BX NR current (one each) 674,070,631 534,628,024 619,932,650 shared, no guard 673,949,686 534,533,467 624,557,387 vm guard, va every call 661,365,466 524,694,895 623,184,743 vm guard + va guard 640,398,586 508,302,607 621,812,459 vm guard, va at end 636,649,933 505,471,702 620,390,021 ... + va guard (kept) 636,268,666 -- 620,008,808 -5.6% on both FDPF flavours. The shared implementation on its own accounts for none of it -- it measures within noise of the old code -- the guard is the whole of it. Wall clock agrees: 6 to 11% off an FDPF solve on case1354 / case2869. The NR pays +25k instructions per solve (out of 206M) for something it did not have: a converged NR now reports an angle in [-pi, pi], like the FDPF. It never wrapped -- it inherited the effect from an atan2 that only fires on a trajectory heading for divergence. 16/16 result dumps bit-identical, both FDPF flavours included. 229 test cases / 590,862 assertions pass in C++17 Release, C++14 Release and Debug (assertions active). Drops the CHANGELOG TODO added last commit: both halves of it are done. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Seventh shape for the same job, from the review: one fused loop reading each element once and writing only the coefficients actually out of form -- Vm < 0, Va outside [-pi, pi] -- with no whole-vector select / cwiseAbs / round and no separate probe. It beats the old code and loses to the guard, consistently, on every grid and both flavours (Ir per three solves, cache reused): case2869 case9241 FDPF_XB current 103,837,561 674,070,631 guarded 98,151,160 636,268,666 coeffwise 100,858,705 654,109,921 FDPF_BX current 124,625,893 534,628,024 guarded 117,670,546 505,090,435 coeffwise 120,985,069 519,022,378 NR current 126,573,030 619,932,650 guarded 126,596,664 620,008,808 coeffwise 127,078,350 623,356,856 ~2.8% behind on both FDPF flavours, and worst of the three on the NR. Writing less is not the operative cost here: a scalar loop does not vectorise, and it reads Vm AND Va where the guard is a vectorised read-only reduction over Vm alone. Wall clock could not separate the two (2.923 vs 2.977 ms on case2869, against 3.104 for current), which is why the deterministic instrument decides. No code change -- src is identical to the previous commit. The rejected variant is deliberately not kept: it looks like an option and is slower than the one next to it. Recorded in the CHANGELOG so it is not re-proposed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Review of #189: both real and complex mismatch in BaseAlgo, filled by the derived classes, so they are readable afterwards (post-processing). Done for the complex per-bus mismatch and the Ybus * V scratch, which are the same concept and the same shape (nb_bus) in both families. NRAlgo owns the NRSystem rather than being one, so it hands the system the addresses of its own BaseAlgo members -- the arrow keeps pointing NRAlgo -> NRSystem and nothing had to learn about the algorithm in the other direction. Three things this turned up: - The two buffers did not hold the same quantity. The FDPF divided the mismatch by Vm IN PLACE, so it ended a solve holding mismatch/Vm where the NR held the raw mismatch: sharing the member as-is would have made get_bus_mismatch() mean two things depending on who ran last. The division is now out of place, and that extra vector is the whole cost of this commit. - Dividing only the rows actually extracted -- mis(pvpq).real() / Vm(pvpq) -- avoids that vector and does strictly less arithmetic. It measured 2.2% SLOWER on both flavours: the gather does not vectorise, the contiguous divide does. Same lesson as the coefficient-wise canonicalisation. - NRSystem is used standalone (test_powerflow_algorithm builds one and calls mismatch() with no NRAlgo), so a null buffer pointer segfaulted. It falls back to its own buffers when unclaimed, and NRAlgo re-points them every compute_pf rather than once, so a copied algorithm cannot write into the original's members. Ir per three solves, case9241pegase, cache reused: FDPF_XB 636,268,666 -> 637,137,124 (+0.14%) FDPF_BX 505,090,435 -> 505,766,473 (+0.13%) NR 620,008,808 -> 620,005,733 (neutral) This is the one commit on the branch that is not faster or free. It buys a mismatch that is exposed and means one thing in both families; the alternative is to keep the in-place divide and accept two meanings behind one accessor. Self-contained, so it reverts on its own if that trade is not wanted. The REAL residual is deliberately not merged: the FDPF's is two vectors (p_ over pvpq, q_ over pq, each solved in place by its own linear solver) while the NR's is one over total_state_variables including custom rows. They differ in shape, not just size, and merging them would cost a copy per iteration. 16/16 result dumps bit-identical. 229 test cases / 590,862 assertions pass in C++17 Release, C++14 Release and Debug (assertions active). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

Four of the five small items from the review of #189: - GenericContainer::_apply_and_track_buses now explains the defence, not just its consequence: an exception from inside a mutation leaves the counts short, that grid must be rebuilt, and the repair is LSGrid::tell_bus_counts_maybe_poisoned() for a mutator or the working-copy protocol in ac_pf / dc_pf for the powerflow -- there is none in the bracket, deliberately, because putting the contribution back on the way out is what the try / catch cost. - GaussSeidelAlgo's diagonal scan breaks at the diagonal. Ybus is column-major and compressed, so row indices ascend within a column; it was reading the whole lower triangle of every column for nothing. - Contrib is constexpr and noexcept throughout -- constructor, structural(), and every accessor including Eigen's triplet protocol. C++14 is enough, confirmed by the C++14 job. - NRAlgo's timer_Va_Vm_ times apply_step and nothing else. It used to start a statement earlier and charge the scaling of the step vector to it; applying the coefficient is the second half of the scaling policy and now counts against timer_scale_. No arithmetic moves, only the attribution. The fifth is NOT done, and the reason is now in the code. Gating refresh_ybus_cache on the change flags -- the O(nnz) transpose it does on every single call -- is unsafe: the batch sweeps mutate Ybus in place between solves (YbusPolicy::Contingency::remove_from_Ybus does `Ybus.coeffRef(i, j) -= value`) and then hand the algorithm a control saying tell_none_changed from the second step on. A gated cache answers from a stale copy against that pattern: reproduced directly -- solve, edit one coefficient, solve again with a "nothing changed" control -- for 0.38 pu of voltage error on case118, silently, where the ungated code gives exactly 0. nnz does not catch it, an in-place `-=` leaves the sparsity pattern alone. The DC family is immune because it holds Ybus internally and is handed each edit through update_internal_Ybus, which is DC-only; gating the AC side would need the same hook. 16/16 result dumps bit-identical. Neutral on instructions (-36 per three solves, case9241pegase). 229 test cases / 590,862 assertions pass in C++17 Release, C++14 Release and Debug (assertions active). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

None of the MATPOWER benchmark cases carries an hvdc line, an SVC or a distributed slack -- case118, case1354pegase, case2869pegase and case9241pegase are all 0 hvdc / 0 SVC / 1 slack generator -- so NRSystem's extensions contribute nothing on them and anything touching that code measures as exactly zero. benchmarks/make_exotic_grid.cpp turns any ".lsb" into the same grid with all of them switched on, as step 1b of the existing profile pipeline. On case9241pegase: 1445 slack generators (weight |target_p|), 5 groups of 4 electrically close generators jointly regulating one bus, 5 voltage-mode SVCs, 4 hvdc lines covering VSC/VSC, VSC/LCC and LCC/LCC with two on angle droop. J goes from 17037x17037 / 129423 nnz to 17087x17087 / 131285, converging in 5 iterations. Two deliberate stress cases: the two droop lines share an end bus, so their dP/dtheta entries collide on one (P row, theta col) position, and the voltage-regulating VSC stations sit on the generator groups' controlled buses. Two constraints found while building it, both enforced by fill_voltage_control_solver_data and both documented in the file: - every controller of a bus must ask for the same voltage; - an SVC must be ALONE in its control group ("not supported in v1"), so an SVC co-regulating a bus with generators is not expressible today. The SVCs are placed on generator-free buses instead. With that grid, the per-feature "needs zeroing" split from the review was measured and is NOT being done: the loop costs 161,550 instructions per three solves out of 699,756,287 -- 0.023% of the program, 0.74% of fill_J -- and realising any of it means reworking the zero -> dS assign -> feature add ordering that makes the shared-bus droop case come out right. The inventory behind it still holds and is in the CHANGELOG: only the hvdc droop's values change between iterations; MultiSlack's weights and every VoltageControl coefficient are fixed for the solve. Benchmark tooling only -- no library code changes. Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

The first version was not convincing: 20 remote-controlled generators where a TSO snapshot has essentially all of them, and 4 hvdc lines where Europe has plenty of cross-border links. Rebuilt around three rules that came out of getting it to solve at all: - every setpoint is read off a reference solve of the UNMODIFIED grid. A setpoint is a property of the bus, not of whoever controls it, so several controllers of one bus agree by construction -- and inventing a target per generator (its own target_vm_pu, belonging to a different bus) diverges: SolverFactor after ~20 iterations, or 300 iterations without converging. - each generator regulates one bus a single branch away, that bus' own base voltage as target. - applied in VERIFIED chunks, last, with the slack, hvdc and SVCs already in place. 895 of the 1445 generator / bus pairs make a system the NR cannot solve, and it is the individual pairs rather than the count: taking 20 controllers at a time along the generator list, seven windows out of eight converge in 6 iterations and one does not. Rather than guess a criterion for the bad ones the tool checks, and a chunk validated against a grid without the other features proves nothing about the grid that gets saved -- which is how the previous version shipped a file that did not solve. case9241pegase now gives 1445 slack generators, 550 remote controllers over 486 controlled buses, 5 SVCs and 20 hvdc lines (10 on droop), J 18330x18330 with 141728 nonzeros against 17037x17037 / 129423 plain, solving in 19 iterations. The tool refuses to save a grid that does not solve. Re-measured, the feature-zeroing split is still not worth doing: old grid this grid program, 3 solves 699,756,287 3,577,650,967 fill_J 21,744,624 86,991,627 zeroing loop 161,550 1,176,480 share of program 0.023% 0.033% share of fill_J 0.74% 1.35% Seven times the absolute cost and the same conclusion. Benchmark tooling only -- no library code changes. Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

benchmark_binary_powerflow.py's usage now points at step 1b, so the grid that exercises NRSystem's extensions is discoverable from the script that would otherwise only ever be given feature-free MATPOWER cases. TODO entry for what building it uncovered: remote voltage control fails on 895 of the 1445 generator / bus pairs of case9241pegase, individually rather than by weight of numbers, on the first iteration, with InifiniteValue or SolverFactor. The base solution satisfies KCL at every bus and each controlled bus is asked for the magnitude it already holds, so a solution exists; feasible setpoints and forbidding control cycles both leave it unchanged. That points at the linear solver or the Jacobian assembly, not at the data. Also worth recording from the same profile, for whoever optimises next: on case9241_exotic KLU is 83.6% of the run (klu_refactor alone 71.6%, klu_scale 5.3%) and lightsim2grid's own code 16.4%; on the plain case it is 73.8% / 26.2%. Docs and CHANGELOG only. Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com>

fix_negative_vm and wrap_va were each called behind an `if` written out at every call site: `if (Vm_.minCoeff() < 0)` in BaseFDPFAlgo::has_converged and in NRSystem::apply_step, `if (va_out_of_range(Va_))` at the end of both compute_pf. Two call sites each, asking the same question, which is how the two families grew two different repairs in the first place -- and the guard, not the shared implementation, is where the whole 5.6% of that change came from, so it should not be something a caller has to remember. Both guards now live inside the function: fix_negative_vm returns immediately when no magnitude is negative, wrap_va when every angle is already in range. canonicalise_vm_va composes the two exactly as before, so the identity test_fdpf_algorithm.cpp checks against the abs()/arg() pair is unchanged (its batch spans negative magnitudes and angles several turns out of range, so both bodies still run there). Behaviour-neutral and free, measured on case9241pegase (3 solves, noreset, callgrind Ir), against the parent commit: FDPF_XB 637,137,124 -> 637,144,075 (+0.001%) FDPF_BX 505,766,473 -> 505,771,525 (+0.001%) NR 620,005,697 -> 620,006,654 (+0.0002%) Element results bit-identical across 20 configurations (case118, case1354pegase, case2869pegase, case9241pegase x NR_KLU, NR_SparseLU, NRSing_KLU, FDPF_XB_KLU, FDPF_BX_KLU). C++ suite: 590862 assertions in 229 test cases, all passing. Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com> Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy

The job died 27 seconds in, at `apt-get update`, before cmake was even installed and long before any of this project's code was compiled. silkeh/clang:11 is Debian buster based, and buster's repositories have moved to archive.debian.org, so the stock sources 404. The two sed lines that fix it were already in the file, commented out since the job was written; compile_gcc_earliest (gcc:8, also buster) has had exactly the same pair active all along and passes. Uncommented here, plus a `cat` of the rewritten sources so the next person can see what apt was pointed at. The third commented line stays commented: it is stretch-specific. Nothing else on this commit is red: compile_gcc_earliest, compile_clang_prev_latest, the C++14 (oldest supported) and C++26 builds, Catch2 + valgrind, the Eigen / libstdc++ assertion build, the standalone cmake builds on all three OSes and every wheel job that has finished are green, and so is DCO. Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com> Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy

get_bus_mismatch was the odd one out: get_Va, get_Vm and get_V, three hundred lines above it in the same class, all return Eigen::Ref<const ...>, as do the read-only vector getters on LSGrid. There was no reason for this one to hand back the storage type instead. Nothing is paid for it. Ref<const CplxVect> built from a contiguous CplxVect lvalue is a pointer and a size, constructed in place; the copy path a Ref<const> can take only fires for a non-matching stride or an expression, which a plain member vector never is. What it buys is that the getter no longer pins the mismatch to being an owned CplxVect. NRAlgo already hands this member's address to the NRSystem it owns, which keeps a fallback buffer of its own; if the mismatch ever becomes a segment or a map of something larger, a const CplxVect& return could only follow it with a copy. No caller exists yet (the getter is new on this branch and is not bound to python), and everything a caller would do -- indexing, size, .array(), binding to a const Eigen::Ref<const CplxVect>& parameter -- is unchanged. The one visible difference is the same one get_V already has: `auto m = get_...()` is now a view rather than a copy. C++ suite: 590862 assertions in 229 test cases, all passing; clang -std=c++14 and -std=c++17 syntax-clean. Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com> Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy

The previous commit repointed this job's apt at archive.debian.org, which was half right. The build log shows the image is bullseye, not buster, and shows where the rewrite lands: archive.debian.org/debian bullseye InRelease 200 (116 kB) archive.debian.org/debian bullseye-updates InRelease 200 archive.debian.org/debian-security bullseye-security InRelease Ign, Release 404 E: The repository 'http://archive.debian.org/debian-security bullseye-security Release' does not have a Release file. bullseye-security is not published on archive.debian.org at all, so rewriting that line only moved the 404. compile_gcc_earliest survives the identical rewrite because buster's security updates ARE there, under `buster/updates`. Dropped instead of repointed. Nothing this job installs comes from the security suite -- cmake, python3, git and pip are all in main, which the rewrite does reach -- and one unreachable entry is enough to fail the whole `apt-get update`. Verified the two seds on a representative bullseye sources.list: what survives is exactly the two suites the log fetched successfully. Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com> Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy

…raced step Repointing the stock sources at archive.debian.org got bullseye/main and bullseye-updates (both 200 in the build log) but not bullseye-security, which is not on the archive; dropping that entry moved the failure somewhere else again, earlier in the job than before. Three failures in, the step-per-command layout is the problem: it says which command exited non-zero and never what apt was actually looking at. So the bootstrap is now one `set -x` block that dumps /etc/os-release, the image's own sources.list and sources.list.d before touching anything, then writes the sources it wants rather than patching whatever is there: deb http://archive.debian.org/debian bullseye main deb http://archive.debian.org/debian bullseye-updates main `main` only, because that is all the image's sources.list had -- the earlier log fetched bullseye/main and nothing else -- and a component the archive does not carry would fail the update on its own. The clang toolchain in /etc/apt/sources.list.d/ is left alone: it is on apt.llvm.org, which is up and answered 200 throughout. Also sets Acquire::Check-Valid-Until false, as apt config rather than an -o flag, so the later `apt` and `apt-get update -y` steps inherit it: an archived suite's Release file is eventually past its Valid-Until, and apt treats that as fatal. `apt-cache policy cmake` closes the block, so the log states whether the first package the job installs is actually reachable. CI-only change; no library code touched. Signed-off-by: DONNOT Benjamin <benjamin.donnot@rte-france.com> Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nEocx4vydkgBNQKNzRfwy

BDonnot and others added 22 commits August 30, 2026 20:07

BDonnot force-pushed the claude/eigen-powerflow-efficiency-wn28tf branch from ad3be14 to b81d10c Compare September 4, 2026 08:31

BDonnot commented Sep 4, 2026

View reviewed changes

BDonnot and others added 2 commits September 4, 2026 15:52

BDonnot force-pushed the claude/eigen-powerflow-efficiency-wn28tf branch from 1184830 to 6f90fb5 Compare September 4, 2026 17:21

BDonnot and others added 9 commits September 4, 2026 18:55

BDonnot merged commit 47af18e into dev_1.0.1 Sep 5, 2026
58 checks passed

BDonnot mentioned this pull request Sep 6, 2026

Derive the voltage-control plan once per powerflow, in one class #191

Merged

Navigation Menu

Uh oh!

Powerflow performance pass, and a cache that cannot outlive a failure - #189

Powerflow performance pass, and a cache that cannot outlive a failure#189
BDonnot merged 35 commits into
dev_1.0.1from
claude/eigen-powerflow-efficiency-wn28tf

BDonnot commented Sep 4, 2026

Uh oh!

BDonnot left a comment

Uh oh!

Uh oh!

BDonnot Sep 4, 2026

Uh oh!

BDonnot Sep 4, 2026

Uh oh!

Uh oh!

Uh oh!

BDonnot Sep 4, 2026

Uh oh!

BDonnot Sep 4, 2026

Uh oh!

BDonnot Sep 4, 2026

Uh oh!

BDonnot Sep 4, 2026

Uh oh!

BDonnot Sep 4, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Uh oh!

Conversation

BDonnot commented Sep 4, 2026

Performance

Internal vs user-facing bounds checks

Correctness

A powerflow that throws can no longer leave a cache claiming to be in sync

The per-bus element counts are rebuilt when the cache cannot be reused

Tests

Notes for review

Uh oh!

BDonnot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

BDonnot Sep 4, 2026

Choose a reason for hiding this comment

Uh oh!

BDonnot Sep 4, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

BDonnot Sep 4, 2026

Choose a reason for hiding this comment

Uh oh!

BDonnot Sep 4, 2026

Choose a reason for hiding this comment

Uh oh!

BDonnot Sep 4, 2026

Choose a reason for hiding this comment

Uh oh!

BDonnot Sep 4, 2026

Choose a reason for hiding this comment

Uh oh!

BDonnot Sep 4, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant