add tikohonov, and subproblem test

fabian-sp · Apr 22, 2024 · ccbd52b · ccbd52b
1 parent bc46c5d
commit ccbd52b
Show file tree

Hide file tree

Showing 4 changed files with 87 additions and 8 deletions.
diff --git a/src/ncopt/sqpgs/cvxpy_subproblem.py b/src/ncopt/sqpgs/cvxpy_subproblem.py
@@ -5,7 +5,7 @@
 
 from ncopt.sqpgs.defaults import DEFAULT_OPTION
 
-CVXPY_SOLVER_DICT = {"osqp": cp.OSQP, "clarabel": cp.CLARABEL}
+CVXPY_SOLVER_DICT = {"osqp-cvxpy": cp.OSQP, "clarabel": cp.CLARABEL}
 
 
 class CVXPYSubproblemSQPGS:

diff --git a/src/ncopt/sqpgs/main.py b/src/ncopt/sqpgs/main.py
@@ -238,11 +238,21 @@ def solve(self):
             ##############################################
             # Subproblem solve
             ##############################################
+            _reg_H = 1e-04
             t1 = time.perf_counter()
             if isinstance(self.SP, OSQPSubproblemSQPGS):
-                self.SP.solve(H, rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k)
+                self.SP.solve(H + _reg_H * np.eye(self.dim), 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)
+                self.SP.solve(
+                    np.linalg.cholesky(H + _reg_H * np.eye(self.dim)),
+                    rho,
+                    D_f,
+                    D_gI,
+                    D_gE,
+                    f_k,
+                    gI_k,
+                    gE_k,
+                )
 
             d_k = self.SP.d.copy()
             t2 = time.perf_counter()

diff --git a/src/ncopt/sqpgs/osqp_subproblem.py b/src/ncopt/sqpgs/osqp_subproblem.py
@@ -4,6 +4,14 @@
 import osqp
 from scipy import sparse
 
+# see: https://osqp.org/docs/interfaces/status_values.html
+OSQP_ALLOWED_STATUS = [
+    "solved",
+    "solved inaccurate",
+    "maximum iterations reached",
+    "run time limit reached",
+]
+
 
 class OSQPSubproblemSQPGS:
     def __init__(
@@ -255,8 +263,7 @@ def solve(
 
         res = problem.solve()
 
-        assert res.info.status not in ["unsolved"]
-        print(res.info.status)
+        assert res.info.status in OSQP_ALLOWED_STATUS, f"OSQP results in status {res.info.status}."
         self._problem = problem
 
         primal_solution = res.x

diff --git a/tests/test_subproblem.py b/tests/test_subproblem.py
@@ -64,7 +64,7 @@ def test_subproblem_eq(subproblem_eq: CVXPYSubproblemSQPGS):
     assert np.isclose(subproblem_eq.objective_val, 1.0)
 
 
-def test_subproblem_consistency():
+def test_subproblem_consistency_ineq():
     dim = 2
     p0 = 2
     pI = np.array([3])
@@ -76,7 +76,7 @@ def test_subproblem_consistency():
     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),
+        L=np.eye(dim, dtype=float),
         rho=0.1,
         D_f=D_f,
         D_gI=D_gI,
@@ -87,7 +87,7 @@ def test_subproblem_consistency():
     )
 
     subproblem2.solve(
-        H=np.eye(2, dtype=float),
+        H=np.eye(dim, dtype=float),
         rho=0.1,
         D_f=D_f,
         D_gI=D_gI,
@@ -101,3 +101,65 @@ def test_subproblem_consistency():
     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)
+
+
+def test_subproblem_consistency_ineq_eq():
+    dim = 4
+    p0 = 2
+    pI = np.array([1])
+    pE = np.array([1])
+    assert_tol = 1e-5
+    subproblem = CVXPYSubproblemSQPGS(dim, p0, pI, pE, assert_tol)
+    subproblem2 = OSQPSubproblemSQPGS(dim, p0, pI, pE, assert_tol)
+    D_f = np.array(
+        [
+            [1.83529234, 2.51893663, -0.87507966, 0.53305111],
+            [0.7042857, -0.19426588, 1.26820232, 0.59255224],
+            [-0.87356341, -0.24994689, -0.82073493, -1.18734854],
+        ]
+    )
+
+    D_gI = [
+        np.array(
+            [
+                [-1.13249659, 0.67854141, -0.0316317, -1.37963152],
+                [-1.64858759, 0.65296873, -0.72343526, 0.60976315],
+            ]
+        )
+    ]
+
+    D_gE = [
+        np.array(
+            [
+                [1.05532136, -0.12589961, 0.49469938, 0.22879848],
+                [2.10668334, -0.00816628, -0.43333072, 0.22656999],
+            ]
+        )
+    ]
+
+    subproblem.solve(
+        L=np.eye(dim, dtype=float),
+        rho=0.1,
+        D_f=D_f,
+        D_gI=D_gI,
+        D_gE=D_gE,
+        f_k=1.0,
+        gI_k=np.array([-1.0]),
+        gE_k=np.array([-1.0]),
+    )
+
+    subproblem2.solve(
+        H=np.eye(dim, dtype=float),
+        rho=0.1,
+        D_f=D_f,
+        D_gI=D_gI,
+        D_gE=D_gE,
+        f_k=1.0,
+        gI_k=np.array([-1.0]),
+        gE_k=np.array([-1.0]),
+    )
+
+    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)