fix: use xarray.Dataset copy instead of constructor

doc/release_notes.rst

@@ -27,6 +27,7 @@ Upcoming Version
 * Add ``format_labels()`` on ``Constraints``/``Variables`` and ``format_infeasibilities()`` on ``Model`` that return strings instead of printing to stdout, allowing usage with logging, storage, or custom output handling. Deprecate ``print_labels()`` and ``print_infeasibilities()``.
 * Add ``fix()``, ``unfix()``, and ``fixed`` to ``Variable`` and ``Variables`` for fixing variables to values via equality constraints. Supports automatic rounding for integer/binary variables.
 * Add ``relax()``, ``unrelax()``, and ``relaxed`` to ``Variable`` and ``Variables`` for LP relaxation of integer/binary variables. Supports partial relaxation via filtered views (e.g. ``m.variables.integers.relax()``). Semi-continuous variables raise ``NotImplementedError``.
+* Add compatibility to latest xarray version.
 
 
 Version 0.6.6

linopy/expressions.py

@@ -284,7 +284,7 @@ def sum(self, use_fallback: bool = False, **kwargs: Any) -> LinearExpression:
                 index.names = [str(col) for col in orig_group.columns]
                 index.name = GROUP_DIM
                 new_coords = Coordinates.from_pandas_multiindex(index, GROUP_DIM)
-                ds = xr.Dataset(ds.assign_coords(new_coords))
+                ds = ds.assign_coords(new_coords).copy()
 
             ds = ds.rename({GROUP_DIM: final_group_name})
             return LinearExpression(ds, self.model)
@@ -391,8 +391,10 @@ def __init__(self, data: Dataset | Any | None, model: Model) -> None:
         coeffs_vars_dict = {str(k): v for k, v in coeffs_vars.items()}
         data = assign_multiindex_safe(data, **coeffs_vars_dict)
 
-        # transpose with new Dataset to really ensure correct order
-        data = Dataset(data.transpose(..., TERM_DIM))
+        # Ensure correct dimension order and materialize a new Dataset.
+        # Wrapping an existing Dataset in Dataset(...) is no longer supported
+        # by newer xarray versions.
+        data = data.transpose(..., TERM_DIM).copy()
 
         # ensure helper dimensions are not set as coordinates
         if drop_dims := set(HELPER_DIMS).intersection(data.coords):
@@ -2098,7 +2100,7 @@ def __init__(self, data: Dataset | None, model: Model) -> None:
             raise ValueError(f"Size of dimension {FACTOR_DIM} must be 2.")
 
         # transpose data to have _term as last dimension and _factor as second last
-        data = xr.Dataset(data.transpose(..., FACTOR_DIM, TERM_DIM))
+        data = data.transpose(..., FACTOR_DIM, TERM_DIM).copy()
         self._data = data
 
     @property
@@ -2384,7 +2386,10 @@ def merge(
     has_quad_expression = any(type(e) is QuadraticExpression for e in exprs)
     has_linear_expression = any(type(e) is LinearExpression for e in exprs)
     if cls is None:
-        cls = QuadraticExpression if has_quad_expression else LinearExpression
+        cls = cast(
+            type[GenericExpression],
+            QuadraticExpression if has_quad_expression else LinearExpression,
+        )
 
     if cls is QuadraticExpression and dim == TERM_DIM and has_linear_expression:
         raise ValueError(

linopy/model.py

@@ -359,7 +359,10 @@ def parameters(self, value: Dataset | Mapping) -> None:
         """
         Set the parameters of the model.
         """
-        self._parameters = Dataset(value)
+        if isinstance(value, Dataset):
+            self._parameters = value.copy()
+        else:
+            self._parameters = Dataset(value)
 
     @property
     def solution(self) -> Dataset:

test/test_optimization.py

@@ -897,8 +897,17 @@ def test_modified_model(
     if solver not in feasible_mip_solvers:
         pytest.skip(f"{solver} does not support MIP")
 
+    # HiGHS 1.14 presolve has a bug where it incorrectly fixes binary variables
+    # when a continuous variable has a lower bound below (but near) the constraint
+    # RHS. Disabling presolve gives the correct optimal solution.
+    # See https://github.com/PyPSA/linopy/issues/648
+    extra_kwargs: dict = {"presolve": "off"} if solver == "highs" else {}
+
     status, condition = modified_model.solve(
-        solver, io_api=io_api, explicit_coordinate_names=explicit_coordinate_names
+        solver,
+        io_api=io_api,
+        explicit_coordinate_names=explicit_coordinate_names,
+        **extra_kwargs,
     )
     assert condition == "optimal"
     assert (modified_model.solution.x == 0).all()