diff --git a/ilpy/impl/solvers/ScipBackend.cpp b/ilpy/impl/solvers/ScipBackend.cpp
index 82fafe6..789c352 100644
--- a/ilpy/impl/solvers/ScipBackend.cpp
+++ b/ilpy/impl/solvers/ScipBackend.cpp
@@ -183,10 +183,12 @@ ScipBackend::addConstraint(const Constraint& constraint) {
 			lhs,
 			rhs));
 
-	_constraints.push_back(c);
 
 	SCIP_CALL_ABORT(SCIPaddCons(_scip, c));
-	SCIP_CALL_ABORT(SCIPreleaseCons(_scip, &c));
+	_constraints.push_back(c);
+	// we do not release the constraint here
+	// so that we can remove the constraints later in freeConstraints()
+	// SCIP_CALL_ABORT(SCIPreleaseCons(_scip, &c));
 }
 
 void
@@ -290,9 +292,21 @@ ScipBackend::freeVariables() {
 }
 
 void
-ScipBackend::freeConstraints() {
+ScipBackend::freeConstraints()
+{
+	// Iterate over all constraints and remove them from the SCIP model
+	for (SCIP_CONS *cons : _constraints)
+	{
+		if (cons != nullptr)
+		{
+			// Remove the constraint from the model
+			SCIP_CALL_ABORT(SCIPdelCons(_scip, cons));
+			// Release the constraint
+			SCIP_CALL_ABORT(SCIPreleaseCons(_scip, &cons));
+		}
+	}
 
-	// SCIPfree should free the constraints for us
+	// Clear the vector of constraint pointers
 	_constraints.clear();
 }
 
diff --git a/tests/test_solvers.py b/tests/test_solvers.py
index 579d7dd..f15aed4 100644
--- a/tests/test_solvers.py
+++ b/tests/test_solvers.py
@@ -199,3 +199,37 @@ def test_solution_indexing() -> None:
     _ = solution[x]  # smoke test
     solution[x] = 2  # can be used to set too
     assert solution[x] == 2
+
+
+@pytest.mark.parametrize("preference", PREFS)
+def test_solve_twice(preference: ilpy.Preference) -> None:
+    solver = ilpy.Solver(2, ilpy.VariableType.Integer, preference=preference)
+    x1 = ilpy.Variable("x1", index=0)
+    x2 = ilpy.Variable("x2", index=1)
+
+    solver.set_objective((x1 + x2).as_objective(ilpy.Maximize))
+
+    c0 = ilpy.Constraints()
+    c0.add(x1 + x2 <= 10)
+    c0.add(2 * x1 + 3 * x2 >= 12)
+    c0.add(x1 - x2 <= 5)
+
+    c1 = ilpy.Constraints()
+    c1.add(x1 + x2 <= 10)
+    c1.add(2 * x1 + 3 * x2 >= 12)
+    c1.add(x1 - x2 >= 5)
+
+    # initial solve
+    solver.set_constraints(c0)
+    solution = solver.solve()
+    assert list(solution) == [7, 3] and solution.get_value() == 10
+
+    # add a constraint and check that the solution has changed
+    solver.set_constraints(c1)
+    solution = solver.solve()
+    assert list(solution) == [18, -8] and solution.get_value() == 10
+
+    # return to initial constraints and check that the solution has changed back
+    solver.set_constraints(c0)
+    solution = solver.solve()
+    assert list(solution) == [7, 3] and solution.get_value() == 10