diff --git a/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py b/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py index 766ef96322a..cedbf430a12 100644 --- a/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py +++ b/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py @@ -289,6 +289,13 @@ class PyomoCyIpoptSolver(object): description="Set the function that will be called each iteration.", ), ) + CONFIG.declare( + "halt_on_evaluation_error", + ConfigValue( + default=None, + description="Whether to halt if a function or derivative evaluation fails", + ), + ) def __init__(self, **kwds): """Create an instance of the CyIpoptSolver. You must @@ -332,7 +339,9 @@ def solve(self, model, **kwds): nlp = pyomo_nlp.PyomoNLP(model) problem = cyipopt_interface.CyIpoptNLP( - nlp, intermediate_callback=config.intermediate_callback + nlp, + intermediate_callback=config.intermediate_callback, + halt_on_evaluation_error=config.halt_on_evaluation_error, ) ng = len(problem.g_lb()) nx = len(problem.x_lb()) diff --git a/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py b/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py index 2a7edb430d4..7ead30117cb 100644 --- a/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py +++ b/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py @@ -24,6 +24,7 @@ if not (numpy_available and scipy_available): raise unittest.SkipTest("Pynumero needs scipy and numpy to run NLP tests") +from pyomo.contrib.pynumero.exceptions import PyNumeroEvaluationError from pyomo.contrib.pynumero.asl import AmplInterface if not AmplInterface.available(): @@ -34,12 +35,18 @@ from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP from pyomo.contrib.pynumero.interfaces.cyipopt_interface import ( + cyipopt, cyipopt_available, CyIpoptNLP, ) from pyomo.contrib.pynumero.algorithms.solvers.cyipopt_solver import CyIpoptSolver +if cyipopt_available: + # We don't raise unittest.SkipTest if not cyipopt_available as there is a + # test below that tests an exception when cyipopt is unavailable. + cyipopt_ge_1_3 = hasattr(cyipopt, "CyIpoptEvaluationError") + def create_model1(): m = pyo.ConcreteModel() @@ -155,6 +162,29 @@ def f(model): return model +def make_hs071_model(): + # This is a model that is mathematically equivalent to the Hock-Schittkowski + # test problem 071, but that will trigger an evaluation error if x[0] goes + # above 1.1. + m = pyo.ConcreteModel() + m.x = pyo.Var([0, 1, 2, 3], bounds=(1.0, 5.0)) + m.x[0] = 1.0 + m.x[1] = 5.0 + m.x[2] = 5.0 + m.x[3] = 1.0 + m.obj = pyo.Objective(expr=m.x[0] * m.x[3] * (m.x[0] + m.x[1] + m.x[2]) + m.x[2]) + # This expression evaluates to zero, but is not well defined when x[0] > 1.1 + trivial_expr_with_eval_error = (pyo.sqrt(1.1 - m.x[0])) ** 2 + m.x[0] - 1.1 + m.ineq1 = pyo.Constraint(expr=m.x[0] * m.x[1] * m.x[2] * m.x[3] >= 25.0) + m.eq1 = pyo.Constraint( + expr=( + m.x[0] ** 2 + m.x[1] ** 2 + m.x[2] ** 2 + m.x[3] ** 2 + == 40.0 + trivial_expr_with_eval_error + ) + ) + return m + + @unittest.skipIf(cyipopt_available, "cyipopt is available") class TestCyIpoptNotAvailable(unittest.TestCase): def test_not_available_exception(self): @@ -257,3 +287,32 @@ def test_options(self): x, info = solver.solve(tee=False) nlp.set_primals(x) self.assertAlmostEqual(nlp.evaluate_objective(), -5.0879028e02, places=5) + + @unittest.skipUnless( + cyipopt_available and cyipopt_ge_1_3, "cyipopt version < 1.3.0" + ) + def test_hs071_evalerror(self): + m = make_hs071_model() + solver = pyo.SolverFactory("cyipopt") + res = solver.solve(m, tee=True) + + x = list(m.x[:].value) + expected_x = np.array([1.0, 4.74299964, 3.82114998, 1.37940829]) + np.testing.assert_allclose(x, expected_x) + + def test_hs071_evalerror_halt(self): + m = make_hs071_model() + solver = pyo.SolverFactory("cyipopt", halt_on_evaluation_error=True) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + res = solver.solve(m, tee=True) + + @unittest.skipIf( + not cyipopt_available or cyipopt_ge_1_3, "cyipopt version >= 1.3.0" + ) + def test_hs071_evalerror_old_cyipopt(self): + m = make_hs071_model() + solver = pyo.SolverFactory("cyipopt") + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + res = solver.solve(m, tee=True) diff --git a/pyomo/contrib/pynumero/interfaces/cyipopt_interface.py b/pyomo/contrib/pynumero/interfaces/cyipopt_interface.py index 19e74625d03..fc9c45c6d1a 100644 --- a/pyomo/contrib/pynumero/interfaces/cyipopt_interface.py +++ b/pyomo/contrib/pynumero/interfaces/cyipopt_interface.py @@ -23,6 +23,7 @@ import abc from pyomo.common.dependencies import attempt_import, numpy as np, numpy_available +from pyomo.contrib.pynumero.exceptions import PyNumeroEvaluationError def _cyipopt_importer(): @@ -252,7 +253,7 @@ def intermediate( class CyIpoptNLP(CyIpoptProblemInterface): - def __init__(self, nlp, intermediate_callback=None): + def __init__(self, nlp, intermediate_callback=None, halt_on_evaluation_error=None): """This class provides a CyIpoptProblemInterface for use with the CyIpoptSolver class that can take in an NLP as long as it provides vectors as numpy ndarrays and @@ -263,6 +264,23 @@ def __init__(self, nlp, intermediate_callback=None): self._nlp = nlp self._intermediate_callback = intermediate_callback + cyipopt_has_eval_error = cyipopt_available and hasattr( + cyipopt, "CyIpoptEvaluationError" + ) + if halt_on_evaluation_error is None: + # If using cyipopt >= 1.3, the default is to continue. + # Otherwise, the default is to halt (because we are forced to). + # + # If CyIpopt is not available, we "halt" (re-raise the original + # exception). + self._halt_on_evaluation_error = not cyipopt_has_eval_error + elif not halt_on_evaluation_error and not cyipopt_has_eval_error: + raise ValueError( + "halt_on_evaluation_error=False is only supported for cyipopt >= 1.3.0" + ) + else: + self._halt_on_evaluation_error = halt_on_evaluation_error + x = nlp.init_primals() y = nlp.init_duals() if np.any(np.isnan(y)): @@ -328,24 +346,54 @@ def scaling_factors(self): return obj_scaling, x_scaling, g_scaling def objective(self, x): - self._set_primals_if_necessary(x) - return self._nlp.evaluate_objective() + try: + self._set_primals_if_necessary(x) + return self._nlp.evaluate_objective() + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError( + "Error in objective function evaluation" + ) def gradient(self, x): - self._set_primals_if_necessary(x) - return self._nlp.evaluate_grad_objective() + try: + self._set_primals_if_necessary(x) + return self._nlp.evaluate_grad_objective() + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError( + "Error in objective gradient evaluation" + ) def constraints(self, x): - self._set_primals_if_necessary(x) - return self._nlp.evaluate_constraints() + try: + self._set_primals_if_necessary(x) + return self._nlp.evaluate_constraints() + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError("Error in constraint evaluation") def jacobianstructure(self): return self._jac_g.row, self._jac_g.col def jacobian(self, x): - self._set_primals_if_necessary(x) - self._nlp.evaluate_jacobian(out=self._jac_g) - return self._jac_g.data + try: + self._set_primals_if_necessary(x) + self._nlp.evaluate_jacobian(out=self._jac_g) + return self._jac_g.data + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError( + "Error in constraint Jacobian evaluation" + ) def hessianstructure(self): if not self._hessian_available: @@ -359,12 +407,20 @@ def hessian(self, x, y, obj_factor): if not self._hessian_available: raise ValueError("Hessian requested, but not supported by the NLP") - self._set_primals_if_necessary(x) - self._set_duals_if_necessary(y) - self._set_obj_factor_if_necessary(obj_factor) - self._nlp.evaluate_hessian_lag(out=self._hess_lag) - data = np.compress(self._hess_lower_mask, self._hess_lag.data) - return data + try: + self._set_primals_if_necessary(x) + self._set_duals_if_necessary(y) + self._set_obj_factor_if_necessary(obj_factor) + self._nlp.evaluate_hessian_lag(out=self._hess_lag) + data = np.compress(self._hess_lower_mask, self._hess_lag.data) + return data + except PyNumeroEvaluationError: + if self._halt_on_evaluation_error: + raise + else: + raise cyipopt.CyIpoptEvaluationError( + "Error in Lagrangian Hessian evaluation" + ) def intermediate( self, diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_cyipopt_interface.py b/pyomo/contrib/pynumero/interfaces/tests/test_cyipopt_interface.py index 2c5d8ff7e4e..f28b7b9b549 100644 --- a/pyomo/contrib/pynumero/interfaces/tests/test_cyipopt_interface.py +++ b/pyomo/contrib/pynumero/interfaces/tests/test_cyipopt_interface.py @@ -10,6 +10,7 @@ # ___________________________________________________________________________ import pyomo.common.unittest as unittest +import pyomo.environ as pyo from pyomo.contrib.pynumero.dependencies import ( numpy as np, @@ -20,20 +21,27 @@ if not (numpy_available and scipy_available): raise unittest.SkipTest("Pynumero needs scipy and numpy to run CyIpopt tests") +from pyomo.contrib.pynumero.exceptions import PyNumeroEvaluationError from pyomo.contrib.pynumero.asl import AmplInterface if not AmplInterface.available(): raise unittest.SkipTest("Pynumero needs the ASL extension to run CyIpopt tests") +from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP from pyomo.contrib.pynumero.interfaces.cyipopt_interface import ( + cyipopt, cyipopt_available, CyIpoptProblemInterface, + CyIpoptNLP, ) if not cyipopt_available: raise unittest.SkipTest("CyIpopt is not available") +cyipopt_ge_1_3 = hasattr(cyipopt, "CyIpoptEvaluationError") + + class TestSubclassCyIpoptInterface(unittest.TestCase): def test_subclass_no_init(self): class MyCyIpoptProblem(CyIpoptProblemInterface): @@ -88,5 +96,129 @@ def hessian(self, x, y, obj_factor): problem.solve(x0) +def _get_model_nlp_interface(halt_on_evaluation_error=None): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], initialize=1.0) + m.obj = pyo.Objective(expr=m.x[1] * pyo.sqrt(m.x[2]) + m.x[1] * m.x[3]) + m.eq1 = pyo.Constraint(expr=m.x[1] * pyo.sqrt(m.x[2]) == 1.0) + nlp = PyomoNLP(m) + interface = CyIpoptNLP(nlp, halt_on_evaluation_error=halt_on_evaluation_error) + bad_primals = np.array([1.0, -2.0, 3.0]) + indices = nlp.get_primal_indices([m.x[1], m.x[2], m.x[3]]) + bad_primals = bad_primals[indices] + return m, nlp, interface, bad_primals + + +class TestCyIpoptVersionDependentConfig(unittest.TestCase): + @unittest.skipIf(cyipopt_ge_1_3, "cyipopt version >= 1.3.0") + def test_config_error(self): + _, nlp, _, _ = _get_model_nlp_interface() + with self.assertRaisesRegex(ValueError, "halt_on_evaluation_error"): + interface = CyIpoptNLP(nlp, halt_on_evaluation_error=False) + + @unittest.skipIf(cyipopt_ge_1_3, "cyipopt version >= 1.3.0") + def test_default_config_with_old_cyipopt(self): + _, nlp, _, bad_x = _get_model_nlp_interface() + interface = CyIpoptNLP(nlp) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.objective(bad_x) + + @unittest.skipIf(not cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_default_config_with_new_cyipopt(self): + _, nlp, _, bad_x = _get_model_nlp_interface() + interface = CyIpoptNLP(nlp) + msg = "Error in objective function" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.objective(bad_x) + + +class TestCyIpoptEvaluationErrors(unittest.TestCase): + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_objective(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in objective function" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.objective(bad_x) + + def test_error_in_objective_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.objective(bad_x) + + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_gradient(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in objective gradient" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.gradient(bad_x) + + def test_error_in_gradient_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.gradient(bad_x) + + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_constraints(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in constraint evaluation" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.constraints(bad_x) + + def test_error_in_constraints_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.constraints(bad_x) + + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_jacobian(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in constraint Jacobian" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.jacobian(bad_x) + + def test_error_in_jacobian_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.jacobian(bad_x) + + @unittest.skipUnless(cyipopt_ge_1_3, "cyipopt version < 1.3.0") + def test_error_in_hessian(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=False + ) + msg = "Error in Lagrangian Hessian" + with self.assertRaisesRegex(cyipopt.CyIpoptEvaluationError, msg): + interface.hessian(bad_x, [1.0], 0.0) + + def test_error_in_hessian_halt(self): + m, nlp, interface, bad_x = _get_model_nlp_interface( + halt_on_evaluation_error=True + ) + msg = "Error in AMPL evaluation" + with self.assertRaisesRegex(PyNumeroEvaluationError, msg): + interface.hessian(bad_x, [1.0], 0.0) + + if __name__ == "__main__": unittest.main() diff --git a/pyomo/core/base/units_container.py b/pyomo/core/base/units_container.py index 1e9685188c7..dd6bb75aec9 100644 --- a/pyomo/core/base/units_container.py +++ b/pyomo/core/base/units_container.py @@ -1514,6 +1514,19 @@ def __getattribute__(self, attr): # present, at which point this instance __class__ will fall back # to PyomoUnitsContainer (where this method is not declared, OR # pint is not available and an ImportError will be raised. + # + # We need special case handling for __class__: gurobipy + # interrogates things by looking at their __class__ during + # python shutdown. Unfortunately, interrogating this + # singleton's __class__ evaluates `pint_available`, which - if + # DASK is installed - imports dask. Importing dask creates + # threading objects. Unfortunately, creating threading objects + # during interpreter shutdown generates a RuntimeError. So, our + # solution is to special-case the resolution of __class__ here + # to avoid accidentally triggering the imports. + if attr == "__class__": + return _DeferredUnitsSingleton + # if pint_available: # If the first thing that is being called is # "units.set_pint_registry(...)", then we will call __init__ diff --git a/pyomo/solvers/plugins/solvers/SCIPAMPL.py b/pyomo/solvers/plugins/solvers/SCIPAMPL.py index 69a24455706..9898b9cdd90 100644 --- a/pyomo/solvers/plugins/solvers/SCIPAMPL.py +++ b/pyomo/solvers/plugins/solvers/SCIPAMPL.py @@ -288,7 +288,7 @@ def _postsolve(self): # UNKNOWN # unknown='unknown' # An uninitialized value - if results.solver.message == "unknown": + if "unknown" in results.solver.message: results.solver.status = SolverStatus.unknown results.solver.termination_condition = TerminationCondition.unknown if len(results.solution) > 0: @@ -296,7 +296,7 @@ def _postsolve(self): # ABORTED # userInterrupt='userInterrupt' # Interrupt signal generated by user - elif results.solver.message == "user interrupt": + elif "user interrupt" in results.solver.message: results.solver.status = SolverStatus.aborted results.solver.termination_condition = TerminationCondition.userInterrupt if len(results.solution) > 0: @@ -304,7 +304,7 @@ def _postsolve(self): # OK # maxEvaluations='maxEvaluations' # Exceeded maximum number of problem evaluations - elif results.solver.message == "node limit reached": + elif "node limit reached" in results.solver.message: results.solver.status = SolverStatus.ok results.solver.termination_condition = TerminationCondition.maxEvaluations if len(results.solution) > 0: @@ -312,7 +312,7 @@ def _postsolve(self): # OK # maxEvaluations='maxEvaluations' # Exceeded maximum number of problem evaluations - elif results.solver.message == "total node limit reached": + elif "total node limit reached" in results.solver.message: results.solver.status = SolverStatus.ok results.solver.termination_condition = TerminationCondition.maxEvaluations if len(results.solution) > 0: @@ -320,7 +320,7 @@ def _postsolve(self): # OK # maxEvaluations='maxEvaluations' # Exceeded maximum number of problem evaluations - elif results.solver.message == "stall node limit reached": + elif "stall node limit reached" in results.solver.message: results.solver.status = SolverStatus.ok results.solver.termination_condition = TerminationCondition.maxEvaluations if len(results.solution) > 0: @@ -328,7 +328,7 @@ def _postsolve(self): # OK # maxTimeLimit='maxTimeLimit' # Exceeded maximum time limited allowed by user but having return a feasible solution - elif results.solver.message == "time limit reached": + elif "time limit reached" in results.solver.message: results.solver.status = SolverStatus.ok results.solver.termination_condition = TerminationCondition.maxTimeLimit if len(results.solution) > 0: @@ -336,7 +336,7 @@ def _postsolve(self): # OK # other='other' # Other, uncategorized normal termination - elif results.solver.message == "memory limit reached": + elif "memory limit reached" in results.solver.message: results.solver.status = SolverStatus.ok results.solver.termination_condition = TerminationCondition.other if len(results.solution) > 0: @@ -344,7 +344,7 @@ def _postsolve(self): # OK # other='other' # Other, uncategorized normal termination - elif results.solver.message == "gap limit reached": + elif "gap limit reached" in results.solver.message: results.solver.status = SolverStatus.ok results.solver.termination_condition = TerminationCondition.other if len(results.solution) > 0: @@ -352,7 +352,7 @@ def _postsolve(self): # OK # other='other' # Other, uncategorized normal termination - elif results.solver.message == "solution limit reached": + elif "solution limit reached" in results.solver.message: results.solver.status = SolverStatus.ok results.solver.termination_condition = TerminationCondition.other if len(results.solution) > 0: @@ -360,7 +360,7 @@ def _postsolve(self): # OK # other='other' # Other, uncategorized normal termination - elif results.solver.message == "solution improvement limit reached": + elif "solution improvement limit reached" in results.solver.message: results.solver.status = SolverStatus.ok results.solver.termination_condition = TerminationCondition.other if len(results.solution) > 0: @@ -368,19 +368,27 @@ def _postsolve(self): # OK # optimal='optimal' # Found an optimal solution - elif results.solver.message == "optimal solution found": + elif "optimal solution" in results.solver.message: results.solver.status = SolverStatus.ok results.solver.termination_condition = TerminationCondition.optimal if len(results.solution) > 0: results.solution(0).status = SolutionStatus.optimal - if results.problem.sense == ProblemSense.minimize: - results.problem.lower_bound = results.solver.primal_bound - else: - results.problem.upper_bound = results.solver.primal_bound + try: + if results.problem.sense == ProblemSense.minimize: + results.problem.lower_bound = results.solver.primal_bound + else: + results.problem.upper_bound = results.solver.primal_bound + except AttributeError: + """ + This may occur if SCIP solves the problem during presolve. In that case, + the log file may not get parsed correctly (self.read_scip_log), and + results.solver.primal_bound will not be populated. + """ + pass # WARNING # infeasible='infeasible' # Demonstrated that the problem is infeasible - elif results.solver.message == "infeasible": + elif "infeasible" in results.solver.message: results.solver.status = SolverStatus.warning results.solver.termination_condition = TerminationCondition.infeasible if len(results.solution) > 0: @@ -388,7 +396,7 @@ def _postsolve(self): # WARNING # unbounded='unbounded' # Demonstrated that problem is unbounded - elif results.solver.message == "unbounded": + elif "unbounded" in results.solver.message: results.solver.status = SolverStatus.warning results.solver.termination_condition = TerminationCondition.unbounded if len(results.solution) > 0: @@ -396,7 +404,7 @@ def _postsolve(self): # WARNING # infeasibleOrUnbounded='infeasibleOrUnbounded' # Problem is either infeasible or unbounded - elif results.solver.message == "infeasible or unbounded": + elif "infeasible or unbounded" in results.solver.message: results.solver.status = SolverStatus.warning results.solver.termination_condition = ( TerminationCondition.infeasibleOrUnbounded