diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index 9c2aa99e..be256a60 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -5,6 +5,7 @@ Upcoming Version
 ----------------
 
 * When writing out an LP file, large variables and constraints are now chunked to avoid memory issues. This is especially useful for large models with constraints with many terms. The chunk size can be set with the `slice_size` argument in the `solve` function.
+* Constraints which of the form `<= infinity` and `>= -infinity` are now automatically filtered out when solving. The `solve` function now has a new argument `sanitize_infinities` to control this feature. Default is set to `True`.
 
 Version 0.3.15
 --------------
diff --git a/linopy/constraints.py b/linopy/constraints.py
index 458fa419..caf7c4ce 100644
--- a/linopy/constraints.py
+++ b/linopy/constraints.py
@@ -51,7 +51,14 @@
     to_polars,
 )
 from linopy.config import options
-from linopy.constants import EQUAL, HELPER_DIMS, TERM_DIM, SIGNS_pretty
+from linopy.constants import (
+    EQUAL,
+    GREATER_EQUAL,
+    HELPER_DIMS,
+    LESS_EQUAL,
+    TERM_DIM,
+    SIGNS_pretty,
+)
 from linopy.types import ConstantLike
 
 if TYPE_CHECKING:
@@ -851,17 +858,17 @@ def equalities(self) -> "Constraints":
         """
         return self[[n for n, s in self.items() if (s.sign == EQUAL).all()]]
 
-    def sanitize_zeros(self):
+    def sanitize_zeros(self) -> None:
         """
         Filter out terms with zero and close-to-zero coefficient.
         """
-        for name in list(self):
+        for name in self:
             not_zero = abs(self[name].coeffs) > 1e-10
             constraint = self[name]
             constraint.vars = self[name].vars.where(not_zero, -1)
             constraint.coeffs = self[name].coeffs.where(not_zero)
 
-    def sanitize_missings(self):
+    def sanitize_missings(self) -> None:
         """
         Set constraints labels to -1 where all variables in the lhs are
         missing.
@@ -872,6 +879,19 @@ def sanitize_missings(self):
                 contains_non_missing, -1
             )
 
+    def sanitize_infinities(self) -> None:
+        """
+        Replace infinite values in the constraints with a large value.
+        """
+        for name in self:
+            constraint = self[name]
+            valid_infinity_values = (
+                (constraint.sign == LESS_EQUAL) & (constraint.rhs == np.inf)
+            ) | ((constraint.sign == GREATER_EQUAL) & (constraint.rhs == -np.inf))
+            self[name].data["labels"] = self[name].labels.where(
+                ~valid_infinity_values, -1
+            )
+
     def get_name_by_label(self, label: Union[int, float]) -> str:
         """
         Get the constraint name of the constraint containing the passed label.
diff --git a/linopy/model.py b/linopy/model.py
index d1f34015..160a16e3 100644
--- a/linopy/model.py
+++ b/linopy/model.py
@@ -32,7 +32,14 @@
     replace_by_map,
     to_path,
 )
-from linopy.constants import HELPER_DIMS, TERM_DIM, ModelStatus, TerminationCondition
+from linopy.constants import (
+    GREATER_EQUAL,
+    HELPER_DIMS,
+    LESS_EQUAL,
+    TERM_DIM,
+    ModelStatus,
+    TerminationCondition,
+)
 from linopy.constraints import AnonymousScalarConstraint, Constraint, Constraints
 from linopy.expressions import (
     LinearExpression,
@@ -583,6 +590,12 @@ def add_constraints(
                 f"Invalid type of `lhs` ({type(lhs)}) or invalid combination of `lhs`, `sign` and `rhs`."
             )
 
+        invalid_infinity_values = (
+            (data.sign == LESS_EQUAL) & (data.rhs == -np.inf)
+        ) | ((data.sign == GREATER_EQUAL) & (data.rhs == np.inf))  # noqa: F821
+        if invalid_infinity_values.any():
+            raise ValueError(f"Constraint {name} contains incorrect infinite values.")
+
         # ensure helper dimensions are not set as coordinates
         if drop_dims := set(HELPER_DIMS).intersection(data.coords):
             # TODO: add a warning here, routines should be safe against this
@@ -953,6 +966,7 @@ def solve(
         keep_files: bool = False,
         env: None = None,
         sanitize_zeros: bool = True,
+        sanitize_infinities: bool = True,
         slice_size: int = 2_000_000,
         remote: None = None,
         **solver_options,
@@ -1003,6 +1017,8 @@ def solve(
             Whether to set terms with zero coefficient as missing.
             This will remove unneeded overhead in the lp file writing.
             The default is True.
+        sanitize_infinities : bool, optional
+            Whether to filter out constraints that are subject to `<= inf` or `>= -inf`.
         slice_size : int, optional
             Size of the slice to use for writing the lp file. The slice size
             is used to split large variables and constraints into smaller
@@ -1083,6 +1099,9 @@ def solve(
         if sanitize_zeros:
             self.constraints.sanitize_zeros()
 
+        if sanitize_infinities:
+            self.constraints.sanitize_infinities()
+
         if self.is_quadratic and solver_name not in quadratic_solvers:
             raise ValueError(
                 f"Solver {solver_name} does not support quadratic problems."
diff --git a/linopy/solvers.py b/linopy/solvers.py
index c9d13065..d478c5ac 100644
--- a/linopy/solvers.py
+++ b/linopy/solvers.py
@@ -1493,7 +1493,7 @@ def get_solver_solution() -> Solution:
                 dual_ = [str(d) for d in m.getConstraint()]
                 dual = pd.Series(m.getDual(dual_), index=dual_, dtype=float)
                 dual = set_int_index(dual)
-            except (xpress.SolverError, SystemError):
+            except (xpress.SolverError, xpress.ModelError, SystemError):
                 logger.warning("Dual values of MILP couldn't be parsed")
                 dual = pd.Series(dtype=float)
 
diff --git a/test/test_constraints.py b/test/test_constraints.py
index e7876642..d13acf8c 100644
--- a/test/test_constraints.py
+++ b/test/test_constraints.py
@@ -181,3 +181,25 @@ def test_constraints_flat():
 
     assert isinstance(m.constraints.flat, pd.DataFrame)
     assert not m.constraints.flat.empty
+
+
+def test_sanitize_infinities():
+    m = Model()
+
+    lower = xr.DataArray(np.zeros((10, 10)), coords=[range(10), range(10)])
+    upper = xr.DataArray(np.ones((10, 10)), coords=[range(10), range(10)])
+    x = m.add_variables(lower, upper, name="x")
+    y = m.add_variables(name="y")
+
+    # Test correct infinities
+    m.add_constraints(x <= np.inf, name="con_inf")
+    m.add_constraints(y >= -np.inf, name="con_neg_inf")
+    m.constraints.sanitize_infinities()
+    assert (m.constraints["con_inf"].labels == -1).all()
+    assert (m.constraints["con_neg_inf"].labels == -1).all()
+
+    # Test incorrect infinities
+    with pytest.raises(ValueError):
+        m.add_constraints(x >= np.inf, name="con_wrong_inf")
+    with pytest.raises(ValueError):
+        m.add_constraints(y <= -np.inf, name="con_wrong_neg_inf")
diff --git a/test/test_optimization.py b/test/test_optimization.py
index c5e040c8..e0f6fc0a 100644
--- a/test/test_optimization.py
+++ b/test/test_optimization.py
@@ -498,9 +498,7 @@ def test_infeasible_model(model, solver, io_api):
             model.compute_infeasibilities()
 
 
-@pytest.mark.parametrize(
-    "solver,io_api", [p for p in params if p[0] not in ["glpk", "cplex", "mindopt"]]
-)
+@pytest.mark.parametrize("solver,io_api", params)
 def test_model_with_inf(model_with_inf, solver, io_api):
     status, condition = model_with_inf.solve(solver, io_api=io_api)
     assert condition == "optimal"