Implicit constraint

pyomo/contrib/incidence_analysis/incidence.py

@@ -6,6 +6,9 @@
 # Solutions of Sandia, LLC, the U.S. Government retains certain rights in this
 # software.  This software is distributed under the 3-clause BSD License.
 # ____________________________________________________________________________________
+#
+#  Additional contributions Copyright (c) 2026 OLI Systems, Inc.
+#  ___________________________________________________________________________________
 """Functionality for identifying variables that participate in expressions"""
 
 from contextlib import nullcontext
@@ -18,6 +21,9 @@
     IncidenceMethod,
     get_config_from_kwds,
 )
+from pyomo.contrib.pynumero.interfaces.external_grey_box_constraint import (
+    EGBConstraintBody,
+)
 
 
 #
@@ -169,7 +175,48 @@ def get_incident_variables(expr, **kwds):
         # Developer error, this should never happen!
         raise RuntimeError("_ampl_repn_visitor must be provided when using ampl_repn")
 
+    return get_variables_incident_to_constraint(
+        expr, method, include_fixed, linear_only, amplrepnvisitor
+    )
+
+
+def get_variables_incident_to_constraint(
+    expr,
+    method: IncidenceMethod,
+    include_fixed: bool = False,
+    linear_only: bool = False,
+    amplrepnvisitor=None,
+):
+    """
+    Helper function to identify variables that are incident on an expression, based on the type of the expression and the method specified.
+
+    Parameters
+    ----------
+    expr: NumericExpression
+        The expression to analyze for incident variables.
+    method: IncidenceMethod
+        The method to use for identifying incident variables (not used with EGBConstraintBody).
+    include_fixed: bool, optional
+        Whether to include fixed variables in the result.
+    linear_only: bool, optional
+        Whether to include only linear variables in the result.
+    amplrepnvisitor: optional
+        The AMPL representation visitor, required if using the AMPL representation method.
+
+    Returns
+    -------
+    list of VarData
+        List containing the variables that are incident on the expression.
+
+    Raises
+    ------
+    ValueError
+        If an unrecognized method is specified for identifying incident variables.
+    """
     # Dispatch to correct method
+    if isinstance(expr, EGBConstraintBody):
+        # If the expression is the body of an implicit constraint, we need to use the get_incident_variables method defined on EGBConstraintBody
+        return expr.get_incident_variables()
     if method is IncidenceMethod.identify_variables:
         return _get_incident_via_identify_variables(expr, include_fixed)
     elif method is IncidenceMethod.standard_repn:

pyomo/contrib/incidence_analysis/interface.py

@@ -6,15 +6,16 @@
 # Solutions of Sandia, LLC, the U.S. Government retains certain rights in this
 # software.  This software is distributed under the 3-clause BSD License.
 # ____________________________________________________________________________________
+#
+#  Additional contributions Copyright (c) 2026 OLI Systems, Inc.
+#  ___________________________________________________________________________________
 """Utility functions and a utility class for interfacing Pyomo components with
 useful graph algorithms.
 
 """
 
-import enum
 import textwrap
 from pyomo.core.base.block import BlockData
-from pyomo.core.base.var import Var
 from pyomo.core.base.constraint import Constraint
 from pyomo.core.base.objective import Objective
 from pyomo.core.expr import EqualityExpression
@@ -29,24 +30,26 @@
 from pyomo.common.deprecation import deprecated, deprecation_warning
 from pyomo.contrib.incidence_analysis.config import get_config_from_kwds
 from pyomo.contrib.incidence_analysis.matching import maximum_matching
-from pyomo.contrib.incidence_analysis.connected import get_independent_submatrices
-from pyomo.contrib.incidence_analysis.triangularize import (
-    get_scc_of_projection,
-    block_triangularize,
-    get_diagonal_blocks,
-    get_blocks_from_maps,
-)
+from pyomo.contrib.incidence_analysis.triangularize import get_scc_of_projection
 from pyomo.contrib.incidence_analysis.dulmage_mendelsohn import (
     dulmage_mendelsohn,
     RowPartition,
     ColPartition,
 )
 from pyomo.contrib.incidence_analysis.incidence import get_incident_variables
 from pyomo.contrib.pynumero.asl import AmplInterface
+from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock
+from pyomo.contrib.pynumero.interfaces.external_grey_box_constraint import (
+    ExternalGreyBoxConstraint,
+)
 
 pyomo_nlp, pyomo_nlp_available = attempt_import(
     "pyomo.contrib.pynumero.interfaces.pyomo_nlp"
 )
+if pyomo_nlp_available:
+    from pyomo.contrib.pynumero.interfaces.pyomo_grey_box_nlp import (
+        PyomoNLPWithGreyBoxBlocks,
+    )
 asl_available = pyomo_nlp_available & AmplInterface.available()
 
 
@@ -283,6 +286,15 @@ def __init__(self, model=None, active=True, include_inequality=True, **kwds):
                 for con in model.component_data_objects(Constraint, active=active)
                 if include_inequality or isinstance(con.expr, EqualityExpression)
             ]
+
+            for egb in model.component_data_objects(
+                ExternalGreyBoxBlock, active=active
+            ):
+                for ic in egb.component_data_objects(
+                    ExternalGreyBoxConstraint, active=active
+                ):
+                    self._constraints.append(ic)
+
             self._variables = list(
                 _generate_variables_in_constraints(self._constraints, **self._config)
             )
@@ -295,7 +307,9 @@ def __init__(self, model=None, active=True, include_inequality=True, **kwds):
             self._incidence_graph = get_bipartite_incidence_graph(
                 self._variables, self._constraints, **self._config
             )
-        elif pyomo_nlp_available and isinstance(model, pyomo_nlp.PyomoNLP):
+        elif pyomo_nlp_available and isinstance(
+            model, (pyomo_nlp.PyomoNLP, PyomoNLPWithGreyBoxBlocks)
+        ):
             if not active:
                 raise ValueError(
                     "Cannot get the Jacobian of inactive constraints from the "

pyomo/contrib/incidence_analysis/tests/test_external_grey_box_integration.py

@@ -0,0 +1,185 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2025
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+#
+#  Additional contributions Copyright (c) 2026 OLI Systems, Inc.
+#  ___________________________________________________________________________
+
+import pyomo.common.unittest as unittest
+import pyomo.environ as pyo
+from pyomo.common.collections import ComponentSet
+
+from pyomo.contrib.incidence_analysis import IncidenceGraphInterface
+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.interfaces.external_grey_box_constraint import (
+    ExternalGreyBoxConstraint,
+)
+
+
+class TestExternalGreyBoxIncidence(unittest.TestCase):
+    def test_pressure_drop_single_output(self):
+        m = pyo.ConcreteModel()
+        m.egb = ExternalGreyBoxBlock()
+        m.egb.set_external_model(ex_models.PressureDropSingleOutput())
+
+        igraph = IncidenceGraphInterface(m, include_inequality=False)
+        var_dm_partition, con_dm_partition = igraph.dulmage_mendelsohn()
+
+        uc_var = var_dm_partition.unmatched + var_dm_partition.underconstrained
+        uc_con = con_dm_partition.underconstrained
+        oc_var = var_dm_partition.overconstrained
+        oc_con = con_dm_partition.overconstrained + con_dm_partition.unmatched
+
+        egb_var = ComponentSet(m.egb.component_data_objects(pyo.Var))
+        self.assertEqual(ComponentSet(uc_var), egb_var)
+
+        uc_cons_set = ComponentSet([m.egb.output_constraints["Pout"]])
+        self.assertEqual(ComponentSet(uc_con), uc_cons_set)
+
+        self.assertEqual(ComponentSet(oc_var), ComponentSet([]))
+        self.assertEqual(ComponentSet(oc_con), ComponentSet([]))
+
+        max_matching = igraph.maximum_matching()
+        self.assertIn(max_matching[m.egb.output_constraints["Pout"]], egb_var)
+
+        cc_vars, cc_cons = igraph.get_connected_components()
+        self.assertEqual(ComponentSet(cc_vars[0]), egb_var)
+        self.assertEqual(cc_cons[0][0].name, "egb.output_constraints[Pout]")
+
+    def test_pressure_drop_single_output_block_triangularization(self):
+        m = pyo.ConcreteModel()
+        m.egb = ExternalGreyBoxBlock()
+        m.egb.set_external_model(ex_models.PressureDropSingleOutput())
+
+        # Add constraints to make model square, then rebuild graph to test block triangularization
+        m.con1 = pyo.Constraint(expr=m.egb.inputs["Pin"] == 1)
+        m.con2 = pyo.Constraint(expr=m.egb.inputs["c"] == 1)
+        m.con3 = pyo.Constraint(expr=m.egb.inputs["F"] == 1)
+        igraph = IncidenceGraphInterface(m, include_inequality=False)
+        bt_vars, bt_cons = igraph.block_triangularize()
+
+        # Expect 4 decomposable sub-sets, one for each linking constraint and one for the grey box
+        self.assertEqual(len(bt_vars), 4)
+        self.assertEqual(len(bt_cons), 4)
+
+        var_set_0 = [m.egb.inputs["Pin"]]
+        var_set_1 = [m.egb.inputs["c"]]
+        var_set_2 = [m.egb.inputs["F"]]
+        var_set_3 = [m.egb.outputs["Pout"]]
+        expected_bt_vars = [var_set_0, var_set_1, var_set_2, var_set_3]
+
+        con_set_0 = [m.con1]
+        con_set_1 = [m.con2]
+        con_set_2 = [m.con3]
+        con_set_3 = [m.egb.output_constraints["Pout"]]
+        expected_bt_cons = [con_set_0, con_set_1, con_set_2, con_set_3]
+
+        self.assertEqual(bt_vars, expected_bt_vars)
+        self.assertEqual(bt_cons, expected_bt_cons)
+
+        self.assertIs(bt_vars[0][0], m.egb.inputs["Pin"])
+        self.assertIs(bt_vars[1][0], m.egb.inputs["c"])
+        self.assertIs(bt_vars[2][0], m.egb.inputs["F"])
+        self.assertIs(bt_vars[3][0], m.egb.outputs["Pout"])
+
+        self.assertIs(bt_cons[0][0], m.con1)
+        self.assertIs(bt_cons[1][0], m.con2)
+        self.assertIs(bt_cons[2][0], m.con3)
+        self.assertIs(bt_cons[3][0], m.egb.output_constraints["Pout"])
+
+        self.assertEqual(
+            ComponentSet(igraph.get_adjacent_to(m.egb.output_constraints["Pout"])),
+            ComponentSet(m.egb.component_data_objects(pyo.Var)),
+        )
+
+    def test_pressure_drop_two_equalities_two_outputs(self):
+        m = pyo.ConcreteModel()
+        m.egb = ExternalGreyBoxBlock()
+        m.egb.set_external_model(ex_models.PressureDropTwoEqualitiesTwoOutputs())
+
+        igraph = IncidenceGraphInterface(m, include_inequality=False)
+        var_dm_partition, con_dm_partition = igraph.dulmage_mendelsohn()
+
+        uc_var = var_dm_partition.unmatched + var_dm_partition.underconstrained
+        uc_con = con_dm_partition.underconstrained
+        oc_var = var_dm_partition.overconstrained
+        oc_con = con_dm_partition.overconstrained + con_dm_partition.unmatched
+
+        uc_var_set = ComponentSet(
+            [
+                m.egb.inputs["F"],
+                m.egb.inputs["P1"],
+                m.egb.inputs["P3"],
+                m.egb.inputs["Pin"],
+                m.egb.inputs["c"],
+                m.egb.outputs["P2"],
+                m.egb.outputs["Pout"],
+            ]
+        )
+        self.assertEqual(ComponentSet(uc_var), uc_var_set)
+
+        uc_con_set = ComponentSet(
+            [
+                m.egb.output_constraints["Pout"],
+                m.egb.output_constraints["P2"],
+                m.egb.eq_constraints["pdrop1"],
+                m.egb.eq_constraints["pdrop3"],
+            ]
+        )
+        self.assertEqual(ComponentSet(uc_con), uc_con_set)
+
+        self.assertEqual(ComponentSet(oc_var), ComponentSet([]))
+        self.assertEqual(ComponentSet(oc_con), ComponentSet([]))
+
+        max_matching = igraph.maximum_matching()
+        egb_var = ComponentSet(m.egb.component_data_objects(pyo.Var))
+        egb_cons = ComponentSet(m.egb.component_data_objects(ExternalGreyBoxConstraint))
+        self.assertIn(max_matching[m.egb.output_constraints["Pout"]], egb_var)
+
+        cc_vars, cc_cons = igraph.get_connected_components()
+        self.assertEqual(ComponentSet(cc_vars[0]), egb_var)
+        self.assertEqual(ComponentSet(cc_cons[0]), egb_cons)
+
+    def test_pressure_drop_two_equalities_two_outputs_block_triangularization(self):
+        m = pyo.ConcreteModel()
+        m.egb = ExternalGreyBoxBlock()
+        m.egb.set_external_model(ex_models.PressureDropTwoEqualitiesTwoOutputs())
+
+        # Add constraints to make model square, then rebuild graph to test block triangularization
+        m.con1 = pyo.Constraint(expr=m.egb.inputs["F"] == 1)
+        m.con2 = pyo.Constraint(expr=m.egb.inputs["Pin"] == 1)
+        m.con3 = pyo.Constraint(expr=m.egb.inputs["c"] == 1)
+        igraph = IncidenceGraphInterface(m, include_inequality=False)
+        bt_vars, bt_cons = igraph.block_triangularize()
+
+        matching = {
+            m.con1: m.egb.inputs["F"],
+            m.con2: m.egb.inputs["Pin"],
+            m.con3: m.egb.inputs["c"],
+            m.egb.eq_constraints["pdrop1"]: m.egb.inputs["P1"],
+            m.egb.eq_constraints["pdrop3"]: m.egb.inputs["P3"],
+            m.egb.output_constraints["P2"]: m.egb.outputs["P2"],
+            m.egb.output_constraints["Pout"]: m.egb.outputs["Pout"],
+        }
+
+        seen = ComponentSet()
+        for vars, cons in zip(bt_vars, bt_cons):
+            self.assertEqual(len(vars), 1)
+            self.assertIs(vars[0], matching[cons[0]])
+            seen.update(vars)
+            # We know that P1 has to come before P2 and P3 in the block triangular form
+            if vars[0] is m.egb.outputs["P2"] or vars[0] is m.egb.inputs["P3"]:
+                self.assertIn(m.egb.inputs["P1"], seen)
+
+        # We know that these constraints have to be in the first three blocks
+        self.assertEqual(
+            set(bt_cons[i][0] for i in range(3)), set([m.con1, m.con2, m.con3])
+        )