Structural analysis: infer the coupling graph from the IR (separable dims, linking constraints) #39

Where this comes from

Prior-art scan of declarative optimisation languages surfaced GBOML (U. Liège, JOSS 2022) — a stand-alone MILP language whose organising abstraction is a hierarchical hypergraph: #NODE blocks own internal variables and local constraints, #HYPEREDGE blocks carry the coupling constraints, and coupling may only reference the external: variables of incident nodes. Their JOSS paper diagnoses a gap we recognise — algebraic modeling languages "rarely provide language constructs to represent the block structure in a model and typically fail to take advantage of it", while framework-style tools (they name PyPSA, Calliope, Balmorel) "lack expressiveness and adding components is often cumbersome". They claim four payoffs from carrying that structure: simpler encoding, model re-use, parallelised model generation, and structure-exploiting solvers.

The observation that motivates this issue: we have the same structural information, and we do not have to ask the user for it. GBOML makes the modeller draw the hypergraph by hand. In our IR the coupling is already marked syntactically, and can be recovered by a static pass over Program.

This issue is only the analysis-and-reporting half. Consuming the structure (block-annotated sinks, decomposition solvers) is explicitly out of scope — see the bottom.

The claim

relational/plan.py Program is a declaration-level object: a handful of ConstraintDeclarations, each with its dims and an affine Expression tree. Nothing is flattened, so the coupling pattern survives — and the expression nodes that couple coords are exactly the ones ARCHITECTURE.md already classifies for executor locality:

From which the rule: a dim that appears in some dims: but never inside a Sum.over, a GroupSum/At, or a Translate.dimension anywhere in the program is separable — the model is block-diagonal along it. This is sound by construction (every coupling route is one of the nodes above) and conservative in the right direction: it can under-report separability, never over-report it.

Caveat worth encoding in the output: Sum and Translate are statically decidable, but GroupSum/At coupling is data-dependent. Whether gen_bus links everything into one component or partitions into five islands lives in the mapping table. So the pass has two tiers — a static verdict from the Program alone, and an optional refinement that runs connected components over the mapping table(s). The refinement is a query over one small parameter, not over the constraint matrix.

What to build

1. A pure function over Program — relational/structure.py, no engine imports, mirroring plan.py's "frozen dataclasses, no execution logic" discipline — returning:
   • separable dims (block-diagonal axes) with the block count where statically known,
   • linking constraint declarations, and which dim each one collapses,
   • time-coupling bandwidth per dim (max |by| over Translate nodes on that dim),
   • the declaration-level bipartite incidence: which constraint families touch which variable families.
2. Optional data-tier refinement for GroupSum: connected components over the mapping parameter, giving the real block count.
3. Surface it. check() is the natural home — it is already the "tell me about this file without building it" entry point, and this needs no data. Overlaps with Diagnostics: a slow build cannot say where the time went #709 (build/solve observability, refiled from Build/solve observability: executor stats #34); structure is a pre-build property, so it should be reportable from check() alone and echoed in the diagnostics when a build happens.

Output should read as a diagnostic a user acts on, e.g. "block-diagonal along region (14 blocks); linking constraints: co2_budget (sums region), transmission_balance (group_sum line_region → region); snapshot coupled with bandwidth 1 (storage_balance)."

Why this is worth it independent of decomposition

Diagnostics alone justify it — knowing your 14 regions are independent except for one CO₂ budget row is a real modelling fact, and today nothing tells you.

But the reason to build it now is that separability is the prerequisite for sequential per-block build-and-solve. If a dim is provably separable, the engine can build and solve one slice at a time, holding a single block in memory. That is an out-of-core scaling story that falls straight out of hard rule 4 rather than straining it, needs no external solver, and is a strictly better use of the structure than exporting it.

Explicitly out of scope

• Block-annotated sinks (SMPS-with-blocks, DSP format). Structurally this is just another sink and the COO batches are already there — but it bets on a small, mostly academic structure-exploiting-solver ecosystem. Separate issue if ever.
• Implementing decomposition (Benders, Dantzig-Wolfe). These are solve algorithms — cut management, convergence, subproblem lifecycle — requiring a Python driver loop around highspy. That is Python-side algorithmic machinery of exactly the kind hard rule 6 exists to keep out. Not a goal.

Definition of done

check() reports the coupling structure of a model with no data bound; tests cover a block-diagonal program, a program with one linking Sum, a Translate-banded program, and a GroupSum program where the static verdict says "coupled" and the data-tier refinement finds N components. ARCHITECTURE.md's locality section notes that the same classification also yields the coupling graph.

🤖 Generated with Claude Code