Text-based (YAML) math definition interface

It would be useful to be able to define all the mathematics of a linopy model in human-readable text format. This would help with transparency and shareability (if that's a word) of math, as well as making it easier to define new math without needing to be so familiar with Python.

I have developed something that I think would be suitable in the Calliope project. It is a vectorised syntax that is documented. The main ideas are:

it works with xarray dataarrays, using a foreach directive to define the array dimensions and a where directive to mask the array so that variables/exprs/constraints can be defined only for a subset.
Math expressions do not reference dimensions unless slicing specific elements in the math. That is, you could define p <= p_nom and it would be applied automatically to all entries in the dimensions defined in foreach, or you could define p[snapshot=2012-01-01] <= p_nom * 0.8 and it would slice p to a given snapshot.
helper functions exist to deal with more complex cases. An obvious use-case is to roll the timeseries dimension forward/backward on one element so that you can get the equivalent of state_of_charge[t] = state_of_charge[t-1] + p_store[t] - p_discharge[t] by calling state_of_charge = roll(state_of_charge, snapshot=1) + p_store - p_discharge. These are designed as plugins so that power users can add and use new ones if the need arises.
It can automatically be rendered as LaTeX math for use in documentation. This means you only document your math once, when defining it in YAML, and it can then be kept up-to-date in the docs automatically.

It may not be fit for all linopy purposes but the fact that it also relies on multi-dim xarrays I think is already a good starting point. I have been exploring how it could also be used for tidy dataframes (if e.g. updating to using pyoframe) and it seems relatively straightforward; I just need to better understand how to apply the where string masking on polars dfs.

To mitigate bad actors injecting code into the expressions, all strings are parsed with pyparsing (ref. where and expression strings).

Examples of use for user-defined math can be seen here.

That would be a wonderful addition especially for making the custom constraints interface of PyPSA declarative (and thus persistable within a saved network) rather than providing only a callback function.

Do you have any thoughts on how to bring in external data (f.ex. an external csv file containing some co2 limits)?

Do you have any thoughts on how to bring in external data (f.ex. an external csv file containing some co2 limits)?

By the data having been loaded in as a parameter that can be called upon or (less ideal) by hardcoding the numbers directly in the YAML.

I really like this.
I think we should probably wait unit we merge #591 and #582, but i sketched sth up here: https://github.com/FBumann/linopy-yaml

Im open for feedback and would like to hear your thoughts

I've opened a dev branch in linopy for this (based on my existing implementation in Calliope): https://github.com/PyPSA/linopy/tree/feature/declarative-yaml-interface

Main things:

Added persistent (linear) expressions that are attached to the model object
Parsers for masking string and equation expression strings. The syntax is opinionated and is open for discussion (do we go for Excel-link nestable functions or pythonic dot notation method stacking?).
Expectation that the user provides a single input dataset with all parameters and groupers or masking arrays (named "lookups")
Can evaluate to a linopy model or to latex math (for auto-documenting an optimisation problem), as we do in Calliope. Latex docs generator not yet implemented here

A working input dataset (based on PyPSA model.energy example) and some YAML math to go along with it: Archive.zip

Things we might want to add, based on offline discussion:

key to flag math components for debugging (would then drop into a debugger in an interactive session)
Enable a dict of input datasets and allow each math component to be based on data in one of the those datasets (e.g. they might each have a different name dimension)

@brynpickering

@brynpickering I argue that having this as a standalone package, that plugs into linopy might be better architecturally.
It allows having breaking changes in the yaml layer until it matures.

Are you interested in discussing your approach compared to mine? (linopy-yaml) https://github.com/FBumann/linopy-yaml

I promped claude Fable 5 to review for architectual blockers of making this a standalone. Thought i should share the findings.

Note

The content below was generated by AI (checked against linopy master).

Does a standalone package block anything? Checked against linopy's actual extension surface — no hard blocks, two clean-integration gaps.

Everything hard about the YAML layer runs on the public API and needs nothing from a standalone package: building variables/constraints/objectives/expressions (add_*), LaTeX rendering (pure), the parser (slicing, sub-expressions), and the helper plugin registry.

Two integration points can't be done cleanly today — neither requires the parser to live in-tree; both are small, generic core hooks that help any extension:

No accessor/plugin registration. linopy has no register_*_accessor / entry_points group, so a model.yaml namespace has to monkeypatch today. Brittle (breaks on refactors, collides across packages). Fix: a tiny public registration point (or an entry_points group).
In-file persistence of the math (the "persistable within a saved network" motivation). linopy persists model-level state via a hardcoded known-key convention in ds.attrs, and read_netcdf drops attrs it doesn't recognize. A standalone package can round-trip by wrapping to_netcdf/read_netcdf — so it's not blocked — but it'll never be part of the format guarantee. Fix: reserve a namespaced attr and have read_netcdf preserve unknown namespaced attrs.

Proposal: standalone package + a small "extensibility" PR to core exposing those two hooks. This yields everything people want from in-tree — a clean model.yaml and guaranteed persistence — while keeping breaking-change freedom for the pre-1.0 YAML schema and keeping ~3–6k lines of parser out of core. Both hooks are defensible on their own merits, independent of this feature.

What this buys core linopy (independent of the YAML layer): mostly it's what core doesn't have to carry.

No SemVer coupling. The YAML schema will churn hard pre-1.0. Behind a plugin boundary that churn — and its breaking changes — never touches linopy's public API or forces a linopy major bump; core releases on its own cadence.
Less maintenance, smaller surface. ~3–6k lines of pyparsing grammar, LaTeX rendering and helper functions live and are tested outside core: no new core dependency, smaller install, less CI, smaller security-review surface, and no support load from monkeypatched installs.
A stable format instead of a lossy one. The reserved-attrs change stops read_netcdf silently dropping third-party metadata, and stops core accreting a bespoke known-key per feature (_relaxed_registry, _piecewise_formulations, …) — one generic mechanism instead of N special cases to keep maintaining.

Net: the hooks are a one-time, generic addition; everything volatile stays outside. Core gets an extension ecosystem (PyPSA custom constraints included) without taking on the maintenance or versioning burden of any of it.

Code paths this is based on (linopy master)

No registration hook: no register_*_accessor, entry_points, or plugin mechanism in linopy/*.py or pyproject.toml.
Persistence convention: to_netcdf (linopy/io.py) writes extra model state as hardcoded ds.attrs keys (_relaxed_registry, _piecewise_formulations, NETCDF_VERSION_ATTR). read_netcdf rebuilds a fresh Model and only restores recognized keys (for k in m.scalar_attrs: if k in ds.attrs, plus explicit _relaxed_registry / _piecewise_formulations handling) — unknown attrs are dropped. So an external _yaml_math attr survives to disk but isn't reattached unless the package also wraps read.

@coroa

Claude has some valid points but I don't agree with the general conclusion. I think the text math functionality would be safer living within the core repository for several reasons:

avoid API drift. Separate repo means hard-linking the linopy version a given YAML math is compatible with since old API and possible breaking new API both would have to be mitigated against.
more visibility & adoption. I'd like to see this as a key point of entry and a separate repo + package always makes it harder to get these things seen and adopted. This includes having a unified docs.
easier maintenance. Likely to be the same team maintaining this as maintaining core linopy and the overhead to ensure the two are aligned and to maintain separate issue triage, CI, etc. doesn't seem worthwhile.

Your separate repo is a lightweight version of what I started implementing in Calliope some years ago and it misses some key functionality and design decisions that would be necessary to make all features I mention above possible. Pyparsing is much more powerful than AST for the mask and expression string parsing, giving us much more control of the syntax we create (we will likely need some funky special character uses like in languages like GAMS to make this whole thing work). And helper functions need to be more strictly managed to work in all cases that they might be evaluated from in the expression or mask strings. I opted for a plugin approach using an ABC as it ensures I have a stable API on all registered helpers that I can then leverage in various other places.

I plan to pass Claude over my branch to tidy it up since it has several rough edges following my porting across from Calliope. I can then open a PR but would also like to continue discussion here with @coroa and @FabianHofmann first, to make a final decision on the best approach for where the implementation should live.

@brynpickering

@brynpickering Just to make this clear: Im perfectly fine discarding my repo. I agree that your approach with years of experience and hardening is for sure superior. Im merely interested in the architectural discussion about this being a standalone package.
For me, having a separate Issue tracker AND separate docs actually seems like an advantage. It keeps each page focused on the layer. And the yaml layer is pretty much its own language. So a separate page seems justified.

Having this in the core linopy docs would probably double the surface, and make for a much harder read.

Having it installable via pip install "linopy[yaml]" should of course be possible.

@brynpickering

I have one case where I think a decision in my approach may be worth adopting @brynpickering. I decided to declare the expected dims of parameters in the yaml. This is then evaluated upfront, and documents the dimensionality of the model in the yaml file. Maybe you have good reasons to not do this, but I think we should consider doing this.

parameters:
   p_max: {dims: [generator]}
   load:  {dims: [snapshot]}

@FBumann

Thanks @FBumann. I expected we would need to iterate over these design decisions so it's good to discuss them here. We'll need @coroa and @FabianHofmann to weigh in on where it lives as we're unlikely to reach consensus just between the two of us!

Not defining the parameter dimensions in advance was a conscious design decision. Frequently, they are dynamically dimensioned in energy system models. For example, load could be static (not indexed over snapshot), efficiency could be static, per-snapshot, per-node, or per-node + per-snapshot. Since the xarray/linopy broadcasting handles it automatically (both dispatch_{n, t} * efficiency_{n} and dispatch_{n, t} * efficiency_{n, t} are fine), it isn't necessary or indeed helpful to specify it in advance.

Saying that, perhaps listing the allowed (or excluded) dimensions could be an option, to at least guide someone reading the math as to what to expect to be working with.

@brynpickering

@brynpickering Yes, of course we'll fait for them.
I think having validation would make the yaml much more robust. But it should be optional.
The validation should have in mind that the parameters are meant for linopy, so subsetting dims is always ok (clean broadcast).

Yeah, having them optional seems reasonable, so that a user can give more context with which linopy can validate. It would be the max dims that a parameter can be indexed over, so that inputs can be a subset of that and linopy is OK with this. If not defined, the full set of dimensions is assumed?

Where I'd like to not be overly constrained (as a modeller) is when I add a new dimension to a problem (e.g. "scenario") and I don't want to have the overhead of adding that dimension to all my parameter dims defs since it's likely there will be a lot of parameters.

Yes, id give the max dims. Subsetting a dim is always possible in linopy. I would also say the parameters should be broadcast to those dims, so all further operations (reductions, selections, ...) can rely on that contract.

Im not sure how we should handle your "scenario" dim though. If it applies to ALL variables, i think a global param outside of the yaml would work as a modifier. If its only for some, we would need to switch of validation or some other workaround.

This comes down to what the yaml files should be for: recipies for certain components/patterns (no dims), or as full model specs (with dims).

I wouldn't broadcast parameters to those dims at the start as that would bloat in-memory footprint. Instead, I would broadcast them just before applying arithmetic so that only those few params that need broadcast exist in that bloated state at any given time. It would mean that we lose that reliable contract, of course, but it will likely be necessary so long as we have xarray as our data layer.

What I mean about my scenario dim is that I might have a base recipe that I want to extend with a partial YAML that extends the foreach of some variables and constraints to include a scenario dimension (e.g. for 2-stage stochastic in energy modelling, it might be only added to the timeseries variables). The number of lists that need extending to include this extra dimension is relatively small, but if I have to also add scenario to all parameters as well it then becomes quite a verbose addition that is more prone to human error. For example, there's 86 params+lookups vs 41 constraints+variables in the Calliope base math.

The way I had imagined it is that you might have a library of recipes that you can apply patches to for a given use-case, e.g.:

base recipe

dimensions: 
  name: ...
  time: ...
parameters:
  foo: {dims: [name], ...}
  bar: {dims: [name, time], ...}
variables:
  capacity:
    foreach: [name]
  dispatch:
    foreach: [name, time]
constraints:
  constrain_dispatch:
    foreach: [name, time]
    ...
  global_constraint:
    foreach: []
    ...
  constrain_capacity:
    foreach: [name]
    ...

scenario patch

dimensions:
  # will be auto-populated with base dims `name` and `time`
  scenario: ...
parameters:
  # will be auto-populated with base param `foo`
  bar: {dims: [name, time, scenario]} # will be auto-populated with other `bar` def key:val pairs
variables:
   # will be auto-populated with base var `capacity`
  dispatch:
    foreach: [name, time, scenario]
constraints:
   # will be auto-populated with base constr `constrain_capacity`
  constrain_dispatch:
    foreach: [name, time, scenario]
  global_constraint:
    ... # might need the math updated to include summing over `scenario`

It's like with snakemake YAML configs, that you can have a file with only the changes relative to default, to avoid duplications. adding the equivalent of bar: {dims: [name, time, scenario]} across all relevant params then might be verbose. But that might well be fine since at least it's explicit!

@brynpickering

@brynpickering about broadcasting: The memory footprint of expanding parameters to their dims will be payed anyway, as soon as they are used in arithmetics or as variable bounds. Its only delayed.

I dont really have an opinion on the recipe/patch stuff yet.

It won't necessarily be paid to the full extent. An example would be that if we have 80 parameters which are all configured to be indexed over a time dim, then 40 vars+constrs also all indexed over a time dim. If you supply all 80 parameters not indexed over time and broadcast them straight away, you will have 120 arrays in memory broadcast over that dimension by the end (params+vars+constrs). If you only broadcast them at the point of applying arithmetic then you might have 2-3 params broadcast over time at any given moment, but not all 80. You'll then end up with only 40 arrays in memory broadcast over that dimension at the end (vars+constrs).

Thats correct, but Im confident that the amount of extra memory is tiny compared to the resulting linopy Objects. Ill get back with a small benchmark of this!

EDIT: Im wrong. See below

@brynpickering

Here's the broadcast-memory benchmark I promised.

Note

The content below was generated by AI (run locally against linopy master).

Same model built both ways; the only difference is whether the resident input params are broadcast to their full dims up front (EAGER) or kept minimal and broadcast by linopy just before arithmetic (LAZY, @brynpickering's position). Peak measured with memray via pytest-benchmem, single run pivoted on the strategy dim:

dev-scripts/test_param_broadcast_mem.py::test_param_broadcast[bigtime-T8760]
            peak (GiB) 
 name              min 
───────────────────────
 (LAZY)     2.46 (1.0) 
 (EAGER)   3.51 (1.42) 

dev-scripts/test_param_broadcast_mem.py::test_param_broadcast[calliope-K2]
             peak (MiB) 
 name               min 
────────────────────────
 (LAZY)     28.08 (1.0) 
 (EAGER)   39.56 (1.41) 

dev-scripts/test_param_broadcast_mem.py::test_param_broadcast[dense-K1]
             peak (MiB) 
 name               min 
────────────────────────
 (LAZY)     18.65 (1.0) 
 (EAGER)   30.13 (1.62) 

dev-scripts/test_param_broadcast_mem.py::test_param_broadcast[fewparams-P20]
             peak (MiB) 
 name               min 
────────────────────────
 (LAZY)     27.80 (1.0) 
 (EAGER)   30.67 (1.10) 

dev-scripts/test_param_broadcast_mem.py::test_param_broadcast[manyterm-K6]
             peak (MiB) 
 name               min 
────────────────────────
 (LAZY)     52.84 (1.0) 
 (EAGER)   64.31 (1.22) 

dev-scripts/test_param_broadcast_mem.py::test_param_broadcast[paramheavy-P160]
             peak (MiB) 
 name               min 
────────────────────────
 (LAZY)     28.45 (1.0) 
 (EAGER)   51.40 (1.81)

P = number of params, K = terms/constraint; all configs have 40 vars+constraints over (name=200, time=T).

The extra peak is exactly the resident broadcast param set (P·N·T·8 bytes): constant across K at fixed P (+11.5 MiB for dense-K1/calliope-K2/manyterm-K6), doubling with P (paramheavy → +22.9 MiB), and scaling with T (bigtime, 8760h → +1.07 GiB). The relative cost only shrinks with K because the vars/constraints in the denominator grow. So it's "tiny" only in the few-params/high-term corner (+10%); in the hourly-resolution regime it's a gigabyte.

Benchmark script (dev-scripts/test_param_broadcast_mem.py)

"""Eager vs lazy parameter broadcasting — peak-memory A/B for the YAML-math layer.

Settles the sub-thread in PyPSA/linopy#561: when the declarative layer stores a
parameter, should it broadcast to the component's full ``dims`` up front (a reliable
dimensional contract) or keep it minimal and let linopy broadcast just before
arithmetic (brynpickering's position)? The disagreement is entirely about **peak**.

Both strategies build the *same* model; only the resident input params differ
(EAGER = broadcast to full dims, LAZY = minimal). A plain pytest-benchmark test —
pytest-benchmem records peak memory onto the same node when run with
``--benchmark-memory``. ``strategy`` and ``config`` are plain ``@parametrize``
params, which pytest-benchmem lifts into analysis **dims** — so the single run
pivots EAGER-vs-LAZY on ``strategy`` with no ``extra_info`` needed.

Run:
    pytest dev-scripts/test_param_broadcast_mem.py \
        --benchmark-only --benchmark-memory \
        --benchmark-json=/tmp/param_bcast.json
    python dev-scripts/test_param_broadcast_mem.py /tmp/param_bcast.json   # report
"""

from __future__ import annotations

import numpy as np
import pytest
import xarray as xr

import linopy


def _params(strategy, P, N, T, name, time):
    """P parameters over `name`; EAGER also broadcasts each to (name, time)."""
    base = [
        xr.DataArray(np.arange(N) + 1.0 + i, dims="name", coords={"name": name})
        for i in range(P)
    ]
    if strategy == "EAGER":
        ones_t = xr.DataArray(np.ones(T), dims="time", coords={"time": time})
        base = [(p * ones_t).transpose("name", "time") for p in base]
    return base


def build(strategy, N, T, P, V, C, K):
    """Build V variables over (name,time) and C constraints of K terms each."""
    name = np.arange(N)
    time = np.arange(T)
    m = linopy.Model()

    params = _params(strategy, P, N, T, name, time)  # resident for the whole build
    variables = [
        m.add_variables(lower=0, coords=[("name", name), ("time", time)], name=f"x{v}")
        for v in range(V)
    ]
    for c in range(C):
        expr = sum(
            params[(c * K + k) % P] * variables[(c * K + k) % V] for k in range(K)
        )
        m.add_constraints(expr <= 10.0, name=f"c{c}")
    return m


# config label -> (N, T, P, V, C, K)
CONFIGS = {
    "dense-K1":        (200,  96,  80, 20, 20, 1),
    "calliope-K2":     (200,  96,  80, 20, 20, 2),
    "manyterm-K6":     (200,  96,  80, 20, 20, 6),
    "paramheavy-P160": (200,  96, 160, 20, 20, 2),
    "fewparams-P20":   (200,  96,  20, 20, 20, 2),
    "bigtime-T8760":   (200, 8760, 80, 20, 20, 2),
}


# `strategy` and `config` are clean scalar params, so pytest-benchmem lifts them
# straight out of the parametrize as analysis dims — no extra_info needed. The run
# then pivots on `strategy` (EAGER vs LAZY), pairing rows on `config`.
@pytest.mark.parametrize("strategy", ["LAZY", "EAGER"])
@pytest.mark.parametrize("config", list(CONFIGS))
def test_param_broadcast(benchmark, strategy, config):
    N, T, P, V, C, K = CONFIGS[config]
    benchmark.pedantic(build, args=(strategy, N, T, P, V, C, K), rounds=1, iterations=1)


def _report(path):
    """Pivot the combined run on `strategy` → paired EAGER-vs-LAZY peak, with Δ."""
    from pytest_benchmem import load_long_df

    df, unit = load_long_df(path, metric="peak", pivot="strategy")
    scale = 2**20 if unit == "B" else 1  # normalise to MiB
    wide = df.pivot_table(index="config", columns="snapshot", values="value") / scale
    order = list(CONFIGS)
    wide = wide.reindex(order)

    print("\n## Eager vs lazy parameter broadcasting — peak memory (MiB, memray)\n")
    print("Same model both ways; only resident input params differ.\n")
    print(f"| {'config':16} | K |   P | V+C |    LAZY |   EAGER |      Δ | Δ % |")
    print("|" + "-" * 63 + "|")
    for label in order:
        lazy, eager = wide.loc[label, "LAZY"], wide.loc[label, "EAGER"]
        N, T, P, V, C, K = CONFIGS[label]
        d = eager - lazy
        print(f"| {label:16} | {K} | {P:3d} | {V + C:3d} | {lazy:7.1f} | {eager:7.1f} | "
              f"{d:6.1f} | {100 * d / lazy:3.0f} |")
    print()


if __name__ == "__main__":
    import sys

    _report(sys.argv[1])

Run: pytest dev-scripts/test_param_broadcast_mem.py --benchmark-only --benchmark-memory --benchmark-json=run.json then benchmem compare run.json --pivot param:strategy --columns peak --stat min.

I'll not be able to review this this week.

As we discussed, a few small points, we'd want to to be able to leverage this math eventually to register streaming constraints. ie. you pass data as a lazy polars long-form frame and the math for a constraint:

the math is parsed into an ast, dimensions and shapes are determined (and maybe sparseness if cheap enough), but no data is evaluated yet. only when the CSR matrix is needed, it is constructed on the fly. and dropped again (unless there was a good reason not to, like persistent solver).

the promise is that we don't need to hold any memory for the constraint anymore.

we also should evaluate where the point is to convert then from the lazy frame to xarray. potentially never. potentially based on specific triggers in the ast.

i plan to read it on Monday.

This lazy construction might become essential when we get further into stochastic space. the mentioned scenario dimension can easily become too large for everything to fit into memory and we'd need to break it down into sub problems and instantiate only parts of anywhere

@FBumann

@FBumann thanks for the benchmarking work. I hadn't considered the fact that each K consists of multiple arrays in a wrapped dataset so I would have probably come round to your position if not for the benchmarking!

@coroa I hadn't thought about streaming constraints so that they never exist as linopy objects for longer than necessary, with the YAML being the source of truth that one returns to if they want to re-construct the constraint at any point. Makes sense when you don't have/want a persistent solver. What you would need in that case, for reproducibility, is the ability to keep the YAML math close to the input data when dumping a model representation to file.

This would say if this is the long term goal for the yaml layer, of course needs to be part of the core package

Navigation Menu

Repository navigation

Text-based (YAML) math definition interface #561

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Text-based (YAML) math definition interface #561

Description

Activity

coroa commented on Jan 30, 2026

brynpickering commented on Mar 6, 2026

FBumann commented on Mar 10, 2026

brynpickering commented on Jul 7, 2026

FBumann commented on Jul 7, 2026

FBumann commented on Jul 7, 2026

brynpickering commented on Jul 8, 2026

FBumann commented on Jul 8, 2026

FBumann commented on Jul 8, 2026

brynpickering commented on Jul 9, 2026

FBumann commented on Jul 9, 2026

brynpickering commented on Jul 9, 2026

FBumann commented on Jul 9, 2026

brynpickering commented on Jul 9, 2026

FBumann commented on Jul 9, 2026

brynpickering commented on Jul 9, 2026

FBumann commented on Jul 9, 2026

FBumann commented on Jul 9, 2026

coroa commented on Jul 9, 2026

coroa commented on Jul 9, 2026

brynpickering commented on Jul 10, 2026

FBumann commented on Jul 10, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions