Add unboundedness detection

pygradflow/iterate.py

@@ -152,6 +152,9 @@ def stat_res(self) -> float:
         r = self.obj_grad + self.cons_jac.T.dot(self.y) + self.bound_duals
         return np.linalg.norm(r, np.inf)
 
+    def is_feasible(self, tol):
+        return (self.cons_violation <= tol) and (self.bound_violation <= tol)
+
     @property
     def total_res(self) -> float:
         return max(self.cons_violation, self.bound_violation, self.stat_res)

pygradflow/params.py

@@ -90,6 +90,8 @@ class Params:
     time_limit: float = np.inf
     display_interval: float = 0.1
 
+    obj_lower_limit: float = -1e10
+
     @property
     def dtype(self):
         return np.float32 if self.precision == Precision.Single else np.float64

pygradflow/solver.py

@@ -23,6 +23,7 @@ class SolverStatus(Enum):
     Converged = (auto(), "Convergence achieved")
     IterationLimit = (auto(), "Reached iteration limit")
     TimeLimit = (auto(), "Reached time limit")
+    Unbounded = (auto(), "Unbounded")
     LocallyInfeasible = (auto(), "Local infeasibility detected")
 
     def __new__(cls, value, description):
@@ -37,13 +38,16 @@ def success(status):
 
 
 class SolverResult:
-    def __init__(self, x, y, d, success, status):
+    def __init__(self, x, y, d, status):
         self.x = x
         self.y = y
         self.d = d
-        self.success = success
         self.status = status
 
+    @property
+    def success(self):
+        return SolverStatus.success(self.status)
+
 
 header_interval = 25
 
@@ -153,8 +157,6 @@ def solve(self, x_0: np.ndarray, y_0: np.ndarray) -> SolverResult:
 
         lamb = params.lamb_init
 
-        success = True
-
         controller = step_controller(problem, params)
 
         self._deriv_check(x, y)
@@ -191,6 +193,12 @@ def solve(self, x_0: np.ndarray, y_0: np.ndarray) -> SolverResult:
                 status = SolverStatus.LocallyInfeasible
                 break
 
+            if (iterate.obj <= params.obj_lower_limit) and \
+               (iterate.is_feasible(params.opt_tol)):
+                logger.debug("Unboundedness detected")
+                status = SolverStatus.Unbounded
+                break
+
             step_result = self.compute_step(controller, iterate, 1.0 / lamb)
 
             x = iterate.x
@@ -252,7 +260,6 @@ def solve(self, x_0: np.ndarray, y_0: np.ndarray) -> SolverResult:
                     break
 
         else:
-            success = False
             status = SolverStatus.IterationLimit
             logger.debug("Iteration limit reached")
 
@@ -274,4 +281,4 @@ def solve(self, x_0: np.ndarray, y_0: np.ndarray) -> SolverResult:
         y = iterate.y
         d = iterate.bound_duals
 
-        return SolverResult(x, y, d, success, status)
+        return SolverResult(x, y, d, status)

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pygradflow"
-version = "0.2.4"
+version = "0.2.5"
 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_solver.py

@@ -1,10 +1,12 @@
 import numpy as np
+import scipy as sp
 
 import pytest
 
 from pygradflow.implicit_func import ImplicitFunc
 from pygradflow.iterate import Iterate
-from pygradflow.solver import Solver
+from pygradflow.problem import Problem
+from pygradflow.solver import Solver, SolverStatus
 from pygradflow.newton import newton_method
 
 from pygradflow.params import (
@@ -251,3 +253,39 @@ 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