diff --git a/src/ncopt/sqpgs/cvxpy_subproblem.py b/src/ncopt/sqpgs/cvxpy_subproblem.py
index 28ca74a..3d0f6a1 100644
--- a/src/ncopt/sqpgs/cvxpy_subproblem.py
+++ b/src/ncopt/sqpgs/cvxpy_subproblem.py
@@ -121,7 +121,7 @@ def solve(
 
         """
 
-        d = self.d
+        d = cp.Variable(self.dim)
         z = cp.Variable()
         if self.has_ineq_constraints:
             r_I = cp.Variable(gI_k.size, nonneg=True)
@@ -150,6 +150,9 @@ def solve(
         assert problem.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE}
         self._problem = problem
 
+        # Extract primal solution
+        self.d = d.value.copy()
+
         # Extract dual variables
         duals = problem.solution.dual_vars
         self.lambda_f = duals[obj_constraint.id]
diff --git a/src/ncopt/sqpgs/main.py b/src/ncopt/sqpgs/main.py
index 0130717..a5a3e37 100644
--- a/src/ncopt/sqpgs/main.py
+++ b/src/ncopt/sqpgs/main.py
@@ -21,6 +21,7 @@
 from ncopt.plot_utils import plot_metrics, plot_timings
 from ncopt.sqpgs.cvxpy_subproblem import CVXPYSubproblemSQPGS
 from ncopt.sqpgs.defaults import DEFAULT_ARG, DEFAULT_OPTION
+from ncopt.sqpgs.osqp_subproblem import OSQPSubproblemSQPGS
 from ncopt.utils import compute_batch_jacobian_vmap, get_logger
 
 
@@ -141,9 +142,12 @@ def solve(self):
         pE = np.repeat(pE_, self.dimE)
         ###############################################################
 
-        self.SP = CVXPYSubproblemSQPGS(
-            self.dim, p0, pI, pE, self.assert_tol, self.options["qp_solver"]
-        )
+        if self.options["qp_solver"] == "osqp":
+            self.SP = OSQPSubproblemSQPGS(self.dim, p0, pI, pE, self.assert_tol)
+        else:
+            self.SP = CVXPYSubproblemSQPGS(
+                self.dim, p0, pI, pE, self.assert_tol, self.options["qp_solver"]
+            )
 
         E_k = np.inf  # for stopping criterion
         x_kmin1 = None  # last iterate
@@ -235,10 +239,14 @@ def solve(self):
             # Subproblem solve
             ##############################################
             t1 = time.perf_counter()
-            self.SP.solve(np.linalg.cholesky(H), rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k)
+            if isinstance(self.SP, OSQPSubproblemSQPGS):
+                self.SP.solve(H, rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k)
+            else:
+                self.SP.solve(np.linalg.cholesky(H), rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k)
+
+            d_k = self.SP.d.copy()
             t2 = time.perf_counter()
 
-            d_k = self.SP.d.value.copy()
             # compute g_k from paper
             g_k = (
                 self.SP.lambda_f @ D_f
diff --git a/src/ncopt/sqpgs/osqp_subproblem.py b/src/ncopt/sqpgs/osqp_subproblem.py
index 3baa9f8..f614082 100644
--- a/src/ncopt/sqpgs/osqp_subproblem.py
+++ b/src/ncopt/sqpgs/osqp_subproblem.py
@@ -56,18 +56,6 @@ def problem(self) -> osqp.interface.OSQP:
     def status(self) -> str:
         return self.problem.status
 
-    @property
-    def objective_val(self) -> float:
-        return self.problem.value
-
-    @property
-    def setup_time(self) -> float:
-        return self.problem.solver_stats.setup_time
-
-    @property
-    def solve_time(self) -> float:
-        return self.problem.solver_stats.solve_time
-
     def _initialize(self):
         """
         The quadratic subrpoblem we solve in every iteration is of the form:
@@ -253,10 +241,22 @@ def solve(
 
         problem = osqp.OSQP()
 
-        problem.setup(P=sparse.csc_matrix(self.P), q=self.q, A=sparse.csc_matrix(A), l=l, u=u)
+        problem.setup(
+            P=sparse.csc_matrix(self.P),
+            q=self.q,
+            A=sparse.csc_matrix(A),
+            l=l,
+            u=u,
+            eps_abs=1e-05,
+            eps_rel=1e-05,
+            verbose=False,
+            polish=1,
+        )
 
         res = problem.solve()
-        # TODO: assert status, where is it stored in OSQP?
+
+        assert res.info.status not in ["unsolved"]
+        print(res.info.status)
         self._problem = problem
 
         primal_solution = res.x
@@ -267,12 +267,6 @@ def solve(
         self.rI = primal_solution[self.dim + 1 : self.dim + 1 + self.nI]
         self.rE = primal_solution[self.dim + 1 + self.nI :]
 
-        # TODO: pipe through assert_tol
-        assert len(self.rE) == self.nE
-        assert len(self.rI) == self.nI
-        # assert np.all(self.rI >= -1e-5) , f"{self.rI}"
-        # assert np.all(self.rE >= -1e-5), f"{self.rE}"
-
         # extract dual variables = KKT multipliers
         dual_solution = res.y
         lambda_f = dual_solution[: self.p0 + 1]
diff --git a/tests/test_subproblem.py b/tests/test_subproblem.py
index bf1e075..5799f7d 100644
--- a/tests/test_subproblem.py
+++ b/tests/test_subproblem.py
@@ -3,6 +3,7 @@
 import pytest
 
 from ncopt.sqpgs.cvxpy_subproblem import CVXPYSubproblemSQPGS
+from ncopt.sqpgs.osqp_subproblem import OSQPSubproblemSQPGS
 
 
 @pytest.fixture
@@ -61,3 +62,42 @@ def test_subproblem_eq(subproblem_eq: CVXPYSubproblemSQPGS):
     )
     assert subproblem_eq.status == cp.OPTIMAL
     assert np.isclose(subproblem_eq.objective_val, 1.0)
+
+
+def test_subproblem_consistency():
+    dim = 2
+    p0 = 2
+    pI = np.array([3])
+    pE = np.array([], dtype=int)
+    assert_tol = 1e-5
+    subproblem = CVXPYSubproblemSQPGS(dim, p0, pI, pE, assert_tol)
+    subproblem2 = OSQPSubproblemSQPGS(dim, p0, pI, pE, assert_tol)
+    D_f = np.array([[-2.0, 1.0], [-2.04236205, -1.0], [-1.92172864, -1.0]])
+    D_gI = [np.array([[0.0, 2.0], [0.0, 2.0], [1.41421356, 0.0], [1.41421356, 0.0]])]
+
+    subproblem.solve(
+        L=np.eye(2, dtype=float),
+        rho=0.1,
+        D_f=D_f,
+        D_gI=D_gI,
+        D_gE=[],
+        f_k=1.0,
+        gI_k=np.array([-1.0]),
+        gE_k=np.array([], dtype=float),
+    )
+
+    subproblem2.solve(
+        H=np.eye(2, dtype=float),
+        rho=0.1,
+        D_f=D_f,
+        D_gI=D_gI,
+        D_gE=[],
+        f_k=1.0,
+        gI_k=np.array([-1.0]),
+        gE_k=np.array([], dtype=float),
+    )
+
+    assert np.allclose(subproblem.d, subproblem2.d)
+    assert np.allclose(subproblem.lambda_f, subproblem2.lambda_f)
+    assert np.allclose(subproblem.lambda_gI, subproblem2.lambda_gI)
+    assert np.allclose(subproblem.lambda_gE, subproblem2.lambda_gE)