FLToptim — Fast Legendre Transform Optim

Fast convex optimisation for energy-storage operation and spatial economic dispatch, built around the Legendre-transform dynamic program for convex piecewise-linear (and quadratic) value functions.

What's inside

• cplfunction / cpqfunction — the convex piecewise-linear/quadratic value-function DP (OptimMargInt). Each function is stored as a map<position, slope-increment>, so the inf-convolutions of the Legendre / dynamic-programming recursion stay near-linear. This is the original dynprogstorage core (Girard, 2013).
• param_simplex — a bespoke parametric right-hand-side dual simplex that traces the exact injection-cost curve φ_n(y) of a small dense network LP in a single pass (one dual pivot per breakpoint), carrying its basis across hours.
• elec / mr_decompose — the spatial-LP ↔ storage-DP decomposition for perfect-foresight annual dispatch: alternate a per-hour spatial network LP with the per-node storage DP, passing each node's convex injection-cost curve (not a scalar price) so the scheme is curvature-damped and converges to the monolithic optimum. Handles reservoirs, batteries (efficiency kink) and coupled multi-energy (electricity + hydrogen) networks, with optional parallelism (hour-chunk curve build; Jacobi across storages).

Interfaces (pommes, PyPSA)

FLToptim's engines consume plain-numpy case dicts, so coupling an external model is a thin adapter. Two ship with the package under FLToptim.interfaces (lazy imports -- the core stays numpy + highspy only):

• pommes (pip install "fltoptim[pommes]", soft dependency: xarray only) — maps a single-year pommes operation Dataset at pinned capacities onto the elec or multi-resource case, including flat-band flexible demand as level-pinned storages:

  from FLToptim.interfaces import pommes_to_elec_case, solve_case_monolith
  from FLToptim.elec import decompose
  
  case, mapping = pommes_to_elec_case(op_params, capacities)  # xr.Datasets
  obj_ref, prices = solve_case_monolith(case, duals=True)     # exact HiGHS LP
  x, hist = decompose(case, sweeps=8)                         # FLT decomposition
  

  The pommes-side Benders driver (pommes-bender) uses this mapping as its FLToptim operation oracle.

• PyPSA / PyPSA-Eur (pip install "fltoptim[pypsa]") — maps a single-carrier electricity Network (Loads, Generators, Links, Lines, StorageUnits) onto the elec case and writes the dispatch back:

  from FLToptim.interfaces import (pypsa_to_elec_case, transport_equivalent,
                                   schedule_to_pypsa, solve_case_monolith)
  
  case, mapping = pypsa_to_elec_case(n)          # n: pypsa.Network, capacities pinned
  obj, prices, sol = solve_case_monolith(case, duals=True, solution=True)
  n2, extras = schedule_to_pypsa(n, case, mapping, sol, prices=prices)
  
  te = transport_equivalent(n)   # lines -> links: PyPSA solves the SAME LP
  te.optimize()                  # -> te.objective == obj (exact reference)
  

  See examples/06_pypsa_toy_dispatch.py and examples/07_pypsa_eur_redispatch.py (PyPSA-Eur networks, e.g. the pre-built ones on Zenodo).

Scope and caveats (details in the module docstrings; out-of-scope inputs raise NotImplementedError):

Install / build

pip install fltoptim          # from PyPI (builds from source on platforms without a wheel)
# or, from a checkout:
pip install -e .              # builds the Cython DP extension + the parametric-simplex shared library
# or
python setup.py build_ext --inplace

Building from source needs a C++17 compiler (Linux/macOS); numpy and highspy are the runtime dependencies. The reference frontal solves in the decomposition use Gurobi when available (optional).

DP value-function storage

The linear storage DP uses the sorted-array backend by default (flat=True, and dp_flat=True in elec.decompose / mr_decompose). It runs the same algorithm and produces schedules bit-for-bit identical to the historical std::map backend, while measured annual workloads are about 4.5--6 times faster in the DP itself. Pass flat=False or dp_flat=False to retain the historical backend for comparison or diagnostics.

With the parametric-simplex engine, dp_fused=True (also the default) passes the batched slope/breakpoint buffers directly into the array DP. This avoids building a temporary map-backed Pycplfunctionvec and then converting it back to arrays. Set dp_fused=False only to compare with that historical two-step path; it does not change the curves, pruning, schedule reconstruction, or feasibility projection.

For the coupled multi-resource engine, fast_cost=True (the default) also reads the exact end-of-sweep system objective from the last committed parametric trace. The trace is accepted only inside its certified domain and when it was not truncated; otherwise the code falls back to direct hourly LP evaluation. The direct monitor uses its own MrSpatial instance and therefore cannot perturb the basis used to seed the next sweep. Set fast_cost=False only for a direct-monitor reference. For an audited run, fast_cost_verify_every=k compares both values every k sweeps and permanently falls back to the direct monitor if their relative difference exceeds fast_cost_verify_rtol (default 1e-9). On the annual electricity--hydrogen benchmark, the traced and direct trajectories are identical and their costs agree to about 1e-15 relatively; avoiding 8,760 monitoring LPs per sweep reduced a 40-sweep local run from 14.09 s to 2.67 s.

The array backend is the recommended choice while intermediate value functions contain the usual handful to few dozen breakpoints. The map can theoretically become preferable when functions retain hundreds or thousands of breakpoints and the recursion performs many insertions in their middle, because array insertion is linear in the number of breakpoints. No such crossover has been observed in the supported electricity or multi-resource annual workloads; benchmark both backends before opting out rather than using a fixed breakpoint threshold.

For the spatial decomposition, stop_rel_cost=tol optionally stops the outer sweeps when two consecutive system costs differ by at most tol relatively, after min_sweeps passes. The default remains a fixed number of sweeps because larger networks can oscillate rather than converge monotonically. A tolerance of 1e-6 reduced the measured 3-node annual case from eight to four passes, but did not trigger within 30 passes on the 10-node case.

reuse_basis=True is retained as an experimental diagnostic: it carries the parametric-simplex basis between outer sweeps and removes repeated HiGHS seeds. Those seeds are already negligible, however, and degenerate optima can send the outer iteration along a different trajectory. It is therefore disabled by default and is not currently recommended as a performance option.

For multi-node cases with one storage direction per node, shared_hourly_bases=True builds one optimal basis per hour at the common current dispatch and uses this atlas to seed every nodal/storage trace. Complete injection-cost curves are still retraced, so this is a seed/factorisation experiment rather than curve caching. Annual measurements reduced trace time by roughly 29% at 3 nodes, 26% at 8 nodes and 38% at 10 nodes. Degenerate hourly optima can nevertheless change the outer Gauss--Seidel trajectory. The option is therefore off by default.

The implementation compares two atlas candidates at each hour: the basis propagated from hour t-1 and the basis of hour t from the preceding outer sweep. hourly_basis_selection exposes five policies: longitudinal, certain (switch only when the previous-hour candidate is infeasible and the previous- sweep candidate is feasible), predict (fewest violated basic variables, then least normalised violation), canonical (the predictor plus a lexicographic, candidate-order-independent signature when both bases are feasible optima), and min_pivots (run both reoptimisations). On the first annual transition, primal admissibility gives a certain zero-pivot choice for 55% / 26% / 20% of hours at 3 / 8 / 10 nodes when the longitudinal candidate is already feasible; the reverse certain choice adds 15% / 14% / 16%. Among the remaining neither-feasible cases, the violation predictor identifies the actual minimum-pivot candidate 72% / 77% / 80% of the time. canonical is the default selection inside the atlas: paired annual runs gave 0.386 s at +0.00258% (3 nodes), 7.84 s at +0.06004% (8 nodes), and 11.17 s at +0.33686% (10 nodes), versus 18.18 s at +0.32719% for the standard 10-node engine. Fewer spatial pivots still do not mathematically guarantee a better outer Gauss--Seidel trajectory; the canonical rule makes the observed degeneracy reproducible rather than perturbing the economic objective.

Keep the atlas disabled for one-node problems, a single curve per hour, ramp-aware or multi-resource decompositions, or whenever constructing it costs more than the repeated reoptimisations it replaces. It is most useful once several nodal storage curves are rebuilt at every outer sweep (the annual 8- and 10-node benchmarks are representative); measure both modes on smaller systems.

Benchmarks

BENCHMARKS.md is the maintained report: every measured compute time of the current version (investment Benders vs frontal, sizing, operation, stochastic), the algorithm inventory, and the history of superseded measurements.

Documentation

Online: https://pages.persee.minesparis.psl.eu/fltoptim-3df8bd/ (rebuilt from main by the pages CI job on every push).

The Sphinx sources (overview, interfaces guide, API reference) live in docs/:

pip install "fltoptim[interfaces,docs]"
sphinx-build -b html docs docs/_build/html    # then open docs/_build/html/index.html

It is published automatically by .github/workflows/docs.yml (GitHub Pages) and by the pages job in .gitlab-ci.yml (GitLab Pages, once the instance has a runner).

The test suite ships with the package:

pip install "fltoptim[test]"
python -m pytest --pyargs FLToptim.tests -m "not slow" -q

Reference

R. Girard, Fast dynamic programming with application to storage planning, 2013. Please cite it if you use this software (see CITATION.cff).