diff --git a/CHANGELOG.md b/CHANGELOG.md
index e70cf08cb..3c4f60820 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -32,6 +32,7 @@ Until here -->
 ## [Unreleased]
 
 ### Added
+- Extra Check for HeatPumpWithSource.COP to be strictly > 1 to avoid division by zero
 
 ### Changed
 - Greatly improved docstrings and documentation of all public classes
@@ -45,6 +46,9 @@ Until here -->
 - Fix error handling in network visualization if networkx is not installed.
 - Fix broken links in docs.
 - Fix missing args in docstrings in `plotting.py`, `solvers.py`, and `core.py`.
+- Fix COP getter and setter of `HeatPumpWithSource` returning and setting wrong conversion factors.
+- Fix custom compression levels in `io.save_dataset_to_netcdf`
+- Fix `total_max` did not work when total min was not used.
 
 ### Known Issues
 
diff --git a/flixopt/features.py b/flixopt/features.py
index eb954944b..14122bb4a 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -616,7 +616,7 @@ def compute_consecutive_hours_in_state(
         if len(hours_per_timestep) < nr_of_indexes_with_consecutive_ones:
             raise ValueError(
                 f'When trying to calculate the consecutive duration, the length of the last duration '
-                f'({len(nr_of_indexes_with_consecutive_ones)}) is longer than the provided hours_per_timestep ({len(hours_per_timestep)}), '
+                f'({nr_of_indexes_with_consecutive_ones}) is longer than the provided hours_per_timestep ({len(hours_per_timestep)}), '
                 f'as {binary_values=}'
             )
 
@@ -945,7 +945,7 @@ def __init__(
 
         # Parameters
         self._shares_are_time_series = shares_are_time_series
-        self._total_max = total_max if total_min is not None else np.inf
+        self._total_max = total_max if total_max is not None else np.inf
         self._total_min = total_min if total_min is not None else -np.inf
         self._max_per_hour = max_per_hour if max_per_hour is not None else np.inf
         self._min_per_hour = min_per_hour if min_per_hour is not None else -np.inf
diff --git a/flixopt/io.py b/flixopt/io.py
index 35d927136..d08c0a1b5 100644
--- a/flixopt/io.py
+++ b/flixopt/io.py
@@ -234,7 +234,7 @@ def save_dataset_to_netcdf(
         path,
         encoding=None
         if not apply_encoding
-        else {data_var: {'zlib': True, 'complevel': 5} for data_var in ds.data_vars},
+        else {data_var: {'zlib': True, 'complevel': compression} for data_var in ds.data_vars},
     )
 
 
diff --git a/flixopt/linear_converters.py b/flixopt/linear_converters.py
index 7f3169c7d..6fdd30cda 100644
--- a/flixopt/linear_converters.py
+++ b/flixopt/linear_converters.py
@@ -570,15 +570,24 @@ def __init__(
         self.Q_ab = Q_ab
         self.Q_th = Q_th
 
+        if np.any(np.asarray(self.COP) <= 1):
+            raise ValueError(f'{self.label_full}.COP must be strictly > 1 for HeatPumpWithSource.')
+
     @property
     def COP(self):  # noqa: N802
-        return self.conversion_factors[0][self.Q_th.label]
+        return self.conversion_factors[0][self.P_el.label]
 
     @COP.setter
     def COP(self, value):  # noqa: N802
         check_bounds(value, 'COP', self.label_full, 1, 20)
-        self.conversion_factors[0][self.Q_th.label] = value
-        self.conversion_factors[1][self.Q_th.label] = value / (value - 1)
+        if np.any(np.asarray(self.COP) <= 1):
+            raise ValueError(f'{self.label_full}.COP must be strictly > 1 for HeatPumpWithSource.')
+        # electricity equation: COP * P_el == 1 * Q_th
+        self.conversion_factors[0][self.P_el.label] = value
+        self.conversion_factors[0][self.Q_th.label] = 1
+        # heat source equation: (COP/(COP-1)) * Q_ab == 1 * Q_th  -> Q_ab = Q_th*(COP-1)/COP
+        self.conversion_factors[1][self.Q_ab.label] = value / (value - 1)
+        self.conversion_factors[1][self.Q_th.label] = 1
 
 
 def check_bounds(
diff --git a/tests/test_on_hours_computation.py b/tests/test_on_hours_computation.py
index 38b340d57..cd0e637d0 100644
--- a/tests/test_on_hours_computation.py
+++ b/tests/test_on_hours_computation.py
@@ -42,7 +42,7 @@ def test_compute_duration(self, binary_values, hours_per_timestep, expected):
     )
     def test_compute_duration_raises_error(self, binary_values, hours_per_timestep):
         """Test error conditions."""
-        with pytest.raises(TypeError):
+        with pytest.raises(ValueError):
             ConsecutiveStateModel.compute_consecutive_hours_in_state(binary_values, hours_per_timestep)
 
 
@@ -67,7 +67,7 @@ class TestComputePreviousOnStates:
             ([np.array([0.1, 0, 0.3]), None, np.array([0, 0, 0])], np.array([1, 0, 1])),
             ([np.array([0, 0, 0]), np.array([0, 1, 0])], np.array([0, 1, 0])),
             ([np.array([0.1, 0, 0]), np.array([0, 0, 0.2])], np.array([1, 0, 1])),
-            # Case 6: Mix of None, 1D and 2D arrays
+            # Case 6: Mix of None and 1D arrays
             ([None, np.array([0, 0, 0]), np.array([0, 1, 0]), np.array([0, 0, 0])], np.array([0, 1, 0])),
             ([np.array([0, 0, 0]), None, np.array([0, 0, 0]), np.array([0, 0, 0])], np.array([0, 0, 0])),
         ],