Merge pull request Pyomo#2994 from Robbybp/cyipopt-evalerror

Handle evaluation errors in CyIpopt solver

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())

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)

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,