diff --git a/pygradflow/iterate.py b/pygradflow/iterate.py index 1ac6b35..6fda583 100644 --- a/pygradflow/iterate.py +++ b/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) diff --git a/pygradflow/params.py b/pygradflow/params.py index 401cfe3..65eca4b 100644 --- a/pygradflow/params.py +++ b/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 diff --git a/pygradflow/solver.py b/pygradflow/solver.py index 64278a7..208875f 100644 --- a/pygradflow/solver.py +++ b/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) diff --git a/pyproject.toml b/pyproject.toml index 2ebb5bd..527be7f 100644 --- a/pyproject.toml +++ b/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 "] readme = "README.md" diff --git a/tests/pygradflow/test_solver.py b/tests/pygradflow/test_solver.py index 0aaf006..3e18e05 100644 --- a/tests/pygradflow/test_solver.py +++ b/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