Merge pull request #19 from chrhansk/feature-local-infeasibility-tests

Add test for (improved) local infeasibility detection

pygradflow/iterate.py

@@ -97,20 +97,22 @@ def aug_lag_deriv_xx(self, rho: float) -> sp.sparse.spmatrix:
     def dist(self, other: "Iterate") -> float:
         return norm_mult(self.x - other.x, self.y - other.y)
 
-    def locally_infeasible(self, tol: float) -> bool:
+    def locally_infeasible(self,
+                           feas_tol: float,
+                           local_infeas_tol: float) -> bool:
         """
         Check if the iterate is locally infeasible. It is
         judged to be locally infeasible if the constraint
         violation is greater than the tolerance and
         optimality conditions for the minimization
         of the constraint violation are (approximately) satisfied.
         """
-        if self.cons_violation <= tol:
+        if self.cons_violation <= feas_tol:
             return False
 
         infeas_opt_res = self.cons_jac.T.dot(self.cons)
 
-        return np.linalg.norm(infeas_opt_res) <= tol
+        return np.linalg.norm(infeas_opt_res, ord=np.inf) <= local_infeas_tol
 
     @functools.cached_property
     def active_set(self) -> ActiveSet:

pygradflow/params.py

@@ -70,6 +70,8 @@ class Params:
     lamb_term: float = 1e-8
     active_tol: float = 1e-8
 
+    local_infeas_tol: float = 1e-8
+
     newton_type: NewtonType = NewtonType.Simplified
     newton_tol: float = 1e-8
 

pygradflow/solver.py

@@ -191,7 +191,8 @@ def solve(self, x_0: np.ndarray, y_0: np.ndarray) -> SolverResult:
                 status = SolverStatus.Converged
                 break
 
-            if iterate.locally_infeasible(params.opt_tol):
+            if iterate.locally_infeasible(params.opt_tol,
+                                          params.local_infeas_tol):
                 logger.debug("Local infeasibility detected")
                 status = SolverStatus.LocallyInfeasible
                 break

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pygradflow"
-version = "0.2.7"
+version = "0.2.8"
 description = "PyGradFlow is a simple implementation of the sequential homotopy method to be used to solve general nonlinear programs."
 authors = ["Christoph Hansknecht <[email protected]>"]
 readme = "README.md"

tests/pygradflow/test_conds.py

@@ -0,0 +1,81 @@
+import numpy as np
+import scipy as sp
+
+from pygradflow.problem import Problem
+from pygradflow.solver import Solver, SolverStatus
+
+
+def test_detect_unbounded():
+    num_vars = 1
+
+    class UnboundedProblem(Problem):
+        def __init__(self):
+            var_lb = np.full(shape=(num_vars,), fill_value=-np.inf)
+            var_ub = np.full(shape=(num_vars,), fill_value=np.inf)
+            super().__init__(var_lb, var_ub)
+
+        def obj(self, x):
+            return x[0]
+
+        def obj_grad(self, x):
+            return np.array([1.] + [0.] * (num_vars - 1))
+
+        def cons(self, x):
+            return np.array([])
+
+        def cons_jac(self, x):
+            return sp.sparse.coo_matrix((0, num_vars))
+
+        def lag_hess(self, x, lag):
+            return sp.sparse.diags([0.]*num_vars)
+
+    problem = UnboundedProblem()
+
+    solver = Solver(problem)
+
+    x0 = np.array([0.0]*num_vars)
+    y0 = np.array([])
+
+    result = solver.solve(x0, y0)
+
+    assert result.status == SolverStatus.Unbounded
+
+
+def test_detect_infeasible():
+    num_vars = 1
+    num_cons = 1
+
+    class InfeasibleProblem(Problem):
+        def __init__(self):
+            var_lb = np.full(shape=(num_vars,), fill_value=-np.inf)
+            var_ub = np.full(shape=(num_vars,), fill_value=np.inf)
+            super().__init__(var_lb, var_ub, num_cons=num_cons)
+
+        def obj(self, x):
+            return x[0]
+
+        def obj_grad(self, x):
+            return np.array([1.] + [0.] * (num_vars - 1))
+
+        def cons(self, x):
+            x = x.item()
+            return np.array([x*x + 1])
+
+        def cons_jac(self, x):
+            x = x.item()
+            jac = np.array([[2*x]])
+            return sp.sparse.coo_matrix(jac)
+
+        def lag_hess(self, x, lag):
+            return sp.sparse.diags([2.]*num_vars)
+
+    problem = InfeasibleProblem()
+
+    solver = Solver(problem)
+
+    x0 = np.array([0.0]*num_vars)
+    y0 = np.array([0.0])
+
+    result = solver.solve(x0, y0)
+
+    assert result.status == SolverStatus.LocallyInfeasible

tests/pygradflow/test_solver.py

@@ -254,39 +254,3 @@ def cons_jac(x):
     invalid_col = jac.col[invalid_index]
 
     assert (e.invalid_indices == [[invalid_row, invalid_col]]).all()
-
-
-def test_detect_unbounded():
-    num_vars = 1
-
-    class UnboundedProblem(Problem):
-        def __init__(self):
-            var_lb = np.full(shape=(num_vars,), fill_value=-np.inf)
-            var_ub = np.full(shape=(num_vars,), fill_value=np.inf)
-            super().__init__(var_lb, var_ub)
-
-        def obj(self, x):
-            return x[0]
-
-        def obj_grad(self, x):
-            return np.array([1.] + [0.] * (num_vars - 1))
-
-        def cons(self, x):
-            return np.array([])
-
-        def cons_jac(self, x):
-            return sp.sparse.coo_matrix((0, num_vars))
-
-        def lag_hess(self, x, lag):
-            return sp.sparse.diags([0.]*num_vars)
-
-    problem = UnboundedProblem()
-
-    solver = Solver(problem)
-
-    x0 = np.array([0.0]*num_vars)
-    y0 = np.array([])
-
-    result = solver.solve(x0, y0)
-
-    assert result.status == SolverStatus.Unbounded