diff --git a/doc/api.rst b/doc/api.rst
index a7612927..28fb0492 100644
--- a/doc/api.rst
+++ b/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
diff --git a/doc/user-guide.rst b/doc/user-guide.rst
index 8b7ee5bd..16df9a15 100644
--- a/doc/user-guide.rst
+++ b/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
 ------------------------------
 
diff --git a/linopy/common.py b/linopy/common.py
index 9ee9777d..16876dc5 100644
--- a/linopy/common.py
+++ b/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.
diff --git a/linopy/constraints.py b/linopy/constraints.py
index 96e2a843..f2402a57 100644
--- a/linopy/constraints.py
+++ b/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,
diff --git a/linopy/io.py b/linopy/io.py
index 1c4714c4..530007c9 100644
--- a/linopy/io.py
+++ b/linopy/io.py
@@ -24,7 +24,11 @@
 from tqdm import tqdm
 
 from linopy import solvers
-from linopy.common import to_polars
+from linopy.common import (
+    constraint_scaling_lookup,
+    to_polars,
+    variable_scaling_lookup,
+)
 from linopy.constants import CONCAT_DIM, SOS_DIM_ATTR, SOS_TYPE_ATTR
 from linopy.objective import Objective
 
@@ -89,6 +93,71 @@ def signed_number(expr: pl.Expr) -> tuple[pl.Expr, pl.Expr]:
     )
 
 
+def _lookup_positive_labels(lookup: np.ndarray, labels: np.ndarray) -> np.ndarray:
+    """Return lookup values for non-negative labels, using 1 for sentinels/nulls."""
+    labels = np.asarray(labels).reshape(-1)
+    values = np.ones(len(labels), dtype=float)
+    valid = ~pd.isna(labels)
+    if not valid.any():
+        return values
+    valid_labels = labels[valid].astype(np.intp, copy=False)
+    positive = valid_labels >= 0
+    valid_positions = np.flatnonzero(valid)
+    if positive.any():
+        values[valid_positions[positive]] = lookup[valid_labels[positive]]
+    return values
+
+
+def _scale_objective_dataframe(
+    df: pl.DataFrame, variable_scaling: np.ndarray, objective_scaling: float
+) -> pl.DataFrame:
+    """Apply column scaling and row-like objective scaling to objective terms."""
+    if df.is_empty():
+        return df
+
+    if "vars" in df.columns:
+        scales = _lookup_positive_labels(variable_scaling, df["vars"].to_numpy())
+    else:
+        scale1 = _lookup_positive_labels(variable_scaling, df["vars1"].to_numpy())
+        scale2 = _lookup_positive_labels(variable_scaling, df["vars2"].to_numpy())
+        scales = scale1 * scale2
+
+    return df.with_columns(
+        (pl.col("coeffs") / pl.Series(scales) * objective_scaling).alias("coeffs")
+    )
+
+
+def _scale_bounds_dataframe(
+    df: pl.DataFrame, variable_scaling: np.ndarray
+) -> pl.DataFrame:
+    """Multiply bounds by solver-side variable scaling factors."""
+    if df.is_empty():
+        return df
+    scales = _lookup_positive_labels(variable_scaling, df["labels"].to_numpy())
+    scale_series = pl.Series(scales)
+    return df.with_columns(
+        (pl.col("lower") * scale_series).alias("lower"),
+        (pl.col("upper") * scale_series).alias("upper"),
+    )
+
+
+def _scale_constraint_dataframe(
+    df: pl.DataFrame,
+    variable_scaling: np.ndarray,
+    constraint_scaling: np.ndarray,
+) -> pl.DataFrame:
+    """Divide by column scaling and multiply by row scaling."""
+    if df.is_empty():
+        return df
+    row_scales = _lookup_positive_labels(constraint_scaling, df["labels"].to_numpy())
+    var_scales = _lookup_positive_labels(variable_scaling, df["vars"].to_numpy())
+    row_scale_series = pl.Series(row_scales)
+    return df.with_columns(
+        (pl.col("coeffs") / pl.Series(var_scales) * row_scale_series).alias("coeffs"),
+        (pl.col("rhs") * row_scale_series).alias("rhs"),
+    )
+
+
 def format_coord(coord: str) -> str:
     from linopy.common import format_coord
 
@@ -213,15 +282,18 @@ def objective_to_file(
     print_variable, _ = get_printers(
         m, explicit_coordinate_names=explicit_coordinate_names
     )
+    variable_scaling = variable_scaling_lookup(m)
 
     sense = m.objective.sense
     f.write(f"{sense}\n\nobj:\n\n".encode())
     df = m.objective.to_polars()
 
     if m.is_linear:
+        df = _scale_objective_dataframe(df, variable_scaling, m.objective.scaling)
         objective_write_linear_terms(f, df, print_variable)
 
     elif m.is_quadratic:
+        df = _scale_objective_dataframe(df, variable_scaling, m.objective.scaling)
         linear_terms = df.filter(pl.col("vars1").eq(-1) | pl.col("vars2").eq(-1))
         linear_terms = linear_terms.with_columns(
             pl.when(pl.col("vars1").eq(-1))
@@ -259,6 +331,7 @@ def bounds_to_file(
     print_variable, _ = get_printers(
         m, explicit_coordinate_names=explicit_coordinate_names
     )
+    variable_scaling = variable_scaling_lookup(m)
 
     f.write(b"\n\nbounds\n\n")
     if progress:
@@ -272,6 +345,7 @@ def bounds_to_file(
         var = m.variables[name]
         for var_slice in var.iterate_slices(slice_size):
             df = var_slice.to_polars()
+            df = _scale_bounds_dataframe(df, variable_scaling)
 
             columns = [
                 *signed_number(pl.col("lower")),
@@ -482,6 +556,7 @@ def indicator_constraints_to_file(
     print_variable_scalar, _ = get_printers_scalar(
         m, explicit_coordinate_names=explicit_coordinate_names
     )
+    variable_scaling = variable_scaling_lookup(m)
 
     for con in m.constraints.indicator.data.values():
         ic_data = con.data
@@ -494,25 +569,31 @@ def indicator_constraints_to_file(
         vars_flat = ic_data.vars.values.reshape(len(labels_flat), -1)
         sign_flat = np.broadcast_to(ic_data.sign.values, labels_flat.shape).flatten()
         rhs_flat = np.broadcast_to(ic_data.rhs.values, labels_flat.shape).flatten()
+        scaling_flat = np.broadcast_to(
+            ic_data.scaling.values, labels_flat.shape
+        ).flatten()
 
         for i in range(len(labels_flat)):
             if labels_flat[i] == -1:
                 continue
 
+            row_scale = float(scaling_flat[i])
             bvar_name = print_variable_scalar(binary_var_flat[i : i + 1])[0]
             valid = vars_flat[i] != -1
             var_names = print_variable_scalar(vars_flat[i][valid])
 
             terms = []
-            for coeff, var_name in zip(coeffs_flat[i][valid], var_names):
-                coeff = float(coeff)
+            for coeff, var_label, var_name in zip(
+                coeffs_flat[i][valid], vars_flat[i][valid], var_names
+            ):
+                coeff = float(coeff) / variable_scaling[int(var_label)] * row_scale
                 prefix = "+" if coeff >= 0 else ""
                 terms.append(f"{prefix}{coeff} {var_name}")
 
             lhs_str = " ".join(terms)
             line = (
                 f"ic{labels_flat[i]}: {bvar_name} = {int(binary_val_flat[i])} -> "
-                f"{lhs_str} {sign_flat[i]} {float(rhs_flat[i])}\n"
+                f"{lhs_str} {sign_flat[i]} {float(rhs_flat[i]) * row_scale}\n"
             )
             f.write(line.encode())
 
@@ -532,6 +613,8 @@ def constraints_to_file(
     print_variable, print_constraint = get_printers(
         m, explicit_coordinate_names=explicit_coordinate_names
     )
+    variable_scaling = variable_scaling_lookup(m)
+    constraint_scaling = constraint_scaling_lookup(m)
 
     f.write(b"\n\ns.t.\n\n")
     names = list(regular)
@@ -548,6 +631,7 @@ def constraints_to_file(
         con = regular[name]
         for con_slice in con.iterate_slices(slice_size):
             df = con_slice.to_polars()
+            df = _scale_constraint_dataframe(df, variable_scaling, constraint_scaling)
 
             if df.height == 0:
                 continue
@@ -962,7 +1046,9 @@ def with_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset:
         for name, con in m.constraints.items()
     ]
     objective = m.objective.data
-    objective = objective.assign_attrs(sense=m.objective.sense)
+    objective = objective.assign_attrs(
+        sense=m.objective.sense, scaling=m.objective.scaling
+    )
     if m.objective.value is not None:
         objective = objective.assign_attrs(value=m.objective.value)
     obj = [with_prefix(objective, "objective")]
@@ -1091,7 +1177,10 @@ def get_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset:
 
     objective = get_prefix(ds, "objective")
     m.objective = Objective(
-        LinearExpression(objective, m), m, objective.attrs.pop("sense")
+        LinearExpression(objective, m),
+        m,
+        objective.attrs.pop("sense"),
+        objective.attrs.pop("scaling", 1),
     )
     m.objective._value = objective.attrs.pop("value", None)
 
@@ -1201,7 +1290,9 @@ def _copy_con_data(con: ConstraintBase) -> xr.Dataset:
     )
 
     obj_expr = LinearExpression(m.objective.expression.data.copy(deep=deep), new_model)
-    new_model._objective = Objective(obj_expr, new_model, m.objective.sense)
+    new_model._objective = Objective(
+        obj_expr, new_model, m.objective.sense, m.objective.scaling
+    )
     new_model._objective._value = (
         float(m.objective.value)
         if (include_solution and m.objective.value is not None)
diff --git a/linopy/matrices.py b/linopy/matrices.py
index 0feba9c9..6a286380 100644
--- a/linopy/matrices.py
+++ b/linopy/matrices.py
@@ -15,6 +15,7 @@
 from numpy import ndarray
 
 from linopy import expressions
+from linopy.common import constraint_scaling_lookup, variable_scaling_lookup
 from linopy.constraints import CSRConstraint
 
 if TYPE_CHECKING:
@@ -58,6 +59,12 @@ def _build_vars(self) -> None:
         m = self._parent
         label_index = m.variables.label_index
         self.vlabels: ndarray = label_index.vlabels
+        var_scaling_by_label = variable_scaling_lookup(m)
+        self.var_scaling: ndarray = (
+            var_scaling_by_label[self.vlabels]
+            if len(self.vlabels)
+            else np.array([], dtype=float)
+        )
 
         lb_list = []
         ub_list = []
@@ -81,8 +88,8 @@ def _build_vars(self) -> None:
             vtypes_list.append(np.full(mask.sum(), vtype))
 
         if lb_list:
-            self.lb: ndarray = np.concatenate(lb_list)
-            self.ub: ndarray = np.concatenate(ub_list)
+            self.lb: ndarray = np.concatenate(lb_list) * self.var_scaling
+            self.ub: ndarray = np.concatenate(ub_list) * self.var_scaling
             self.vtypes: ndarray = np.concatenate(vtypes_list)
         else:
             self.lb = np.array([])
@@ -93,13 +100,34 @@ def _build_cons(self) -> None:
         m = self._parent
         label_index = m.variables.label_index
         label_to_pos = label_index.label_to_pos
+        con_scaling_by_label = constraint_scaling_lookup(m)
+
+        def scale_rows_and_cols(
+            csr: scipy.sparse.csr_array, con_labels: np.ndarray, b: np.ndarray
+        ) -> tuple[scipy.sparse.csr_array, np.ndarray]:
+            if csr.shape[0] == 0:
+                return csr, b
+            row_scaling = con_scaling_by_label[con_labels]
+            # With solver variables y = Scol * x, constraints A x = b become
+            # Srow * A * Scol^-1 * y = Srow * b.
+            csr = cast(
+                scipy.sparse.csr_array,
+                csr.multiply(row_scaling[:, np.newaxis]).tocsr(),
+            )
+            if csr.shape[1] and len(self.var_scaling):
+                csr = cast(
+                    scipy.sparse.csr_array,
+                    csr.multiply(1 / self.var_scaling[np.newaxis, :]).tocsr(),
+                )
+            return csr, b * row_scaling
 
         reg_csrs, reg_b, reg_sense = [], [], []
         ind_csrs, ind_b, ind_sense, ind_binvar, ind_binval = [], [], [], [], []
         for c in m.constraints.data.values():
             if c.is_indicator:
                 cc = c if isinstance(c, CSRConstraint) else c.freeze()
-                csr, _, b, sense = cc.to_matrix_with_rhs(label_index)
+                csr, con_labels, b, sense = cc.to_matrix_with_rhs(label_index)
+                csr, b = scale_rows_and_cols(csr, con_labels, b)
                 ind_csrs.append(csr)
                 ind_b.append(b)
                 ind_sense.append(sense)
@@ -107,7 +135,8 @@ def _build_cons(self) -> None:
                 binval = cast("int | np.ndarray", cc._binval)
                 ind_binval.append(_binval_per_row(binval, len(b)))
             else:
-                csr, _, b, sense = c.to_matrix_with_rhs(label_index)
+                csr, con_labels, b, sense = c.to_matrix_with_rhs(label_index)
+                csr, b = scale_rows_and_cols(csr, con_labels, b)
                 reg_csrs.append(csr)
                 reg_b.append(b)
                 reg_sense.append(sense)
@@ -144,7 +173,10 @@ def c(self) -> ndarray:
             coeffs = expr.data.coeffs.values.ravel()
 
         mask = var_labels != -1
-        np.add.at(result, label_to_pos[var_labels[mask]], coeffs[mask])
+        positions = label_to_pos[var_labels[mask]]
+        scaled_coeffs = coeffs[mask] / self.var_scaling[positions]
+        scaled_coeffs = scaled_coeffs * m.objective.scaling
+        np.add.at(result, positions, scaled_coeffs)
         return result
 
     @cached_property
@@ -154,7 +186,13 @@ def Q(self) -> scipy.sparse.csc_matrix | None:
         expr = m.objective.expression
         if not isinstance(expr, expressions.QuadraticExpression):
             return None
-        return expr.to_matrix()[self.vlabels][:, self.vlabels]
+        q = expr.to_matrix()[self.vlabels][:, self.vlabels]
+        if q.nnz:
+            # Quadratic coefficients get one inverse column scaling per factor.
+            q = q.multiply(1 / self.var_scaling[:, np.newaxis])
+            q = q.multiply(1 / self.var_scaling[np.newaxis, :])
+            q = q * self._parent.objective.scaling
+        return q.tocsc()
 
     @cached_property
     def sol(self) -> ndarray:
diff --git a/linopy/model.py b/linopy/model.py
index 884d59db..69cd734b 100644
--- a/linopy/model.py
+++ b/linopy/model.py
@@ -34,6 +34,8 @@
     maybe_replace_signs,
     replace_by_map,
     to_path,
+    validate_scaling,
+    variable_solver_scaling,
 )
 from linopy.constants import (
     GREATER_EQUAL,
@@ -589,11 +591,12 @@ def _check_valid_dim_names(self, ds: DataArray | Dataset) -> None:
         -------
         None.
         """
-        unsupported_dim_names = ["labels", "coeffs", "vars", "sign", "rhs"]
+        unsupported_dim_names = ["labels", "coeffs", "vars", "sign", "rhs", "scaling"]
         if any(dim in unsupported_dim_names for dim in ds.dims):
             raise ValueError(
                 "Added data contains unsupported dimension names. "
-                "Dimensions cannot be named 'labels', 'coeffs', 'vars', 'sign' or 'rhs'."
+                "Dimensions cannot be named 'labels', 'coeffs', 'vars', 'sign', "
+                "'rhs' or 'scaling'."
             )
 
     def add_variables(
@@ -606,6 +609,7 @@ def add_variables(
         binary: bool = False,
         integer: bool = False,
         semi_continuous: bool = False,
+        scaling: Any = 1,
         **kwargs: Any,
     ) -> Variable:
         """
@@ -653,6 +657,13 @@ def add_variables(
             Boolean mask with False values for variables which are skipped.
             The shape of the mask has to match the shape the added variables.
             Default is None.
+        scaling : float/array_like, optional
+            Positive finite scaling factor(s) used when exporting the model to
+            a solver. Continuous and semi-continuous solver variables are
+            represented as ``scaling * variable``, so coefficients in those
+            columns are divided by ``scaling`` and bounds are multiplied by it.
+            Binary and integer variables keep their ordinary discrete solver
+            columns. The default is 1.
         binary : bool
             Whether the new variable is a binary variable which are used for
             Mixed-Integer problems.
@@ -786,10 +797,14 @@ def add_variables(
 
         lower_da = broadcast_to_coords(lower, coords, label="lower bound", **kwargs)
         upper_da = broadcast_to_coords(upper, coords, label="upper bound", **kwargs)
+        scaling_da = broadcast_to_coords(
+            scaling, coords, label="variable scaling", **kwargs
+        )
         data = Dataset(
             {
                 "lower": lower_da,
                 "upper": upper_da,
+                "scaling": validate_scaling(scaling_da, "variable scaling"),
                 "labels": -1,
             }
         )
@@ -985,6 +1000,7 @@ def add_constraints(
         coords: Sequence[Sequence | pd.Index] | Mapping | None = ...,
         mask: MaskLike | None = ...,
         freeze: Literal[False] = ...,
+        scaling: ConstantLike = ...,
     ) -> Constraint: ...
 
     @overload
@@ -1001,6 +1017,7 @@ def add_constraints(
         coords: Sequence[Sequence | pd.Index] | Mapping | None = ...,
         mask: MaskLike | None = ...,
         freeze: Literal[True] = ...,
+        scaling: ConstantLike = ...,
     ) -> CSRConstraint: ...
 
     def add_constraints(
@@ -1016,6 +1033,7 @@ def add_constraints(
         coords: Sequence[Sequence | pd.Index] | Mapping | None = None,
         mask: MaskLike | None = None,
         freeze: bool | None = None,
+        scaling: ConstantLike = 1,
     ) -> ConstraintBase:
         """
         Assign a new, possibly multi-dimensional array of constraints to the
@@ -1054,6 +1072,10 @@ def add_constraints(
             If True, convert the constraint to an immutable CSR-backed CSRConstraint
             for better memory efficiency. If None, uses the model default
             ``Model.freeze_constraints`` setting (default False).
+        scaling : float/array_like, optional
+            Positive finite scaling factor(s) for constraint rows. Solver-side
+            left-hand-side coefficients and right-hand-side values are multiplied
+            by this factor. The default is 1.
 
         Returns
         -------
@@ -1098,6 +1120,14 @@ def add_constraints(
 
         data["labels"] = -1
         (data,) = xr.broadcast(data, exclude=[TERM_DIM])
+        row_coords = data.labels.coords
+        data = assign_multiindex_safe(
+            data,
+            scaling=validate_scaling(
+                broadcast_to_coords(scaling, row_coords, label="constraint scaling"),
+                "constraint scaling",
+            ),
+        )
 
         if mask is not None:
             mask = broadcast_to_coords(mask, data.coords, label="mask").astype(bool)
@@ -1143,6 +1173,7 @@ def add_indicator_constraints(
         sign: SignLike | None = None,
         rhs: ConstantLike | None = None,
         name: str | None = None,
+        scaling: ConstantLike = 1,
     ) -> ConstraintBase:
         """
         Add indicator constraints to the model.
@@ -1170,6 +1201,10 @@ def add_indicator_constraints(
             Right-hand side. Required when ``lhs`` is an expression.
         name : str, optional
             Name for the indicator constraint group.
+        scaling : float/array_like, optional
+            Positive finite scaling factor(s) for indicator rows. Solver-side
+            left-hand-side coefficients and right-hand-side values are multiplied
+            by this factor. The default is 1.
 
         Returns
         -------
@@ -1196,6 +1231,16 @@ def add_indicator_constraints(
 
         data["labels"] = -1
         (data,) = xr.broadcast(data, exclude=[TERM_DIM])
+        row_coords = data.labels.coords
+        data = assign_multiindex_safe(
+            data,
+            scaling=validate_scaling(
+                broadcast_to_coords(
+                    scaling, row_coords, label="indicator constraint scaling"
+                ),
+                "indicator constraint scaling",
+            ),
+        )
 
         data = self._allocate_constraint_labels(data, name)
 
@@ -1217,6 +1262,7 @@ def add_objective(
         | Sequence[tuple[ConstantLike, VariableLike]],
         overwrite: bool = False,
         sense: str = "min",
+        scaling: ConstantLike = 1,
     ) -> None:
         """
         Add an objective function to the model.
@@ -1227,6 +1273,12 @@ def add_objective(
             Expression describing the objective function.
         overwrite : False, optional
             Whether to overwrite the existing objective. The default is False.
+        sense : "min" or "max", optional
+            Objective sense. The default is "min".
+        scaling : float, optional
+            Positive finite scaling factor for the objective. Solver-side
+            objective coefficients are multiplied by this factor and reported
+            objective values are divided back. The default is 1.
 
         Returns
         -------
@@ -1242,6 +1294,7 @@ def add_objective(
             expr = 1 * expr
         self.objective.expression = expr
         self.objective.sense = sense
+        self.objective.scaling = scaling
 
     def remove_variables(self, name: str) -> None:
         """
@@ -1969,7 +2022,7 @@ def assign_result(
         result.info()
 
         if result.solution is not None:
-            self.objective._value = result.solution.objective
+            self.objective._value = result.solution.objective / self.objective.scaling
 
         status_value = result.status.status.value
         termination_condition = result.status.termination_condition.value
@@ -1985,9 +2038,11 @@ def assign_result(
         primal = result.solution.primal
         for _, var in self.variables.items():
             start, end = var.range
-            var.solution = xr.DataArray(
-                primal[start:end].reshape(var.shape), var.coords
+            values = (
+                primal[start:end].reshape(var.shape)
+                / variable_solver_scaling(var).values
             )
+            var.solution = xr.DataArray(values, var.coords)
 
         if len(result.solution.dual):
             dual = result.solution.dual
@@ -1996,9 +2051,12 @@ def assign_result(
                     continue
                 start, end = con.range
                 coords = {dim: con.coords[dim] for dim in con.coord_dims}
-                con.dual = xr.DataArray(
-                    dual[start:end].reshape(con.shape), coords, dims=con.coord_dims
+                values = (
+                    dual[start:end].reshape(con.shape)
+                    * con.scaling.values
+                    / self.objective.scaling
                 )
+                con.dual = xr.DataArray(values, coords, dims=con.coord_dims)
 
         return status_value, termination_condition
 
diff --git a/linopy/objective.py b/linopy/objective.py
index a51b2207..2be6f5ff 100644
--- a/linopy/objective.py
+++ b/linopy/objective.py
@@ -38,7 +38,10 @@ def _objwrap(obj: Objective, *args: Any, **kwargs: Any) -> Objective:
         for k, v in new_default_kwargs.items():
             kwargs.setdefault(k, v)
         return obj.__class__(
-            method(obj.expression, *default_args, *args, **kwargs), obj.model, obj.sense
+            method(obj.expression, *default_args, *args, **kwargs),
+            obj.model,
+            obj.sense,
+            obj.scaling,
         )
 
     _objwrap.__doc__ = (
@@ -55,7 +58,7 @@ class Objective:
     An objective expression containing all relevant information.
     """
 
-    __slots__ = ("_expression", "_model", "_sense", "_value")
+    __slots__ = ("_expression", "_model", "_sense", "_scaling", "_value")
     __array_ufunc__ = None
     __array_priority__ = 10000
     __pandas_priority__ = 10000
@@ -67,11 +70,13 @@ def __init__(
         expression: expressions.LinearExpression | expressions.QuadraticExpression,
         model: Model,
         sense: str = "min",
+        scaling: ConstantLike = 1,
     ) -> None:
         self._model: Model = model
         self._value: float | None = None
 
         self.sense: str = sense
+        self.scaling = scaling
         self.expression: (
             expressions.LinearExpression | expressions.QuadraticExpression
         ) = expression
@@ -217,6 +222,22 @@ def sense(self, sense: str) -> None:
             raise ValueError("Invalid sense. Must be 'min' or 'max'.")
         self._sense = sense
 
+    @property
+    def scaling(self) -> float:
+        """
+        Returns the positive objective scaling factor.
+        """
+        return self._scaling
+
+    @scaling.setter
+    def scaling(self, scaling: ConstantLike) -> None:
+        if not isinstance(scaling, int | float | np.floating | np.integer):
+            raise TypeError("Objective scaling must be a numeric scalar.")
+        scaling = float(scaling)
+        if not np.isfinite(scaling) or scaling <= 0:
+            raise ValueError("Objective scaling must be finite and positive.")
+        self._scaling = scaling
+
     @property
     def value(self) -> float64 | float | None:
         """
@@ -256,24 +277,26 @@ def __add__(
         if isinstance(expr, Objective):
             expr = expr.expression
 
-        return Objective(self.expression + expr, self.model, self.sense)
+        return Objective(self.expression + expr, self.model, self.sense, self.scaling)
 
     def __sub__(self, expr: LinearExpression | Objective) -> Objective:
         if not isinstance(expr, Objective):
             expr = Objective(expr, self.model)
-        return Objective(self.expression - expr.expression, self.model, self.sense)
+        return Objective(
+            self.expression - expr.expression, self.model, self.sense, self.scaling
+        )
 
     def __mul__(self, expr: ConstantLike) -> Objective:
         # only allow scalar multiplication
         if not isinstance(expr, int | float | np.floating | np.integer):
             raise ValueError("Invalid type for multiplication.")
-        return Objective(self.expression * expr, self.model, self.sense)
+        return Objective(self.expression * expr, self.model, self.sense, self.scaling)
 
     def __neg__(self) -> Objective:
-        return Objective(-self.expression, self.model, self.sense)
+        return Objective(-self.expression, self.model, self.sense, self.scaling)
 
     def __truediv__(self, expr: ConstantLike) -> Objective:
         # only allow scalar division
         if not isinstance(expr, int | float | np.floating | np.integer):
             raise ValueError("Invalid type for division.")
-        return Objective(self.expression / expr, self.model, self.sense)
+        return Objective(self.expression / expr, self.model, self.sense, self.scaling)
diff --git a/linopy/variables.py b/linopy/variables.py
index 0eacfdc4..4063aa3f 100644
--- a/linopy/variables.py
+++ b/linopy/variables.py
@@ -53,6 +53,7 @@
     set_int_index,
     to_dataframe,
     to_polars,
+    validate_scaling,
 )
 from linopy.config import options
 from linopy.constants import (
@@ -85,7 +86,7 @@
 
 logger = logging.getLogger(__name__)
 
-FILL_VALUE = {"labels": -1, "lower": np.nan, "upper": np.nan}
+FILL_VALUE = {"labels": -1, "lower": np.nan, "upper": np.nan, "scaling": 1.0}
 
 
 def varwrap(
@@ -197,6 +198,14 @@ def __init__(
         data = data.assign_attrs(name=name)
         if not skip_broadcast:
             (data,) = broadcast(data)
+        if "scaling" not in data:
+            data = assign_multiindex_safe(
+                data, scaling=DataArray(1.0).broadcast_like(data.labels)
+            )
+        data = assign_multiindex_safe(
+            data,
+            scaling=validate_scaling(data.scaling, f"scaling for variable '{name}'"),
+        )
         for attr in ("lower", "upper"):
             # convert to float, important for  operations like "shift"
             if not issubdtype(data[attr].dtype, floating):
@@ -311,6 +320,24 @@ def to_pandas(self) -> pd.Series:
         """
         return self.labels.to_pandas()
 
+    @property
+    def scaling(self) -> DataArray:
+        """Scaling factor supplied for this variable."""
+        return self.data.scaling
+
+    @scaling.setter
+    def scaling(self, value: ConstantLike) -> None:
+        value = broadcast_to_coords(
+            value,
+            coords=self.coords,
+            dims=self.dims,
+            label=f"scaling for variable '{self.name}'",
+            strict=False,
+        )
+        value = value.reindex_like(self.labels, fill_value=1.0)
+        value = validate_scaling(value, f"scaling for variable '{self.name}'")
+        self._data = assign_multiindex_safe(self.data, scaling=value)
+
     def to_linexpr(
         self,
         coefficient: ConstantLike = 1,
@@ -1009,7 +1036,8 @@ def flat(self) -> DataFrame:
         Convert the variable to a pandas DataFrame.
 
         The resulting DataFrame represents a long table format of the variable
-        with columns `labels`, `lower`, `upper` which are not masked.
+        with columns ``labels``, ``lower``, ``upper`` and ``scaling`` which are
+        not masked.
 
         Returns
         -------
@@ -1028,8 +1056,8 @@ def to_polars(self) -> pl.DataFrame:
         """
         Convert all variables to a single polars DataFrame.
 
-        The resulting dataframe is a long format of the variables
-        with columns `labels`, `lower`, 'upper` and `mask`.
+        The resulting dataframe is a long format of the variables with columns
+        ``labels``, ``lower``, ``upper`` and ``scaling``.
 
         Returns
         -------
@@ -1135,9 +1163,24 @@ def where(
         elif isinstance(other, Variable):
             _other = other.data
         elif isinstance(other, ScalarVariable):
-            _other = {"labels": other.label, "lower": other.lower, "upper": other.upper}
-        elif isinstance(other, dict | Dataset):
+            _other = {
+                "labels": other.label,
+                "lower": other.lower,
+                "upper": other.upper,
+                "scaling": other.scaling,
+            }
+        elif isinstance(other, dict):
+            _other = dict(other)
+            _other.setdefault("scaling", self._fill_value["scaling"])
+        elif isinstance(other, Dataset):
             _other = other
+            if "scaling" not in _other:
+                _other = assign_multiindex_safe(
+                    _other,
+                    scaling=DataArray(self._fill_value["scaling"]).broadcast_like(
+                        _other.labels
+                    ),
+                )
         else:
             raise ValueError(
                 f"other must be a Variable, ScalarVariable, dict or Dataset, got {type(other)}"
@@ -1479,8 +1522,8 @@ class Variables:
     _label_position_index: LabelPositionIndex | None = None
     _variable_label_index: VariableLabelIndex | None = None
 
-    dataset_attrs = ["labels", "lower", "upper"]
-    dataset_names = ["Labels", "Lower bounds", "Upper bounds"]
+    dataset_attrs = ["labels", "lower", "upper", "scaling"]
+    dataset_names = ["Labels", "Lower bounds", "Upper bounds", "Scaling factors"]
 
     def _formatted_names(self) -> dict[str, str]:
         """
@@ -1657,6 +1700,13 @@ def upper(self) -> Dataset:
         """
         return save_join(*[v.upper.rename(k) for k, v in self.items()])
 
+    @property
+    def scaling(self) -> Dataset:
+        """
+        Get the scaling factors of all variables.
+        """
+        return save_join(*[v.scaling.rename(k) for k, v in self.items()])
+
     @property
     def nvars(self) -> int:
         """
@@ -1941,8 +1991,8 @@ def flat(self) -> pd.DataFrame:
         """
         Convert all variables to a single pandas Dataframe.
 
-        The resulting dataframe is a long format of the variables
-        with columns `labels`, `lower`, 'upper` and `mask`.
+        The resulting dataframe is a long format of the variables with columns
+        ``labels``, ``lower``, ``upper``, ``scaling`` and ``key``.
 
         Returns
         -------
@@ -2029,6 +2079,14 @@ def upper(self) -> float:
         name, position = self.model.variables.get_label_position(self.label)
         return self.model.variables[name].upper.sel(position).item()
 
+    @property
+    def scaling(self) -> float:
+        """
+        Get the scaling factor of the variable.
+        """
+        name, position = self.model.variables.get_label_position(self.label)
+        return self.model.variables[name].scaling.sel(position).item()
+
     @property
     def model(self) -> Model:
         """
diff --git a/test/test_scaling.py b/test/test_scaling.py
new file mode 100644
index 00000000..14687170
--- /dev/null
+++ b/test/test_scaling.py
@@ -0,0 +1,194 @@
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+import pytest
+import xarray as xr
+
+from linopy import Model, read_netcdf
+from linopy.constants import Result, Solution, Status
+from linopy.solvers import available_solvers
+
+
+def test_variable_constraint_and_objective_scaling_in_matrices() -> None:
+    m = Model()
+    i = pd.Index(["a", "b"], name="i")
+    scaling = xr.DataArray([10.0, 100.0], coords=[i])
+
+    x = m.add_variables(
+        lower=xr.DataArray([1.0, 2.0], coords=[i]),
+        upper=xr.DataArray([3.0, 4.0], coords=[i]),
+        name="x",
+        scaling=scaling,
+    )
+    b = m.add_variables(binary=True, name="b", scaling=50.0)
+
+    row_scaling = xr.DataArray([2.0, 4.0], coords=[i])
+    m.add_constraints(
+        2 * x + 3 * b,
+        ">=",
+        xr.DataArray([20.0, 40.0], coords=[i]),
+        name="c",
+        scaling=row_scaling,
+    )
+    m.add_objective((5 * x).sum() + 7 * b, scaling=10.0)
+
+    matrices = m.matrices
+
+    np.testing.assert_allclose(matrices.lb, [10.0, 200.0, 0.0])
+    np.testing.assert_allclose(matrices.ub, [30.0, 400.0, 1.0])
+    np.testing.assert_allclose(matrices.b, [40.0, 160.0])
+    np.testing.assert_allclose(
+        matrices.A.toarray(),
+        [
+            [2 / 10 * 2, 0.0, 3 * 2],
+            [0.0, 2 / 100 * 4, 3 * 4],
+        ],
+    )
+    np.testing.assert_allclose(matrices.c, [5 / 10 * 10, 5 / 100 * 10, 7 * 10])
+
+    assert float(b.scaling) == 50.0
+
+
+def test_quadratic_objective_scaling_in_matrix() -> None:
+    m = Model()
+    i = pd.Index([0, 1], name="i")
+    x = m.add_variables(coords=[i], name="x", scaling=[2.0, 4.0])
+
+    m.add_objective((x * x).sum(), scaling=10.0)
+
+    np.testing.assert_allclose(
+        m.matrices.Q.toarray(),
+        np.diag([2 / 2 / 2 * 10, 2 / 4 / 4 * 10]),
+    )
+
+
+def test_indicator_constraint_scaling_in_matrix() -> None:
+    m = Model()
+    i = pd.Index(["a", "b"], name="i")
+    x = m.add_variables(coords=[i], name="x", scaling=[10.0, 100.0])
+    b = m.add_variables(binary=True, name="b")
+
+    rhs = xr.DataArray([20.0, 40.0], coords=[i])
+    scaling = xr.DataArray([2.0, 4.0], coords=[i])
+    con = m.add_indicator_constraints(
+        b,
+        1,
+        2 * x,
+        "<=",
+        rhs,
+        name="ic",
+        scaling=scaling,
+    )
+
+    np.testing.assert_allclose(con.scaling.values, [2.0, 4.0])
+    np.testing.assert_allclose(m.matrices.indicator_b, [40.0, 160.0])
+    np.testing.assert_allclose(
+        m.matrices.indicator_A.toarray(),
+        [
+            [2 / 10 * 2, 0.0, 0.0],
+            [0.0, 2 / 100 * 4, 0.0],
+        ],
+    )
+
+
+def test_assign_result_unscales_solution_objective_and_dual() -> None:
+    m = Model()
+    i = pd.Index(["a", "b"], name="i")
+    x = m.add_variables(coords=[i], name="x", scaling=[10.0, 100.0])
+    b = m.add_variables(binary=True, name="b", scaling=50.0)
+
+    m.add_constraints(x + b >= 1, name="c", scaling=[2.0, 4.0])
+    m.add_objective(x.sum() + b, scaling=10.0)
+
+    primal = np.full(m._xCounter, np.nan)
+    primal[x.labels.values.ravel()] = [20.0, 300.0]
+    primal[int(b.labels)] = 1.0
+    dual = np.full(m._cCounter, np.nan)
+    dual[m.constraints["c"].labels.values.ravel()] = [4.0, 8.0]
+    result = Result(
+        status=Status.from_termination_condition("optimal"),
+        solution=Solution(primal=primal, dual=dual, objective=12.3),
+        solver_name="mock",
+    )
+
+    m.assign_result(result)
+
+    np.testing.assert_allclose(x.solution.values, [2.0, 3.0])
+    assert float(b.solution) == 1.0
+    assert m.objective.value == pytest.approx(1.23)
+    np.testing.assert_allclose(m.constraints["c"].dual.values, [0.8, 3.2])
+
+
+def test_scaling_is_preserved_in_netcdf(tmp_path) -> None:
+    m = Model()
+    i = pd.Index(["a", "b"], name="i")
+    x = m.add_variables(coords=[i], name="x", scaling=[10.0, 100.0])
+    m.add_constraints(x >= 1, name="c", scaling=[2.0, 4.0], freeze=True)
+    m.add_objective(x.sum(), scaling=5.0)
+
+    path = tmp_path / "scaled.nc"
+    m.to_netcdf(path)
+    restored = read_netcdf(path)
+
+    np.testing.assert_allclose(restored.variables["x"].scaling.values, [10.0, 100.0])
+    np.testing.assert_allclose(restored.constraints["c"].scaling.values, [2.0, 4.0])
+    assert restored.objective.scaling == 5.0
+
+
+def test_lp_export_uses_scaled_values(tmp_path) -> None:
+    m = Model()
+    x = m.add_variables(lower=0, upper=10, name="x", scaling=10.0)
+    m.add_constraints(2 * x >= 20, name="c", scaling=2.0)
+    m.add_objective(5 * x, scaling=5.0)
+
+    path = tmp_path / "scaled.lp"
+    m.to_file(path, io_api="lp", progress=False)
+    text = path.read_text()
+
+    assert "2.5 x0" in text
+    assert ">= 40.0" in text
+    assert "<= +100.0" in text
+
+
+def test_scaling_validation() -> None:
+    m = Model()
+    with pytest.raises(ValueError, match="finite positive"):
+        m.add_variables(name="x", scaling=0)
+
+    x = m.add_variables(name="x")
+    with pytest.raises(ValueError, match="finite positive"):
+        m.add_constraints(x >= 1, scaling=-1)
+
+    with pytest.raises(ValueError, match="finite and positive"):
+        m.add_objective(x, scaling=np.inf)
+
+
+def test_constraint_scaling_setter_broadcasts_to_rows() -> None:
+    m = Model()
+    i = pd.Index(["a", "b"], name="i")
+    x = m.add_variables(coords=[i], name="x")
+    con = m.add_constraints(x >= 1, name="c")
+
+    con.scaling = 2.0
+
+    np.testing.assert_allclose(con.scaling.values, [2.0, 2.0])
+    np.testing.assert_allclose(m.matrices.b, [2.0, 2.0])
+
+
+@pytest.mark.skipif("highs" not in available_solvers, reason="HiGHS is not installed")
+def test_scaled_solve_preserves_user_units() -> None:
+    reference = Model()
+    rx = reference.add_variables(lower=0, name="x")
+    reference.add_constraints(rx >= 3)
+    reference.add_objective(rx)
+    reference.solve("highs", io_api="direct", output_flag=False)
+
+    scaled = Model()
+    sx = scaled.add_variables(lower=0, name="x", scaling=10.0)
+    scaled.add_constraints(sx >= 3, scaling=2.0)
+    scaled.add_objective(sx, scaling=5.0)
+    scaled.solve("highs", io_api="direct", output_flag=False)
+
+    assert float(sx.solution) == pytest.approx(float(rx.solution))
+    assert scaled.objective.value == pytest.approx(reference.objective.value)