diff --git a/README.md b/README.md
index 3cd19df8..381978de 100644
--- a/README.md
+++ b/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
diff --git a/benchmark/scripts/benchmark_scaling.py b/benchmark/scripts/benchmark_scaling.py
new file mode 100644
index 00000000..b013278f
--- /dev/null
+++ b/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()
diff --git a/doc/index.rst b/doc/index.rst
index f05f89f3..9449edba 100644
--- a/doc/index.rst
+++ b/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
diff --git a/doc/scaling.rst b/doc/scaling.rst
new file mode 100644
index 00000000..c4b5695c
--- /dev/null
+++ b/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.
diff --git a/linopy/io.py b/linopy/io.py
index 7065adbb..ba4a4753 100644
--- a/linopy/io.py
+++ b/linopy/io.py
@@ -25,11 +25,14 @@
 from linopy import solvers
 from linopy.constants import CONCAT_DIM
 from linopy.objective import Objective
+from linopy.scaling import ScalingContext
 
 if TYPE_CHECKING:
     from highspy.highs import Highs
 
+    from linopy.matrices import MatrixAccessor
     from linopy.model import Model
+    from linopy.scaling import ScaledMatrices
 
 
 logger = logging.getLogger(__name__)
@@ -162,6 +165,7 @@ def objective_to_file(
     f: BufferedWriter,
     progress: bool = False,
     explicit_coordinate_names: bool = False,
+    scaling: ScalingContext | None = None,
 ) -> None:
     """
     Write out the objective of a model to a lp file.
@@ -176,6 +180,16 @@ def objective_to_file(
     sense = m.objective.sense
     f.write(f"{sense}\n\nobj:\n\n".encode())
     df = m.objective.to_polars()
+    if scaling is not None and getattr(scaling, "col_inv", None) is not None:
+        scale_map = pl.DataFrame(
+            {"vars": pl.Series(scaling.matrices.vlabels), "col_inv": scaling.col_inv}
+        )
+        df = (
+            df.join(scale_map, on="vars", how="left")
+            .with_columns(pl.col("col_inv").fill_null(1.0))
+            .with_columns(pl.col("coeffs") * pl.col("col_inv"))
+            .drop("col_inv")
+        )
 
     if m.is_linear:
         objective_write_linear_terms(f, df, print_variable)
@@ -200,6 +214,7 @@ def bounds_to_file(
     progress: bool = False,
     slice_size: int = 2_000_000,
     explicit_coordinate_names: bool = False,
+    scaling: ScalingContext | None = None,
 ) -> None:
     """
     Write out variables of a model to a lp file.
@@ -224,6 +239,24 @@ def bounds_to_file(
         var = m.variables[name]
         for var_slice in var.iterate_slices(slice_size):
             df = var_slice.to_polars()
+            if scaling is not None:
+                scale_map = pl.DataFrame(
+                    {
+                        "labels": pl.Series(scaling.matrices.vlabels),
+                        "col_scale": scaling.col_scale,
+                    }
+                )
+                df = (
+                    df.join(scale_map, on="labels", how="left")
+                    .with_columns(pl.col("col_scale").fill_null(1.0))
+                    .with_columns(
+                        [
+                            pl.col("lower") * pl.col("col_scale"),
+                            pl.col("upper") * pl.col("col_scale"),
+                        ]
+                    )
+                    .drop("col_scale")
+                )
 
             columns = [
                 pl.when(pl.col("lower") >= 0).then(pl.lit("+")).otherwise(pl.lit("")),
@@ -334,6 +367,7 @@ def constraints_to_file(
     lazy: bool = False,
     slice_size: int = 2_000_000,
     explicit_coordinate_names: bool = False,
+    scaling: ScalingContext | None = None,
 ) -> None:
     if not len(m.constraints):
         return
@@ -358,6 +392,37 @@ def constraints_to_file(
         for con_slice in con.iterate_slices(slice_size):
             df = con_slice.to_polars()
 
+            if scaling is not None:
+                row_map = pl.DataFrame(
+                    {
+                        "labels": pl.Series(scaling.matrices.clabels),
+                        "row_scale": scaling.row_scale,
+                    }
+                )
+                col_map = pl.DataFrame(
+                    {
+                        "vars": pl.Series(scaling.matrices.vlabels),
+                        "col_inv": scaling.col_inv,
+                    }
+                )
+                df = (
+                    df.join(row_map, on="labels", how="left")
+                    .join(col_map, on="vars", how="left")
+                    .with_columns(
+                        [
+                            pl.col("row_scale").fill_null(1.0),
+                            pl.col("col_inv").fill_null(1.0),
+                        ]
+                    )
+                    .with_columns(
+                        [
+                            pl.col("coeffs") * pl.col("row_scale") * pl.col("col_inv"),
+                            pl.col("rhs") * pl.col("row_scale"),
+                        ]
+                    )
+                    .drop(["row_scale", "col_inv"])
+                )
+
             if df.height == 0:
                 continue
 
@@ -428,12 +493,17 @@ def to_lp_file(
     slice_size: int = 2_000_000,
     progress: bool = True,
     explicit_coordinate_names: bool = False,
+    scaling: ScalingContext | None = None,
 ) -> None:
     with open(fn, mode="wb") as f:
         start = time.time()
 
         objective_to_file(
-            m, f, progress=progress, explicit_coordinate_names=explicit_coordinate_names
+            m,
+            f,
+            progress=progress,
+            explicit_coordinate_names=explicit_coordinate_names,
+            scaling=scaling,
         )
         constraints_to_file(
             m,
@@ -441,6 +511,7 @@ def to_lp_file(
             progress=progress,
             slice_size=slice_size,
             explicit_coordinate_names=explicit_coordinate_names,
+            scaling=scaling,
         )
         bounds_to_file(
             m,
@@ -448,6 +519,7 @@ def to_lp_file(
             progress=progress,
             slice_size=slice_size,
             explicit_coordinate_names=explicit_coordinate_names,
+            scaling=scaling,
         )
         binaries_to_file(
             m,
@@ -477,6 +549,7 @@ def to_file(
     slice_size: int = 2_000_000,
     progress: bool | None = None,
     explicit_coordinate_names: bool = False,
+    scaling: ScalingContext | None = None,
 ) -> Path:
     """
     Write out a model to a lp or mps file.
@@ -502,6 +575,7 @@ def to_file(
             slice_size=slice_size,
             progress=progress,
             explicit_coordinate_names=explicit_coordinate_names,
+            scaling=scaling,
         )
 
     elif io_api == "mps":
@@ -512,7 +586,10 @@ def to_file(
 
         # Use very fast highspy implementation
         # Might be replaced by custom writer, however needs C/Rust bindings for performance
-        h = m.to_highspy(explicit_coordinate_names=explicit_coordinate_names)
+        h = m.to_highspy(
+            explicit_coordinate_names=explicit_coordinate_names,
+            matrices=scaling.matrices if scaling else None,
+        )
         h.writeModel(str(fn))
     else:
         raise ValueError(
@@ -523,7 +600,10 @@ def to_file(
 
 
 def to_mosek(
-    m: Model, task: Any | None = None, explicit_coordinate_names: bool = False
+    m: Model,
+    task: Any | None = None,
+    explicit_coordinate_names: bool = False,
+    matrices: MatrixAccessor | ScaledMatrices | None = None,
 ) -> Any:
     """
     Export model to MOSEK.
@@ -552,7 +632,7 @@ def to_mosek(
     task.appendvars(m.nvars)
     task.appendcons(m.ncons)
 
-    M = m.matrices
+    M = m.matrices if matrices is None else matrices
     # for j, n in enumerate(("x" + M.vlabels.astype(str).astype(object))):
     #    task.putvarname(j, n)
 
@@ -635,7 +715,10 @@ def to_mosek(
 
 
 def to_gurobipy(
-    m: Model, env: Any | None = None, explicit_coordinate_names: bool = False
+    m: Model,
+    env: Any | None = None,
+    explicit_coordinate_names: bool = False,
+    matrices: MatrixAccessor | ScaledMatrices | None = None,
 ) -> Any:
     """
     Export the model to gurobipy.
@@ -662,7 +745,7 @@ def to_gurobipy(
     m.constraints.sanitize_missings()
     model = gurobipy.Model(env=env)
 
-    M = m.matrices
+    M = m.matrices if matrices is None else matrices
 
     names = np.vectorize(print_variable)(M.vlabels).astype(object)
     kwargs = {}
@@ -687,7 +770,11 @@ def to_gurobipy(
     return model
 
 
-def to_highspy(m: Model, explicit_coordinate_names: bool = False) -> Highs:
+def to_highspy(
+    m: Model,
+    explicit_coordinate_names: bool = False,
+    matrices: MatrixAccessor | ScaledMatrices | None = None,
+) -> Highs:
     """
     Export the model to highspy.
 
@@ -710,7 +797,7 @@ def to_highspy(m: Model, explicit_coordinate_names: bool = False) -> Highs:
         m, explicit_coordinate_names=explicit_coordinate_names
     )
 
-    M = m.matrices
+    M = m.matrices if matrices is None else matrices
     h = highspy.Highs()
     h.addVars(len(M.vlabels), M.lb, M.ub)
     if len(m.binaries) + len(m.integers):
diff --git a/linopy/model.py b/linopy/model.py
index 3982b84d..e99f5f9c 100644
--- a/linopy/model.py
+++ b/linopy/model.py
@@ -59,6 +59,7 @@
 from linopy.matrices import MatrixAccessor
 from linopy.objective import Objective
 from linopy.remote import OetcHandler, RemoteHandler
+from linopy.scaling import ScaleOptions, ScalingContext, resolve_options
 from linopy.solvers import (
     IO_APIS,
     NO_SOLUTION_FILE_SOLVERS,
@@ -1059,6 +1060,7 @@ def solve(
         slice_size: int = 2_000_000,
         remote: RemoteHandler | OetcHandler = None,  # type: ignore
         progress: bool | None = None,
+        scale: bool | str | ScaleOptions | None = None,
         **solver_options: Any,
     ) -> tuple[str, str]:
         """
@@ -1126,6 +1128,10 @@ def solve(
             Whether to show a progress bar of writing the lp file. The default is
             None, which means that the progress bar is shown if the model has more
             than 10000 variables and constraints.
+        scale : bool | str | ScaleOptions, optional
+            Enable pre-solve scaling. Use True for default row scaling,
+            a string ('row-max', 'row-l2') to choose the norm, or pass a
+            ScaleOptions instance for fine-grained control.
         **solver_options : kwargs
             Options passed to the solver.
 
@@ -1137,6 +1143,8 @@ def solve(
         """
         # clear cached matrix properties potentially present from previous solve commands
         self.matrices.clean_cached_properties()
+        scale_options = resolve_options(scale)
+        scaling: ScalingContext | None = None
 
         # check io_api
         if io_api is not None and io_api not in IO_APIS:
@@ -1145,6 +1153,10 @@ def solve(
             )
 
         if remote is not None:
+            if scale_options.enabled:
+                raise NotImplementedError(
+                    "Scaling is not supported for remote solving yet."
+                )
             if isinstance(remote, OetcHandler):
                 solved = remote.solve_on_oetc(self)
             else:
@@ -1208,6 +1220,10 @@ def solve(
         if sanitize_infinities:
             self.constraints.sanitize_infinities()
 
+        if scale_options.enabled:
+            scaling = ScalingContext.from_model(self.matrices, scale_options)
+        matrices_override = scaling.matrices if scaling is not None else None
+
         if self.is_quadratic and solver_name not in quadratic_solvers:
             raise ValueError(
                 f"Solver {solver_name} does not support quadratic problems."
@@ -1229,6 +1245,7 @@ def solve(
                     basis_fn=to_path(basis_fn),
                     env=env,
                     explicit_coordinate_names=explicit_coordinate_names,
+                    matrices=matrices_override,
                 )
             else:
                 if solver_name in ["glpk", "cbc"] and explicit_coordinate_names:
@@ -1242,6 +1259,7 @@ def solve(
                     explicit_coordinate_names=explicit_coordinate_names,
                     slice_size=slice_size,
                     progress=progress,
+                    scaling=scaling,
                 )
                 result = solver.solve_problem_from_file(
                     problem_fn=to_path(problem_fn),
@@ -1257,6 +1275,9 @@ def solve(
                 if fn is not None and (os.path.exists(fn) and not keep_files):
                     os.remove(fn)
 
+        if scaling is not None:
+            result = scaling.unscale_solution(result)
+
         result.info()
 
         self.objective._value = result.solution.objective
diff --git a/linopy/scaling.py b/linopy/scaling.py
new file mode 100644
index 00000000..622cfc54
--- /dev/null
+++ b/linopy/scaling.py
@@ -0,0 +1,262 @@
+"""
+Utilities for scaling linear problems before handing them to a solver.
+
+The scaling performed here is intentionally conservative: it supports
+row (constraint) scaling for all problem types and optional column
+(variable) scaling for continuous variables. Scaling is applied to the
+data structures used by the solver interfaces (`matrices`), leaving the
+original model untouched. Primal and dual values can be mapped back to
+the unscaled space via the provided helpers.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import Literal
+
+import numpy as np
+import pandas as pd
+from numpy import ndarray
+from scipy.sparse import csc_matrix, diags
+
+from linopy.common import set_int_index
+from linopy.constants import Result, Solution
+from linopy.matrices import MatrixAccessor
+
+ScaleMethod = Literal["row-max", "row-l2"]
+
+
+@dataclass
+class ScaleOptions:
+    """
+    Configuration for model scaling.
+
+    Attributes
+    ----------
+    enabled : bool
+        Whether scaling is enabled.
+    method : {"row-max", "row-l2"}
+        Norm used to calculate row/column magnitudes.
+    variable_scaling : bool
+        Enable column/variable scaling in addition to row scaling.
+    scale_integer_variables : bool
+        Whether integer/binary variables are scaled. Default False keeps
+        integrality-safe scaling by leaving them untouched.
+    scale_bounds : bool
+        Whether bounds are scaled along with variables. Required for
+        consistent variable scaling; ignored when `variable_scaling` is
+        False.
+    target : float
+        Target magnitude per row/column after scaling.
+    zero_floor : float
+        Lower bound for norms to avoid division by (near) zero.
+    """
+
+    enabled: bool = True
+    method: ScaleMethod = "row-max"
+    variable_scaling: bool = False
+    scale_integer_variables: bool = False
+    scale_bounds: bool = True
+    target: float = 1.0
+    zero_floor: float = 1e-12
+
+
+@dataclass
+class ScaledMatrices:
+    """
+    Container for scaled matrix data passed to solver frontends.
+    """
+
+    A: csc_matrix | None
+    b: ndarray
+    c: ndarray
+    lb: ndarray
+    ub: ndarray
+    sense: ndarray
+    vtypes: ndarray
+    vlabels: ndarray
+    clabels: ndarray
+    Q: csc_matrix | None
+
+
+def _row_norms(A: csc_matrix, method: ScaleMethod) -> ndarray:
+    """
+    Compute per-row magnitudes for a sparse matrix.
+    """
+    A_csr = A.tocsr()
+    if method == "row-l2":
+        norms = np.sqrt(np.array(A_csr.power(2).mean(axis=1)).ravel(), dtype=float)
+    else:
+        A_abs = A_csr.copy()
+        A_abs.data = np.abs(A_abs.data)
+        norms = np.array(A_abs.max(axis=1).toarray()).ravel().astype(float)
+
+    # rows without entries yield 0 or nan; keep them unscaled
+    norms = np.where(np.isnan(norms) | (norms == 0), 1.0, norms)
+    return norms
+
+
+def _col_norms(A: csc_matrix, method: ScaleMethod) -> ndarray:
+    """
+    Compute per-column magnitudes for a sparse matrix.
+    """
+    A_csc = A.tocsc()
+    if method == "row-l2":
+        norms = np.sqrt(np.array(A_csc.power(2).mean(axis=0)).ravel(), dtype=float)
+    else:
+        A_abs = A_csc.copy()
+        A_abs.data = np.abs(A_abs.data)
+        norms = np.array(A_abs.max(axis=0).toarray()).ravel().astype(float)
+
+    norms = np.where(np.isnan(norms) | (norms == 0), 1.0, norms)
+    return norms
+
+
+def _safe_scale(norms: ndarray, target: float, floor: float) -> ndarray:
+    norms = np.maximum(norms, floor)
+    return target / norms
+
+
+def resolve_options(scale: bool | str | ScaleOptions | None) -> ScaleOptions:
+    """
+    Normalize scale input into a ScaleOptions instance.
+    """
+    if scale is None or scale is False:
+        return ScaleOptions(enabled=False)
+    if scale is True:
+        return ScaleOptions()
+    if isinstance(scale, str):
+        if scale not in ("row-max", "row-l2"):
+            msg = f"Invalid scale method: {scale!r}. Must be 'row-max' or 'row-l2'."
+            raise ValueError(msg)
+        return ScaleOptions(method=scale)  # type: ignore[arg-type]
+    if isinstance(scale, ScaleOptions):
+        return scale
+
+    msg = "scale must be one of None, bool, str {'row-max','row-l2'} or ScaleOptions."
+    raise ValueError(msg)
+
+
+class ScalingContext:
+    """
+    Compute and apply scaling factors for a model's matrices.
+    """
+
+    def __init__(
+        self,
+        options: ScaleOptions,
+        row_scale: ndarray,
+        col_scale: ndarray,
+        matrices: ScaledMatrices,
+    ) -> None:
+        self.options = options
+        self.row_scale = row_scale
+        self.col_scale = col_scale
+        self.col_inv = np.reciprocal(col_scale)
+        self.matrices = matrices
+
+    @classmethod
+    def from_model(
+        cls,
+        matrices: MatrixAccessor,
+        options: ScaleOptions,
+    ) -> ScalingContext:
+        """
+        Build a scaling context from a model's matrix accessor.
+        """
+        if not options.enabled:
+            raise ValueError("ScalingContext.from_model called with scaling disabled.")
+
+        A = matrices.A
+
+        if A is None:
+            row_scale = np.ones_like(matrices.clabels, dtype=float)
+        else:
+            row_norm = _row_norms(A, options.method)
+            row_scale = _safe_scale(row_norm, options.target, options.zero_floor)
+
+        b = matrices.b * row_scale if len(matrices.b) else matrices.b
+
+        Q = matrices.Q
+        if options.variable_scaling and not options.scale_bounds:
+            raise ValueError(
+                "Variable scaling requires scaling bounds for consistency. "
+                "Set `scale_bounds=True` or disable `variable_scaling`."
+            )
+
+        if options.variable_scaling and A is not None:
+            A_rows = diags(row_scale) @ A
+            col_norm = _col_norms(A_rows, options.method)
+            col_scale = _safe_scale(col_norm, options.target, options.zero_floor)
+        else:
+            col_scale = np.ones_like(matrices.vlabels, dtype=float)
+
+        if not options.scale_integer_variables:
+            vtypes = matrices.vtypes
+            mask_int = (vtypes == "I") | (vtypes == "B")
+            col_scale[mask_int] = 1.0
+
+        col_inv = np.reciprocal(np.maximum(col_scale, options.zero_floor))
+
+        if A is not None:
+            A_scaled = diags(row_scale) @ A @ diags(col_inv)
+        else:
+            A_scaled = None
+
+        c = matrices.c * col_inv
+        lb = matrices.lb.copy()
+        ub = matrices.ub.copy()
+        if options.variable_scaling and options.scale_bounds:
+            lb = lb * col_scale
+            ub = ub * col_scale
+
+        if Q is not None:
+            D_inv = diags(col_inv)
+            Q = D_inv @ Q @ D_inv
+
+        scaled = ScaledMatrices(
+            A=A_scaled,
+            b=b,
+            c=c,
+            lb=lb,
+            ub=ub,
+            sense=matrices.sense,
+            vtypes=matrices.vtypes,
+            vlabels=matrices.vlabels,
+            clabels=matrices.clabels,
+            Q=Q,
+        )
+        return cls(
+            options=options, row_scale=row_scale, col_scale=col_scale, matrices=scaled
+        )
+
+    def unscale_primal(self, primal: pd.Series) -> pd.Series:
+        """
+        Map primal values from scaled to original space.
+        """
+        primal = set_int_index(primal)
+        factors = pd.Series(self.col_scale, index=self.matrices.vlabels)
+        return primal.div(factors, fill_value=np.nan)
+
+    def unscale_dual(self, dual: pd.Series) -> pd.Series:
+        """
+        Map dual values from scaled to original space.
+        """
+        dual = set_int_index(dual)
+        factors = pd.Series(self.row_scale, index=self.matrices.clabels)
+        return dual.mul(factors, fill_value=np.nan)
+
+    def unscale_solution(self, result: Result | None) -> Result | None:
+        """
+        Apply unscaling to primal and dual values stored in a solver result.
+        """
+        if result is None or result.solution is None:
+            return result
+
+        sol: Solution = result.solution
+        if not sol.primal.empty:
+            sol.primal = self.unscale_primal(sol.primal)
+        if not sol.dual.empty:
+            sol.dual = self.unscale_dual(sol.dual)
+        result.solution = sol
+        return result
diff --git a/linopy/solvers.py b/linopy/solvers.py
index a121a2b5..e16fc101 100644
--- a/linopy/solvers.py
+++ b/linopy/solvers.py
@@ -35,7 +35,9 @@
 if TYPE_CHECKING:
     import gurobipy
 
+    from linopy.matrices import MatrixAccessor
     from linopy.model import Model
+    from linopy.scaling import ScaledMatrices
 
 EnvType = TypeVar("EnvType")
 
@@ -269,6 +271,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: EnvType | None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         """
         Abstract method to solve a linear problem from a model.
@@ -308,6 +311,7 @@ def solve_problem(
         basis_fn: Path | None = None,
         env: EnvType | None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         """
         Solve a linear problem either from a model or a problem file.
@@ -328,6 +332,7 @@ def solve_problem(
                 basis_fn=basis_fn,
                 env=env,
                 explicit_coordinate_names=explicit_coordinate_names,
+                matrices=matrices,
             )
         elif problem_fn is not None:
             return self.solve_problem_from_file(
@@ -372,6 +377,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         msg = "Direct API not implemented for CBC"
         raise NotImplementedError(msg)
@@ -558,6 +564,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         msg = "Direct API not implemented for GLPK"
         raise NotImplementedError(msg)
@@ -738,6 +745,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         """
         Solve a linear problem directly from a linopy model using the Highs solver.
@@ -780,7 +788,9 @@ def solve_problem_from_model(
                 "Drop the solver option or use 'choose' to enable quadratic terms / integrality."
             )
 
-        h = model.to_highspy(explicit_coordinate_names=explicit_coordinate_names)
+        h = model.to_highspy(
+            explicit_coordinate_names=explicit_coordinate_names, matrices=matrices
+        )
         self._set_solver_params(h, log_fn)
 
         return self._solve(
@@ -979,6 +989,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: gurobipy.Env | dict[str, Any] | None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         """
         Solve a linear problem directly from a linopy model using the Gurobi solver.
@@ -1015,7 +1026,9 @@ def solve_problem_from_model(
                 env_ = env
 
             m = model.to_gurobipy(
-                env=env_, explicit_coordinate_names=explicit_coordinate_names
+                env=env_,
+                explicit_coordinate_names=explicit_coordinate_names,
+                matrices=matrices,
             )
 
             return self._solve(
@@ -1218,6 +1231,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         msg = "Direct API not implemented for Cplex"
         raise NotImplementedError(msg)
@@ -1370,6 +1384,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         msg = "Direct API not implemented for SCIP"
         raise NotImplementedError(msg)
@@ -1526,6 +1541,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         msg = "Direct API not implemented for Xpress"
         raise NotImplementedError(msg)
@@ -1675,6 +1691,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         """
         Solve a linear problem directly from a linopy model using the MOSEK solver.
@@ -1711,7 +1728,11 @@ def solve_problem_from_model(
                 stacklevel=2,
             )
         with mosek.Task() as m:
-            m = model.to_mosek(m, explicit_coordinate_names=explicit_coordinate_names)
+            m = model.to_mosek(
+                m,
+                explicit_coordinate_names=explicit_coordinate_names,
+                matrices=matrices,
+            )
 
             return self._solve(
                 m,
@@ -2009,6 +2030,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         msg = "Direct API not implemented for COPT"
         raise NotImplementedError(msg)
@@ -2150,6 +2172,7 @@ def solve_problem_from_model(
         basis_fn: Path | None = None,
         env: None = None,
         explicit_coordinate_names: bool = False,
+        matrices: MatrixAccessor | ScaledMatrices | None = None,
     ) -> Result:
         msg = "Direct API not implemented for MindOpt"
         raise NotImplementedError(msg)
diff --git a/test/test_optimization.py b/test/test_optimization.py
index da1050fc..90fe0349 100644
--- a/test/test_optimization.py
+++ b/test/test_optimization.py
@@ -20,6 +20,7 @@
 
 from linopy import GREATER_EQUAL, LESS_EQUAL, Model, solvers
 from linopy.common import to_path
+from linopy.scaling import ScaleOptions
 from linopy.solvers import _new_highspy_mps_layout, available_solvers, quadratic_solvers
 
 logger = logging.getLogger(__name__)
@@ -532,6 +533,99 @@ def test_non_aligned_variables(
             assert np.issubdtype(dtype, np.floating)
 
 
+def _build_scaling_model() -> Model:
+    m = Model()
+    x = m.add_variables(lower=0, name="x")
+    y = m.add_variables(lower=0, name="y")
+    m.add_constraints(1000 * x + 2 * y >= 10, name="c0")
+    m.add_constraints(0.1 * x + 0.5 * y >= 1, name="c1")
+    m.objective = 10 * x + y
+    return m
+
+
+@pytest.mark.skipif("highs" not in available_solvers, reason="Highs not installed")
+def test_scaling_integration_row_only() -> None:
+    base = _build_scaling_model()
+    status, _ = base.solve("highs", io_api="direct")
+    assert status == "ok"
+    base_solution = base.solution.to_pandas()
+    base_obj = base.objective.value or 0.0
+
+    scaled = _build_scaling_model()
+    status, _ = scaled.solve("highs", io_api="direct", scale=True)
+    assert status == "ok"
+    scaled_solution = scaled.solution.to_pandas()
+    scaled_obj = scaled.objective.value or 0.0
+
+    assert np.allclose(base_solution.values, scaled_solution.values)
+    assert np.isclose(base_obj, scaled_obj)
+
+
+@pytest.mark.skipif("highs" not in available_solvers, reason="Highs not installed")
+def test_scaling_integration_row_and_column() -> None:
+    base = _build_scaling_model()
+    status, _ = base.solve("highs", io_api="direct")
+    assert status == "ok"
+    base_solution = base.solution.to_pandas()
+    base_obj = base.objective.value or 0.0
+
+    scaled = _build_scaling_model()
+    status, _ = scaled.solve(
+        "highs",
+        io_api="direct",
+        scale=ScaleOptions(
+            enabled=True, variable_scaling=True, scale_integer_variables=False
+        ),
+    )
+    assert status == "ok"
+    scaled_solution = scaled.solution.to_pandas()
+    scaled_obj = scaled.objective.value or 0.0
+
+    assert np.allclose(base_solution.values, scaled_solution.values)
+    assert np.isclose(base_obj, scaled_obj)
+
+
+@pytest.mark.skipif("highs" not in available_solvers, reason="Highs not installed")
+def test_scaling_integration_lp_file_io() -> None:
+    """Test scaling with LP file-based IO to cover io.py scaling paths."""
+    base = _build_scaling_model()
+    status, _ = base.solve("highs", io_api="lp-polars")
+    assert status == "ok"
+    base_solution = base.solution.to_pandas()
+    base_obj = base.objective.value or 0.0
+
+    scaled = _build_scaling_model()
+    status, _ = scaled.solve(
+        "highs",
+        io_api="lp-polars",
+        scale=ScaleOptions(
+            enabled=True, variable_scaling=True, scale_integer_variables=False
+        ),
+    )
+    assert status == "ok"
+    scaled_solution = scaled.solution.to_pandas()
+    scaled_obj = scaled.objective.value or 0.0
+
+    assert np.allclose(base_solution.values, scaled_solution.values, rtol=1e-4)
+    assert np.isclose(base_obj, scaled_obj, rtol=1e-4)
+
+
+@pytest.mark.skipif("highs" not in available_solvers, reason="Highs not installed")
+def test_scaling_with_l2_norm() -> None:
+    """Test scaling with row-l2 method."""
+    base = _build_scaling_model()
+    status, _ = base.solve("highs", io_api="direct")
+    assert status == "ok"
+    base_obj = base.objective.value or 0.0
+
+    scaled = _build_scaling_model()
+    status, _ = scaled.solve("highs", io_api="direct", scale="row-l2")
+    assert status == "ok"
+    scaled_obj = scaled.objective.value or 0.0
+
+    assert np.isclose(base_obj, scaled_obj, rtol=1e-4)
+
+
 @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params)
 def test_set_files(
     tmp_path: Any,
diff --git a/test/test_scaling.py b/test/test_scaling.py
new file mode 100644
index 00000000..93fc3acc
--- /dev/null
+++ b/test/test_scaling.py
@@ -0,0 +1,185 @@
+import numpy as np
+import pandas as pd
+import pytest
+
+from linopy.constants import (
+    Result,
+    Solution,
+    SolverStatus,
+    Status,
+    TerminationCondition,
+)
+from linopy.model import Model
+from linopy.scaling import ScaleOptions, ScalingContext, resolve_options
+
+
+def _build_simple_model() -> Model:
+    m = Model()
+    x = m.add_variables(lower=0, name="x", coords=[range(2)])
+    # Row 0: max coeff 4 -> scale factor 0.25
+    m.add_constraints(2 * x.isel(dim_0=0) + 4 * x.isel(dim_0=1) == 8, name="row0")
+    # Row 1: max coeff 0.5 -> scale factor 2.0
+    m.add_constraints(0.5 * x.isel(dim_0=0) + 0.25 * x.isel(dim_0=1) >= 1, name="row1")
+    return m
+
+
+def test_row_scaling_and_dual_unscale() -> None:
+    m = _build_simple_model()
+    opts = ScaleOptions(enabled=True, variable_scaling=False, method="row-max")
+    ctx = ScalingContext.from_model(m.matrices, opts)
+    M = ctx.matrices
+
+    assert np.isclose(ctx.row_scale[0], 0.25)
+    assert np.isclose(ctx.row_scale[1], 2.0)
+    assert M.A is not None
+    A = np.asarray(M.A.todense())
+    # Row 0 scaled
+    assert np.isclose(A[0, 0], 0.5)
+    assert np.isclose(A[0, 1], 1.0)
+    # RHS scaled with the same factor
+    assert np.isclose(M.b[0], 2.0)
+
+    dual = pd.Series([1.0, 1.5], index=m.matrices.clabels, dtype=float)
+    unscaled_dual = ctx.unscale_dual(dual)
+    assert np.isclose(unscaled_dual.loc[m.matrices.clabels[0]], 0.25)
+    assert np.isclose(unscaled_dual.loc[m.matrices.clabels[1]], 3.0)
+
+
+def test_column_scaling_continuous_only() -> None:
+    m = Model()
+    x = m.add_variables(lower=0, name="x")
+    y = m.add_variables(lower=0, name="y")
+    b = m.add_variables(name="b", binary=True)
+
+    m.add_constraints(1000 * x + 2 * y + b >= 10, name="row0")
+
+    opts = ScaleOptions(
+        enabled=True,
+        variable_scaling=True,
+        scale_integer_variables=False,
+        method="row-max",
+    )
+    ctx = ScalingContext.from_model(m.matrices, opts)
+
+    # After row scaling: row factor = 1/1000; col norms become [1, 0.002, 0.001]
+    # Only continuous vars should be scaled.
+    assert np.isclose(ctx.col_scale[0], 1.0)
+    assert np.isclose(ctx.col_scale[1], 500.0)
+    # binary/int columns remain unscaled
+    assert np.isclose(ctx.col_scale[2], 1.0)
+
+    scaled = pd.Series([1.0, 5.0, 1.0], index=m.matrices.vlabels, dtype=float)
+    unscaled = ctx.unscale_primal(scaled)
+    assert np.isclose(unscaled.loc[m.matrices.vlabels[0]], 1.0)
+    assert np.isclose(unscaled.loc[m.matrices.vlabels[1]], 0.01)
+    assert np.isclose(unscaled.loc[m.matrices.vlabels[2]], 1.0)
+
+
+def test_row_l2_scaling() -> None:
+    m = _build_simple_model()
+    opts = ScaleOptions(enabled=True, variable_scaling=False, method="row-l2")
+    ctx = ScalingContext.from_model(m.matrices, opts)
+
+    # RMS of row 0: sqrt(mean([2^2, 4^2])) = sqrt(10) ≈ 3.162
+    # scale factor = 1 / sqrt(10)
+    assert np.isclose(ctx.row_scale[0], 1 / np.sqrt(10), rtol=1e-5)
+
+
+def test_column_scaling_with_l2_norm() -> None:
+    m = Model()
+    x = m.add_variables(lower=0, name="x")
+    y = m.add_variables(lower=0, name="y")
+    m.add_constraints(100 * x + 10 * y >= 10, name="row0")
+
+    opts = ScaleOptions(enabled=True, variable_scaling=True, method="row-l2")
+    ctx = ScalingContext.from_model(m.matrices, opts)
+
+    # After row scaling, column norms are computed with L2
+    assert ctx.col_scale[0] != 1.0
+    assert ctx.col_scale[1] != 1.0
+
+
+def test_resolve_options_variations() -> None:
+    # None -> disabled
+    opts = resolve_options(None)
+    assert not opts.enabled
+
+    # False -> disabled
+    opts = resolve_options(False)
+    assert not opts.enabled
+
+    # True -> enabled with defaults
+    opts = resolve_options(True)
+    assert opts.enabled
+    assert opts.method == "row-max"
+
+    # String -> enabled with method
+    opts = resolve_options("row-l2")
+    assert opts.enabled
+    assert opts.method == "row-l2"
+
+    # ScaleOptions passthrough
+    custom = ScaleOptions(enabled=True, variable_scaling=True)
+    opts = resolve_options(custom)
+    assert opts is custom
+
+
+def test_resolve_options_invalid_string() -> None:
+    with pytest.raises(ValueError, match="Invalid scale method"):
+        resolve_options("invalid-method")
+
+
+def test_resolve_options_invalid_type() -> None:
+    with pytest.raises(ValueError, match="scale must be one of"):
+        resolve_options(123)  # type: ignore[arg-type]
+
+
+def test_scaling_disabled_error() -> None:
+    m = _build_simple_model()
+    opts = ScaleOptions(enabled=False)
+    with pytest.raises(ValueError, match="scaling disabled"):
+        ScalingContext.from_model(m.matrices, opts)
+
+
+def test_variable_scaling_without_bounds_error() -> None:
+    m = _build_simple_model()
+    opts = ScaleOptions(enabled=True, variable_scaling=True, scale_bounds=False)
+    with pytest.raises(ValueError, match="scaling bounds"):
+        ScalingContext.from_model(m.matrices, opts)
+
+
+def test_unscale_solution_none_result() -> None:
+    m = _build_simple_model()
+    opts = ScaleOptions(enabled=True)
+    ctx = ScalingContext.from_model(m.matrices, opts)
+
+    # None result returns None
+    assert ctx.unscale_solution(None) is None
+
+
+def test_unscale_solution_none_solution() -> None:
+    m = _build_simple_model()
+    opts = ScaleOptions(enabled=True)
+    ctx = ScalingContext.from_model(m.matrices, opts)
+
+    # Result with None solution returns unchanged
+    result = Result(
+        Status(SolverStatus.warning, TerminationCondition.infeasible), solution=None
+    )
+    unscaled = ctx.unscale_solution(result)
+    assert unscaled is result
+    assert unscaled is not None and unscaled.solution is None
+
+
+def test_unscale_solution_empty_primal_dual() -> None:
+    m = _build_simple_model()
+    opts = ScaleOptions(enabled=True)
+    ctx = ScalingContext.from_model(m.matrices, opts)
+
+    # Empty primal/dual should not fail
+    sol = Solution(primal=pd.Series(dtype=float), dual=pd.Series(dtype=float))
+    result = Result(Status(SolverStatus.ok, TerminationCondition.optimal), solution=sol)
+    unscaled = ctx.unscale_solution(result)
+    assert unscaled is not None and unscaled.solution is not None
+    assert unscaled.solution.primal.empty
+    assert unscaled.solution.dual.empty