Add pre-solve scaling support

README.md

@@ -38,6 +38,7 @@
 * Use **lazy operations** for large linear programs  with [dask](https://dask.org/)
 * Choose from **different commercial and non-commercial solvers**
 * Fast **import and export** a linear model using xarray's netcdf IO
+* Optional **pre-solve scaling** of constraints (and continuous variables) to improve numerical robustness
 
 
 ## Installation

benchmark/scripts/benchmark_scaling.py

@@ -0,0 +1,100 @@
+#!/usr/bin/env python3
+"""
+Ad-hoc benchmark to compare solve times with and without pre-solve scaling.
+
+This is intentionally lightweight and meant for local experimentation.
+It relies on HiGHS (highspy) being installed. Adjust sizes or iterations
+via CLI flags if you want to stress test further.
+"""
+
+from __future__ import annotations
+
+import argparse
+import time
+from collections.abc import Iterable
+
+import numpy as np
+
+from linopy import Model
+from linopy.scaling import ScaleOptions
+from linopy.solvers import available_solvers
+
+
+def build_model(n_vars: int, n_cons: int, density: float) -> Model:
+    rng = np.random.default_rng(123)
+    m = Model()
+    x = m.add_variables(lower=0, name="x", coords=[range(n_vars)])
+
+    data = rng.lognormal(mean=0.0, sigma=2.0, size=int(n_vars * n_cons * density))
+    rows = rng.integers(0, n_cons, size=data.size)
+    cols = rng.integers(0, n_vars, size=data.size)
+
+    # accumulate entries per row
+    for i in range(n_cons):
+        mask = rows == i
+        if not mask.any():
+            continue
+        coeffs = data[mask]
+        vars_idx = cols[mask]
+        lhs = sum(coeff * x.isel(dim_0=idx) for coeff, idx in zip(coeffs, vars_idx))
+        rhs = abs(coeffs).sum() * 0.1
+        m.add_constraints(lhs >= rhs, name=f"c{i}")
+
+    obj_coeffs = rng.uniform(0.1, 1.0, size=n_vars)
+    m.objective = (obj_coeffs * x).sum()
+    return m
+
+
+def time_solve(m: Model, scale: bool | ScaleOptions, repeats: int) -> Iterable[float]:
+    for _ in range(repeats):
+        start = time.perf_counter()
+        status, _ = m.solve("highs", io_api="direct", scale=scale)
+        end = time.perf_counter()
+        if status != "ok":
+            raise RuntimeError(f"Solve failed with status {status}")
+        yield end - start
+
+
+def run_benchmark(
+    n_vars: int, n_cons: int, density: float, repeats: int
+) -> tuple[np.ndarray, np.ndarray]:
+    base_model = build_model(n_vars, n_cons, density)
+    scaled_model = build_model(n_vars, n_cons, density)
+
+    base_times = np.fromiter(time_solve(base_model, False, repeats), dtype=float)
+    scaled_times = np.fromiter(time_solve(scaled_model, True, repeats), dtype=float)
+    return base_times, scaled_times
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--vars", type=int, default=400, help="Number of variables.")
+    parser.add_argument("--cons", type=int, default=300, help="Number of constraints.")
+    parser.add_argument(
+        "--density",
+        type=float,
+        default=0.01,
+        help="Constraint density (0-1) for random coefficients.",
+    )
+    parser.add_argument(
+        "--repeats", type=int, default=3, help="Number of solve repetitions."
+    )
+    args = parser.parse_args()
+
+    if "highs" not in available_solvers:
+        raise RuntimeError("HiGHS (highspy) is required for this benchmark.")
+
+    base_times, scaled_times = run_benchmark(
+        n_vars=args.vars, n_cons=args.cons, density=args.density, repeats=args.repeats
+    )
+
+    print(f"Solve times without scaling: {base_times}")
+    print(f"Solve times with scaling   : {scaled_times}")
+    print(
+        f"Median speedup: {np.median(base_times) / np.median(scaled_times):.2f}x "
+        f"(lower is better for scaled)"
+    )
+
+
+if __name__ == "__main__":
+    main()

doc/index.rst

@@ -110,6 +110,7 @@ This package is published under MIT license.
    creating-expressions
    creating-constraints
    manipulating-models
+   scaling
    testing-framework
    transport-tutorial
    infeasible-model

doc/scaling.rst

@@ -0,0 +1,72 @@
+.. _scaling:
+
+Pre-solve Scaling
+=================
+
+Linopy can rescale your model before handing it to a solver. Scaling can
+improve numerical robustness when constraint coefficients or bounds span
+different orders of magnitude.
+
+How it works
+------------
+
+- **Row scaling** (default) rescales constraint rows by their largest (or
+  RMS) coefficient and applies the same factor to the RHS and duals.
+- **Column scaling** (optional) rescales continuous variables so column
+  norms are close to 1. Bounds, objective coefficients, and primal values
+  are adjusted accordingly. Integer/binary variables stay unscaled by default
+  to preserve integrality.
+- Scaling is undone automatically on primal/dual values and the objective
+  before they are stored on the model.
+
+API
+---
+
+Enable scaling via ``Model.solve(scale=...)``:
+
+.. code-block:: python
+
+    import pandas as pd
+    import linopy
+    from linopy.scaling import ScaleOptions
+
+    m = linopy.Model()
+
+    hours = pd.RangeIndex(24, name="h")
+    p = m.add_variables(lower=0, name="prod", coords=[hours])
+
+    # Coefficients with very different magnitudes
+    m.add_constraints(1e3 * p <= 5e4, name="capacity")
+    m.add_constraints(0.01 * p.sum() >= 10, name="energy")
+
+    # Default: row scaling using max-norm
+    m.solve(scale=True)
+
+    # Custom: row+column scaling on continuous variables, RMS norm
+    m.solve(
+        scale=ScaleOptions(
+            enabled=True,
+            method="row-l2",
+            variable_scaling=True,
+            scale_integer_variables=False,
+        )
+    )
+
+Options
+-------
+
+``scale`` accepts:
+
+- ``False`` / ``None``: disable scaling (default behavior prior to this feature)
+- ``True``: enable row scaling with max-norm
+- ``"row-max"`` or ``"row-l2"``: select norm for row scaling
+- ``ScaleOptions``: full control (row/column scaling, integer handling,
+  target magnitude, zero-floor)
+
+Notes
+-----
+
+- Remote execution currently disables scaling.
+- Many solvers also perform internal scaling; Linopy’s scaling is optional
+  and aims to help when you need deterministic pre-conditioning on the
+  problem you pass to the solver.