Opt-in allocator arena release after Model.solve() to bound RSS across sequential solves

[MaykThewessen]

Summary

Sequential PyPSA-Eur-style chunked solves (weekly LPs, identical structure, N iterations in one process) exhibit monotonic RSS growth until OOM, even after del model + gc.collect() between iterations. Root cause is platform allocator arena retention (glibc on Linux, libSystem on macOS), not Python-level leaks. Proposing a small opt-in helper that releases arena pages back to the kernel after the existing solve path completes.

This is orthogonal to:

• perf: matrix accessor rewrite #630 (CSR freeze=True) — reduces peak within a single solve.
• Persistent solver: in-place model updates via ModelDiff #699 (persistent solver / ModelDiff) — avoids rebuilds between solves.
• Optionally release memory during solving process #219 — NetCDF round-trip to lower peak during a single solve (has an unresolved SOS-undo round-trip gap, see discussion there).

Neither addresses the allocator-arena pattern: glibc / libSystem keep freed pages cached, so RSS grows monotonically across the chunk loop even when the Python-side working set is bounded.

Proposal

linopy/_memory.py (~30 LOC):

def release_allocator_pages() -> None:
    """gc.collect() then ask platform allocator to return arenas to kernel.

    Linux:   libc.malloc_trim(0)
    macOS:   libSystem.malloc_zone_pressure_relief(NULL, 0)
    Other:   gc only, silent no-op.
    """

Wired into Model.solve(release_memory: bool = False). Default unchanged. Also exposed at top level (linopy.release_allocator_pages()) so callers in chunked-solve loops can trim mid-pipeline without going through Model.solve.

Benchmark (macOS arm64, linopy master 37af4ba)

5 sequential solves of an identical 5k var × 2k con dense-block LP via HiGHS, single Python process:

metric	baseline	release_memory=True	Δ
RSS post-cleanup, iter 4 (GB)	4.41	2.02	−54 %
ru_maxrss (GB)	5.30	4.32	−19 %
wall-clock total (s)	30.97	32.09	+3.6 % (noise)

Baseline RSS climbs 3.71 → 4.41 GB monotonically across the 5 iterations. With trim, post-cleanup RSS drops back to 1.5–2.0 GB after each solve. Linux numbers via malloc_trim typically more dramatic — happy to add a Linux benchmark before opening the PR.

Tests

• Platform-mocked dispatch (Linux / Darwin / Windows / FreeBSD).
• Kwarg default backward-compat smoke test.
• Repeated-call stability (no crashes if libc / libSystem lookup misses).
• Swallows OSError / AttributeError, debug-logs on miss.

Open questions

1. Naming: release_memory=True on Model.solve vs. top-level linopy.release_allocator_pages() vs. both. Leaning toward both — kwarg for the common single-solve case, public helper for chunked loops where users want to trim between solves without re-entering Model.solve.
2. Interaction with Persistent solver: in-place model updates via ModelDiff #699: once persistent solver lands, native solver model + factor stay resident. Python-side arena cache still grows; trim hook stays useful. Agree?
3. Should the helper also clear m._xCounter / m._cCounter style state on Model.solve(release_memory=True)? Leaning no — out of scope for an allocator helper, separate concern.

Will open PR once direction is confirmed.

Refs: #219, #630, #699

[FabianHofmann]

@MaykThewessen, let's think about this. I am not fully sure which direction that should go in the end. could you help me out with an example workflow where the memory should be released? are we talking about creating linopy models in a for-loop?

[MaykThewessen]

Yes, it's a for-loop that builds a fresh Model per chunk, solves, extracts, discards, all in one long-lived process. Three cases where this matters to me:

• Rolling-horizon dispatch: a full-year LP split into ~52 weekly chunks, solved sequentially with storage state carried forward.
• Myopic pathway: sequential planning-year solves (2030 → 2035 → …), each seeded with the previous year's capacities.
• Scenario sweeps: same network, hundreds of parameter draws.

In all three the Python working set is bounded but RSS isn't, and the natural trim point is right after solve() returns. So the two-surface proposal (Model.solve(release_memory=False) + linopy.release_allocator_pages()) fits. I'm happy to open the PR if you're OK with the direction.

I also ran a reproducer to check whether the hook is actually needed vs. just tuning the allocator. Marked + collapsed below.

Note

Benchmark below produced with an AI agent (script + runs). Numbers are real; prose in this block is generated.

Reproducer + results: RSS across sequential solves (macOS libSystem & Linux glibc)

Identical-structure dense-block LPs (5000 var × 2000 con), built/solved sequentially, del + gc.collect() between iterations. --vary oscillates the per-week constraint count (the realistic case: active constraint count differs per chunk → fragmentation). Each mode runs in its own fresh process so the trim run doesn't inherit the baseline run's arenas.

Post-cleanup RSS per iteration, --vary, 8 solves:

platform	baseline trajectory (GB)	+ release_memory (GB)
macOS arm64 (libSystem)	3.81 → 4.31 → 5.04 → 5.26 → 6.41 → 7.42 → 7.64 (climbs, no plateau)	~1.4 → 2.4 (bounded)
Linux aarch64 (glibc)	~3.0, flat from iter 2 (plateaus)	~1.5 (−47%)

So unbounded growth is mainly a macOS/libSystem behaviour; glibc reaches steady state but still holds ~2× the trimmed RSS. Wall-clock cost of the hook is ~0% on both.

glibc env tuning is not a substitute (Linux, same workload):

config	baseline RSS (GB)	+ hook (GB)	baseline wall (s)
default	~2.9	1.60	47
MALLOC_ARENA_MAX=1	~2.1	1.69	49
MALLOC_TRIM_THRESHOLD_=0	~2.5	1.45	hook run 56 (+18%)

ARENA_MAX=1 helps but stays ~35% above the hook; TRIM_THRESHOLD_=0 matches RSS but trims on every free() (+18% wall). A single post-solve malloc_trim(0) gets bounded RSS at ~0% runtime cost, which supports keeping it opt-in (release_memory=False).

benchmark script (also happy to drop in dev-scripts/ with the PR)

#!/usr/bin/env python
"""Reproducer for linopy#720: allocator-arena RSS across sequential solves."""
from __future__ import annotations
import argparse, ctypes, ctypes.util, gc, json, platform, resource, subprocess, sys, time
import numpy as np, pandas as pd, psutil, xarray as xr
import linopy

GB = 1024**3
PROC = psutil.Process()

def release_allocator_pages() -> None:
    gc.collect()
    system = platform.system()
    try:
        if system == "Linux":
            ctypes.CDLL("libc.so.6", use_errno=True).malloc_trim(0)
        elif system == "Darwin":
            lib = ctypes.CDLL(ctypes.util.find_library("System"))
            lib.malloc_zone_pressure_relief.argtypes = [ctypes.c_void_p, ctypes.c_size_t]
            lib.malloc_zone_pressure_relief.restype = ctypes.c_size_t
            lib.malloc_zone_pressure_relief(None, 0)
    except (OSError, AttributeError):
        pass

def rss_gb() -> float:
    return PROC.memory_info().rss / GB

def build_and_solve(week: int, n_var: int, n_con: int) -> float:
    rng = np.random.default_rng(week)
    A = rng.uniform(0.1, 1.0, size=(n_con, n_var))
    b = A.sum(axis=1) * rng.uniform(0.4, 0.6, size=n_con)
    cost = rng.uniform(1.0, 50.0, size=n_var)
    m = linopy.Model()
    x = m.add_variables(lower=0.0, upper=10.0, coords=[pd.RangeIndex(n_var, name="i")], name="x")
    A_da = xr.DataArray(A, dims=["c", "i"], coords={"c": np.arange(n_con)})
    b_da = xr.DataArray(b, dims=["c"], coords={"c": np.arange(n_con)})
    cost_da = xr.DataArray(cost, dims=["i"], coords={"i": np.arange(n_var)})
    m.add_constraints((A_da * x).sum("i") <= b_da, name="dense")
    m.add_objective((cost_da * x).sum())
    m.solve(solver_name="highs", output_flag=False)
    obj = float(m.objective.value)
    del m, x, A, b, cost, A_da, b_da, cost_da
    return obj

def week_n_con(week: int, n_con: int, vary: bool) -> int:
    if not vary:
        return n_con
    factors = np.random.default_rng(0).uniform(0.5, 1.6, size=max(week + 1, 64))
    return int(n_con * factors[week])

def run_loop(weeks, n_var, n_con, trim, vary):
    trace = []
    for week in range(weeks):
        build_and_solve(week, n_var, week_n_con(week, n_con, vary))
        release_allocator_pages() if trim else gc.collect()
        trace.append(rss_gb())
    return trace

def run_single_mode(args):
    trim = args.mode == "release_memory"
    t0 = time.perf_counter()
    trace = run_loop(args.weeks, args.n_var, args.n_con, trim, args.vary)
    wall = time.perf_counter() - t0
    ru = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
    maxrss = ru / GB if platform.system() == "Darwin" else ru / (1024**2)
    return dict(trace=trace, wall=wall, maxrss=maxrss)

def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--weeks", type=int, default=5)
    ap.add_argument("--n-var", type=int, default=5000)
    ap.add_argument("--n-con", type=int, default=2000)
    ap.add_argument("--vary", action="store_true")
    ap.add_argument("--mode", choices=["baseline", "release_memory"], default=None)
    args = ap.parse_args()
    if args.mode is not None:
        print(json.dumps(run_single_mode(args)))
        return
    print(f"platform={platform.system()} {platform.machine()} | linopy={linopy.__version__} | "
          f"weeks={args.weeks} | LP={args.n_var}x{args.n_con} (dense)\n")
    for mode in ("baseline", "release_memory"):
        cmd = [sys.executable, __file__, "--mode", mode, "--weeks", str(args.weeks),
               "--n-var", str(args.n_var), "--n-con", str(args.n_con)] + (["--vary"] if args.vary else [])
        r = json.loads(subprocess.run(cmd, capture_output=True, text=True, check=True).stdout.splitlines()[-1])
        print(f"[{mode}] RSS/iter (GB): " + ", ".join(f"{v:.2f}" for v in r["trace"]))
        print(f"[{mode}] wall={r['wall']:.1f}s  ru_maxrss={r['maxrss']:.2f} GB\n")

if __name__ == "__main__":
    main()

[FBumann]

@MaykThewessen Please read our contribution guidelines about AI usage and try to keep volume down, or collapse verbose stuff. See Agent.md