Greybox Objectives

pyomo/contrib/pynumero/examples/external_grey_box/external_with_objective.py

@@ -0,0 +1,151 @@
+import math
+import numpy as np
+from scipy.sparse import coo_matrix
+import pyomo.environ as pe
+from pyomo.contrib.pynumero.interfaces.external_grey_box import (
+    ExternalGreyBoxModel,
+    ExternalGreyBoxBlock,
+)
+
+
+class Unconstrained(ExternalGreyBoxModel):
+    """
+    min (x+2)**2 + (y-2)**2
+    """
+
+    def input_names(self):
+        return ['x', 'y']
+
+    def set_input_values(self, input_values):
+        self._input_values = list(input_values)
+
+    def has_objective(self):
+        return True
+
+    def evaluate_objective(self):
+        x = self._input_values[0]
+        y = self._input_values[1]
+        return (x + 2) ** 2 + (y - 2) ** 2
+
+    def evaluate_grad_objective(self):
+        x = self._input_values[0]
+        y = self._input_values[1]
+        return np.asarray([2 * (x + 2), 2 * (y - 2)], dtype=float)
+
+
+class Constrained(ExternalGreyBoxModel):
+    """
+    min x**2 + y**2
+    s.t. 0 == y - exp(x)
+    """
+
+    def input_names(self):
+        return ['x', 'y']
+
+    def set_input_values(self, input_values):
+        self._input_values = list(input_values)
+
+    def has_objective(self):
+        return True
+
+    def evaluate_objective(self):
+        x = self._input_values[0]
+        y = self._input_values[1]
+        return x**2 + y**2
+
+    def evaluate_grad_objective(self):
+        x = self._input_values[0]
+        y = self._input_values[1]
+        return np.asarray([2 * x, 2 * y], dtype=float)
+
+    def equality_constraint_names(self):
+        return ['c1']
+
+    def evaluate_equality_constraints(self):
+        x = self._input_values[0]
+        y = self._input_values[1]
+        return np.asarray([y - math.exp(x)], dtype=float)
+
+    def evaluate_jacobian_equality_constraints(self):
+        x = self._input_values[0]
+        row = [0, 0]
+        col = [0, 1]
+        data = [-math.exp(x), 1]
+        jac = coo_matrix((data, (row, col)), shape=(1, 2))
+        return jac
+
+
+class ConstrainedWithHessian(Constrained):
+    def evaluate_hessian_objective(self):
+        row = [0, 1]
+        col = [0, 1]
+        data = [2, 2]
+        hess = coo_matrix((data, (row, col)), shape=(2, 2))
+        return hess
+
+    def set_equality_constraint_multipliers(self, eq_con_multiplier_values):
+        self._dual = eq_con_multiplier_values[0]
+
+    def evaluate_hessian_equality_constraints(self):
+        x = self._input_values[0]
+        row = [0]
+        col = [0]
+        data = [-math.exp(x) * self._dual]
+        hess = coo_matrix((data, (row, col)), shape=(2, 2))
+        return hess
+
+
+def solve_unconstrained():
+    m = pe.ConcreteModel()
+    m.z = pe.Var()
+    m.grey_box = ExternalGreyBoxBlock(external_model=Unconstrained())
+    m.c = pe.Constraint(expr=m.z == m.grey_box.inputs['x'] + 1)
+
+    opt = pe.SolverFactory('cyipopt')
+    opt.config.options['hessian_approximation'] = 'limited-memory'
+    res = opt.solve(m, tee=True)
+    pe.assert_optimal_termination(res)
+    x = m.grey_box.inputs['x'].value
+    y = m.grey_box.inputs['y'].value
+    assert math.isclose(x, -2)
+    assert math.isclose(y, 2)
+    return m
+
+
+def solve_constrained():
+    m = pe.ConcreteModel()
+    m.z = pe.Var()
+    m.grey_box = ExternalGreyBoxBlock(external_model=Constrained())
+    m.c2 = pe.Constraint(expr=m.z == m.grey_box.inputs['x'] + 1)
+
+    opt = pe.SolverFactory('cyipopt')
+    opt.config.options['hessian_approximation'] = 'limited-memory'
+    res = opt.solve(m, tee=True)
+    pe.assert_optimal_termination(res)
+    x = m.grey_box.inputs['x'].value
+    y = m.grey_box.inputs['y'].value
+    assert math.isclose(x, -0.4263027509962655)
+    assert math.isclose(y, 0.6529186403960969)
+    return m
+
+
+def solve_constrained_with_hessian():
+    m = pe.ConcreteModel()
+    m.z = pe.Var()
+    m.grey_box = ExternalGreyBoxBlock(external_model=ConstrainedWithHessian())
+    m.c2 = pe.Constraint(expr=m.z == m.grey_box.inputs['x'] + 1)
+
+    opt = pe.SolverFactory('cyipopt')
+    res = opt.solve(m, tee=True)
+    pe.assert_optimal_termination(res)
+    x = m.grey_box.inputs['x'].value
+    y = m.grey_box.inputs['y'].value
+    assert math.isclose(x, -0.4263027509962655)
+    assert math.isclose(y, 0.6529186403960969)
+    return m
+
+
+if __name__ == '__main__':
+    m = solve_constrained_with_hessian()
+    print(f"x: {m.grey_box.inputs['x'].value}")
+    print(f"y: {m.grey_box.inputs['y'].value}")

pyomo/contrib/pynumero/interfaces/external_grey_box.py

@@ -11,6 +11,7 @@
 
 import abc
 import logging
+import numpy as np
 from scipy.sparse import coo_matrix
 from pyomo.common.dependencies import numpy as np
 
@@ -45,6 +46,10 @@
       we have a specialized interface built on PyNumero that provides
       an interface to the CyIpopt solver.
 
+constraints: c(x) = 0
+outputs: y = c(x)
+
+
 To use this interface:
    * Create a class that is derived from ExternalGreyBoxModel and
      implement the necessary methods. This derived class must provide
@@ -126,6 +131,9 @@ def evaluate_hessian_outputs(self):
         the input variables. This method must return
         H_o^k = sum_i (y_o^k)_i * grad^2_{uu} w_o(u^k)
 
+    def evaluate_hessian_objective(self):
+        Compute the hessian of the objective
+
     Examples that show Hessian support are also found in:
     pyomo/contrib/pynumero/examples/external_grey_box/react-example/
 
@@ -315,6 +323,29 @@ def evaluate_jacobian_outputs(self):
     # def evaluate_hessian_outputs(self):
     #
 
+    # Support for objectives
+    def has_objective(self):
+        return False
+
+    def evaluate_objective(self) -> float:
+        """
+        Compute the objective from the  values set in
+        input_values
+        """
+        raise NotImplementedError(
+            'evaluate_objective called but not ' 'implemented in the derived class.'
+        )
+
+    def evaluate_grad_objective(self, out=None):
+        """
+        Compute the gradient of the objective from the
+        values set in input_values
+        """
+        raise NotImplementedError(
+            'evaluate_grad_objective called but not '
+            'implemented in the derived class.'
+        )
+
 
 class ExternalGreyBoxBlockData(BlockData):
     def set_external_model(self, external_grey_box_model, inputs=None, outputs=None):