constraints: sanitize infinity values (#370)

* constraints: sanitize infinity values

* constraints: move check of invalid rhs to constraint assignment

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
 --------------

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.

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."

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)
 

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")

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"