Use Hessian approximation in CyIpopt when grey-box without Hessian is present

pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py

@@ -350,6 +350,7 @@ def solve(self, model, **kwds):
             intermediate_callback=config.intermediate_callback,
             halt_on_evaluation_error=config.halt_on_evaluation_error,
         )
+
         ng = len(problem.g_lb())
         nx = len(problem.x_lb())
         cyipopt_solver = problem
@@ -371,6 +372,20 @@ def solve(self, model, **kwds):
                 set_scaling = cyipopt_solver.setProblemScaling
             set_scaling(obj_scaling, x_scaling, g_scaling)
 
+        # has_hessian_support only exists for PyomoNLPWithGreyBoxBlocks, but if
+        # we have grey box blocks, we know our NLP is of this class.
+        if grey_box_blocks and not nlp.has_hessian_support():
+            # Note that `config` is a copy of the instance-level self.config so
+            # we don't have to reset the option after the solve.
+            if config.options.get("hessian_approximation", None) == "exact":
+                logger.warning(
+                    "'hessian_approximation' option is set to 'exact', but at"
+                    " least one grey box model does not support Hessians."
+                    " Overriding this option with"
+                    " hessian_approximation='limited-memory'."
+                )
+            config.options["hessian_approximation"] = "limited-memory"
+
         # add options
         try:
             add_option = cyipopt_solver.add_option

pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py

@@ -39,6 +39,11 @@
 )
 
 from pyomo.contrib.pynumero.algorithms.solvers.cyipopt_solver import CyIpoptSolver
+from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock
+import pyomo.contrib.pynumero.interfaces.tests.external_grey_box_models as ex_models
+from pyomo.contrib.pynumero.examples.external_grey_box.react_example.reactor_model_outputs import (
+    ReactorConcentrationsOutputModel,
+)
 
 if cyipopt_available:
     # We don't raise unittest.SkipTest if not cyipopt_available as there is a
@@ -207,7 +212,7 @@ def test_model1(self):
         nlp.set_primals(x)
         nlp.set_duals(y_sol)
         self.assertAlmostEqual(nlp.evaluate_objective(), -428.6362455416348, places=5)
-        self.assertTrue(np.allclose(info['mult_g'], y_sol, rtol=1e-4))
+        self.assertTrue(np.allclose(info["mult_g"], y_sol, rtol=1e-4))
 
     def test_model1_with_scaling(self):
         m = create_model1()
@@ -219,35 +224,35 @@ def test_model1_with_scaling(self):
 
         with TempfileManager.new_context() as temp:
             cynlp = CyIpoptNLP(PyomoNLP(m))
-            logfile = temp.create_tempfile('_cyipopt-scaling.log')
+            logfile = temp.create_tempfile("_cyipopt-scaling.log")
             options = {
-                'nlp_scaling_method': 'user-scaling',
-                'output_file': logfile,
-                'file_print_level': 10,
-                'max_iter': 0,
+                "nlp_scaling_method": "user-scaling",
+                "output_file": logfile,
+                "file_print_level": 10,
+                "max_iter": 0,
             }
             solver = CyIpoptSolver(cynlp, options=options)
             x, info = solver.solve()
             cynlp.close()
 
-            with open(logfile, 'r') as fd:
+            with open(logfile, "r") as fd:
                 solver_trace = fd.read()
 
         # check for the following strings in the log
-        self.assertIn('nlp_scaling_method = user-scaling', solver_trace)
+        self.assertIn("nlp_scaling_method = user-scaling", solver_trace)
         self.assertIn(f"output_file = {logfile}", solver_trace)
-        self.assertIn('objective scaling factor = 1e-06', solver_trace)
-        self.assertIn('x scaling provided', solver_trace)
-        self.assertIn('c scaling provided', solver_trace)
-        self.assertIn('d scaling provided', solver_trace)
+        self.assertIn("objective scaling factor = 1e-06", solver_trace)
+        self.assertIn("x scaling provided", solver_trace)
+        self.assertIn("c scaling provided", solver_trace)
+        self.assertIn("d scaling provided", solver_trace)
         self.assertIn('DenseVector "x scaling vector" with 3 elements:', solver_trace)
-        self.assertIn('x scaling vector[    1]= 1.0000000000000000e+00', solver_trace)
-        self.assertIn('x scaling vector[    2]= 1.0000000000000000e+00', solver_trace)
-        self.assertIn('x scaling vector[    3]= 4.0000000000000000e+00', solver_trace)
+        self.assertIn("x scaling vector[    1]= 1.0000000000000000e+00", solver_trace)
+        self.assertIn("x scaling vector[    2]= 1.0000000000000000e+00", solver_trace)
+        self.assertIn("x scaling vector[    3]= 4.0000000000000000e+00", solver_trace)
         self.assertIn('DenseVector "c scaling vector" with 1 elements:', solver_trace)
-        self.assertIn('c scaling vector[    1]= 2.0000000000000000e+00', solver_trace)
+        self.assertIn("c scaling vector[    1]= 2.0000000000000000e+00", solver_trace)
         self.assertIn('DenseVector "d scaling vector" with 1 elements:', solver_trace)
-        self.assertIn('d scaling vector[    1]= 3.0000000000000000e+00', solver_trace)
+        self.assertIn("d scaling vector[    1]= 3.0000000000000000e+00", solver_trace)
 
     def test_model2(self):
         model = create_model2()
@@ -260,7 +265,7 @@ def test_model2(self):
         nlp.set_primals(x)
         nlp.set_duals(y_sol)
         self.assertAlmostEqual(nlp.evaluate_objective(), -31.000000057167462, places=5)
-        self.assertTrue(np.allclose(info['mult_g'], y_sol, rtol=1e-4))
+        self.assertTrue(np.allclose(info["mult_g"], y_sol, rtol=1e-4))
 
     def test_model3(self):
         G = np.array([[6, 2, 1], [2, 5, 2], [1, 2, 4]])
@@ -278,12 +283,12 @@ def test_model3(self):
         nlp.set_primals(x)
         nlp.set_duals(y_sol)
         self.assertAlmostEqual(nlp.evaluate_objective(), -3.5, places=5)
-        self.assertTrue(np.allclose(info['mult_g'], y_sol, rtol=1e-4))
+        self.assertTrue(np.allclose(info["mult_g"], y_sol, rtol=1e-4))
 
     def test_options(self):
         model = create_model1()
         nlp = PyomoNLP(model)
-        solver = CyIpoptSolver(CyIpoptNLP(nlp), options={'max_iter': 1})
+        solver = CyIpoptSolver(CyIpoptNLP(nlp), options={"max_iter": 1})
         x, info = solver.solve(tee=False)
         nlp.set_primals(x)
         self.assertAlmostEqual(nlp.evaluate_objective(), -5.0879028e02, places=5)
@@ -414,3 +419,67 @@ def intermediate(
         x, y = iterate_data[-1]
         self.assertTrue(np.allclose(x_sol, x))
         self.assertTrue(np.allclose(y_sol, y))
+
+
+@unittest.skipUnless(cyipopt_available, "cyipopt is not available")
+class TestCyIpoptGreyBox(unittest.TestCase):
+    """Most of the grey-box functionality is tested elsewhere. Here we just
+    need to test that we correctly override the hessian_approximation option
+    when Hessians aren't supported.
+    """
+
+    def test_solve_hessian_not_supported(self):
+        m = pyo.ConcreteModel()
+        m.reactor = ExternalGreyBoxBlock(
+            external_model=ReactorConcentrationsOutputModel()
+        )
+        m.k1con = pyo.Constraint(expr=m.reactor.inputs["k1"] == 5 / 6)
+        m.k2con = pyo.Constraint(expr=m.reactor.inputs["k2"] == 5 / 3)
+        m.k3con = pyo.Constraint(expr=m.reactor.inputs["k3"] == 1 / 6000)
+        m.cafcon = pyo.Constraint(expr=m.reactor.inputs["caf"] == 10000)
+        m.obj = pyo.Objective(expr=m.reactor.outputs["cb"], sense=pyo.maximize)
+        solver = pyo.SolverFactory("cyipopt")
+        results = solver.solve(m, tee=True)
+        pyo.assert_optimal_termination(results)
+        self.assertNotIn("hessian_approximation", solver.config.options)
+        self.assertAlmostEqual(pyo.value(m.reactor.inputs["sv"]), 1.34381, places=3)
+        self.assertAlmostEqual(pyo.value(m.reactor.outputs["cb"]), 1072.4372, places=2)
+
+    def test_solve_hessian_not_supported_override_user_option(self):
+        """Make sure we do this even when the user says to use an exact Hessian"""
+        m = pyo.ConcreteModel()
+        m.reactor = ExternalGreyBoxBlock(
+            external_model=ReactorConcentrationsOutputModel()
+        )
+        m.k1con = pyo.Constraint(expr=m.reactor.inputs["k1"] == 5 / 6)
+        m.k2con = pyo.Constraint(expr=m.reactor.inputs["k2"] == 5 / 3)
+        m.k3con = pyo.Constraint(expr=m.reactor.inputs["k3"] == 1 / 6000)
+        m.cafcon = pyo.Constraint(expr=m.reactor.inputs["caf"] == 10000)
+        m.obj = pyo.Objective(expr=m.reactor.outputs["cb"], sense=pyo.maximize)
+        solver = pyo.SolverFactory("cyipopt")
+        solver.config.options["hessian_approximation"] = "exact"
+        results = solver.solve(m, tee=True)
+        pyo.assert_optimal_termination(results)
+        self.assertEqual(solver.config.options["hessian_approximation"], "exact")
+        self.assertAlmostEqual(pyo.value(m.reactor.inputs["sv"]), 1.34381, places=3)
+        self.assertAlmostEqual(pyo.value(m.reactor.outputs["cb"]), 1072.4372, places=2)
+
+    def test_hessian_supported_no_override(self):
+        """Make sure we didn't accidentally override the hessian_approximation
+        option when Hessian *is* supported.
+        """
+        m = pyo.ConcreteModel()
+        external_model = ex_models.PressureDropSingleOutputWithHessian()
+        m.x = pyo.Var(range(external_model.n_inputs()))
+        m.eq = pyo.Constraint(expr=sum(m.x.values()) == 1)
+        solver = pyo.SolverFactory("cyipopt")
+        with TempfileManager.new_context() as temp:
+            logfile = temp.create_tempfile("hessian-no-override.log")
+            options = dict(output_file=logfile, max_iter=0, print_user_options="yes")
+            solver.config.options.update(options)
+            _, cynlp = solver.solve(m, tee=True, return_nlp=True)
+            # IIRC this may be necessary to avoid file IO race condition
+            # cynlp.close()
+            with open(logfile, "r") as fd:
+                solver_trace = fd.read()
+        self.assertNotIn("hessian_approximation", solver_trace)

pyomo/contrib/pynumero/interfaces/pyomo_grey_box_nlp.py

@@ -503,6 +503,9 @@ def load_state_into_pyomo(self, bound_multipliers=None):
                     zip(self._pyomo_model_var_datas, -obj_sign * bound_multipliers[1])
                 )
 
+    def has_hessian_support(self):
+        return self._has_hessian_support
+
 
 def _default_if_none(value, default):
     if value is None:

pyomo/contrib/pynumero/interfaces/tests/test_pyomo_grey_box_nlp.py

@@ -2586,6 +2586,28 @@ def test_duals_after_solve(self):
         )
         self.assertAlmostEqual(m.dual[m.egb]['egb.u2_con'], 62.5, places=3)
 
+    def test_has_hessian_support_false(self):
+        external_model = ex_models.PressureDropSingleOutput()
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(range(external_model.n_inputs()))
+        m.gb = ExternalGreyBoxBlock()
+        m.gb.set_external_model(external_model, inputs=list(m.x.values()))
+        # Some random constraint to let us construct the NLP...
+        m.eq = pyo.Constraint(expr=sum(m.x.values()) == 1)
+        nlp = PyomoNLPWithGreyBoxBlocks(m)
+        self.assertFalse(nlp.has_hessian_support())
+
+    def test_has_hessian_support_true(self):
+        external_model = ex_models.PressureDropSingleOutputWithHessian()
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(range(external_model.n_inputs()))
+        m.gb = ExternalGreyBoxBlock()
+        m.gb.set_external_model(external_model, inputs=list(m.x.values()))
+        # Some random constraint to let us construct the NLP...
+        m.eq = pyo.Constraint(expr=sum(m.x.values()) == 1)
+        nlp = PyomoNLPWithGreyBoxBlocks(m)
+        self.assertTrue(nlp.has_hessian_support())
+
 
 class TestGreyBoxObjectives(unittest.TestCase):
     @unittest.skipIf(not cyipopt_available, "CyIpopt needed to run tests with solve")