diff --git a/src/sage/numerical/backends/generic_backend.pyx b/src/sage/numerical/backends/generic_backend.pyx
index c5c7dadc4c1..d6eafc65429 100644
--- a/src/sage/numerical/backends/generic_backend.pyx
+++ b/src/sage/numerical/backends/generic_backend.pyx
@@ -434,7 +434,9 @@ cdef class GenericBackend:
             sage: p.add_linear_constraint([(0, 3), (1, 2)], None, 6)
             sage: p.remove_constraints([0, 1])
         """
-        if isinstance(constraints, int): self.remove_constraint(constraints)
+        if isinstance(constraints, int):
+            self.remove_constraint(constraints)
+            return
 
         cdef int last = self.nrows() + 1
 
diff --git a/src/sage/numerical/backends/scip_backend.pyx b/src/sage/numerical/backends/scip_backend.pyx
index 606b16190a8..1ae54f483aa 100644
--- a/src/sage/numerical/backends/scip_backend.pyx
+++ b/src/sage/numerical/backends/scip_backend.pyx
@@ -61,6 +61,20 @@ cdef class SCIPBackend(GenericBackend):
     def get_constraints(self):
         """
         Get all constraints of the problem.
+
+        EXAMPLES::
+
+            sage: from sage.numerical.backends.generic_backend import get_solver
+            sage: lp = get_solver(solver="SCIP")
+            sage: lp.add_variables(3)
+            2
+            sage: lp.add_linear_constraint(zip([0, 1, 2], [8, 6, 1]), None, 48)
+            sage: lp.add_linear_constraint(zip([0, 1, 2], [2, 1.5, 0.5]), None, 8)
+
+            sage: lp.get_constraints()
+            [c1, c2]
+            sage: lp.row(1)                                                     # indirect doctest
+            ([0, 1, 2], [2.0, 1.5, 0.5])
         """
         if self.constraints is None:
             self.constraints = self.model.getConss()
@@ -71,10 +85,11 @@ cdef class SCIPBackend(GenericBackend):
         Get the model as a pyscipopt Model.
 
         EXAMPLES::
-        sage: from sage.numerical.backends.generic_backend import get_solver
-        sage: p = get_solver(solver = "SCIP")
-        sage: p._get_model()
-        <pyscipopt.scip.Model object at ...
+
+            sage: from sage.numerical.backends.generic_backend import get_solver
+            sage: p = get_solver(solver = "SCIP")
+            sage: p._get_model()
+            <pyscipopt.scip.Model object at ...
         """
         return self.model
 
@@ -369,6 +384,42 @@ cdef class SCIPBackend(GenericBackend):
         self.model.delCons(self.get_constraints()[i])
         self.constraints = None
 
+    cpdef remove_constraints(self, constraints) noexcept:
+        r"""
+        Remove several constraints.
+
+        INPUT:
+
+        - ``constraints`` -- an iterable containing the indices of the rows to remove
+
+        EXAMPLES::
+
+            sage: from sage.numerical.backends.generic_backend import get_solver
+            sage: p = get_solver(solver='SCIP')
+            sage: p.add_variables(2)
+            1
+            sage: p.add_linear_constraint([(0, 2), (1, 3)], None, 6)
+            sage: p.add_linear_constraint([(0, 3), (1, 2)], None, 6)
+            sage: p.row(0)
+            ([0, 1], [2.0, 3.0])
+            sage: p.remove_constraints([0, 1])
+            sage: p.nrows()
+            0
+        """
+        if isinstance(constraints, int):
+            self.remove_constraint(constraints)
+            return
+
+        if self.model.getStatus() != 'unknown':
+            self.model.freeTransform()
+            self.constraints = None
+
+        all_constraints = self.get_constraints()
+        to_remove = [all_constraints[i] for i in constraints]
+        for constraint in to_remove:
+            self.model.delCons(constraint)
+        self.constraints = None
+
     cpdef add_linear_constraint(self, coefficients, lower_bound, upper_bound, name=None) noexcept:
         """
         Add a linear constraint.
@@ -528,7 +579,7 @@ cdef class SCIPBackend(GenericBackend):
 
         INPUT:
 
-        - ``indices`` (list of integers) -- this list constains the
+        - ``indices`` (list of integers) -- this list contains the
           indices of the constraints in which the variable's
           coefficient is nonzero
 
@@ -692,6 +743,20 @@ cdef class SCIPBackend(GenericBackend):
 
         EXAMPLES::
 
+            sage: # needs sage.graphs
+            sage: g = graphs.CubeGraph(9)
+            sage: p = MixedIntegerLinearProgram(solver="SCIP")
+            sage: p.solver_parameter("limits/gap", 100)
+            sage: b = p.new_variable(binary=True)
+            sage: p.set_objective(p.sum(b[v] for v in g))
+            sage: for v in g:
+            ....:     p.add_constraint(b[v]+p.sum(b[u] for u in g.neighbors(v)) <= 1)
+            sage: p.add_constraint(b[v] == 1) # Force an easy non-0 solution
+            sage: p.solve() # rel tol 100
+            1.0
+            sage: backend = p.get_backend()
+            sage: backend.best_known_objective_bound() # random
+            31.0
         """
         return self.model.getPrimalbound()
 
@@ -713,12 +778,28 @@ cdef class SCIPBackend(GenericBackend):
 
         EXAMPLES::
 
+            sage: # needs sage.graphs
+            sage: g = graphs.CubeGraph(9)
+            sage: p = MixedIntegerLinearProgram(solver="SCIP")
+            sage: p.solver_parameter("limits/gap", 100)
+            sage: b = p.new_variable(binary=True)
+            sage: p.set_objective(p.sum(b[v] for v in g))
+            sage: for v in g:
+            ....:     p.add_constraint(b[v]+p.sum(b[u] for u in g.neighbors(v)) <= 1)
+            sage: p.add_constraint(b[v] == 1) # Force an easy non-0 solution
+            sage: p.solve() # rel tol 100
+            1.0
+            sage: backend = p.get_backend()
+            sage: backend.get_relative_objective_gap() # random
+            46.99999999999999
 
         TESTS:
 
         Just make sure that the variable *has* been defined, and is not just
         undefined::
 
+            sage: backend.get_relative_objective_gap() > 1                              # needs sage.graphs
+            True
         """
         return self.model.getGap()
 
@@ -811,7 +892,26 @@ cdef class SCIPBackend(GenericBackend):
             sage: p.add_linear_constraints(2, 2, None)
             sage: p.nrows()
             2
+
+        TESTS::
+
+        After calling :meth:`remove_constraints` we know that
+        `self.constraints is None`.  `SCIP` keeps track of the number
+        of constraints, so we can do the optimization::
+
+            sage: from sage.numerical.backends.generic_backend import get_solver
+            sage: p = get_solver(solver='SCIP')
+            sage: p.add_variables(2)
+            1
+            sage: p.add_linear_constraint([(0, 2), (1, 3)], None, 6)
+            sage: p.row(0)
+            ([0, 1], [2.0, 3.0])
+            sage: p.remove_constraints([0])
+            sage: p.nrows()
+            0
         """
+        if self.constraints is None:
+            return self.model.getNConss()
         return len(self.get_constraints())
 
     cpdef col_name(self, int index) noexcept:
@@ -870,7 +970,6 @@ cdef class SCIPBackend(GenericBackend):
             sage: p.set_variable_type(0,0)
             sage: p.is_variable_binary(0)
             True
-
         """
         return self.variables[index].vtype() == 'BINARY'
 
@@ -1155,27 +1254,10 @@ cdef class SCIPBackend(GenericBackend):
             6.0
         """
         cdef SCIPBackend cp = type(self)(maximization=self.is_maximization())
+        cp.model = Model(sourceModel=self.model, origcopy=True)
         cp.problem_name(self.problem_name())
-        for i, v in enumerate(self.variables):
-            vtype = v.vtype()
-            cp.add_variable(self.variable_lower_bound(i),
-                            self.variable_upper_bound(i),
-                            binary=vtype == 'BINARY',
-                            continuous=vtype == 'CONTINUOUS',
-                            integer=vtype == 'INTEGER',
-                            obj=self.objective_coefficient(i),
-                            name=self.col_name(i))
-        assert self.ncols() == cp.ncols()
-
-        for i in range(self.nrows()):
-            coefficients = zip(*self.row(i))
-            lower_bound, upper_bound = self.row_bounds(i)
-            name = self.row_name(i)
-            cp.add_linear_constraint(coefficients,
-                                     lower_bound,
-                                     upper_bound,
-                                     name=name)
-        assert self.nrows() == cp.nrows()
+        cp.obj_constant_term = self.obj_constant_term
+        cp.variables = cp.model.getVars()
         return cp
 
     cpdef solver_parameter(self, name, value=None) noexcept:
@@ -1189,6 +1271,13 @@ cdef class SCIPBackend(GenericBackend):
         - ``value`` -- the parameter's value if it is to be defined,
           or ``None`` (default) to obtain its current value.
 
+        EXAMPLES:
+
+            sage: from sage.numerical.backends.generic_backend import get_solver
+            sage: lp = get_solver(solver="SCIP")
+            sage: p.solver_parameter("limits/time", 1)
+            sage: p.solver_parameter("limits/time")
+            1.0
         """
         if value is not None:
             if name.lower() == 'timelimit':