diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index d79eb96c..1512ee1d 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -54,6 +54,8 @@ Upcoming Version
 **Bug fixes**
 
 * A multi-key ``groupby`` now returns its groups sorted by key tuple, like the single-key path. The key combinations were numbered by iterating a ``set``, so the group order was arbitrary and changed between processes with ``PYTHONHASHSEED``.
+* ``Model.copy`` (and the ``copy.copy``/``copy.deepcopy`` protocols) no longer downgrades a quadratic objective to a linear one. The objective was rebuilt as a ``LinearExpression`` regardless of its type, so the copy silently solved a different problem. ``linopy.testing.assert_model_equal`` now compares the objective by expression type as well, which it could not do before. (`#903 <https://github.com/PyPSA/linopy/issues/903>`__)
+* ``Model.to_netcdf``/``linopy.read_netcdf`` downgraded a quadratic objective the same way; the expression type is now stored alongside the objective and restored on read. Files written by earlier versions are read as before. (`#903 <https://github.com/PyPSA/linopy/issues/903>`__)
 * ``Solver.close()`` no longer leaves dangling native handles behind. The solver model is now dropped before the environment that owns it, instead of after. And the COPT and MindOpt file interfaces no longer hand back a model they already disposed: after a file-based COPT or MindOpt solve, ``model.solver_model`` is ``None`` rather than a handle into freed memory. (`#899 <https://github.com/PyPSA/linopy/pull/899>`__)
 
 **Breaking Changes**
diff --git a/linopy/io.py b/linopy/io.py
index e6196d47..b5ea8bc1 100644
--- a/linopy/io.py
+++ b/linopy/io.py
@@ -40,6 +40,7 @@
 logger = logging.getLogger(__name__)
 
 NETCDF_VERSION_ATTR = "_linopy_version"
+EXPR_TYPE_ATTR = "_linopy_expr_type"
 
 
 ufunc_kwargs = dict(vectorize=True)
@@ -986,13 +987,15 @@ def with_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset:
     ]
     exprs = [
         with_prefix(
-            expr.data.assign_attrs(name=name, _linopy_expr_type=expr.type),
+            expr.data.assign_attrs(name=name, **{EXPR_TYPE_ATTR: expr.type}),
             f"expressions-{name}",
         )
         for name, expr in m.expressions.items()
     ]
     objective = m.objective.data
-    objective = objective.assign_attrs(sense=m.objective.sense)
+    objective = objective.assign_attrs(
+        sense=m.objective.sense, **{EXPR_TYPE_ATTR: m.objective.expression.type}
+    )
     if m.objective.value is not None:
         objective = objective.assign_attrs(value=m.objective.value)
     obj = [with_prefix(objective, "objective")]
@@ -1114,7 +1117,7 @@ def get_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset:
     for k in sorted(expr_names):
         name = remove_prefix(k, "expressions")
         expr_ds = get_prefix(ds, k)
-        expr_type = expr_ds.attrs.pop("_linopy_expr_type", None)
+        expr_type = expr_ds.attrs.pop(EXPR_TYPE_ATTR, None)
         expr_ds.attrs.pop("name", None)  # re-attached below, after construction
         expr: LinearExpression | QuadraticExpression
         if expr_type == "QuadraticExpression" or (
@@ -1141,9 +1144,14 @@ def get_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset:
     m._constraints = Constraints(constraints, m)
 
     objective = get_prefix(ds, "objective")
-    m.objective = Objective(
-        LinearExpression(objective, m), m, objective.attrs.pop("sense")
+    obj_type = objective.attrs.pop(EXPR_TYPE_ATTR, None)
+    obj_cls: type[LinearExpression] | type[QuadraticExpression] = (
+        QuadraticExpression
+        if obj_type == "QuadraticExpression"
+        or (obj_type is None and FACTOR_DIM in objective.dims)
+        else LinearExpression
     )
+    m.objective = Objective(obj_cls(objective, m), m, objective.attrs.pop("sense"))
     m.objective._value = objective.attrs.pop("value", None)
 
     m.parameters = get_prefix(ds, "parameters")
@@ -1267,7 +1275,9 @@ def _copy_con_data(con: ConstraintBase) -> xr.Dataset:
         new_model,
     )
 
-    obj_expr = LinearExpression(m.objective.expression.data.copy(deep=deep), new_model)
+    obj_expr = type(m.objective.expression)(
+        m.objective.expression.data.copy(deep=deep), new_model
+    )
     new_model._objective = Objective(obj_expr, new_model, m.objective.sense)
     new_model._objective._value = (
         float(m.objective.value)
diff --git a/linopy/testing.py b/linopy/testing.py
index d9c67f7b..5fd16778 100644
--- a/linopy/testing.py
+++ b/linopy/testing.py
@@ -130,7 +130,7 @@ def assert_model_equal(a: Model, b: Model) -> None:
     for e in a.expressions:
         assert_exprequal(a.expressions[e], b.expressions[e])
 
-    assert_linequal(a.objective.expression, b.objective.expression)
+    assert_exprequal(a.objective.expression, b.objective.expression, check_name=False)
     assert a.objective.sense == b.objective.sense
     assert a.objective.value == b.objective.value
 
diff --git a/test/test_io.py b/test/test_io.py
index 7c37a7dd..50f53b76 100644
--- a/test/test_io.py
+++ b/test/test_io.py
@@ -277,6 +277,20 @@ def test_model_to_netcdf_quadratic_expression(
     assert_exprequal(m.expressions["quad"], p.expressions["quad"])
 
 
+def test_model_to_netcdf_quadratic_objective(tmp_path: Path) -> None:
+    """A quadratic objective survives the round-trip as a QP, not an LP."""
+    m = Model()
+    x = m.add_variables(4, pd.Series([8, 10]), name="x")
+    m.add_objective((x * x + 2 * x).sum())
+
+    fn = tmp_path / "test.nc"
+    m.to_netcdf(fn)
+    p = read_netcdf(fn)
+
+    assert isinstance(p.objective.expression, QuadraticExpression)
+    assert_model_equal(m, p)
+
+
 def test_model_to_netcdf_expression_dash_name(
     model_with_expressions: Model, tmp_path: Path
 ) -> None:
diff --git a/test/test_model.py b/test/test_model.py
index f55f280c..de0b3cb6 100644
--- a/test/test_model.py
+++ b/test/test_model.py
@@ -7,6 +7,7 @@
 
 import copy as pycopy
 import weakref
+from collections.abc import Callable
 from pathlib import Path
 from tempfile import gettempdir
 
@@ -394,6 +395,26 @@ def test_copy_model_with_expressions(
     assert m.expressions["lin"].coeffs.values.flat[0] == original_coeff
 
 
+@pytest.mark.parametrize(
+    "copy_fn",
+    [lambda m: m.copy(), pycopy.copy, pycopy.deepcopy],
+    ids=["copy", "shallowcopy", "deepcopy"],
+)
+def test_model_copy_preserves_quadratic_objective(
+    copy_fn: Callable[[Model], Model],
+) -> None:
+    """Copying a model keeps its objective quadratic instead of linearizing it."""
+    m: Model = Model()
+    x = m.add_variables(lower=-5, upper=5, name="x")
+    m.add_objective(x * x + 2 * x)
+
+    c = copy_fn(m)
+
+    assert type(c.objective.expression) is type(m.objective.expression)
+    assert c.type == m.type == "QP"
+    assert_model_equal(m, c)
+
+
 @pytest.mark.skipif(not available_solvers, reason="No solver installed")
 class TestModelCopySolved:
     def test_model_deepcopy_protocol_excludes_solution(