Support indicator, SOS and semi-continuous constraints in SCIP

doc/release_notes.rst

@@ -4,6 +4,11 @@ Release Notes
 Upcoming Version
 ----------------
 
+* The SCIP solver now accepts indicator constraints, SOS constraints and
+  semi-continuous variables. SCIP already read all three from the problem
+  file; it rewrites them into rows and variables of its own, which the
+  solution reader now leaves out of the returned solution and duals.
+
 
 *Strict "v1" arithmetic semantics (opt-in)*
 

linopy/solvers.py

@@ -2500,6 +2500,9 @@ class SCIP(Solver[None]):
             SolverFeature.LP_FILE_NAMES,
             SolverFeature.READ_MODEL_FROM_FILE,
             SolverFeature.SOLUTION_FILE_NOT_NEEDED,
+            SolverFeature.SOS_CONSTRAINTS,
+            SolverFeature.INDICATOR_CONSTRAINTS,
+            SolverFeature.SEMI_CONTINUOUS_VARIABLES,
         }
     )
 
@@ -2583,9 +2586,25 @@ def _run_file(
         def get_solver_solution() -> Solution:
             objective = m.getObjVal()
             vars_to_ignore = {"quadobjvar", "qmatrixvar", "quadobj", "qmatrix"}
+            # SCIP reformulates indicator, SOS and semi-continuous
+            # declarations into rows and variables of its own, none of which
+            # the model knows about.
+            aux_var_prefixes = ("indslack_",)
+            aux_con_prefixes = ("indlin_",)
+            aux_con_handlers = {
+                "indicator",
+                "SOS1",
+                "SOS2",
+                "bounddisjunction",
+            }
 
             s = m.getSols()[0]
-            kept_vars = [v for v in m.getVars() if v.name not in vars_to_ignore]
+            kept_vars = [
+                v
+                for v in m.getVars()
+                if v.name not in vars_to_ignore
+                and not v.name.startswith(aux_var_prefixes)
+            ]
             sol = _solution_from_names(
                 np.array([s[v] for v in kept_vars], dtype=float),
                 [v.name for v in kept_vars],
@@ -2594,7 +2613,13 @@ def get_solver_solution() -> Solution:
 
             cons = m.getConss(False)
             if len(cons) != 0:
-                kept_cons = [c for c in cons if c.name not in vars_to_ignore]
+                kept_cons = [
+                    c
+                    for c in cons
+                    if c.name not in vars_to_ignore
+                    and not c.name.startswith(aux_con_prefixes)
+                    and c.getConshdlrName() not in aux_con_handlers
+                ]
                 dual = _solution_from_names(
                     np.array([m.getDualSolVal(c) for c in kept_cons], dtype=float),
                     [c.name for c in kept_cons],

test/test_indicator_constraints.py

@@ -15,6 +15,10 @@
     "gurobi" not in available_solvers, reason="Gurobi not installed"
 )
 
+requires_scip = pytest.mark.skipif(
+    "scip" not in available_solvers, reason="SCIP not installed"
+)
+
 
 @pytest.fixture
 def mbx() -> tuple[Model, Variable, Variable]:
@@ -423,3 +427,43 @@ def test_mip_has_no_duals(self, mbx: tuple[Model, Variable, Variable]) -> None:
         m.solve(solver_name="gurobi")
         with pytest.raises(AttributeError, match="dual"):
             _ = m.matrices.dual
+
+
+class TestSCIP:
+    """SCIP reads indicator constraints from the problem file."""
+
+    @requires_scip
+    def test_indicator_is_enforced(self, mbx: tuple[Model, Variable, Variable]) -> None:
+        m, b, x = mbx
+        m.add_constraints(b >= 1, name="fix_b")
+        m.add_indicator_constraints(b, 1, x, "<=", 3, name="ic0")
+        m.add_objective(x, sense="max")
+        m.solve(solver_name="scip")
+        assert m.objective.value is not None
+        assert np.isclose(m.objective.value, 3, atol=1e-6)
+
+    @requires_scip
+    def test_indicator_is_slack_when_binary_is_off(
+        self, mbx: tuple[Model, Variable, Variable]
+    ) -> None:
+        """With the binary held at zero the bound does not apply."""
+        m, b, x = mbx
+        m.add_constraints(b <= 0, name="fix_b")
+        m.add_indicator_constraints(b, 1, x, "<=", 3, name="ic0")
+        m.add_objective(x, sense="max")
+        m.solve(solver_name="scip")
+        assert m.objective.value is not None
+        assert np.isclose(m.objective.value, 10, atol=1e-6)
+
+    @requires_scip
+    def test_solution_excludes_reformulation_variables(
+        self, mbx: tuple[Model, Variable, Variable]
+    ) -> None:
+        """SCIP adds a slack variable per indicator that is not part of the model."""
+        m, b, x = mbx
+        m.add_constraints(b >= 1, name="fix_b")
+        m.add_indicator_constraints(b, 1, x, "<=", 3, name="ic0")
+        m.add_objective(x, sense="max")
+        m.solve(solver_name="scip")
+        assert set(m.variables) == {"b", "x"}
+        assert np.isclose(float(m.variables["x"].solution), 3, atol=1e-6)

test/test_semi_continuous.py

@@ -203,3 +203,14 @@ def test_semi_continuous_solve_cuopt() -> None:
     m.solve(solver_name="cuopt", io_api="direct", log_to_console=False)
     assert m.objective.value is not None
     assert np.isclose(m.objective.value, 0, atol=1e-6)
+
+
+@pytest.mark.skipif("scip" not in available_solvers, reason="SCIP not installed")
+def test_semi_continuous_scip() -> None:
+    """A semi-continuous variable takes zero when its lower bound is out of reach."""
+    m = Model()
+    x = m.add_variables(lower=3, upper=10, name="x", semi_continuous=True)
+    m.add_constraints(x <= 2, name="cap")
+    m.add_objective(x, sense="max")
+    m.solve(solver_name="scip")
+    assert np.isclose(float(x.solution), 0, atol=1e-6)

test/test_sos_constraints.py

@@ -253,3 +253,30 @@ def test_to_highspy_raises_when_sos_present() -> None:
 
     with pytest.raises(ValueError, match="does not support SOS constraints"):
         m.to_highspy()
+
+
+@pytest.mark.skipif("scip" not in available_solvers, reason="SCIP not installed")
+def test_sos1_scip() -> None:
+    """SCIP reads SOS1 sets from the problem file and enforces them."""
+    m = Model()
+    locations = pd.Index([0, 1, 2], name="locations")
+    build = m.add_variables(coords=[locations], name="build", binary=True)
+    m.add_sos_constraints(build, sos_type=1, sos_dim="locations")
+    m.add_objective(build * np.array([1, 2, 3]), sense="max")
+    m.solve(solver_name="scip")
+    assert np.isclose(build.solution.values, [0, 0, 1]).all()
+    assert m.objective.value is not None
+    assert np.isclose(m.objective.value, 3)
+
+
+@pytest.mark.skipif("scip" not in available_solvers, reason="SCIP not installed")
+def test_sos2_scip() -> None:
+    """An SOS2 set allows two adjacent members to be nonzero."""
+    m = Model()
+    segments = pd.Index([0, 1, 2], name="seg")
+    var = m.add_variables(coords=[segments], lower=0, upper=1, name="lambda")
+    m.add_sos_constraints(var, sos_type=2, sos_dim="seg")
+    m.add_objective(var.sum(), sense="max")
+    m.solve(solver_name="scip")
+    assert m.objective.value is not None
+    assert np.isclose(m.objective.value, 2)