Add scaling

doc/api.rst

@@ -130,6 +130,7 @@ Attributes
 
    variables.Variable.lower
    variables.Variable.upper
+   variables.Variable.scaling
    variables.Variable.type
    variables.Variable.solution
 
@@ -182,6 +183,7 @@ Attributes
 
    variables.Variables.lower
    variables.Variables.upper
+   variables.Variables.scaling
    variables.Variables.solution
 
 Modification
@@ -329,6 +331,7 @@ Structure
    constraints.Constraint.lhs
    constraints.Constraint.sign
    constraints.Constraint.rhs
+   constraints.Constraint.scaling
    constraints.Constraint.coeffs
    constraints.Constraint.vars
 
@@ -372,6 +375,7 @@ Structure
    constraints.CSRConstraint.vars
    constraints.CSRConstraint.sign
    constraints.CSRConstraint.rhs
+   constraints.CSRConstraint.scaling
    constraints.CSRConstraint.ncons
    constraints.CSRConstraint.nterm
 
@@ -422,6 +426,7 @@ Aggregate access
    constraints.Constraints.vars
    constraints.Constraints.sign
    constraints.Constraints.rhs
+   constraints.Constraints.scaling
    constraints.Constraints.dual
 
 Conversion
@@ -445,6 +450,7 @@ Wraps the objective expression on a model. Accessed via
    objective.Objective
    objective.Objective.expression
    objective.Objective.sense
+   objective.Objective.scaling
    objective.Objective.value
    objective.Objective.is_linear
    objective.Objective.is_quadratic

doc/user-guide.rst

@@ -33,6 +33,42 @@ you're unsure what a particular operator or argument does.
 After these four you can build any LP/MIP/QP linopy supports.
 
 
+Numerical scaling
+-----------------
+
+Large or very small coefficients can make mathematical programs harder
+for solvers to process. Linopy therefore accepts optional positive,
+finite ``scaling`` factors when adding variables, constraints, and the
+objective.
+
+For continuous and semi-continuous variables, ``scaling=s`` means the
+solver-side column represents ``s * x``. Bounds are multiplied by ``s``,
+and coefficients that reference ``x`` are divided by ``s``. Binary and
+integer variable scaling is stored and round-tripped, but their solver
+columns remain ordinary discrete columns.
+
+For constraints, ``scaling=s`` multiplies both the left-hand-side
+coefficients and the right-hand side by ``s``. Objective scaling is
+row-like: ``scaling=s`` multiplies linear and quadratic objective
+coefficients by ``s``. Primal values, dual values, and objective values
+are transformed back to the original user units after solving.
+
+Equivalently, for constraint matrix ``A``, right-hand side ``b``,
+linear objective coefficients ``c``, variable scaling matrix ``Sx``,
+constraint scaling matrix ``Sc``, and scalar objective scaling ``so``,
+linopy exports ``y = Sx x`` and
+``A_solver = Sc A Sx^-1``, ``b_solver = Sc b``, and
+``c_solver = so c Sx^-1``. Quadratic objective coefficients receive the
+inverse variable scaling once for each variable factor and the scalar
+objective scaling once globally.
+
+.. code-block:: python
+
+    x = m.add_variables(lower=0, scaling=1e3, name="x")
+    m.add_constraints(2 * x >= 10, scaling=10, name="demand")
+    m.add_objective(5 * x, scaling=100)
+
+
 Working with an existing model
 ------------------------------
 

linopy/common.py

@@ -40,6 +40,7 @@
 if TYPE_CHECKING:
     from linopy.constraints import ConstraintBase
     from linopy.expressions import LinearExpression, QuadraticExpression
+    from linopy.model import Model
     from linopy.variables import Variable
 
 
@@ -97,6 +98,61 @@ def maybe_replace_signs(sign: DataArray) -> DataArray:
     return apply_ufunc(func, sign, dask="parallelized", output_dtypes=[sign.dtype])
 
 
+def validate_scaling(scaling: DataArray, label: str = "scaling") -> DataArray:
+    """
+    Validate and normalize a scaling array.
+
+    Scaling values are positive numeric factors used during solver export.
+    Column scaling converts variable units, while row and objective scaling
+    multiply exported coefficients. They must be finite and strictly positive.
+    """
+    scaling = scaling.astype(float)
+    values = np.asarray(scaling.values, dtype=float)
+    if not np.isfinite(values).all() or (values <= 0).any():
+        raise ValueError(f"{label} must contain only finite positive values.")
+    return scaling
+
+
+def variable_solver_scaling(variable: Variable) -> DataArray:
+    """
+    Return solver-side scaling for a variable.
+
+    Discrete variables keep ordinary solver columns, so their supplied scaling
+    is metadata only.
+    """
+    if variable.attrs.get("binary", False) or variable.attrs.get("integer", False):
+        return DataArray(1.0).broadcast_like(variable.labels)
+    return variable.scaling
+
+
+def variable_scaling_lookup(model: Model) -> np.ndarray:
+    """Return solver-side variable scaling indexed by raw variable label."""
+    scaling = np.ones(model._xCounter, dtype=float)
+    for var in model.variables.data.values():
+        labels = var.labels.values.ravel()
+        mask = labels != -1
+        if mask.any():
+            scaling[labels[mask]] = variable_solver_scaling(var).values.ravel()[mask]
+    return scaling
+
+
+def constraint_scaling_lookup(model: Model) -> np.ndarray:
+    """Return constraint scaling indexed by raw constraint label."""
+    scaling = np.ones(model._cCounter, dtype=float)
+    for con in model.constraints.data.values():
+        if hasattr(con, "_con_labels") and hasattr(con, "_scaling"):
+            labels = con._con_labels
+            if len(labels):
+                scaling[labels] = con._scaling
+            continue
+
+        labels = con.labels.values.ravel()
+        mask = labels != -1
+        if mask.any():
+            scaling[labels[mask]] = con.scaling.values.ravel()[mask]
+    return scaling
+
+
 def format_string_as_variable_name(name: Hashable) -> str:
     """
     Format a string to a valid python variable name.

linopy/constraints.py

@@ -31,6 +31,7 @@
 from xarray.core.utils import Frozen
 
 from linopy import expressions, variables
+from linopy.alignment import broadcast_to_coords
 from linopy.common import (
     ConstraintLabelIndex,
     LabelPositionIndex,
@@ -59,6 +60,7 @@
     save_join,
     to_dataframe,
     to_polars,
+    validate_scaling,
 )
 from linopy.config import options
 from linopy.constants import (
@@ -82,7 +84,14 @@
     from linopy.model import Model
 
 
-FILL_VALUE = {"labels": -1, "rhs": np.nan, "coeffs": 0, "vars": -1, "sign": "="}
+FILL_VALUE = {
+    "labels": -1,
+    "rhs": np.nan,
+    "coeffs": 0,
+    "vars": -1,
+    "sign": "=",
+    "scaling": 1.0,
+}
 
 
 def conwrap(
@@ -169,6 +178,11 @@ def sign(self) -> DataArray:
     def rhs(self) -> DataArray:
         """Get the RHS DataArray."""
 
+    @property
+    @abstractmethod
+    def scaling(self) -> DataArray:
+        """Get the row scaling DataArray."""
+
     @property
     @abstractmethod
     def dual(self) -> DataArray:
@@ -525,6 +539,7 @@ class CSRConstraint(ConstraintBase):
         "_con_labels",
         "_rhs",
         "_sign",
+        "_scaling",
         "_coords",
         "_model",
         "_name",
@@ -540,6 +555,7 @@ def __init__(
         con_labels: np.ndarray,
         rhs: np.ndarray,
         sign: str | np.ndarray,
+        scaling: np.ndarray | None,
         coords: list[pd.Index],
         model: Model,
         name: str = "",
@@ -552,6 +568,11 @@ def __init__(
         self._con_labels = con_labels
         self._rhs = rhs
         self._sign = sign
+        self._scaling = (
+            np.asarray(scaling, dtype=float)
+            if scaling is not None
+            else np.ones_like(rhs, dtype=float)
+        )
         self._coords = coords
         self._model = model
         self._name = name
@@ -685,6 +706,11 @@ def rhs(self, value: ConstantLike) -> None:
             "CSRConstraint.rhs is read-only; call .mutable() to modify."
         )
 
+    @property
+    def scaling(self) -> DataArray:
+        """Get row scaling DataArray, shape (*coord_dims)."""
+        return self._active_to_dataarray(self._scaling, fill=1.0)
+
     @property
     def lhs(self) -> expressions.LinearExpression:
         """Get LHS as LinearExpression (triggers Dataset reconstruction)."""
@@ -784,7 +810,11 @@ def _to_dataset(self, nterm: int) -> Dataset:
     def data(self) -> Dataset:
         """Reconstruct the xarray Dataset from the CSR representation."""
         ds = self._to_dataset(self.nterm)
-        extra: dict[str, Any] = {"sign": self.sign, "rhs": self.rhs}
+        extra: dict[str, Any] = {
+            "sign": self.sign,
+            "rhs": self.rhs,
+            "scaling": self.scaling,
+        }
         if self._dual is not None:
             extra["dual"] = self._active_to_dataarray(self._dual, fill=np.nan)
         if self._binvar_labels is not None:
@@ -866,6 +896,7 @@ def to_netcdf_ds(self) -> Dataset:
             "indices": DataArray(csr.indices, dims=["_nnz"]),
             "data": DataArray(csr.data, dims=["_nnz"]),
             "rhs": DataArray(self._rhs, dims=["_flat"]),
+            "scaling": DataArray(self._scaling, dims=["_flat"]),
             "_con_labels": DataArray(self._con_labels, dims=["_flat"]),
         }
         if isinstance(self._sign, np.ndarray):
@@ -903,6 +934,7 @@ def from_netcdf_ds(cls, ds: Dataset, model: Model, name: str) -> CSRConstraint:
             shape=shape,
         )
         rhs = ds["rhs"].values
+        scaling = ds["scaling"].values if "scaling" in ds else None
         sign: str | np.ndarray = ds["_sign"].values if "_sign" in ds else attrs["sign"]
         _cindex_raw = int(attrs["cindex"])
         cindex: int | None = _cindex_raw if _cindex_raw >= 0 else None
@@ -927,6 +959,7 @@ def from_netcdf_ds(cls, ds: Dataset, model: Model, name: str) -> CSRConstraint:
             con_labels,
             rhs,
             sign,
+            scaling,
             coords,
             model,
             name,
@@ -984,6 +1017,7 @@ def sanitize_infinities(self) -> CSRConstraint:
         self._csr = self._csr[keep]
         self._con_labels = self._con_labels[keep]
         self._rhs = self._rhs[keep]
+        self._scaling = self._scaling[keep]
         if not isinstance(self._sign, str):
             self._sign = self._sign[keep]
         return self
@@ -1054,6 +1088,7 @@ def iterate_slices(
                 con_labels=self._con_labels[rows],
                 rhs=self._rhs[rows],
                 sign=sign,
+                scaling=self._scaling[rows],
                 coords=[],
                 model=self._model,
                 name=self._name,
@@ -1084,6 +1119,7 @@ def from_mutable(
         vars_flat = con.vars.values.reshape(len(labels_flat), -1)
         active_mask = (labels_flat != -1) & (vars_flat != -1).any(axis=1)
         rhs = con.rhs.values.ravel()[active_mask]
+        scaling = con.scaling.values.ravel()[active_mask]
         sign_vals = con.sign.values.ravel()
         active_signs = sign_vals[active_mask]
         unique_signs = np.unique(active_signs)
@@ -1111,6 +1147,7 @@ def from_mutable(
             con_labels,
             rhs,
             sign,
+            scaling,
             coords,
             con.model,
             con.name,
@@ -1153,6 +1190,14 @@ def __init__(
 
         if not skip_broadcast:
             (data,) = xr.broadcast(data, exclude=[TERM_DIM])
+        if "scaling" not in data:
+            data = assign_multiindex_safe(
+                data, scaling=DataArray(1.0).broadcast_like(data.rhs)
+            )
+        data = assign_multiindex_safe(
+            data,
+            scaling=validate_scaling(data.scaling, f"scaling for constraint '{name}'"),
+        )
 
         self._assigned = "labels" in data
         self._data = data
@@ -1242,6 +1287,21 @@ def rhs(self, value: ExpressionLike) -> None:
         self.lhs = self.lhs - value.reset_const()
         self._data = assign_multiindex_safe(self.data, rhs=value.const)
 
+    @property
+    def scaling(self) -> DataArray:
+        return self.data.scaling
+
+    @scaling.setter
+    def scaling(self, value: ConstantLike) -> None:
+        scaling = broadcast_to_coords(
+            value,
+            coords=self.labels.coords,
+            dims=self.coord_dims,
+            label=f"scaling for constraint '{self.name}'",
+        )
+        scaling = validate_scaling(scaling, f"scaling for constraint '{self.name}'")
+        self._data = assign_multiindex_safe(self.data, scaling=scaling)
+
     @property
     def is_indicator(self) -> bool:
         return "binary_var" in self._data
@@ -1557,13 +1617,14 @@ class Constraints:
     _label_position_index: LabelPositionIndex | None = None
     _constraint_label_index: ConstraintLabelIndex | None = None
 
-    dataset_attrs = ["labels", "coeffs", "vars", "sign", "rhs"]
+    dataset_attrs = ["labels", "coeffs", "vars", "sign", "rhs", "scaling"]
     dataset_names = [
         "Labels",
         "Left-hand-side coefficients",
         "Left-hand-side variables",
         "Signs",
         "Right-hand-side constants",
+        "Scaling factors",
     ]
 
     def _formatted_names(self) -> dict[str, str]:
@@ -1732,6 +1793,13 @@ def rhs(self) -> Dataset:
         """
         return save_join(*[v.rhs.rename(k) for k, v in self.items()])
 
+    @property
+    def scaling(self) -> Dataset:
+        """
+        Get the scaling factors of all constraints.
+        """
+        return save_join(*[v.scaling.rename(k) for k, v in self.items()])
+
     @property
     def dual(self) -> Dataset:
         """
@@ -1999,6 +2067,7 @@ def reset_dual(self) -> None:
                         c._con_labels,
                         c._rhs,
                         c._sign,
+                        c._scaling,
                         c._coords,
                         c._model,
                         c._name,