diff --git a/pyomo/contrib/incidence_analysis/incidence.py b/pyomo/contrib/incidence_analysis/incidence.py index 9e0000ee2a4..08d8ab749a3 100644 --- a/pyomo/contrib/incidence_analysis/incidence.py +++ b/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: diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 86da622cc13..9869b8b643d 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/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,13 +30,7 @@ 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, @@ -43,10 +38,18 @@ ) 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 " diff --git a/pyomo/contrib/incidence_analysis/tests/test_external_grey_box_integration.py b/pyomo/contrib/incidence_analysis/tests/test_external_grey_box_integration.py new file mode 100644 index 00000000000..3af44ee080a --- /dev/null +++ b/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]) + ) diff --git a/pyomo/contrib/incidence_analysis/tests/test_incidence.py b/pyomo/contrib/incidence_analysis/tests/test_incidence.py index 36bcb783218..5b930028e54 100644 --- a/pyomo/contrib/incidence_analysis/tests/test_incidence.py +++ b/pyomo/contrib/incidence_analysis/tests/test_incidence.py @@ -14,7 +14,11 @@ from pyomo.contrib.incidence_analysis.incidence import ( IncidenceMethod, get_incident_variables, + get_variables_incident_to_constraint, ) +from pyomo.contrib.incidence_analysis.config import get_config_from_kwds +from pyomo.repn.ampl import AMPLRepnVisitor +from unittest.mock import MagicMock class TestAssumedBehavior(unittest.TestCase): @@ -256,5 +260,247 @@ def _get_incident_variables(self, expr, **kwds): return get_incident_variables(expr, method=method, **kwds) +class TestGetVariablesIncidentToConstraint(unittest.TestCase): + """ + Unit tests for the get_variables_incident_to_constraint helper function. + This function was extracted to handle the logic of dispatching to different + methods for identifying incident variables. + """ + + def test_identify_variables_method(self): + """Test using identify_variables method""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + expr = m.x[1] + m.x[1] * m.x[2] + m.x[1] * pyo.exp(m.x[3]) + + variables = get_variables_incident_to_constraint( + expr, method=IncidenceMethod.identify_variables, include_fixed=False + ) + self.assertEqual(ComponentSet(variables), ComponentSet(m.x[:])) + + def test_identify_variables_with_fixed(self): + """Test identify_variables method with include_fixed=True""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], initialize=1.0) + expr = m.x[1] + m.x[1] * m.x[2] + m.x[1] * pyo.exp(m.x[3]) + m.x[2].fix() + + # With include_fixed=True, should include all variables + variables = get_variables_incident_to_constraint( + expr, method=IncidenceMethod.identify_variables, include_fixed=True + ) + self.assertEqual(ComponentSet(variables), ComponentSet(m.x[:])) + + # With include_fixed=False, should exclude fixed variables + variables = get_variables_incident_to_constraint( + expr, method=IncidenceMethod.identify_variables, include_fixed=False + ) + self.assertEqual(ComponentSet(variables), ComponentSet([m.x[1], m.x[3]])) + + def test_standard_repn_method(self): + """Test using standard_repn method""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + expr = m.x[1] + 2 * m.x[2] + 3 * m.x[3] ** 2 + + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.standard_repn, + include_fixed=False, + linear_only=False, + ) + self.assertEqual(ComponentSet(variables), ComponentSet(m.x[:])) + + def test_standard_repn_linear_only(self): + """Test standard_repn method with linear_only=True""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + expr = 2 * m.x[1] + 2 * m.x[2] * m.x[3] + 3 * m.x[2] + + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.standard_repn, + include_fixed=False, + linear_only=True, + ) + # Only x[1] is purely linear + self.assertEqual(ComponentSet(variables), ComponentSet([m.x[1]])) + + def test_standard_repn_compute_values_method(self): + """Test using standard_repn_compute_values method""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + m.p = pyo.Param([1, 2], mutable=True, initialize=1.0) + expr = m.p[1] * m.x[1] + m.p[2] * m.x[2] + m.x[3] ** 2 + + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.standard_repn_compute_values, + include_fixed=False, + linear_only=False, + ) + self.assertEqual(ComponentSet(variables), ComponentSet(m.x[:])) + + def test_standard_repn_compute_values_zero_coefficient(self): + """Test standard_repn_compute_values filters zero coefficients""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + m.p = pyo.Param([1, 2], mutable=True, initialize=1.0) + m.p[1].set_value(0) + expr = m.p[1] * m.x[1] + m.p[2] * m.x[2] + m.x[3] ** 2 + + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.standard_repn_compute_values, + include_fixed=False, + linear_only=False, + ) + # x[1] should be filtered out due to zero coefficient + self.assertEqual(ComponentSet(variables), ComponentSet([m.x[2], m.x[3]])) + + def test_ampl_repn_method(self): + """Test using ampl_repn method""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + expr = m.x[1] + 2 * m.x[2] + 3 * m.x[3] ** 2 + + # Create AMPLRepnVisitor + config = get_config_from_kwds(method=IncidenceMethod.ampl_repn) + visitor = config._ampl_repn_visitor + + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.ampl_repn, + linear_only=False, + amplrepnvisitor=visitor, + ) + self.assertEqual(ComponentSet(variables), ComponentSet(m.x[:])) + + def test_ampl_repn_linear_only(self): + """Test ampl_repn method with linear_only=True""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + expr = 2 * m.x[1] + 2 * m.x[2] * m.x[3] + 3 * m.x[2] + + # Create AMPLRepnVisitor + config = get_config_from_kwds(method=IncidenceMethod.ampl_repn) + visitor = config._ampl_repn_visitor + + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.ampl_repn, + linear_only=True, + amplrepnvisitor=visitor, + ) + # Only x[1] is purely linear + self.assertEqual(ComponentSet(variables), ComponentSet([m.x[1]])) + + def test_egb_constraint_body(self): + """Test special handling for EGBConstraintBody expressions""" + # Create a mock EGBConstraintBody with a get_incident_variables method + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + + mock_egb_body = MagicMock() + mock_egb_body.get_incident_variables.return_value = [m.x[1], m.x[2]] + + # Import EGBConstraintBody to use isinstance check + from pyomo.contrib.pynumero.interfaces.external_grey_box_constraint import ( + EGBConstraintBody, + ) + + # Make the mock an instance of EGBConstraintBody + mock_egb_body.__class__ = EGBConstraintBody + + variables = get_variables_incident_to_constraint( + mock_egb_body, + method=IncidenceMethod.standard_repn, # Method is ignored for EGBConstraintBody + include_fixed=False, + linear_only=False, + ) + + # Should call get_incident_variables on the EGBConstraintBody object + mock_egb_body.get_incident_variables.assert_called_once() + self.assertEqual(variables, [m.x[1], m.x[2]]) + + def test_invalid_method_raises_error(self): + """Test that invalid method raises ValueError""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + expr = m.x[1] + m.x[2] + + # Create an invalid method (not a real IncidenceMethod) + invalid_method = "not_a_real_method" + + with self.assertRaises(ValueError) as cm: + get_variables_incident_to_constraint( + expr, method=invalid_method, include_fixed=False, linear_only=False + ) + + self.assertIn("Unrecognized value", str(cm.exception)) + self.assertIn("for the method used to identify incident", str(cm.exception)) + + def test_standard_repn_with_fixed_variables(self): + """Test standard_repn with include_fixed parameter""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], initialize=1.0) + expr = m.x[1] + m.x[2] + m.x[3] ** 2 + m.x[2].fix() + + # With include_fixed=True, temporarily unfix variables + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.standard_repn, + include_fixed=True, + linear_only=False, + ) + self.assertEqual(ComponentSet(variables), ComponentSet(m.x[:])) + # Variable should still be fixed after the call + self.assertTrue(m.x[2].fixed) + + # With include_fixed=False, exclude fixed variables + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.standard_repn, + include_fixed=False, + linear_only=False, + ) + self.assertEqual(ComponentSet(variables), ComponentSet([m.x[1], m.x[3]])) + + def test_zero_coefficient_filtering(self): + """Test that variables with zero coefficients are properly filtered""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + + # Expression where x[1] cancels out + expr = m.x[1] + m.x[2] * m.x[3] - m.x[1] + + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.standard_repn, + include_fixed=False, + linear_only=False, + ) + # x[1] should be filtered out + self.assertEqual(ComponentSet(variables), ComponentSet([m.x[2], m.x[3]])) + + def test_nonlinear_variables_filtering(self): + """Test filtering of variables that only appear nonlinearly when linear_only=True""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + + # x[1] is linear, x[2] is both linear and nonlinear, x[3] is only nonlinear + expr = 2 * m.x[1] + 3 * m.x[2] + m.x[2] * m.x[3] + m.x[3] ** 2 + + variables = get_variables_incident_to_constraint( + expr, + method=IncidenceMethod.standard_repn, + include_fixed=False, + linear_only=True, + ) + # Only x[1] should be included (purely linear) + self.assertEqual(ComponentSet(variables), ComponentSet([m.x[1]])) + + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/mindtpy/tests/test_mindtpy_grey_box.py b/pyomo/contrib/mindtpy/tests/test_mindtpy_grey_box.py index 4c1ed2f3621..7083efac124 100644 --- a/pyomo/contrib/mindtpy/tests/test_mindtpy_grey_box.py +++ b/pyomo/contrib/mindtpy/tests/test_mindtpy_grey_box.py @@ -14,6 +14,9 @@ from pyomo.environ import SolverFactory, value, maximize from pyomo.opt import TerminationCondition from pyomo.common.dependencies import numpy_available, scipy_available + +if not (numpy_available and scipy_available): + raise unittest.SkipTest("Pynumero needs scipy and numpy to run NLP tests") from pyomo.contrib.mindtpy.tests.MINLP_simple import SimpleMINLP as SimpleMINLP model_list = [SimpleMINLP(grey_box=True)] diff --git a/pyomo/contrib/pynumero/interfaces/external_grey_box.py b/pyomo/contrib/pynumero/interfaces/external_grey_box.py index 8674be34cd7..ce824dc2e68 100644 --- a/pyomo/contrib/pynumero/interfaces/external_grey_box.py +++ b/pyomo/contrib/pynumero/interfaces/external_grey_box.py @@ -6,24 +6,25 @@ # 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 abc import logging -import numpy as np -from scipy.sparse import coo_matrix -from pyomo.common.dependencies import numpy as np from pyomo.common.deprecation import RenamedClass from pyomo.common.log import is_debug_set from pyomo.common.timing import ConstructionTimer -from pyomo.core.base import Var, Set, Constraint, value -from pyomo.core.base.block import BlockData, Block, declare_custom_block +from pyomo.core.base import Var, Set +from pyomo.core.base.block import BlockData, Block from pyomo.core.base.global_set import UnindexedComponent_index from pyomo.core.base.initializer import Initializer from pyomo.core.base.set import UnindexedComponent_set from pyomo.core.base.reference import Reference -from ..sparse.block_matrix import BlockMatrix +from pyomo.contrib.pynumero.interfaces.external_grey_box_constraint import ( + ExternalGreyBoxConstraint, +) logger = logging.getLogger('pyomo.contrib.pynumero') @@ -407,9 +408,26 @@ def set_external_model(self, external_grey_box_model, inputs=None, outputs=None) # call the callback so the model can set initialization, bounds, etc. external_grey_box_model.finalize_block_construction(self) + # Construct the ExternalGreyBoxConstraint objects + self._construct_implicit_constraints() + def get_external_model(self): return self._ex_model + def _construct_implicit_constraints(self): + """ + Construct the implicit constraints for this block. This should be + called by the solver interface before solving to ensure that the + implicit constraints are constructed and available on the block. + """ + # Let the EGBConstraints infer names from the indexing sets + self._equality_constraint_set = Set( + initialize=self._equality_constraint_names, ordered=True + ) + + self.eq_constraints = ExternalGreyBoxConstraint(self._equality_constraint_set) + self.output_constraints = ExternalGreyBoxConstraint(self._output_names_set) + class ExternalGreyBoxBlock(Block): _ComponentDataClass = ExternalGreyBoxBlockData diff --git a/pyomo/contrib/pynumero/interfaces/external_grey_box_constraint.py b/pyomo/contrib/pynumero/interfaces/external_grey_box_constraint.py new file mode 100644 index 00000000000..6827425aa58 --- /dev/null +++ b/pyomo/contrib/pynumero/interfaces/external_grey_box_constraint.py @@ -0,0 +1,774 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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. +# ___________________________________________________________________________________ +""" +This module implements the ExternalGreyBoxConstraint class, which is used to +represent implicit constraints defined by external grey-box models within Pyomo. +""" + +from __future__ import annotations +from operator import index +import sys +import logging +from collections.abc import Mapping +from typing import Union, Type +from weakref import ref as weakref_ref + +from pyomo.common.pyomo_typing import overload +from pyomo.common.formatting import tabular_writer +from pyomo.common.modeling import NOTSET + +from pyomo.core.expr.numvalue import value +from pyomo.core.base.component import ComponentData, ModelComponentFactory +from pyomo.core.base.global_set import UnindexedComponent_index +from pyomo.core.base.indexed_component import IndexedComponent, UnindexedComponent_set +from pyomo.core.base.disable_methods import disable_methods + +logger = logging.getLogger('pyomo.contrib.pynumero') + + +JAC_ZERO_TOLERANCE = 1e-8 + + +class EGBConstraintBody: + """ + This class creates a representation of the "body" of an implicit constraint in + an ExternalGreyBox model. + + Currently, this supports: + * evaluation of the residual of the implicit constraint + * identification of incident variables in the implicit constraint + """ + + def __init__(self, parent_model, implicit_constraint_id): + self._parent_model = parent_model + self._implicit_constraint_id = implicit_constraint_id + + self._ext_output_idx = None + self._ext_eq_cons_idx = None + + ext_model = parent_model.get_external_model() + + if self._implicit_constraint_id in ext_model.equality_constraint_names(): + self._ext_eq_cons_idx = ext_model.equality_constraint_names().index( + self._implicit_constraint_id + ) + elif self._implicit_constraint_id in ext_model.output_names(): + self._ext_output_idx = ext_model.output_names().index( + self._implicit_constraint_id + ) + else: + raise ValueError( + f"Implicit_constraint_id '{self._implicit_constraint_id}' is not a valid identifier in " + f"the external model." + ) + + def get_output_var(self): + """ + If this EGBConstraintBody corresponds to an output variable, return the corresponding Pyomo VarData object. + Otherwise, return None. + """ + if self._ext_output_idx is not None: + out_var = list(self._parent_model.outputs.values())[self._ext_output_idx] + return out_var + return None + + @property + def is_numeric_type(self): + """ + Returns True if the body of this constraint is a numeric type (i.e., it can be evaluated to a number). + """ + return True + + def __call__(self, exception=NOTSET): + """Compute the value of the body of this constraint.""" + if self._ext_eq_cons_idx is not None: + # For an implicit constraint, return the residual + try: + return self._parent_model.get_external_model().evaluate_equality_constraints()[ + self._ext_eq_cons_idx + ] + except Exception as e: + raise RuntimeError( + f"Error evaluating implicit equality constraint '{self._implicit_constraint_id}' " + "in external model. Have the external model inputs been set?" + ) from e + # For an output, the ExternalGreyBox will always return the value + # of the output as a function of the inputs. + evaluated_value = self._parent_model.get_external_model().evaluate_outputs()[ + self._ext_output_idx + ] + var_value = value(self.get_output_var(), exception=exception) + return var_value - evaluated_value + + def get_incident_variables(self): + """ + Get the variables that are incident on this implicit constraint. + + Returns + ------- + list of VarData + List containing the variables that participate in the expression + + """ + # There are (at least) three ways incident variables could be defined for an implicit constraint: + # 1) We consider all inputs to the external model to be incident on the implicit constraint + # 2) We trust the shape of the Jacobian returned by the ExternalGreyBoxModel to contain all elements + # that could possibly be non-zero (i.e. incident) + # 3) We consider only the inputs that have a non-zero Jacobian entry for the implicit constraint + # to be incident on the implicit constraint + # For now, we will go with option 2 as it is the most general. + ext_model = self._parent_model.get_external_model() + incident_variables = [] + + # First, if this implicit constraint corresponds to an output variable, we need to include that output + # variable as an incident variable since it is part of the expression for the implicit constraint. + out_var = self.get_output_var() + if out_var is not None: + incident_variables.append(out_var) + + # Next, get the Jacobian for the external model + try: + if self._ext_eq_cons_idx is not None: + jac = ext_model.evaluate_jacobian_equality_constraints().tocsr() + con_idx = self._ext_eq_cons_idx + else: + jac = ext_model.evaluate_jacobian_outputs().tocsr() + con_idx = self._ext_output_idx + except Exception as e: + raise RuntimeError( + f"Error evaluating Jacobian for external model when getting incident variables for implicit constraint " + f"'{self._implicit_constraint_id}'. Original error message: {str(e)}" + ) from e + + # Get all variables with entries for this constraint + # We do not check value, as we assume entries indicate potential incident variables, + # however they may currently have zero derivative values. + var_indices = jac.getrow(con_idx).indices + incident_variables.extend( + list(self._parent_model.inputs.values())[j] for j in var_indices + ) + + return incident_variables + + +class ExternalGreyBoxConstraintData(ComponentData): + """This class defines the data for a single algebraic constraint. + + Parameters + ---------- + expr : ExpressionBase + The Pyomo expression stored in this constraint. + + component : ExternalGreyBoxConstraint + The ExternalGreyBoxConstraint object that owns this data. + + """ + + __slots__ = ('_implicit_constraint_id', '_body', '_index') + + def __init__(self, implicit_constraint_id=None, component=None, index=NOTSET): + # + # These lines represent in-lining of the + # following constructors: + # - ExternalGreyBoxConstraintData + # - ComponentData + self._component = weakref_ref(component) if (component is not None) else None + self._index = index + self._implicit_constraint_id = implicit_constraint_id + + # Placeholder for body + self._body = None + + def __call__(self, exception=NOTSET): + """Compute the value of the body of this constraint.""" + body = value(self.body, exception=exception) + return body + + def to_bounded_expression(self, *args, **kwargs): + """Duck-type method from ConstraintData. + + Raises + ------ + + TypeError + Always. ExternalGreyBoxConstraints do not have an explicit expression. + + """ + raise TypeError( + "ExternalGreyBoxConstraints do not have an explicit expression." + ) + + @property + def body(self): + """Value (residual) of the implicit ExternalGreyBoxConstraint.""" + if self._body is None: + # Create the EGBConstraintBody object + self._body = EGBConstraintBody( + parent_model=self.parent_block(), + implicit_constraint_id=self._implicit_constraint_id, + ) + return self._body + + @property + def lower(self): + """The lower bound of a ExternalGreyBoxConstraint. + + Implicit constraints always have a lower bound of 0. + + """ + return 0.0 + + @property + def upper(self): + """Access the upper bound of a ExternalGreyBoxConstraint. + + Implicit constraints always have an upper bound of 0. + + """ + return 0.0 + + @property + def lb(self): + """float : the value of the lower bound of a ExternalGreyBoxConstraint expression. + + Implicit constraints always have a lower bound of 0. + """ + return 0.0 + + @property + def ub(self): + """float : the value of the upper bound of a ExternalGreyBoxConstraint expression. + + Implicit constraints always have an upper bound of 0. + """ + return 0.0 + + @property + def equality(self): + """bool : True. ExternalGreyBoxConstraints are always equalities.""" + return True + + @property + def strict_lower(self): + """bool : True if this ExternalGreyBoxConstraint has a strict lower bound.""" + return False + + @property + def strict_upper(self): + """bool : True if this ExternalGreyBoxConstraint has a strict upper bound.""" + return False + + def has_lb(self): + """Returns :const:`True`. Implicit constraints always have a lower bound.""" + return True + + def has_ub(self): + """Returns :const:`True`. Implicit constraints always have an upper bound.""" + return True + + @property + def expr(self): + """Return the expression associated with this ExternalGreyBoxConstraint. + + Raises: + TypeError + Always. ExternalGreyBoxConstraints do not have an explicit expression. + """ + raise TypeError( + "ExternalGreyBoxConstraints do not have an explicit expression." + ) + + def get_value(self): + """Get the expression on this ExternalGreyBoxConstraint. + + Raises: + TypeError + Always. ExternalGreyBoxConstraints do not have an explicit expression. + """ + return self.expr + + def set_value(self, expr): + """Set the expression on this ExternalGreyBoxConstraint. + + Raises: + TypeError + Always. ExternalGreyBoxConstraints do not have an explicit expression. + """ + raise TypeError( + "ExternalGreyBoxConstraints do not have an explicit expression." + ) + + def lslack(self): + """ + Returns the value of f(x)-L for ExternalGreyBoxConstraints of the form: + L <= f(x) (<= U) + (U >=) f(x) >= L + """ + return value(self.body) + + def uslack(self): + """ + Returns the value of U-f(x) for ExternalGreyBoxConstraints of the form: + (L <=) f(x) <= U + U >= f(x) (>= L) + """ + return -value(self.body) + + def slack(self): + """ + Returns the smaller of lslack and uslack values + """ + return -abs(value(self.body)) + + # Duck-typing a few common Constraint methods and properties + @property + def active(self): + """bool : True if this ExternalGreyBoxConstraint is active.""" + return self.parent_block().active + + def activate(self): + """Raise a TypeError, as ExternalGreyBoxConstraints cannot be activated or deactivated.""" + raise TypeError( + "ExternalGreyBoxConstraints cannot be activated or deactivated individually. " + "Activate or deactivate the parent ExternalGreyBoxBlock instead." + ) + + def deactivate(self): + """Raise a TypeError, as ExternalGreyBoxConstraints cannot be activated or deactivated.""" + # Refer back to the activate method to ensure message consistency. + self.activate() + + +@ModelComponentFactory.register("General ExternalGreyBoxConstraint expressions.") +class ExternalGreyBoxConstraint(IndexedComponent): + """ + This modeling component defines a ExternalGreyBoxConstraint for either an + implicit equality constraint or output variable in an ExternalGreyBox model. + + Constructor arguments: + implicit_constraint_id + The identifier for this implicit constraint or output variable + name + A name for this component + doc + A text string describing this component + + Public class attributes: + doc + A text string describing this component + name + A name for this component + active + A boolean that is true if this component will be used to + construct a model instance + implicit_constraint_id + The identifier for this implicit constraint or output variable + + Private class attributes: + _constructed + A boolean that is true if this component has been constructed + _data + A dictionary from the index set to component data objects + _index + The set of valid indices + _model + A weakref to the model that owns this component + _parent + A weakref to the parent block that owns this component + _type + The class type for the derived subclass + """ + + _ComponentDataClass = ExternalGreyBoxConstraintData + + @overload + def __new__( + cls: Type[ScalarExternalGreyBoxConstraint], *args, **kwds + ) -> ScalarExternalGreyBoxConstraint: ... + + @overload + def __new__( + cls: Type[IndexedExternalGreyBoxConstraint], *args, **kwds + ) -> IndexedExternalGreyBoxConstraint: ... + + @overload + def __new__( + cls: Type[ExternalGreyBoxConstraint], *args, **kwds + ) -> Union[ScalarExternalGreyBoxConstraint, IndexedExternalGreyBoxConstraint]: ... + + def __new__(cls, *args, **kwds): + if cls != ExternalGreyBoxConstraint: + return super().__new__(cls) + if not args or (args[0] is UnindexedComponent_set and len(args) == 1): + return super().__new__(AbstractScalarExternalGreyBoxConstraint) + return super().__new__(IndexedExternalGreyBoxConstraint) + + @overload + def __init__( + self, + *indexes, + expr=None, + rule=None, + implicit_constraint_id: str | dict[tuple, str] | None = None, + name=None, + doc=None, + ): ... + + def __init__(self, *args, **kwargs): + # Get id of the implicit constraint (either the equality_constraint_name or output_name) + self._implicit_constraint_id = kwargs.pop('implicit_constraint_id', None) + + # Check for normal Constraint arguments, and raise a TypeError if found + rule = kwargs.pop('rule', None) + expr = kwargs.pop('expr', None) + + if rule is not None: + raise TypeError( + "The 'rule' argument is not supported by ExternalGreyBoxConstraint. " + "Use the 'implicit_constraint_id' argument instead." + ) + if expr is not None: + raise TypeError( + "The 'expr' argument is not supported by ExternalGreyBoxConstraint. " + "ExternalGreyBoxConstraints do not have explicit expressions." + ) + + kwargs.setdefault('ctype', ExternalGreyBoxConstraint) + IndexedComponent.__init__(self, *args, **kwargs) + + # Validate implicit_constraint_id + if not self.is_indexed(): + # Scalar EGBConstraint - ensure that implicit_constraint_id is provided + if self._implicit_constraint_id is None: + raise ValueError( + "The 'implicit_constraint_id' argument must be provided for non-indexed ExternalGreyBoxConstraints." + ) + elif isinstance(self._implicit_constraint_id, Mapping): + # Indexed EGBConstraint with implicit_constraint_id mapping - ensure that keys match index set + valid_indices = set(self.index_set()) + provided_indices = set(self._implicit_constraint_id.keys()) + + missing = valid_indices - provided_indices + extra = provided_indices - valid_indices + + if missing or extra: + msg = ( + "For indexed ExternalGreyBoxConstraints, the " + "'implicit_constraint_id' mapping keys must exactly match " + "the component index set." + ) + if missing: + msg += f" Missing keys: {sorted(missing, key=str)}." + if extra: + msg += f" Invalid keys: {sorted(extra, key=str)}." + raise ValueError(msg) + elif self._implicit_constraint_id is not None: + # Indexed EGBConstraint with unexpected implicit_constraint_id type - raise an error + raise TypeError( + "For indexed ExternalGreyBoxConstraints, the 'implicit_constraint_id' argument must be a mapping from index to identifier or None." + ) + + def construct(self, data=None): + """ + Construct the ExternalGreyBoxConstraint. + """ + if self._constructed: + return + + # First, check that the parent_block is an ExternalGreyBoxBlock + if self.parent_block() is None or not hasattr( + self.parent_block(), "get_external_model" + ): + raise ValueError( + "ExternalGreyBoxConstraint components must be " + "added to an ExternalGreyBoxBlock." + ) + + super().construct(data=data) + + self._constructed = True + + def _pprint(self): + """ + Return data that will be printed for this component. + """ + return ( + [ + ("Size", len(self)), + ("Index", self._index_set if self.is_indexed() else None), + ("Active", self.active), + ], + self.items, + ("Lower", "Body", "Upper", "Active"), + lambda k, v: [ + "-Inf" if v.lower is None else v.lower, + v.body, + "+Inf" if v.upper is None else v.upper, + v.active, + ], + ) + + @property + def implicit_constraint_id(self): + """ + Identifier for this implicit constraint or output variable in the external model. + """ + return self._implicit_constraint_id + + def display(self, prefix="", ostream=None): + """ + Print component state information + + This duplicates logic in Component.pprint() + """ + if not self.active: + return + if ostream is None: + ostream = sys.stdout + tab = " " + ostream.write(prefix + self.local_name + " : ") + ostream.write("Size=" + str(len(self))) + + ostream.write("\n") + tabular_writer( + ostream, + prefix + tab, + ((k, v) for k, v in self._data.items() if v.active), + ("Lower", "Body", "Upper"), + lambda k, v: [ + value(v.lower, exception=False), + value(v.body, exception=False), + value(v.upper, exception=False), + ], + ) + + +class ScalarExternalGreyBoxConstraint( + ExternalGreyBoxConstraintData, ExternalGreyBoxConstraint +): + """ + ScalarExternalGreyBoxConstraint is the implementation representing a single, + non-indexed ExternalGreyBoxConstraint. + """ + + def __init__(self, *args, **kwds): + ExternalGreyBoxConstraintData.__init__(self, component=self) + ExternalGreyBoxConstraint.__init__(self, *args, **kwds) + self._index = UnindexedComponent_index + + # Set _data here, as it isn't getting set elsewhere + self._data[None] = self + + def construct(self, data=None): + if self._constructed: + return + + ExternalGreyBoxConstraint.construct(self, data=data) + # Validate implicit_constraint_id for this scalar constraint + _validate_implicit_constraint_id(self, self._implicit_constraint_id) + + # + # Singleton ExternalGreyBoxConstraints are strange in that we want them to be + # both be constructed but have len() == 0 when not initialized with + # anything (at least according to the unit tests that are + # currently in place). So during initialization only, we will + # treat them as "indexed" objects where things like + # Constraint.Skip are managed. But after that they will behave + # like ExternalGreyBoxConstraintData objects where set_value does not handle + # Constraint.Skip but expects a valid expression or None. + # + @property + def body(self): + """The body (variable portion) of a ExternalGreyBoxConstraint expression.""" + if not self._data: + raise ValueError( + f"Accessing the body of ScalarExternalGreyBoxConstraint " + f"'{self.name}' before the ExternalGreyBoxConstraint has been assigned. " + "There is currently nothing to access." + ) + return ExternalGreyBoxConstraintData.body.fget(self) + + @property + def lower(self): + """The lower bound of a ExternalGreyBoxConstraint expression. + + This is the fixed lower bound of a ExternalGreyBoxConstraint as a Pyomo + expression. This may contain potentially variable terms + that are currently fixed. If there is no lower bound, this will + return `None`. + + """ + if not self._data: + raise ValueError( + f"Accessing the lower bound of ScalarExternalGreyBoxConstraint " + f"'{self.name}' before the ExternalGreyBoxConstraint has been assigned. " + "There is currently nothing to access." + ) + return ExternalGreyBoxConstraintData.lower.fget(self) + + @property + def upper(self): + """Access the upper bound of a ExternalGreyBoxConstraint expression. + + This is the fixed upper bound of a ExternalGreyBoxConstraint as a Pyomo + expression. This may contain potentially variable terms + that are currently fixed. If there is no upper bound, this will + return `None`. + + """ + if not self._data: + raise ValueError( + f"Accessing the upper bound of ScalarExternalGreyBoxConstraint " + f"'{self.name}' before the ExternalGreyBoxConstraint has been assigned. " + "There is currently nothing to access." + ) + return ExternalGreyBoxConstraintData.upper.fget(self) + + @property + def equality(self): + """bool : True if this is an equality ExternalGreyBoxConstraint.""" + if not self._data: + raise ValueError( + f"Accessing the equality flag of ScalarExternalGreyBoxConstraint " + f"'{self.name}' before the ExternalGreyBoxConstraint has been assigned. " + "There is currently nothing to access." + ) + return ExternalGreyBoxConstraintData.equality.fget(self) + + @property + def strict_lower(self): + """bool : True if this ExternalGreyBoxConstraint has a strict lower bound.""" + if not self._data: + raise ValueError( + f"Accessing the strict_lower flag of ScalarExternalGreyBoxConstraint " + f"'{self.name}' before the ExternalGreyBoxConstraint has been assigned. " + "There is currently nothing to access." + ) + return ExternalGreyBoxConstraintData.strict_lower.fget(self) + + @property + def strict_upper(self): + """bool : True if this ExternalGreyBoxConstraint has a strict upper bound.""" + if not self._data: + raise ValueError( + f"Accessing the strict_upper flag of ScalarExternalGreyBoxConstraint " + f"'{self.name}' before the ExternalGreyBoxConstraint has been assigned. " + "There is currently nothing to access." + ) + return ExternalGreyBoxConstraintData.strict_upper.fget(self) + + def clear(self): + self._data = {} + + def set_value(self, expr): + """Set the expression on this ExternalGreyBoxConstraint.""" + if not self._data: + self._data[None] = self + return super().set_value(expr) + + # + # Leaving this method for backward compatibility reasons. + # (probably should be removed) + # + def add(self, index, expr): + """Add a ExternalGreyBoxConstraint with a given index.""" + if index is not None: + raise ValueError( + f"ScalarExternalGreyBoxConstraint object '{self.name}' does not accept " + f"index values other than None. Invalid value: {index}" + ) + self.set_value(expr) + return self + + +@disable_methods( + { + '__call__', + 'add', + 'set_value', + 'to_bounded_expression', + 'expr', + 'body', + 'lower', + 'upper', + 'equality', + 'strict_lower', + 'strict_upper', + } +) +class AbstractScalarExternalGreyBoxConstraint(ScalarExternalGreyBoxConstraint): + """ + Implementation of abstract ExternalGreyBoxConstraints. + """ + + +class IndexedExternalGreyBoxConstraint(ExternalGreyBoxConstraint): + """ + Implementation of indexed ExternalGreyBoxConstraints. + """ + + # + # Leaving this method for backward compatibility reasons + # + # Note: Beginning after Pyomo 5.2 this method will now validate that + # the index is in the underlying index set (through 5.2 the index + # was not checked). + # + def add(self, index, expr): + """Add a ExternalGreyBoxConstraint with a given index.""" + return self.__setitem__(index, expr) + + def construct(self, data=None): + super().construct(data=data) + + if not self._data: + for idx in self.index_set(): + # Get implicit_constraint_id for this index + if self._implicit_constraint_id is not None: + # If implicit_constraint_id is provided, get the corresponding id for this index + implicit_constraint_id = self._implicit_constraint_id[idx] + else: + # Infer implicit_constraint_id from index (assumes index is a single value that can be converted to a string) + implicit_constraint_id = idx + + # Check that implicit_constraint_id is a valid constraint in the external model + _validate_implicit_constraint_id(self, implicit_constraint_id) + + self._data[idx] = self._ComponentDataClass( + component=self, + implicit_constraint_id=implicit_constraint_id, + index=idx, + ) + + @overload + def __getitem__(self, index) -> ExternalGreyBoxConstraintData: ... + + __getitem__ = IndexedComponent.__getitem__ # type: ignore + + +def _validate_implicit_constraint_id(obj, implicit_constraint_id): + if not isinstance(implicit_constraint_id, str): + raise TypeError( + "ExternalGreyBoxConstraint implicit_constraint_id values must be strings. " + f"Invalid value: {implicit_constraint_id!r}" + ) + + external_model = obj.parent_block().get_external_model() + + if not ( + implicit_constraint_id in external_model.equality_constraint_names() + or implicit_constraint_id in external_model.output_names() + ): + raise ValueError( + f"implicit_constraint_id '{implicit_constraint_id}' does not exist in the " + f"external model associated with ExternalGreyBoxBlock '{obj.parent_block().name}'." + ) diff --git a/pyomo/contrib/pynumero/interfaces/pyomo_grey_box_nlp.py b/pyomo/contrib/pynumero/interfaces/pyomo_grey_box_nlp.py index bddc8167cc1..6426f0b0ddf 100644 --- a/pyomo/contrib/pynumero/interfaces/pyomo_grey_box_nlp.py +++ b/pyomo/contrib/pynumero/interfaces/pyomo_grey_box_nlp.py @@ -6,19 +6,19 @@ # 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. +# ___________________________________________________________________________________ """ This module defines the classes that provide an NLP interface based on the Ampl Solver Library (ASL) implementation """ -import os import numpy as np import logging from scipy.sparse import coo_matrix, identity -from pyomo.common.deprecation import deprecated import pyomo.core.base as pyo -from pyomo.common.collections import ComponentMap from pyomo.common.modeling import unique_component_name from pyomo.contrib.pynumero.sparse.block_matrix import BlockMatrix from pyomo.contrib.pynumero.sparse.block_vector import BlockVector @@ -31,6 +31,9 @@ from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock from pyomo.contrib.pynumero.interfaces.nlp_projections import ProjectedNLP from pyomo.core.base.suffix import SuffixFinder +from pyomo.contrib.pynumero.interfaces.external_grey_box_constraint import ( + ExternalGreyBoxConstraint, +) # Todo: make some of the numpy arrays not writable from __init__ @@ -77,18 +80,34 @@ def __init__(self, pyomo_model): if number_of_objectives == 0: pyomo_model.del_component(objective) + # With Vars we need to account for Vars that are not part of the Block + # but appear within Constraints + # First, get all the Vars from the model self._pyomo_model_var_names_to_datas = { v.getname(fully_qualified=True): v for v in pyomo_model.component_data_objects( ctype=pyo.Var, descend_into=True ) } + # Next, check the PyomoNLP for any Vars that are missing + # This can occur if a constraint in the model references a Var that is not part of the model + for v in self._pyomo_nlp.get_pyomo_variables(): + self._pyomo_model_var_names_to_datas[ + v.getname(fully_qualified=True) + ] = v + self._pyomo_model_constraint_names_to_datas = { c.getname(fully_qualified=True): c for c in pyomo_model.component_data_objects( - ctype=pyo.Constraint, descend_into=True + ctype=pyo.Constraint, descend_into=True, active=True ) } + # Check for ExternalGreyBoxConstraint objects and add + # them too + for c in pyomo_model.component_data_objects( + ExternalGreyBoxConstraint, active=True, descend_into=True + ): + self._pyomo_model_constraint_names_to_datas[c.name] = c finally: # Restore the ctypes of the ExternalGreyBoxBlock components @@ -506,6 +525,36 @@ def load_state_into_pyomo(self, bound_multipliers=None): def has_hessian_support(self): return self._has_hessian_support + # Compatibility API for PyomoNLP - this is only a partial implementation + def get_pyomo_variables(self): + return self._pyomo_model_var_datas + + def get_pyomo_constraints(self): + return list(self._pyomo_model_constraint_names_to_datas.values()) + + def get_pyomo_equality_constraints(self): + return [c for c in self.get_pyomo_constraints() if c.equality] + + def get_pyomo_inequality_constraints(self): + return [c for c in self.get_pyomo_constraints() if not c.equality] + + def get_primal_indices(self, var): + # get the name of the variable + var_name = var.getname(fully_qualified=True) + # get the index of the variable in the primals + try: + return self._primals_names.index(var_name) + except ValueError: + raise ValueError(f'Variable {var_name} not found in primals.') + + def get_constraint_indices(self, constraint): + constraint_name = constraint.getname(fully_qualified=True) + # get the index of the constraint in the constraints + try: + return self._constraint_names.index(constraint_name) + except ValueError: + raise ValueError(f'Constraint {constraint_name} not found in constraints.') + def _default_if_none(value, default): if value is None: @@ -528,7 +577,6 @@ def __init__(self, external_grey_box_block): self._obj_factor = 1.0 n_inputs = len(self._block.inputs) assert n_inputs == self._ex_model.n_inputs() - n_eq_constraints = self._ex_model.n_equality_constraints() n_outputs = len(self._block.outputs) assert n_outputs == self._ex_model.n_outputs() @@ -542,23 +590,15 @@ def __init__(self, external_grey_box_block): self._block.outputs[k].getname(fully_qualified=True) for k in self._block.outputs ) - n_primals = len(self._primals_names) - prefix = self._block.getname(fully_qualified=True) - self._constraint_names = [ - '{}.{}'.format(prefix, nm) - for nm in self._ex_model.equality_constraint_names() - ] - output_var_names = [ - self._block.outputs[k].getname(fully_qualified=False) - for k in self._block.outputs - ] - self._constraint_names.extend( - [ - '{}.output_constraints[{}]'.format(prefix, nm) - for nm in self._ex_model.output_names() - ] - ) + # Collect implicit constraints. + self._constraint_names = [] + self._constraint_datas = [] + for c in self._block.component_data_objects( + ExternalGreyBoxConstraint, active=True, descend_into=False + ): + self._constraint_names.append(c.getname(fully_qualified=True)) + self._constraint_datas.append(c) # create the numpy arrays of bounds on the primals self._primals_lb = BlockVector(2) @@ -642,6 +682,9 @@ def n_constraints(self): def constraint_names(self): return list(self._constraint_names) + def constraint_datas(self): + return list(self._constraint_datas) + def nnz_jacobian(self): if self._nnz_jacobian is None: J = self.evaluate_jacobian() diff --git a/pyomo/contrib/pynumero/interfaces/tests/external_grey_box_models.py b/pyomo/contrib/pynumero/interfaces/tests/external_grey_box_models.py index 19daf1b0602..6c20702b85d 100644 --- a/pyomo/contrib/pynumero/interfaces/tests/external_grey_box_models.py +++ b/pyomo/contrib/pynumero/interfaces/tests/external_grey_box_models.py @@ -14,6 +14,7 @@ scipy_available, ) from pyomo.common.dependencies.scipy import sparse as spa +import pyomo.common.unittest as unittest if not (numpy_available and scipy_available): raise unittest.SkipTest("Pynumero needs scipy and numpy to run NLP tests") diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_dynamic_model.py b/pyomo/contrib/pynumero/interfaces/tests/test_dynamic_model.py index 9c97feaaa0a..1d288e9bd6a 100644 --- a/pyomo/contrib/pynumero/interfaces/tests/test_dynamic_model.py +++ b/pyomo/contrib/pynumero/interfaces/tests/test_dynamic_model.py @@ -556,16 +556,16 @@ def test_compare_evaluations(self): 'h20', ] c2list = [ - 'egb.h1bal_1', - 'egb.h1bal_2', - 'egb.h1bal_3', - 'egb.h1bal_4', - 'egb.h1bal_5', - 'egb.h2bal_1', - 'egb.h2bal_2', - 'egb.h2bal_3', - 'egb.h2bal_4', - 'egb.h2bal_5', + 'egb.eq_constraints[h1bal_1]', + 'egb.eq_constraints[h1bal_2]', + 'egb.eq_constraints[h1bal_3]', + 'egb.eq_constraints[h1bal_4]', + 'egb.eq_constraints[h1bal_5]', + 'egb.eq_constraints[h2bal_1]', + 'egb.eq_constraints[h2bal_2]', + 'egb.eq_constraints[h2bal_3]', + 'egb.eq_constraints[h2bal_4]', + 'egb.eq_constraints[h2bal_5]', 'egb.output_constraints[F12_0]', 'egb.output_constraints[F12_1]', 'egb.output_constraints[F12_2]', diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_constraint.py b/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_constraint.py new file mode 100644 index 00000000000..3d195f87954 --- /dev/null +++ b/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_constraint.py @@ -0,0 +1,1769 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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 unittest.mock import patch + +from pyomo.contrib.pynumero.dependencies import numpy as np + +from pyomo.contrib.pynumero.interfaces.external_grey_box import ( + ExternalGreyBoxBlock, + ExternalGreyBoxBlockData, +) +from pyomo.contrib.pynumero.interfaces.external_grey_box_constraint import ( + ExternalGreyBoxConstraint, + ExternalGreyBoxConstraintData, + ScalarExternalGreyBoxConstraint, + EGBConstraintBody, +) +import pyomo.contrib.pynumero.interfaces.tests.external_grey_box_models as ex_models + + +def _no_op_construct_implicit_constraints(self): + """ + No-op implementation of _construct_implicit_constraints. + + This prevents automatic construction of ExternalGreyBoxConstraints, + allowing tests to manually construct them to test the constraint + construction logic itself. + """ + pass + + +# Decorator to patch _construct_implicit_constraints with a no-op for test classes +def skip_implicit_constraint_construction(test_class): + """ + Decorator that patches _construct_implicit_constraints to prevent + automatic construction of implicit constraints in tests. + + This is scoped to the decorated test class only and won't affect other tests. + """ + return patch.object( + ExternalGreyBoxBlockData, + '_construct_implicit_constraints', + _no_op_construct_implicit_constraints, + )(test_class) + + +@skip_implicit_constraint_construction +class TestExternalGreyBoxConstraintConstruction(unittest.TestCase): + """Test construction and initialization of ExternalGreyBoxConstraint.""" + + def test_construction_without_implicit_constraint_id_raises(self): + """Test that constructing without implicit_constraint_id raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + with self.assertRaises(ValueError) as context: + m.egb.c = ExternalGreyBoxConstraint() + self.assertIn("implicit_constraint_id", str(context.exception)) + + def test_construction_with_rule_raises(self): + """Test that passing 'rule' argument raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + with self.assertRaises(TypeError) as context: + m.egb.c = ExternalGreyBoxConstraint( + implicit_constraint_id='pdrop', rule=lambda m: None + ) + self.assertIn("rule", str(context.exception)) + + def test_construction_with_expr_raises(self): + """Test that passing 'expr' argument raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + with self.assertRaises(TypeError) as context: + m.egb.c = ExternalGreyBoxConstraint( + implicit_constraint_id='pdrop', expr=pyo.Constraint.Skip + ) + self.assertIn("expr", str(context.exception)) + + def test_construction_with_invalid_implicit_constraint_id_raises(self): + """Test that invalid implicit_constraint_id raises ValueError on construct.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + # Construction happens automatically when adding to block + with self.assertRaises(ValueError) as context: + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='nonexistent') + self.assertIn("does not exist", str(context.exception)) + + def test_construction_not_in_external_grey_box_block_raises(self): + """Test that construction outside ExternalGreyBoxBlock raises an error.""" + m = pyo.ConcreteModel() + + # Construction happens automatically when added to block + # This should raise either ValueError or AttributeError depending on validation order + with self.assertRaises((ValueError, AttributeError)) as context: + m.c = ExternalGreyBoxConstraint(implicit_constraint_id='test') + # Check that error message indicates the problem is related to ExternalGreyBoxBlock + self.assertTrue( + "ExternalGreyBoxBlock" in str(context.exception) + or "get_external_model" in str(context.exception) + ) + + def test_scalar_construction_with_equality_constraint(self): + """Test scalar constraint construction with equality constraint.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertIsInstance(m.egb.c, ScalarExternalGreyBoxConstraint) + self.assertEqual(m.egb.c.implicit_constraint_id, 'pdrop') + self.assertTrue(m.egb.c.active) + + def test_scalar_construction_with_output(self): + """Test scalar constraint construction with output variable.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleOutput() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + self.assertIsInstance(m.egb.c, ScalarExternalGreyBoxConstraint) + self.assertEqual(m.egb.c.implicit_constraint_id, 'Pout') + + +class TestExternalGreyBoxConstraintProperties(unittest.TestCase): + """Test properties of ExternalGreyBoxConstraint.""" + + def test_body_with_equality_constraint(self): + """Test body property returns EGBConstraintBody that evaluates to residual.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + # Set input values directly on external model: Pin=100, c=2, F=3, Pout=50 + # Expected residual: Pout - (Pin - 4*c*F^2) = 50 - (100 - 4*2*9) = 50 - 28 = 22 + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + + body_value = pyo.value(m.egb.c.body) + self.assertAlmostEqual(body_value, 22.0, places=6) + + def test_body_with_output(self): + """Test body property returns EGBConstraintBody that evaluates residual for outputs.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleOutput() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + # Set input values directly on external model + # Pin=100, c=2, F=3 => Pout_evaluated = 100 - 4*2*3^2 = 100 - 72 = 28 + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + + # Set output variable to match the evaluated value + m.egb.outputs['Pout'].set_value(28.0) + + # For outputs, when variable matches evaluated value, residual should be 0.0 + body_value = pyo.value(m.egb.c.body) + self.assertAlmostEqual(body_value, 0.0, places=6) + + def test_body_with_invalid_constraint_id_raises(self): + """Test body property raises ValueError for invalid constraint ID.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + # Create constraint with valid id, then manually change it + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + m.egb.c._implicit_constraint_id = 'invalid_id' + + with self.assertRaises(ValueError) as context: + _ = m.egb.c.body + self.assertIn("invalid_id", str(context.exception)) + + def test_body_without_inputs_set_evaluates_with_defaults(self): + """Test body property evaluates with default zero inputs when not explicitly set.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + # External model initializes with zeros, so evaluation should work + # Expected: 0 - (0 - 4*0*0) = 0 + body_value = pyo.value(m.egb.c.body) + self.assertAlmostEqual(body_value, 0.0, places=6) + + def test_lower_property(self): + """Test lower bound is always 0.0.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertEqual(m.egb.c.lower, 0.0) + + def test_upper_property(self): + """Test upper bound is always 0.0.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertEqual(m.egb.c.upper, 0.0) + + def test_lb_property(self): + """Test lb (lower bound value) is always 0.0.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertEqual(m.egb.c.lb, 0.0) + + def test_ub_property(self): + """Test ub (upper bound value) is always 0.0.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertEqual(m.egb.c.ub, 0.0) + + def test_equality_property(self): + """Test equality property is always True.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertTrue(m.egb.c.equality) + + def test_strict_lower_property(self): + """Test strict_lower is always False.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertFalse(m.egb.c.strict_lower) + + def test_strict_upper_property(self): + """Test strict_upper is always False.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertFalse(m.egb.c.strict_upper) + + def test_has_lb_method(self): + """Test has_lb() returns True.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertTrue(m.egb.c.has_lb()) + + def test_has_ub_method(self): + """Test has_ub() returns True.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertTrue(m.egb.c.has_ub()) + + def test_expr_property_raises(self): + """Test expr property raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + with self.assertRaises(TypeError) as context: + _ = m.egb.c.expr + self.assertIn("do not have an explicit expression", str(context.exception)) + + def test_get_value_raises(self): + """Test get_value() raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + with self.assertRaises(TypeError) as context: + m.egb.c.get_value() + self.assertIn("do not have an explicit expression", str(context.exception)) + + def test_set_value_raises(self): + """Test set_value() raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + with self.assertRaises(TypeError) as context: + m.egb.c.set_value(None) + self.assertIn("do not have an explicit expression", str(context.exception)) + + def test_to_bounded_expression_raises(self): + """Test to_bounded_expression() raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + with self.assertRaises(TypeError) as context: + m.egb.c.to_bounded_expression() + self.assertIn("do not have an explicit expression", str(context.exception)) + + +@skip_implicit_constraint_construction +class TestEGBConstraintBody(unittest.TestCase): + """Test the EGBConstraintBody object returned by the body property.""" + + def test_body_returns_egb_constraint_body_object(self): + """Test that body property returns an EGBConstraintBody object.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + body_obj = m.egb.c.body + self.assertIsInstance(body_obj, EGBConstraintBody) + + def test_body_object_is_numeric_type(self): + """Test that EGBConstraintBody object has is_numeric_type property set to True.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + body_obj = m.egb.c.body + self.assertTrue(body_obj.is_numeric_type) + + def test_body_object_can_be_called(self): + """Test that EGBConstraintBody object can be called directly to get residual.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + + body_obj = m.egb.c.body + body_value = body_obj() + self.assertAlmostEqual(body_value, 22.0, places=6) + + def test_body_object_with_pyo_value(self): + """Test that pyo.value() works with EGBConstraintBody object.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + + body_obj = m.egb.c.body + body_value = pyo.value(body_obj) + self.assertAlmostEqual(body_value, 22.0, places=6) + + def test_body_object_caching(self): + """Test that body property returns the same object on repeated access.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + body_obj1 = m.egb.c.body + body_obj2 = m.egb.c.body + self.assertIs(body_obj1, body_obj2) + + def test_body_object_evaluates_with_different_inputs(self): + """Test that body object evaluates correctly with different external model inputs.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + body_obj = m.egb.c.body + + # First evaluation + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + body_value1 = pyo.value(body_obj) + self.assertAlmostEqual(body_value1, 22.0, places=6) + + # Second evaluation with different inputs + external_model.set_input_values(np.asarray([100, 2, 3, 28], dtype=np.float64)) + body_value2 = pyo.value(body_obj) + self.assertAlmostEqual(body_value2, 0.0, places=6) + + def test_body_object_for_output_constraint(self): + """Test EGBConstraintBody object for output-based constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleOutput() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + # Pin=100, c=2, F=3 => Pout_evaluated = 100 - 4*2*3^2 = 28 + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + + # Set output variable to match evaluated value + m.egb.outputs['Pout'].set_value(28.0) + + body_obj = m.egb.c.body + # For outputs, when variable matches evaluated value, residual should be 0 + self.assertAlmostEqual(pyo.value(body_obj), 0.0, places=6) + + def test_output_constraint_residual_when_variable_matches_evaluated(self): + """Test that residual is zero when output variable value matches evaluated value.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleOutput() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + # Set inputs: Pin=100, c=2, F=3 => Pout_evaluated = 100 - 4*2*3^2 = 28 + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + + # Set the output variable to match the evaluated value + evaluated_value = external_model.evaluate_outputs()[0] + m.egb.outputs['Pout'].set_value(evaluated_value) + + # Residual should be: var_value - evaluated_value = 28 - 28 = 0 + residual = pyo.value(m.egb.c.body) + self.assertAlmostEqual(residual, 0.0, places=6) + + def test_output_constraint_residual_when_variable_does_not_match_evaluated(self): + """Test that residual is non-zero when output variable value does not match evaluated value.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleOutput() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + # Set inputs: Pin=100, c=2, F=3 => Pout_evaluated = 100 - 4*2*3^2 = 28 + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + + # Set the output variable to a different value than evaluated + m.egb.outputs['Pout'].set_value(50.0) + + # Residual should be: var_value - evaluated_value = 50 - 28 = 22 + residual = pyo.value(m.egb.c.body) + self.assertAlmostEqual(residual, 22.0, places=6) + + # Test with another non-matching value + m.egb.outputs['Pout'].set_value(20.0) + + # Residual should be: var_value - evaluated_value = 20 - 28 = -8 + residual = pyo.value(m.egb.c.body) + self.assertAlmostEqual(residual, -8.0, places=6) + + def test_output_constraint_residual_updates_with_inputs_and_variable(self): + """Test that residual updates correctly when inputs or output variable changes.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleOutput() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + # Initial inputs: Pin=100, c=2, F=3 => Pout_evaluated = 100 - 4*2*9 = 28 + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + m.egb.outputs['Pout'].set_value(30.0) + + # Residual should be: 30 - 28 = 2 + residual = pyo.value(m.egb.c.body) + self.assertAlmostEqual(residual, 2.0, places=6) + + # Change inputs: Pin=100, c=2, F=2 => Pout_evaluated = 100 - 4*2*4 = 68 + external_model.set_input_values(np.asarray([100, 2, 2], dtype=np.float64)) + + # Residual should be: 30 - 68 = -38 + residual = pyo.value(m.egb.c.body) + self.assertAlmostEqual(residual, -38.0, places=6) + + # Now set variable to match the new evaluated value + m.egb.outputs['Pout'].set_value(68.0) + + # Residual should be: 68 - 68 = 0 + residual = pyo.value(m.egb.c.body) + self.assertAlmostEqual(residual, 0.0, places=6) + + def test_body_object_invalid_constraint_id_raises(self): + """Test that body object raises error for invalid constraint ID.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + # Directly create body object with invalid constraint ID + with self.assertRaises(ValueError) as context: + body_obj = EGBConstraintBody(m.egb, 'invalid_constraint_id') + self.assertIn("invalid_constraint_id", str(context.exception)) + + def test_get_incident_variables_without_jacobian(self): + """Test get_incident_variables returns all input variables.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + body_obj = m.egb.c.body + incident_vars = body_obj.get_incident_variables() + + # Should return all 4 input variables + self.assertEqual(len(incident_vars), 4) + expected_names = ['Pin', 'c', 'F', 'Pout'] + actual_names = [var.name for var in incident_vars] + self.assertEqual( + actual_names, [f'egb.inputs[{name}]' for name in expected_names] + ) + + def test_get_incident_variables_with_output_constraint(self): + """Test get_incident_variables for output-based constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + body_obj = m.egb.c.body + + incident_vars = body_obj.get_incident_variables() + self.assertEqual(len(incident_vars), 4) + expected_names = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.outputs[Pout]', + ] + assert all(var.name in expected_names for var in incident_vars) + + def test_get_incident_variables_multiple_outputs(self): + """Test get_incident_variables for different output constraints in same model.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.egb.c1 = ExternalGreyBoxConstraint(implicit_constraint_id='P2') + m.egb.c2 = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + + body_obj1 = m.egb.c1.body + body_obj2 = m.egb.c2.body + + incident_vars1 = body_obj1.get_incident_variables() + incident_vars2 = body_obj2.get_incident_variables() + + self.assertEqual(len(incident_vars1), 4) + self.assertEqual(len(incident_vars2), 4) + + # Compare variable names + for v in incident_vars1: + expected1 = [ + "egb.inputs[Pin]", + "egb.inputs[c]", + "egb.inputs[F]", + "egb.outputs[P2]", + ] + assert v.name in expected1 + for v in incident_vars2: + expected2 = [ + "egb.inputs[Pin]", + "egb.inputs[c]", + "egb.inputs[F]", + "egb.outputs[Pout]", + ] + assert v.name in expected2 + + def test_get_incident_variables_multiple_constraints_and_outputs(self): + """Test get_incident_variables for different implicit constraints in same model.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + # Manually create the implicit constraint objects + m.egb.pdrop1 = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop1') + m.egb.pdrop3 = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop3') + m.egb.P2_constraint = ExternalGreyBoxConstraint(implicit_constraint_id='P2') + m.egb.Pout_constraint = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + # Implicit constraint: 'pdrop1' + body_obj1 = m.egb.pdrop1.body + incident_vars1 = body_obj1.get_incident_variables() + self.assertEqual(len(incident_vars1), 4) + expected_names = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + ] + for v in incident_vars1: + self.assertIn(v.name, expected_names) + + # Implicit constraint: 'pdrop3' + body_obj1 = m.egb.pdrop3.body + incident_vars1 = body_obj1.get_incident_variables() + self.assertEqual(len(incident_vars1), 4) + expected_names = [ + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + ] + for v in incident_vars1: + self.assertIn(v.name, expected_names) + + # Implicit constraint: 'P2_constraint' + body_obj1 = m.egb.P2_constraint.body + incident_vars1 = body_obj1.get_incident_variables() + self.assertEqual(len(incident_vars1), 4) + expected_names = [ + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.outputs[P2]', + ] + for v in incident_vars1: + self.assertIn(v.name, expected_names) + + # Implicit constraint: 'Pout_constraint' + body_obj1 = m.egb.Pout_constraint.body + incident_vars1 = body_obj1.get_incident_variables() + self.assertEqual(len(incident_vars1), 4) + expected_names = [ + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[Pin]', + 'egb.outputs[Pout]', + ] + for v in incident_vars1: + self.assertIn(v.name, expected_names) + + def test_get_incident_variables_default_parameters(self): + """Test get_incident_variables with default parameters (use_jacobian=False).""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + body_obj = m.egb.c.body + # Call without parameters (should default to use_jacobian=False) + incident_vars = body_obj.get_incident_variables() + + # Should return all 4 input variables + self.assertEqual(len(incident_vars), 4) + + def test_get_incident_variables_returns_var_data_objects(self): + """Test that get_incident_variables returns actual variable data objects.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + body_obj = m.egb.c.body + incident_vars = body_obj.get_incident_variables() + + # Verify that each element is a Pyomo VarData object + for var in incident_vars: + self.assertTrue(hasattr(var, 'value')) + self.assertTrue(hasattr(var, 'fixed')) + self.assertTrue(hasattr(var, 'lb')) + self.assertTrue(hasattr(var, 'ub')) + + +@skip_implicit_constraint_construction +class TestExternalGreyBoxConstraintSlack(unittest.TestCase): + """Test slack methods of ExternalGreyBoxConstraint.""" + + def test_lslack(self): + """Test lslack() returns body value.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + + lslack_value = m.egb.c.lslack() + body_value = pyo.value(m.egb.c.body) + self.assertAlmostEqual(lslack_value, body_value, places=6) + + def test_uslack(self): + """Test uslack() returns negative body value.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + + uslack_value = m.egb.c.uslack() + body_value = pyo.value(m.egb.c.body) + self.assertAlmostEqual(uslack_value, -body_value, places=6) + + def test_slack(self): + """Test slack() returns negative absolute value of body.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + + slack_value = m.egb.c.slack() + body_value = pyo.value(m.egb.c.body) + self.assertAlmostEqual(slack_value, -abs(body_value), places=6) + + def test_call_method(self): + """Test __call__() method returns body value.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + + call_value = m.egb.c() + body_value = pyo.value(m.egb.c.body) + self.assertAlmostEqual(call_value, body_value, places=6) + + +@skip_implicit_constraint_construction +class TestExternalGreyBoxConstraintActive(unittest.TestCase): + """Test active status methods of ExternalGreyBoxConstraint.""" + + def test_active_property(self): + """Test active property follows parent block.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + self.assertTrue(m.egb.c.active) + + m.egb.deactivate() + self.assertFalse(m.egb.c.active) + + m.egb.activate() + self.assertTrue(m.egb.c.active) + + def test_activate_raises(self): + """Test activate() raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + with self.assertRaises(TypeError) as context: + m.egb.c.activate() + self.assertIn("cannot be activated or deactivated", str(context.exception)) + + def test_deactivate_raises(self): + """Test deactivate() raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + with self.assertRaises(TypeError) as context: + m.egb.c.deactivate() + self.assertIn("cannot be activated or deactivated", str(context.exception)) + + +@skip_implicit_constraint_construction +class TestScalarExternalGreyBoxConstraint(unittest.TestCase): + """Test ScalarExternalGreyBoxConstraint specific functionality.""" + + def test_scalar_body_before_assignment_raises(self): + """Test accessing body before assignment raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + m.egb.c.clear() # Clear the data + + with self.assertRaises(ValueError) as context: + _ = m.egb.c.body + self.assertIn( + "before the ExternalGreyBoxConstraint has been assigned", + str(context.exception), + ) + + def test_scalar_lower_before_assignment_raises(self): + """Test accessing lower before assignment raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + m.egb.c.clear() + + with self.assertRaises(ValueError) as context: + _ = m.egb.c.lower + self.assertIn( + "before the ExternalGreyBoxConstraint has been assigned", + str(context.exception), + ) + + def test_scalar_upper_before_assignment_raises(self): + """Test accessing upper before assignment raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + m.egb.c.clear() + + with self.assertRaises(ValueError) as context: + _ = m.egb.c.upper + self.assertIn( + "before the ExternalGreyBoxConstraint has been assigned", + str(context.exception), + ) + + def test_scalar_equality_before_assignment_raises(self): + """Test accessing equality before assignment raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + m.egb.c.clear() + + with self.assertRaises(ValueError) as context: + _ = m.egb.c.equality + self.assertIn( + "before the ExternalGreyBoxConstraint has been assigned", + str(context.exception), + ) + + def test_scalar_strict_lower_before_assignment_raises(self): + """Test accessing strict_lower before assignment raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + m.egb.c.clear() + + with self.assertRaises(ValueError) as context: + _ = m.egb.c.strict_lower + self.assertIn( + "before the ExternalGreyBoxConstraint has been assigned", + str(context.exception), + ) + + def test_scalar_strict_upper_before_assignment_raises(self): + """Test accessing strict_upper before assignment raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + m.egb.c.clear() + + with self.assertRaises(ValueError) as context: + _ = m.egb.c.strict_upper + self.assertIn( + "before the ExternalGreyBoxConstraint has been assigned", + str(context.exception), + ) + + def test_scalar_add_with_invalid_index_raises(self): + """Test add() with non-None index raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + with self.assertRaises(ValueError) as context: + m.egb.c.add(1, None) + self.assertIn( + "does not accept index values other than None", str(context.exception) + ) + + +@skip_implicit_constraint_construction +class TestExternalGreyBoxConstraintMultipleConstraints(unittest.TestCase): + """Test with models having multiple equality constraints.""" + + def test_two_equality_constraints(self): + """Test with model having two equality constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualities() + m.egb.set_external_model(external_model) + + m.egb.c1 = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop2') + m.egb.c2 = ExternalGreyBoxConstraint(implicit_constraint_id='pdropout') + + # Set input values directly on external model: Pin=100, c=2, F=3, P2=82, Pout=64 + external_model.set_input_values( + np.asarray([100, 2, 3, 82, 64], dtype=np.float64) + ) + + # Expected residual for pdrop2: P2 - (Pin - 2*c*F^2) = 82 - (100 - 2*2*9) = 82 - 64 = 18 + body1 = pyo.value(m.egb.c1.body) + self.assertAlmostEqual(body1, 18.0, places=6) + + # Expected residual for pdropout: Pout - (P2 - 2*c*F^2) = 64 - (82 - 2*2*9) = 64 - 46 = 18 + body2 = pyo.value(m.egb.c2.body) + self.assertAlmostEqual(body2, 18.0, places=6) + + def test_two_outputs(self): + """Test with model having two outputs.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.egb.c1 = ExternalGreyBoxConstraint(implicit_constraint_id='P2') + m.egb.c2 = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + # Set input values directly on external model + # Pin=100, c=2, F=3 + # P2_evaluated = 100 - 2*2*9 = 64 + # Pout_evaluated = 100 - 4*2*9 = 28 + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + + # Set output variables to match evaluated values + m.egb.outputs['P2'].set_value(64.0) + m.egb.outputs['Pout'].set_value(28.0) + + # For outputs, when variables match evaluated values, residuals should be 0.0 + self.assertAlmostEqual(pyo.value(m.egb.c1.body), 0.0, places=6) + self.assertAlmostEqual(pyo.value(m.egb.c2.body), 0.0, places=6) + + +@skip_implicit_constraint_construction +class TestExternalGreyBoxConstraintDisplay(unittest.TestCase): + """Test display and printing methods.""" + + def test_display_method(self): + """Test display() method executes without error.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + + # Should not raise an exception + import io + + output = io.StringIO() + m.egb.c.display(ostream=output) + result = output.getvalue() + + # Check that output contains expected elements + self.assertIn('Lower', result) + self.assertIn('Body', result) + self.assertIn('Upper', result) + + def test_display_inactive_does_nothing(self): + """Test display() on inactive component produces no output.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + m.egb.deactivate() + + import io + + output = io.StringIO() + m.egb.c.display(ostream=output) + result = output.getvalue() + + # Should produce empty or minimal output + self.assertEqual(result, '') + + def test_pprint_method(self): + """Test _pprint() returns expected data structure.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + external_model.set_input_values(np.asarray([100, 2, 3, 50], dtype=np.float64)) + + data = m.egb.c._pprint() + + # Check structure + self.assertEqual(len(data), 4) + headers, items_method, columns, formatter = data + + # Check headers + self.assertIn(("Size", 1), headers) + self.assertIn(("Active", True), headers) + + # Check columns + self.assertEqual(columns, ("Lower", "Body", "Upper", "Active")) + + +@skip_implicit_constraint_construction +class TestExternalGreyBoxConstraintImplicitConstraintId(unittest.TestCase): + """Test implicit_constraint_id property.""" + + def test_implicit_constraint_id_property(self): + """Test implicit_constraint_id property returns correct value.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + constraint_id = 'pdrop' + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id=constraint_id) + + self.assertEqual(m.egb.c.implicit_constraint_id, constraint_id) + + def test_implicit_constraint_id_stored_correctly(self): + """Test implicit_constraint_id is stored in _implicit_constraint_id.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualities() + m.egb.set_external_model(external_model) + + m.egb.c1 = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop2') + m.egb.c2 = ExternalGreyBoxConstraint(implicit_constraint_id='pdropout') + + self.assertEqual(m.egb.c1._implicit_constraint_id, 'pdrop2') + self.assertEqual(m.egb.c2._implicit_constraint_id, 'pdropout') + + +@skip_implicit_constraint_construction +class TestExternalGreyBoxConstraintIntegration(unittest.TestCase): + """Integration tests with various external models.""" + + def test_with_hessian_model_equality(self): + """Test with model that supports Hessian evaluation.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEqualityWithHessian() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + external_model.set_input_values(np.asarray([100, 2, 3, 28], dtype=np.float64)) + + # At correct solution, residual should be 0 + body_value = pyo.value(m.egb.c.body) + self.assertAlmostEqual(body_value, 0.0, places=6) + + def test_with_hessian_model_output(self): + """Test with output model that supports Hessian evaluation.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleOutputWithHessian() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + # Pin=100, c=2, F=3 => Pout_evaluated = 100 - 4*2*9 = 28 + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + + # Set output variable to match evaluated value + m.egb.outputs['Pout'].set_value(28.0) + + # For outputs, when variable matches evaluated value, residual should be 0 + self.assertAlmostEqual(pyo.value(m.egb.c.body), 0.0, places=6) + + def test_constraint_in_different_blocks(self): + """Test constraints in multiple ExternalGreyBoxBlocks.""" + m = pyo.ConcreteModel() + + m.egb1 = ExternalGreyBoxBlock() + external_model1 = ex_models.PressureDropSingleEquality() + m.egb1.set_external_model(external_model1) + m.egb1.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + m.egb2 = ExternalGreyBoxBlock() + external_model2 = ex_models.PressureDropSingleOutput() + m.egb2.set_external_model(external_model2) + m.egb2.c = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + # Set inputs for first block + external_model1.set_input_values(np.asarray([100, 2, 3, 28], dtype=np.float64)) + + # Set inputs for second block + # Pin=50, c=1, F=2 => Pout_evaluated = 50 - 4*1*4 = 34 + external_model2.set_input_values(np.asarray([50, 1, 2], dtype=np.float64)) + + # Set output variable for second block to match evaluated value + m.egb2.outputs['Pout'].set_value(34.0) + + # Check both constraints work independently + self.assertAlmostEqual(pyo.value(m.egb1.c.body), 0.0, places=6) + self.assertAlmostEqual(pyo.value(m.egb2.c.body), 0.0, places=6) + + +@skip_implicit_constraint_construction +class TestExternalGreyBoxConstraintEdgeCases(unittest.TestCase): + """Test edge cases and boundary conditions.""" + + def test_constraint_with_zero_inputs(self): + """Test constraint evaluation with zero input values.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + # Set all inputs to zero directly on external model + external_model.set_input_values(np.asarray([0, 0, 0, 0], dtype=np.float64)) + + # Expected residual: 0 - (0 - 4*0*0) = 0 + body_value = pyo.value(m.egb.c.body) + self.assertAlmostEqual(body_value, 0.0, places=6) + + +@skip_implicit_constraint_construction +class TestIndexedExternalGreyBoxConstraint(unittest.TestCase): + """Test indexed ExternalGreyBoxConstraint functionality.""" + + def test_indexed_with_explicit_mapping(self): + """Test indexed constraint with explicit implicit_constraint_id mapping.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout', 'pdrop1', 'pdrop3']) + + # Create indexed constraint with explicit mapping + m.egb.c = ExternalGreyBoxConstraint( + m.set, implicit_constraint_id={i: i for i in m.set} + ) + + # Verify construction + self.assertTrue(m.egb.c.is_indexed()) + self.assertEqual(len(m.egb.c), 4) + + # Verify each index has correct implicit_constraint_id + for idx in m.set: + self.assertEqual(m.egb.c[idx]._implicit_constraint_id, idx) + + def test_indexed_with_implicit_ids(self): + """Test indexed constraint with inferred implicit_constraint_id from index.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout', 'pdrop1', 'pdrop3']) + + # Create indexed constraint without explicit mapping (ids inferred from index) + m.egb.c = ExternalGreyBoxConstraint(m.set) + + # Verify construction + self.assertTrue(m.egb.c.is_indexed()) + self.assertEqual(len(m.egb.c), 4) + + # Verify each index has implicit_constraint_id equal to the index + for idx in m.set: + self.assertEqual(m.egb.c[idx]._implicit_constraint_id, idx) + + def test_indexed_body_evaluation(self): + """Test body evaluation for indexed constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + # Set inputs: Pin=100, c=2, F=3 + # P2_evaluated = 100 - 2*2*9 = 64 + # Pout_evaluated = 100 - 4*2*9 = 28 + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + + # Set output variables + m.egb.outputs['P2'].set_value(70.0) + m.egb.outputs['Pout'].set_value(30.0) + + # Check residuals + # P2: 70 - 64 = 6 + self.assertAlmostEqual(pyo.value(m.egb.c['P2'].body), 6.0, places=6) + # Pout: 30 - 28 = 2 + self.assertAlmostEqual(pyo.value(m.egb.c['Pout'].body), 2.0, places=6) + + def test_indexed_iteration(self): + """Test iterating over indexed constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout', 'pdrop1', 'pdrop3']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + # Test iteration + count = 0 + indices = [] + for idx in m.egb.c: + count += 1 + indices.append(idx) + + self.assertEqual(count, 4) + self.assertEqual(set(indices), set(m.set)) + + def test_indexed_items_method(self): + """Test items() method for indexed constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + # Test items() method + items = list(m.egb.c.items()) + self.assertEqual(len(items), 2) + + for idx, constraint_data in items: + self.assertIn(idx, m.set) + self.assertEqual(constraint_data._implicit_constraint_id, idx) + + def test_indexed_properties(self): + """Test properties work correctly for indexed constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + for idx in m.set: + # Test bounds + self.assertEqual(m.egb.c[idx].lower, 0.0) + self.assertEqual(m.egb.c[idx].upper, 0.0) + self.assertEqual(m.egb.c[idx].lb, 0.0) + self.assertEqual(m.egb.c[idx].ub, 0.0) + + # Test equality + self.assertTrue(m.egb.c[idx].equality) + + # Test strict bounds + self.assertFalse(m.egb.c[idx].strict_lower) + self.assertFalse(m.egb.c[idx].strict_upper) + + # Test has_lb/has_ub + self.assertTrue(m.egb.c[idx].has_lb()) + self.assertTrue(m.egb.c[idx].has_ub()) + + def test_indexed_with_tuple_index(self): + """Test indexed constraint with tuple indices.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set( + initialize=[(1, 'P2'), (1, 'Pout'), (2, 'pdrop1'), (2, 'pdrop3')] + ) + + # Create mapping from tuple indices to string constraint ids + id_map = { + (1, 'P2'): 'P2', + (1, 'Pout'): 'Pout', + (2, 'pdrop1'): 'pdrop1', + (2, 'pdrop3'): 'pdrop3', + } + + m.egb.c = ExternalGreyBoxConstraint(m.set, implicit_constraint_id=id_map) + + # Verify construction + self.assertEqual(len(m.egb.c), 4) + + for idx in m.set: + self.assertEqual(m.egb.c[idx]._implicit_constraint_id, id_map[idx]) + + +@skip_implicit_constraint_construction +class TestIndexedExternalGreyBoxConstraintValidation(unittest.TestCase): + """Test validation errors for indexed ExternalGreyBoxConstraint.""" + + def test_indexed_missing_keys_raises(self): + """Test that missing keys in implicit_constraint_id mapping raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout', 'pdrop1', 'pdrop3']) + + # Create mapping missing some keys + incomplete_map = {'P2': 'P2', 'Pout': 'Pout'} # Missing 'pdrop1' and 'pdrop3' + + with self.assertRaises(ValueError) as context: + m.egb.c = ExternalGreyBoxConstraint( + m.set, implicit_constraint_id=incomplete_map + ) + + self.assertIn("Missing keys", str(context.exception)) + self.assertIn("pdrop1", str(context.exception)) + self.assertIn("pdrop3", str(context.exception)) + + def test_indexed_extra_keys_raises(self): + """Test that extra keys in implicit_constraint_id mapping raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + # Create mapping with extra keys + mapping_with_extras = { + 'P2': 'P2', + 'Pout': 'Pout', + 'extra1': 'pdrop1', + 'extra2': 'pdrop3', + } + + with self.assertRaises(ValueError) as context: + m.egb.c = ExternalGreyBoxConstraint( + m.set, implicit_constraint_id=mapping_with_extras + ) + + self.assertIn("Invalid keys", str(context.exception)) + self.assertIn("extra1", str(context.exception)) + self.assertIn("extra2", str(context.exception)) + + def test_indexed_missing_and_extra_keys_raises(self): + """Test error message when both missing and extra keys exist.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout', 'pdrop1']) + + # Create mapping with missing and extra keys + bad_map = { + 'P2': 'P2', + 'extra': 'Pout', + } # Missing 'Pout' and 'pdrop1', extra 'extra' + + with self.assertRaises(ValueError) as context: + m.egb.c = ExternalGreyBoxConstraint(m.set, implicit_constraint_id=bad_map) + + error_msg = str(context.exception) + self.assertIn("Missing keys", error_msg) + self.assertIn("Invalid keys", error_msg) + + def test_indexed_invalid_type_raises(self): + """Test that invalid type for implicit_constraint_id raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + # Pass a string instead of mapping (invalid for indexed) + with self.assertRaises(TypeError) as context: + m.egb.c = ExternalGreyBoxConstraint(m.set, implicit_constraint_id='P2') + + self.assertIn("must be a mapping", str(context.exception)) + + def test_indexed_invalid_constraint_id_value_raises(self): + """Test that invalid constraint id value in mapping raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + # Create mapping with invalid constraint id + bad_map = { + 'P2': 'P2', + 'Pout': 'invalid_constraint', # This constraint doesn't exist + } + + with self.assertRaises(ValueError) as context: + m.egb.c = ExternalGreyBoxConstraint(m.set, implicit_constraint_id=bad_map) + + self.assertIn("invalid_constraint", str(context.exception)) + self.assertIn("does not exist", str(context.exception)) + + def test_indexed_non_string_constraint_id_raises(self): + """Test that non-string implicit_constraint_id value raises TypeError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + # Create mapping with non-string value + bad_map = {'P2': 123, 'Pout': 'Pout'} # Not a string + + with self.assertRaises(TypeError) as context: + m.egb.c = ExternalGreyBoxConstraint(m.set, implicit_constraint_id=bad_map) + + self.assertIn("must be strings", str(context.exception)) + + def test_indexed_with_inferred_id_invalid_raises(self): + """Test that inferred implicit_constraint_id that doesn't exist raises ValueError.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + # Use indices that don't match any constraint/output names + m.set = pyo.Set(initialize=['invalid1', 'invalid2']) + + with self.assertRaises(ValueError) as context: + m.egb.c = ExternalGreyBoxConstraint(m.set) + + self.assertIn("does not exist", str(context.exception)) + + +@skip_implicit_constraint_construction +class TestIndexedExternalGreyBoxConstraintAdvanced(unittest.TestCase): + """Advanced tests for indexed ExternalGreyBoxConstraint.""" + + def test_indexed_add_method(self): + """Test add() method for indexed constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + # The add method should work the same as __setitem__ + # Test that we can access via indexing + self.assertIsNotNone(m.egb.c['P2']) + self.assertIsNotNone(m.egb.c['Pout']) + + def test_indexed_getitem(self): + """Test __getitem__ for indexed constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + # Test getitem + c_p2 = m.egb.c['P2'] + c_pout = m.egb.c['Pout'] + + self.assertEqual(c_p2._implicit_constraint_id, 'P2') + self.assertEqual(c_pout._implicit_constraint_id, 'Pout') + + def test_indexed_with_different_constraint_types(self): + """Test indexed constraints mixing equality constraints and outputs.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + # Mix of outputs and equality constraints + m.set = pyo.Set(initialize=['P2', 'pdrop1']) + + id_map = {'P2': 'P2', 'pdrop1': 'pdrop1'} # output # equality constraint + + m.egb.c = ExternalGreyBoxConstraint(m.set, implicit_constraint_id=id_map) + + # Set inputs + # Using 6 inputs: Pin, c, F, P1, P3, and one for the missing input of the equality constraint + # Inputs: Pin=100, c=2, F=3, P1=82, P3=46, plus empty output placeholders + external_model.set_input_values( + np.asarray([100, 2, 3, 82, 46], dtype=np.float64) + ) + + # For output P2, set value + m.egb.outputs['P2'].set_value(64.0) + + # Both constraints should be accessible and evaluatable + self.assertIsNotNone(m.egb.c['P2'].body) + self.assertIsNotNone(m.egb.c['pdrop1'].body) + + def test_indexed_component_property(self): + """Test that indexed constraint data has correct component reference.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + for idx in m.set: + constraint_data = m.egb.c[idx] + # Verify parent component reference + self.assertIs(constraint_data.parent_component(), m.egb.c) + + def test_indexed_active_status(self): + """Test active status for indexed constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + # All should be active + for idx in m.set: + self.assertTrue(m.egb.c[idx].active) + + # Deactivate parent block + m.egb.deactivate() + + # All should now be inactive + for idx in m.set: + self.assertFalse(m.egb.c[idx].active) + + def test_indexed_slack_methods(self): + """Test slack methods for indexed constraints.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + # Set inputs + external_model.set_input_values(np.asarray([100, 2, 3], dtype=np.float64)) + m.egb.outputs['P2'].set_value(70.0) + m.egb.outputs['Pout'].set_value(30.0) + + for idx in m.set: + body_val = pyo.value(m.egb.c[idx].body) + + # Test lslack + self.assertAlmostEqual(m.egb.c[idx].lslack(), body_val, places=6) + + # Test uslack + self.assertAlmostEqual(m.egb.c[idx].uslack(), -body_val, places=6) + + # Test slack + self.assertAlmostEqual(m.egb.c[idx].slack(), -abs(body_val), places=6) + + def test_constraint_with_negative_inputs(self): + """Test constraint evaluation with negative input values.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + # Set negative inputs directly on external model + external_model.set_input_values( + np.asarray([-100, -2, -3, -50], dtype=np.float64) + ) + + # Should evaluate without error + body_value = pyo.value(m.egb.c.body) + # Expected: -50 - (-100 - 4*(-2)*(-3)^2) = -50 - (-100 - 4*(-2)*9) = -50 - (-100 + 72) = -50 + 28 = -22 + self.assertAlmostEqual(body_value, -22.0, places=6) + + def test_constraint_with_large_inputs(self): + """Test constraint evaluation with large input values.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropSingleEquality() + m.egb.set_external_model(external_model) + + m.egb.c = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop') + + # Set large inputs directly on external model + external_model.set_input_values( + np.asarray([1e6, 1e3, 1e2, 1e5], dtype=np.float64) + ) + + # Should evaluate without error + body_value = pyo.value(m.egb.c.body) + self.assertIsInstance(body_value, (float, np.floating)) + + +@skip_implicit_constraint_construction +def test_component_data_objects_with_EGBC(): + """Test that ExternalGreyBoxConstraints can be iterated over using component_data_objects.""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputsWithHessian() + m.egb.set_external_model(external_model) + + # Manually create the implicit constraint objects + m.egb.pdrop1 = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop1') + m.egb.pdrop3 = ExternalGreyBoxConstraint(implicit_constraint_id='pdrop3') + m.egb.P2_constraint = ExternalGreyBoxConstraint(implicit_constraint_id='P2') + m.egb.Pout_constraint = ExternalGreyBoxConstraint(implicit_constraint_id='Pout') + + count = 0 + for c in m.egb.component_data_objects( + ctype=ExternalGreyBoxConstraint, descend_into=False + ): + assert isinstance(c, ExternalGreyBoxConstraintData) + assert c.local_name in ['P2_constraint', 'Pout_constraint', 'pdrop1', 'pdrop3'] + count += 1 + assert count == 4 + + +@skip_implicit_constraint_construction +def test_indexed_egbc_no_implicit_constraint_id(): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout', 'pdrop1', 'pdrop3']) + + m.egb.c = ExternalGreyBoxConstraint(m.set) + + assert m.egb.c["P2"]._implicit_constraint_id == "P2" + assert m.egb.c["Pout"]._implicit_constraint_id == "Pout" + assert m.egb.c["pdrop1"]._implicit_constraint_id == "pdrop1" + assert m.egb.c["pdrop3"]._implicit_constraint_id == "pdrop3" + + +@skip_implicit_constraint_construction +def test_indexed_egbc_implicit_constraint_id_mapping(): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputs() + m.egb.set_external_model(external_model) + + m.set = pyo.Set(initialize=['P2', 'Pout', 'pdrop1', 'pdrop3']) + + m.egb.c = ExternalGreyBoxConstraint( + m.set, implicit_constraint_id={i: i for i in m.set} + ) + + assert m.egb.c["P2"]._implicit_constraint_id == "P2" + assert m.egb.c["Pout"]._implicit_constraint_id == "Pout" + assert m.egb.c["pdrop1"]._implicit_constraint_id == "pdrop1" + assert m.egb.c["pdrop3"]._implicit_constraint_id == "pdrop3" + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_model.py b/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_model.py index be7bf311f9f..2676a7357d0 100644 --- a/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_model.py +++ b/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_model.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. +# ___________________________________________________________________________________ import pyomo.common.unittest as unittest import pyomo.environ as pyo diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_model_with_constraints.py b/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_model_with_constraints.py new file mode 100644 index 00000000000..c9910d26be3 --- /dev/null +++ b/pyomo/contrib/pynumero/interfaces/tests/test_external_grey_box_model_with_constraints.py @@ -0,0 +1,876 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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.pynumero.dependencies import ( + numpy as np, + numpy_available, + scipy_available, +) + +if not (numpy_available and scipy_available): + raise unittest.SkipTest("Pynumero needs scipy and numpy to run NLP tests") + +from scipy.sparse import coo_matrix + +from pyomo.contrib.pynumero.asl import AmplInterface + +if not AmplInterface.available(): + raise unittest.SkipTest("ASL interface is not available") + +from pyomo.contrib.pynumero.interfaces.external_grey_box import ( + ExternalGreyBoxBlock, + ExternalGreyBoxModel, +) +from pyomo.contrib.pynumero.interfaces.external_grey_box_constraint import ( + ExternalGreyBoxConstraint, + ExternalGreyBoxConstraintData, +) +from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoGreyBoxNLP +from pyomo.contrib.pynumero.interfaces.tests.compare_utils import ( + check_vectors_specific_order, + check_sparse_matrix_specific_order, +) +import pyomo.contrib.pynumero.interfaces.tests.external_grey_box_models as ex_models +from pyomo.contrib.incidence_analysis import IncidenceGraphInterface + + +class TestExternalGreyBoxModelWithConstraints(unittest.TestCase): + """Tests for ExternalGreyBoxBlock with implicit_constraint_objects""" + + def test_pressure_drop_single_output_constraint_creation(self): + """Test that constraint objects are created for outputs""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.PressureDropSingleOutput()) + + # Check that the constraint object was created for the output + self.assertTrue(hasattr(m.egb, 'eq_constraints')) + self.assertTrue(len(m.egb.eq_constraints) == 0) + self.assertTrue(hasattr(m.egb, 'output_constraints')) + self.assertIsInstance( + m.egb.output_constraints['Pout'], ExternalGreyBoxConstraintData + ) + + # Check that no equality constraint objects were created (no equality constraints) + egbm = m.egb.get_external_model() + eq_con_names = egbm.equality_constraint_names() + self.assertEqual(eq_con_names, []) + + def test_pressure_drop_single_equality_constraint_creation(self): + """Test that constraint objects are created for equality constraints""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.PressureDropSingleEquality()) + + # Check that the constraint object was created for the equality constraint + self.assertTrue(hasattr(m.egb, 'eq_constraints')) + self.assertIsInstance( + m.egb.eq_constraints['pdrop'], ExternalGreyBoxConstraintData + ) + self.assertTrue(hasattr(m.egb, 'output_constraints')) + self.assertTrue(len(m.egb.output_constraints) == 0) + + # Check that no output constraint objects were created (no outputs) + egbm = m.egb.get_external_model() + output_names = egbm.output_names() + self.assertEqual(output_names, []) + + def test_pressure_drop_two_outputs_constraint_creation(self): + """Test that constraint objects are created for multiple outputs""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.PressureDropTwoOutputs()) + + # Check that constraint objects were created for both outputs + self.assertTrue(hasattr(m.egb, 'output_constraints')) + self.assertIsInstance( + m.egb.output_constraints['P2'], ExternalGreyBoxConstraintData + ) + self.assertIsInstance( + m.egb.output_constraints['Pout'], ExternalGreyBoxConstraintData + ) + + def test_pressure_drop_two_equalities_constraint_creation(self): + """Test that constraint objects are created for multiple equality constraints""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.PressureDropTwoEqualities()) + + # Check that constraint objects were created for both equality constraints + self.assertTrue(hasattr(m.egb, 'eq_constraints')) + self.assertIsInstance( + m.egb.eq_constraints['pdrop2'], ExternalGreyBoxConstraintData + ) + self.assertIsInstance( + m.egb.eq_constraints['pdropout'], ExternalGreyBoxConstraintData + ) + + def test_pressure_drop_two_equalities_two_outputs_constraint_creation(self): + """Test that constraint objects are created for both equality constraints and outputs""" + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.PressureDropTwoEqualitiesTwoOutputs()) + + # Check that constraint objects were created for equality constraints + self.assertTrue(hasattr(m.egb, 'eq_constraints')) + self.assertIsInstance( + m.egb.eq_constraints['pdrop1'], ExternalGreyBoxConstraintData + ) + self.assertIsInstance( + m.egb.eq_constraints['pdrop3'], ExternalGreyBoxConstraintData + ) + + # Check that constraint objects were created for outputs + self.assertTrue(hasattr(m.egb, 'output_constraints')) + self.assertIsInstance( + m.egb.output_constraints['P2'], ExternalGreyBoxConstraintData + ) + self.assertIsInstance( + m.egb.output_constraints['Pout'], ExternalGreyBoxConstraintData + ) + + def test_pressure_drop_single_equality_with_constraints(self): + """Test PyomoGreyBoxNLP with single equality constraint and constraint objects""" + self._test_pressure_drop_single_equality( + ex_models.PressureDropSingleEquality(), False + ) + self._test_pressure_drop_single_equality( + ex_models.PressureDropSingleEqualityWithHessian(), True + ) + + def _test_pressure_drop_single_equality(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['Pout'].value = 50 + m.egb.inputs['Pout'].setlb(0) + m.egb.inputs['Pout'].setub(100) + m.obj = pyo.Objective(expr=(m.egb.inputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoGreyBoxNLP(m) + + self.assertEqual(4, pyomo_nlp.n_primals()) + self.assertEqual(1, pyomo_nlp.n_constraints()) + self.assertEqual(4, pyomo_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(3, pyomo_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[Pout]', + ] + x_order = pyomo_nlp.variable_names() + comparison_c_order = ['egb.pdrop'] + c_order = pyomo_nlp.constraint_names() + + xlb = pyomo_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = pyomo_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = pyomo_nlp.constraints_lb() + comparison_clb = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = pyomo_nlp.constraints_ub() + comparison_cub = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = pyomo_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = pyomo_nlp.init_duals() + comparison_duals_init = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(4, len(pyomo_nlp.create_new_vector('primals'))) + self.assertEqual(1, len(pyomo_nlp.create_new_vector('constraints'))) + self.assertEqual(1, len(pyomo_nlp.create_new_vector('duals'))) + self.assertEqual(1, len(pyomo_nlp.create_new_vector('eq_constraints'))) + self.assertEqual(0, len(pyomo_nlp.create_new_vector('ineq_constraints'))) + self.assertEqual(1, len(pyomo_nlp.create_new_vector('duals_eq'))) + self.assertEqual(0, len(pyomo_nlp.create_new_vector('duals_ineq'))) + + pyomo_nlp.set_primals(np.asarray([1, 2, 3, 4], dtype=np.float64)) + x = pyomo_nlp.get_primals() + self.assertTrue(np.array_equal(x, np.asarray([1, 2, 3, 4], dtype=np.float64))) + pyomo_nlp.set_primals(pyomo_nlp.init_primals()) + + pyomo_nlp.set_duals(np.asarray([42], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42], dtype=np.float64))) + pyomo_nlp.set_duals(np.asarray([21], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21], dtype=np.float64))) + + fac = pyomo_nlp.get_obj_factor() + self.assertEqual(fac, 1) + pyomo_nlp.set_obj_factor(42) + self.assertEqual(pyomo_nlp.get_obj_factor(), 42) + pyomo_nlp.set_obj_factor(1) + + f = pyomo_nlp.evaluate_objective() + self.assertEqual(f, 900) + + gradf = pyomo_nlp.evaluate_grad_objective() + comparison_gradf = np.asarray([0, 0, 0, 60], dtype=np.float64) + check_vectors_specific_order( + self, gradf, x_order, comparison_gradf, comparison_x_order + ) + c = pyomo_nlp.evaluate_constraints() + comparison_c = np.asarray([22], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(1) + pyomo_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = pyomo_nlp.evaluate_jacobian() + comparison_j = np.asarray([[-1, 36, 48, 1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + pyomo_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = pyomo_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (4, 4)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0], + [0, 0, 0, 0], + [0, 8 * 3 * 21, 8 * 2 * 21, 0], + [0, 0, 0, 2 * 1], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(AttributeError): + h = pyomo_nlp.evaluate_hessian_lag() + + def test_pressure_drop_two_equalities_with_constraints(self): + """Test PyomoGreyBoxNLP with two equality constraints and constraint objects""" + self._test_pressure_drop_two_equalities( + ex_models.PressureDropTwoEqualities(), False + ) + self._test_pressure_drop_two_equalities( + ex_models.PressureDropTwoEqualitiesWithHessian(), True + ) + + def _test_pressure_drop_two_equalities(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P2'].value = 80 + m.egb.inputs['P2'].setlb(10) + m.egb.inputs['P2'].setub(90) + m.egb.inputs['Pout'].value = 50 + m.egb.inputs['Pout'].setlb(0) + m.egb.inputs['Pout'].setub(100) + m.obj = pyo.Objective(expr=(m.egb.inputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoGreyBoxNLP(m) + + self.assertEqual(5, pyomo_nlp.n_primals()) + self.assertEqual(2, pyomo_nlp.n_constraints()) + self.assertEqual(8, pyomo_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(3, pyomo_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P2]', + 'egb.inputs[Pout]', + ] + x_order = pyomo_nlp.variable_names() + comparison_c_order = ['egb.pdrop2', 'egb.pdropout'] + c_order = pyomo_nlp.constraint_names() + + xlb = pyomo_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 10, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = pyomo_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 90, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = pyomo_nlp.constraints_lb() + comparison_clb = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = pyomo_nlp.constraints_ub() + comparison_cub = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = pyomo_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 80, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = pyomo_nlp.init_duals() + comparison_duals_init = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(5, len(pyomo_nlp.create_new_vector('primals'))) + self.assertEqual(2, len(pyomo_nlp.create_new_vector('constraints'))) + self.assertEqual(2, len(pyomo_nlp.create_new_vector('duals'))) + self.assertEqual(2, len(pyomo_nlp.create_new_vector('eq_constraints'))) + self.assertEqual(0, len(pyomo_nlp.create_new_vector('ineq_constraints'))) + self.assertEqual(2, len(pyomo_nlp.create_new_vector('duals_eq'))) + self.assertEqual(0, len(pyomo_nlp.create_new_vector('duals_ineq'))) + + pyomo_nlp.set_primals(np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + x = pyomo_nlp.get_primals() + self.assertTrue( + np.array_equal(x, np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + ) + pyomo_nlp.set_primals(pyomo_nlp.init_primals()) + + pyomo_nlp.set_duals(np.asarray([42, 10], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42, 10], dtype=np.float64))) + pyomo_nlp.set_duals(np.asarray([21, 5], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21, 5], dtype=np.float64))) + + fac = pyomo_nlp.get_obj_factor() + self.assertEqual(fac, 1) + pyomo_nlp.set_obj_factor(42) + self.assertEqual(pyomo_nlp.get_obj_factor(), 42) + pyomo_nlp.set_obj_factor(1) + + f = pyomo_nlp.evaluate_objective() + self.assertEqual(f, 900) + + gradf = pyomo_nlp.evaluate_grad_objective() + comparison_gradf = np.asarray([0, 0, 0, 0, 60], dtype=np.float64) + check_vectors_specific_order( + self, gradf, x_order, comparison_gradf, comparison_x_order + ) + c = pyomo_nlp.evaluate_constraints() + comparison_c = np.asarray([16, 6], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(2) + pyomo_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = pyomo_nlp.evaluate_jacobian() + comparison_j = np.asarray([[-1, 18, 24, 1, 0], [0, 18, 24, -1, 1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + pyomo_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = pyomo_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (5, 5)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, (4 * 3 * 21) + (4 * 3 * 5), (4 * 2 * 21) + (4 * 2 * 5), 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 2 * 1], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(AttributeError): + h = pyomo_nlp.evaluate_hessian_lag() + + +class TestExternalGreyBoxModelWithIncidenceAnalysis(unittest.TestCase): + """Tests for integration of ExternalGreyBoxBlock with incidence analysis""" + + def build_model(self): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputsWithHessian() + m.egb.set_external_model(external_model) + + return m + + def build_model_with_pyomo_components(self): + m = self.build_model() + + # Add Vars and linking constraints to m + m.Pin = pyo.Var() + m.c = pyo.Var() + m.F = pyo.Var() + m.P1 = pyo.Var() + m.P3 = pyo.Var() + m.P2 = pyo.Var() + m.Pout = pyo.Var() + + m.link_Pin = pyo.Constraint(expr=m.Pin == m.egb.inputs['Pin']) + m.link_c = pyo.Constraint(expr=m.c == m.egb.inputs['c']) + m.link_F = pyo.Constraint(expr=m.F == m.egb.inputs['F']) + m.link_P1 = pyo.Constraint(expr=m.P1 == m.egb.inputs['P1']) + m.link_P3 = pyo.Constraint(expr=m.P3 == m.egb.inputs['P3']) + m.link_P2 = pyo.Constraint(expr=m.P2 == m.egb.outputs['P2']) + m.link_Pout = pyo.Constraint(expr=m.Pout == m.egb.outputs['Pout']) + + return m + + def test_grey_box_only(self): + """ + Test that the incidence analysis correctly determines the DM partition for + a grey box model with two equality constraints and two outputs + """ + m = self.build_model() + + # Check that the get_incident_variables method on the implicit constraint body returns the correct variables + # Implicit constraint: 'pdrop1' + body_obj1 = m.egb.eq_constraints['pdrop1'].body + incident_vars1 = body_obj1.get_incident_variables() + + incident_var_set = ComponentSet( + [ + m.egb.inputs['Pin'], + m.egb.inputs['c'], + m.egb.inputs['F'], + m.egb.inputs['P1'], + ] + ) + self.assertEqual(ComponentSet(incident_vars1), incident_var_set) + + # Implicit constraint: 'pdrop3' + body_obj1 = m.egb.eq_constraints['pdrop3'].body + incident_vars1 = body_obj1.get_incident_variables() + incident_var_set = ComponentSet( + [ + m.egb.inputs['c'], + m.egb.inputs['F'], + m.egb.inputs['P1'], + m.egb.inputs['P3'], + ] + ) + self.assertEqual(ComponentSet(incident_vars1), incident_var_set) + + # Implicit constraint: 'P2_constraint' + body_obj1 = m.egb.output_constraints['P2'].body + incident_vars1 = body_obj1.get_incident_variables() + incident_var_set = ComponentSet( + [ + m.egb.inputs['c'], + m.egb.inputs['F'], + m.egb.inputs['P1'], + m.egb.outputs['P2'], + ] + ) + self.assertEqual(ComponentSet(incident_vars1), incident_var_set) + + # Implicit constraint: 'Pout_constraint' + body_obj1 = m.egb.output_constraints['Pout'].body + incident_vars1 = body_obj1.get_incident_variables() + incident_var_set = ComponentSet( + [ + m.egb.inputs['Pin'], + m.egb.inputs['c'], + m.egb.inputs['F'], + m.egb.outputs['Pout'], + ] + ) + self.assertEqual(ComponentSet(incident_vars1), incident_var_set) + + # Check Dulmage-Mendelsohn partitioning of the incidence graph + igraph = IncidenceGraphInterface(m, include_inequality=False) + var_dm_partition, con_dm_partition = igraph.dulmage_mendelsohn() + + # In this case, as we have not fixed any variables, we expect the system to be under-constrained. + # All variables will be in the under-constrained or unmatched sets + # All constraints will be in the under-constrained set + self.assertEqual(var_dm_partition.overconstrained, []) + self.assertEqual(var_dm_partition.square, []) + self.assertEqual(con_dm_partition.unmatched, []) + self.assertEqual(con_dm_partition.overconstrained, []) + self.assertEqual(con_dm_partition.square, []) + + self.assertEqual(len(var_dm_partition.underconstrained), 4) + self.assertEqual(len(var_dm_partition.unmatched), 3) + + for v in var_dm_partition.underconstrained: + # output variables should be in the under-constrained set + # The other two variables should be drawn from the inputs, but we cannot guarantee which ones + self.assertIn( + v.name, + [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + ], + ) + + for v in var_dm_partition.unmatched: + # Unmatched set will have the remaining 3 input variables, but again we cannot guarantee which ones + # We will instead check that the name is one of the inputs and that it is not in the under-constrained set + self.assertNotIn( + v.name, [u.name for u in var_dm_partition.underconstrained] + ) + self.assertIn( + v.name, + [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + ], + ) + + self.assertEqual(len(con_dm_partition.underconstrained), 4) + con_names = [c.name for c in con_dm_partition.underconstrained] + for c in con_names: + self.assertIn( + c, + [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + ], + ) + + def test_grey_box_w_pyomo_components(self): + """ + Test that the incidence analysis correctly determines the DM partition for + a model containing both grey box and other components + """ + m = self.build_model_with_pyomo_components() + + # Check Dulmage-Mendelsohn partitioning of the incidence graph + igraph = IncidenceGraphInterface(m, include_inequality=False) + var_dm_partition, con_dm_partition = igraph.dulmage_mendelsohn() + + # In this case, as we have not fixed any variables, we expect the system to be under-constrained. + # All variables will be in the under-constrained or unmatched sets + # All constraints will be in the under-constrained set + self.assertEqual(var_dm_partition.overconstrained, []) + self.assertEqual(var_dm_partition.square, []) + self.assertEqual(con_dm_partition.unmatched, []) + self.assertEqual(con_dm_partition.overconstrained, []) + self.assertEqual(con_dm_partition.square, []) + + self.assertEqual(len(var_dm_partition.underconstrained), 11) + self.assertEqual(len(var_dm_partition.unmatched), 3) + var_names = [ + v.name + for v in var_dm_partition.underconstrained + var_dm_partition.unmatched + ] + for v in var_names: + self.assertIn( + v, + [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + 'Pin', + 'c', + 'F', + 'P1', + 'P3', + 'P2', + 'Pout', + ], + ) + + self.assertEqual(len(con_dm_partition.underconstrained), 11) + con_names = [c.name for c in con_dm_partition.underconstrained] + for c in con_names: + self.assertIn( + c, + [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + 'link_Pin', + 'link_c', + 'link_F', + 'link_P1', + 'link_P3', + 'link_P2', + 'link_Pout', + ], + ) + + def test_grey_box_w_pyomo_components_square(self): + """ + Test that the incidence analysis correctly determines the DM partition for + a model containing both grey box and other components + """ + m = self.build_model_with_pyomo_components() + + # Fix 3 inputs + # Note that we have 2 implicit constraints that cross-link inputs + m.Pin.fix(1) + m.c.fix(1) + m.F.fix(1) + + # Check Dulmage-Mendelsohn partitioning of the incidence graph + igraph = IncidenceGraphInterface(m, include_inequality=False) + var_dm_partition, con_dm_partition = igraph.dulmage_mendelsohn() + + # In this case, as we have fixed all input variables, we expect the system to be square. + # All variables and constraints will be in the square sets + self.assertEqual(var_dm_partition.unmatched, []) + self.assertEqual(var_dm_partition.overconstrained, []) + self.assertEqual(var_dm_partition.underconstrained, []) + self.assertEqual(con_dm_partition.unmatched, []) + self.assertEqual(con_dm_partition.overconstrained, []) + self.assertEqual(con_dm_partition.underconstrained, []) + + self.assertEqual(len(var_dm_partition.square), 11) + var_names = [v.name for v in var_dm_partition.square] + for v in var_names: + self.assertIn( + v, + [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + # These three re fixed, so do not appear by default + # 'Pin', + # 'c', + # 'F', + 'P1', + 'P3', + 'P2', + 'Pout', + ], + ) + + self.assertEqual(len(con_dm_partition.square), 11) + con_names = [c.name for c in con_dm_partition.square] + for c in con_names: + self.assertIn( + c, + [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + 'link_Pin', + 'link_c', + 'link_F', + 'link_P1', + 'link_P3', + 'link_P2', + 'link_Pout', + ], + ) + + +class MyGreyBox(ExternalGreyBoxModel): + def n_inputs(self): + return 4 + + def input_names(self): + return [f"x[{i}]" for i in range(1, self.n_inputs() + 1)] + + def n_equality_constraints(self): + return 3 + + def equality_constraint_names(self): + return [f"eq[{i}]" for i in range(1, self.n_equality_constraints() + 1)] + + def set_input_values(self, input_values): + if len(input_values) != self.n_inputs(): + raise ValueError( + f"Expected {self.n_inputs()} inputs, got {len(input_values)}." + ) + self._inputs = np.asarray(input_values, dtype=float) + + def evaluate_equality_constraints(self): + x1, x2, x3, x4 = self._inputs + return np.array( + [x2 * x3**1.5 * x4 - 5.0, x1 - x2 + x4 - 4.0, x1 * np.exp(-x3) - 1.0], + dtype=float, + ) + + def evaluate_jacobian_equality_constraints(self): + x1, x2, x3, x4 = self._inputs + + # Rows correspond to eq[1], eq[2], eq[3] + # Cols correspond to x[1], x[2], x[3], x[4] + rows = np.array([0, 0, 0, 1, 1, 1, 2, 2], dtype=int) + cols = np.array([1, 2, 3, 0, 1, 3, 0, 2], dtype=int) + data = np.array( + [ + x3**1.5 * x4, # d(eq1)/d(x2) + 1.5 * x2 * x3**0.5 * x4, # d(eq1)/d(x3) + x2 * x3**1.5, # d(eq1)/d(x4) + 1.0, # d(eq2)/d(x1) + -1.0, # d(eq2)/d(x2) + 1.0, # d(eq2)/d(x4) + np.exp(-x3), # d(eq3)/d(x1) + -x1 * np.exp(-x3), # d(eq3)/d(x3) + ], + dtype=float, + ) + return coo_matrix( + (data, (rows, cols)), shape=(self.n_equality_constraints(), self.n_inputs()) + ) + + def finalize_block_construction(self, block): + self._inputs = np.full(self.n_inputs(), 1.0, dtype=float) + for v in block.inputs.values(): + v.set_value(1.0) + + +def test_with_custom_input_names(): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3, 4], bounds=(0.0, None), initialize=1.0) + m.objective = pyo.Objective( + expr=m.x[1] ** 2 + 2 * m.x[2] ** 2 + 3 * m.x[3] ** 2 + 4 * m.x[4] ** 2 + ) + m.grey_box = ExternalGreyBoxBlock() + m.grey_box.set_external_model(MyGreyBox(), inputs=[m.x[i] for i in range(1, 5)]) + + igraph = IncidenceGraphInterface(m) + matching = igraph.maximum_matching() + + # Minimal check on results, as we really only want to make sure the code runs + # when given custom input names + assert len(matching) == 3 + + +def test_custom_input_and_output_names(): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3, 4, 5], bounds=(0.0, None), initialize=1.0) + m.y = pyo.Var([1, 2], bounds=(0.0, None), initialize=1.0) + + m.egb = ExternalGreyBoxBlock() + external_model = ex_models.PressureDropTwoEqualitiesTwoOutputsWithHessian() + m.egb.set_external_model( + external_model, + inputs=[m.x[i] for i in range(1, 6)], + outputs=[m.y[i] for i in range(1, 3)], + ) + + igraph = IncidenceGraphInterface(m) + matching = igraph.maximum_matching() + + # Minimal check on results, as we really only want to make sure the code runs + # when given custom input names + assert len(matching) == 4 + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_external_pyomo_block.py b/pyomo/contrib/pynumero/interfaces/tests/test_external_pyomo_block.py index e49346a338b..729e89f6751 100644 --- a/pyomo/contrib/pynumero/interfaces/tests/test_external_pyomo_block.py +++ b/pyomo/contrib/pynumero/interfaces/tests/test_external_pyomo_block.py @@ -792,8 +792,8 @@ def test_set_and_evaluate(self): "linking_constraint[0]", "linking_constraint[1]", "linking_constraint[2]", - "ex_block.residual_0", - "ex_block.residual_1", + "ex_block.eq_constraints[residual_0]", + "ex_block.eq_constraints[residual_1]", ] self.assertEqual(constraint_names, nlp.constraint_names()) residuals = np.array( diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_pyomo_grey_box_nlp.py b/pyomo/contrib/pynumero/interfaces/tests/test_pyomo_grey_box_nlp.py index 5dab96d3489..c2fd8f5ad77 100644 --- a/pyomo/contrib/pynumero/interfaces/tests/test_pyomo_grey_box_nlp.py +++ b/pyomo/contrib/pynumero/interfaces/tests/test_pyomo_grey_box_nlp.py @@ -235,7 +235,7 @@ def _test_pressure_drop_single_equality(self, ex_model, hessian_support): 'egb.inputs[Pout]', ] x_order = egb_nlp.primals_names() - comparison_c_order = ['egb.pdrop'] + comparison_c_order = ['egb.eq_constraints[pdrop]'] c_order = egb_nlp.constraint_names() xlb = egb_nlp.primals_lb() @@ -547,7 +547,10 @@ def _test_pressure_drop_two_equalities(self, ex_model, hessian_support): 'egb.inputs[Pout]', ] x_order = egb_nlp.primals_names() - comparison_c_order = ['egb.pdrop2', 'egb.pdropout'] + comparison_c_order = [ + 'egb.eq_constraints[pdrop2]', + 'egb.eq_constraints[pdropout]', + ] c_order = egb_nlp.constraint_names() xlb = egb_nlp.primals_lb() @@ -713,8 +716,8 @@ def _test_pressure_drop_two_equalities_two_outputs(self, ex_model, hessian_suppo ] x_order = egb_nlp.primals_names() comparison_c_order = [ - 'egb.pdrop1', - 'egb.pdrop3', + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]', ] @@ -934,8 +937,8 @@ def test_scaling_greybox_only(self): egb_nlp = _ExternalGreyBoxAsNLP(m.egb) comparison_c_order = [ - 'egb.pdrop1', - 'egb.pdrop3', + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]', ] @@ -1204,7 +1207,7 @@ def _test_pressure_drop_single_equality(self, ex_model, hessian_support): 'egb.inputs[Pout]', ] x_order = pyomo_nlp.primals_names() - comparison_c_order = ['egb.pdrop'] + comparison_c_order = ['egb.eq_constraints[pdrop]'] c_order = pyomo_nlp.constraint_names() xlb = pyomo_nlp.primals_lb() @@ -1544,7 +1547,10 @@ def _test_pressure_drop_two_equalities(self, ex_model, hessian_support): 'egb.inputs[Pout]', ] x_order = pyomo_nlp.primals_names() - comparison_c_order = ['egb.pdrop2', 'egb.pdropout'] + comparison_c_order = [ + 'egb.eq_constraints[pdrop2]', + 'egb.eq_constraints[pdropout]', + ] c_order = pyomo_nlp.constraint_names() xlb = pyomo_nlp.primals_lb() @@ -1721,8 +1727,8 @@ def _test_pressure_drop_two_equalities_two_outputs(self, ex_model, hessian_suppo ] x_order = pyomo_nlp.primals_names() comparison_c_order = [ - 'egb.pdrop1', - 'egb.pdrop3', + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]', ] @@ -1935,8 +1941,8 @@ def _test_external_additional_constraints_vars(self, ex_model, hessian_support): ] x_order = pyomo_nlp.primals_names() comparison_c_order = [ - 'egb.pdrop1', - 'egb.pdrop3', + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]', 'incon', @@ -2242,8 +2248,8 @@ def test_scaling_pyomo_model_only(self): ] x_order = pyomo_nlp.primals_names() comparison_c_order = [ - 'egb.pdrop1', - 'egb.pdrop3', + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]', 'incon', @@ -2288,8 +2294,8 @@ def test_scaling_greybox_only(self): ] x_order = pyomo_nlp.primals_names() comparison_c_order = [ - 'egb.pdrop1', - 'egb.pdrop3', + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]', 'incon', @@ -2367,8 +2373,8 @@ def test_scaling_pyomo_model_and_greybox(self): ] x_order = pyomo_nlp.primals_names() comparison_c_order = [ - 'egb.pdrop1', - 'egb.pdrop3', + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]', 'incon', @@ -2465,7 +2471,7 @@ def test_external_greybox_solve_scaling(self): self.assertIn('c scaling provided', solver_trace) self.assertIn('d scaling provided', solver_trace) # x_order: ['egb.inputs[F]', 'egb.inputs[P1]', 'egb.inputs[P3]', 'egb.inputs[Pin]', 'egb.inputs[c]', 'egb.outputs[P2]', 'egb.outputs[Pout]', 'mu'] - # c_order: ['ccon', 'pcon', 'pincon', 'egb.pdrop1', 'egb.pdrop3', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]'] + # c_order: ['ccon', 'pcon', 'pincon', 'egb.eq_constraints[pdrop1]', 'egb.eq_constraints[pdrop3]', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]'] self.assertIn('DenseVector "x scaling vector" with 8 elements:', solver_trace) self.assertIn( 'x scaling vector[ 1]= 1.3000000000000000e+00', solver_trace @@ -2584,7 +2590,9 @@ def test_duals_after_solve(self): self.assertAlmostEqual( m.dual[m.egb]['egb.output_constraints[o]'], -25.0, places=3 ) - self.assertAlmostEqual(m.dual[m.egb]['egb.u2_con'], 62.5, places=3) + self.assertAlmostEqual( + m.dual[m.egb]['egb.eq_constraints[u2_con]'], 62.5, places=3 + ) def test_has_hessian_support_false(self): external_model = ex_models.PressureDropSingleOutput() diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_pyomo_grey_box_nlp_with_constraints.py b/pyomo/contrib/pynumero/interfaces/tests/test_pyomo_grey_box_nlp_with_constraints.py new file mode 100644 index 00000000000..2dff96bf4ee --- /dev/null +++ b/pyomo/contrib/pynumero/interfaces/tests/test_pyomo_grey_box_nlp_with_constraints.py @@ -0,0 +1,2752 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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.tempfiles import TempfileManager + +from pyomo.contrib.pynumero.dependencies import ( + numpy as np, + numpy_available, + scipy, + scipy_available, +) +from pyomo.common.dependencies.scipy import sparse as spa + +if not (numpy_available and scipy_available): + raise unittest.SkipTest("Pynumero needs scipy and numpy to run NLP tests") + +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.algorithms.solvers.cyipopt_solver import cyipopt_available + +from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock +from pyomo.contrib.pynumero.interfaces.pyomo_grey_box_nlp import ( + _ExternalGreyBoxAsNLP, + PyomoNLPWithGreyBoxBlocks, +) +from pyomo.contrib.pynumero.interfaces.tests.compare_utils import ( + check_vectors_specific_order, + check_sparse_matrix_specific_order, +) +import pyomo.contrib.pynumero.interfaces.tests.external_grey_box_models as ex_models +from pyomo.contrib.pynumero.examples.external_grey_box.external_with_objective import ( + solve_unconstrained, + solve_constrained, + solve_constrained_with_hessian, +) + + +class TestExternalGreyBoxAsNLP(unittest.TestCase): + def test_pressure_drop_single_output(self): + self._test_pressure_drop_single_output( + ex_models.PressureDropSingleOutput(), False + ) + self._test_pressure_drop_single_output( + ex_models.PressureDropSingleOutputWithHessian(), True + ) + + def _test_pressure_drop_single_output(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(0) + m.egb.outputs['Pout'].setub(100) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + + self.assertEqual(4, egb_nlp.n_primals()) + self.assertEqual(1, egb_nlp.n_constraints()) + self.assertEqual(4, egb_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(3, egb_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.outputs[Pout]', + ] + x_order = egb_nlp.primals_names() + comparison_c_order = ['egb.output_constraints[Pout]'] + c_order = egb_nlp.constraint_names() + + xlb = egb_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = egb_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = egb_nlp.constraints_lb() + comparison_clb = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = egb_nlp.constraints_ub() + comparison_cub = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = egb_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = egb_nlp.init_duals() + comparison_duals_init = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(4, len(egb_nlp.create_new_vector('primals'))) + self.assertEqual(1, len(egb_nlp.create_new_vector('constraints'))) + self.assertEqual(1, len(egb_nlp.create_new_vector('duals'))) + + egb_nlp.set_primals(np.asarray([1, 2, 3, 4], dtype=np.float64)) + x = egb_nlp.get_primals() + self.assertTrue(np.array_equal(x, np.asarray([1, 2, 3, 4], dtype=np.float64))) + egb_nlp.set_primals(egb_nlp.init_primals()) + + egb_nlp.set_duals(np.asarray([42], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42], dtype=np.float64))) + egb_nlp.set_duals(np.asarray([21], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21], dtype=np.float64))) + + c = egb_nlp.evaluate_constraints() + comparison_c = np.asarray([-22], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(1, dtype=np.float64) + egb_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = egb_nlp.evaluate_jacobian() + comparison_j = np.asarray([[1, -36, -48, -1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + egb_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = egb_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (4, 4)) + # hessian should be "full", not lower or upper triangular + comparison_h = np.asarray( + [ + [0, 0, 0, 0], + [0, 0, -8 * 3 * 21, 0], + [0, -8 * 3 * 21, -8 * 2 * 21, 0], + [0, 0, 0, 0], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = egb_nlp.evaluate_hessian_lag() + + def test_pressure_drop_single_equality(self): + self._test_pressure_drop_single_equality( + ex_models.PressureDropSingleEquality(), False + ) + self._test_pressure_drop_single_equality( + ex_models.PressureDropSingleEqualityWithHessian(), True + ) + + def _test_pressure_drop_single_equality(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['Pout'].value = 50 + m.egb.inputs['Pout'].setlb(0) + m.egb.inputs['Pout'].setub(100) + m.obj = pyo.Objective(expr=(m.egb.inputs['Pout'] - 20) ** 2) + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + + self.assertEqual(4, egb_nlp.n_primals()) + self.assertEqual(1, egb_nlp.n_constraints()) + self.assertEqual(4, egb_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(3, egb_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[Pout]', + ] + x_order = egb_nlp.primals_names() + comparison_c_order = ['egb.eq_constraints[pdrop]'] + c_order = egb_nlp.constraint_names() + + xlb = egb_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = egb_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = egb_nlp.constraints_lb() + comparison_clb = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = egb_nlp.constraints_ub() + comparison_cub = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = egb_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = egb_nlp.init_duals() + comparison_duals_init = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(4, len(egb_nlp.create_new_vector('primals'))) + self.assertEqual(1, len(egb_nlp.create_new_vector('constraints'))) + self.assertEqual(1, len(egb_nlp.create_new_vector('duals'))) + + egb_nlp.set_primals(np.asarray([1, 2, 3, 4], dtype=np.float64)) + x = egb_nlp.get_primals() + self.assertTrue(np.array_equal(x, np.asarray([1, 2, 3, 4], dtype=np.float64))) + egb_nlp.set_primals(egb_nlp.init_primals()) + + egb_nlp.set_duals(np.asarray([42], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42], dtype=np.float64))) + egb_nlp.set_duals(np.asarray([21], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21], dtype=np.float64))) + + c = egb_nlp.evaluate_constraints() + comparison_c = np.asarray([22], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(1, dtype=np.float64) + egb_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = egb_nlp.evaluate_jacobian() + comparison_j = np.asarray([[-1, 36, 48, 1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + egb_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = egb_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (4, 4)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0], + [0, 0, 8 * 3 * 21, 0], + [0, 8 * 3 * 21, 8 * 2 * 21, 0], + [0, 0, 0, 0], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = egb_nlp.evaluate_hessian_lag() + + def test_pressure_drop_two_outputs(self): + self._test_pressure_drop_two_outputs(ex_models.PressureDropTwoOutputs(), False) + self._test_pressure_drop_two_outputs( + ex_models.PressureDropTwoOutputsWithHessian(), True + ) + + def _test_pressure_drop_two_outputs(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.outputs['P2'].value = 80 + m.egb.outputs['P2'].setlb(10) + m.egb.outputs['P2'].setub(90) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(0) + m.egb.outputs['Pout'].setub(100) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + + self.assertEqual(5, egb_nlp.n_primals()) + self.assertEqual(2, egb_nlp.n_constraints()) + self.assertEqual(8, egb_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(3, egb_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + ] + x_order = egb_nlp.primals_names() + comparison_c_order = [ + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + ] + c_order = egb_nlp.constraint_names() + + xlb = egb_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 10, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = egb_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 90, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = egb_nlp.constraints_lb() + comparison_clb = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = egb_nlp.constraints_ub() + comparison_cub = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = egb_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 80, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = egb_nlp.init_duals() + comparison_duals_init = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(5, len(egb_nlp.create_new_vector('primals'))) + self.assertEqual(2, len(egb_nlp.create_new_vector('constraints'))) + self.assertEqual(2, len(egb_nlp.create_new_vector('duals'))) + + egb_nlp.set_primals(np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + x = egb_nlp.get_primals() + self.assertTrue( + np.array_equal(x, np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + ) + egb_nlp.set_primals(egb_nlp.init_primals()) + + egb_nlp.set_duals(np.asarray([42, 10], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42, 10], dtype=np.float64))) + egb_nlp.set_duals(np.asarray([21, 5], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21, 5], dtype=np.float64))) + + c = egb_nlp.evaluate_constraints() + comparison_c = np.asarray([-16, -22], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(2) + egb_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = egb_nlp.evaluate_jacobian() + comparison_j = np.asarray([[1, -18, -24, -1, 0], [1, -36, -48, 0, -1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + egb_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = egb_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (5, 5)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0, 0], + [0, 0, (-4 * 3 * 21) + (-8 * 3 * 5), 0, 0], + [ + 0, + (-4 * 3 * 21) + (-8 * 3 * 5), + (-4 * 2 * 21) + (-8 * 2 * 5), + 0, + 0, + ], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = egb_nlp.evaluate_hessian_lag() + + def test_pressure_drop_two_equalities(self): + self._test_pressure_drop_two_equalities( + ex_models.PressureDropTwoEqualities(), False + ) + self._test_pressure_drop_two_equalities( + ex_models.PressureDropTwoEqualitiesWithHessian(), True + ) + + def _test_pressure_drop_two_equalities(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P2'].value = 80 + m.egb.inputs['P2'].setlb(10) + m.egb.inputs['P2'].setub(90) + m.egb.inputs['Pout'].value = 50 + m.egb.inputs['Pout'].setlb(0) + m.egb.inputs['Pout'].setub(100) + m.obj = pyo.Objective(expr=(m.egb.inputs['Pout'] - 20) ** 2) + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + + self.assertEqual(5, egb_nlp.n_primals()) + self.assertEqual(2, egb_nlp.n_constraints()) + self.assertEqual(8, egb_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(3, egb_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P2]', + 'egb.inputs[Pout]', + ] + x_order = egb_nlp.primals_names() + comparison_c_order = [ + 'egb.eq_constraints[pdrop2]', + 'egb.eq_constraints[pdropout]', + ] + c_order = egb_nlp.constraint_names() + + xlb = egb_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 10, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = egb_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 90, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = egb_nlp.constraints_lb() + comparison_clb = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = egb_nlp.constraints_ub() + comparison_cub = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = egb_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 80, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = egb_nlp.init_duals() + comparison_duals_init = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(5, len(egb_nlp.create_new_vector('primals'))) + self.assertEqual(2, len(egb_nlp.create_new_vector('constraints'))) + self.assertEqual(2, len(egb_nlp.create_new_vector('duals'))) + + egb_nlp.set_primals(np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + x = egb_nlp.get_primals() + self.assertTrue( + np.array_equal(x, np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + ) + egb_nlp.set_primals(egb_nlp.init_primals()) + + egb_nlp.set_duals(np.asarray([42, 10], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42, 10], dtype=np.float64))) + egb_nlp.set_duals(np.asarray([21, 5], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21, 5], dtype=np.float64))) + + c = egb_nlp.evaluate_constraints() + comparison_c = np.asarray([16, 6], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(2) + egb_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = egb_nlp.evaluate_jacobian() + comparison_j = np.asarray([[-1, 18, 24, 1, 0], [0, 18, 24, -1, 1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + egb_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = egb_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (5, 5)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0, 0], + [0, 0, (4 * 3 * 21) + (4 * 3 * 5), 0, 0], + [0, (4 * 3 * 21) + (4 * 3 * 5), (4 * 2 * 21) + (4 * 2 * 5), 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = egb_nlp.evaluate_hessian_lag() + + def test_pressure_drop_two_equalities_two_outputs(self): + self._test_pressure_drop_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputs(), False + ) + self._test_pressure_drop_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputsWithHessian(), True + ) + + def _test_pressure_drop_two_equalities_two_outputs(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P1'].value = 80 + m.egb.inputs['P1'].setlb(10) + m.egb.inputs['P1'].setub(90) + m.egb.inputs['P3'].value = 70 + m.egb.inputs['P3'].setlb(20) + m.egb.inputs['P3'].setub(80) + m.egb.outputs['P2'].value = 75 + m.egb.outputs['P2'].setlb(15) + m.egb.outputs['P2'].setub(85) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(30) + m.egb.outputs['Pout'].setub(70) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + + self.assertEqual(7, egb_nlp.n_primals()) + self.assertEqual(4, egb_nlp.n_constraints()) + self.assertEqual(16, egb_nlp.nnz_jacobian()) + if hessian_support: + # this number is larger than expected because the nnz for the + # hessian of equality and output constraints are concatenated + # even if they occur in the same place + self.assertEqual(6, egb_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + ] + x_order = egb_nlp.primals_names() + comparison_c_order = [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + ] + c_order = egb_nlp.constraint_names() + + xlb = egb_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 10, 20, 15, 30], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = egb_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 90, 80, 85, 70], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = egb_nlp.constraints_lb() + comparison_clb = np.asarray([0, 0, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = egb_nlp.constraints_ub() + comparison_cub = np.asarray([0, 0, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = egb_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 80, 70, 75, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = egb_nlp.init_duals() + comparison_duals_init = np.asarray([0, 0, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(7, len(egb_nlp.create_new_vector('primals'))) + self.assertEqual(4, len(egb_nlp.create_new_vector('constraints'))) + self.assertEqual(4, len(egb_nlp.create_new_vector('duals'))) + + egb_nlp.set_primals(np.asarray([1, 2, 3, 4, 5, 6, 7], dtype=np.float64)) + x = egb_nlp.get_primals() + self.assertTrue( + np.array_equal(x, np.asarray([1, 2, 3, 4, 5, 6, 7], dtype=np.float64)) + ) + egb_nlp.set_primals(egb_nlp.init_primals()) + + egb_nlp.set_duals(np.asarray([42, 10, 11, 12], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue( + np.array_equal(y, np.asarray([42, 10, 11, 12], dtype=np.float64)) + ) + egb_nlp.set_duals(np.asarray([21, 5, 6, 7], dtype=np.float64)) + y = egb_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21, 5, 6, 7], dtype=np.float64))) + + c = egb_nlp.evaluate_constraints() + comparison_c = np.asarray([-2, 26, -13, -22], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(4) + egb_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = egb_nlp.evaluate_jacobian() + comparison_j = np.asarray( + [ + [-1, 9, 12, 1, 0, 0, 0], + [0, 18, 24, -1, 1, 0, 0], + [0, -9, -12, 1, 0, -1, 0], + [1, -36, -48, 0, 0, 0, -1], + ] + ) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + egb_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = egb_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (7, 7)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0], + [ + 0, + 0, + (2 * 3 * 21) + (4 * 3 * 5) + (-2 * 3 * 6) + (-8 * 3 * 7), + 0, + 0, + 0, + 0, + ], + [ + 0, + (2 * 3 * 21) + (4 * 3 * 5) + (-2 * 3 * 6) + (-8 * 3 * 7), + (2 * 2 * 21) + (4 * 2 * 5) + (-2 * 2 * 6) + (-8 * 2 * 7), + 0, + 0, + 0, + 0, + ], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = egb_nlp.evaluate_hessian_lag() + + def create_model_two_equalities_two_outputs(self, external_model): + m = pyo.ConcreteModel() + m.hin = pyo.Var(bounds=(0, None), initialize=10) + m.hout = pyo.Var(bounds=(0, None)) + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(external_model) + m.incon = pyo.Constraint(expr=0 <= m.egb.inputs['Pin'] - 10 * m.hin) + m.outcon = pyo.Constraint(expr=0 == m.egb.outputs['Pout'] - 10 * m.hout) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P1'].value = 80 + m.egb.inputs['P1'].setlb(10) + m.egb.inputs['P1'].setub(90) + m.egb.inputs['P3'].value = 70 + m.egb.inputs['P3'].setlb(20) + m.egb.inputs['P3'].setub(80) + m.egb.outputs['P2'].value = 75 + m.egb.outputs['P2'].setlb(15) + m.egb.outputs['P2'].setub(85) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(30) + m.egb.outputs['Pout'].setub(70) + return m + + def test_scaling_all_missing(self): + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputs() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + with self.assertRaises(NotImplementedError): + fs = egb_nlp.get_obj_scaling() + with self.assertRaises(NotImplementedError): + xs = egb_nlp.get_primals_scaling() + cs = egb_nlp.get_constraints_scaling() + self.assertIsNone(cs) + + def test_scaling_pyomo_model_only(self): + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputs() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + # m.scaling_factor[m.obj] = 0.1 # scale the objective + m.scaling_factor[m.egb.inputs['Pin']] = 1.1 # scale the variable + m.scaling_factor[m.egb.inputs['c']] = 1.2 # scale the variable + m.scaling_factor[m.egb.inputs['F']] = 1.3 # scale the variable + # m.scaling_factor[m.egb.inputs['P1']] = 1.4 # scale the variable + m.scaling_factor[m.egb.inputs['P3']] = 1.5 # scale the variable + m.scaling_factor[m.egb.outputs['P2']] = 1.6 # scale the variable + m.scaling_factor[m.egb.outputs['Pout']] = 1.7 # scale the variable + # m.scaling_factor[m.hin] = 1.8 + m.scaling_factor[m.hout] = 1.9 + # m.scaling_factor[m.incon] = 2.1 + m.scaling_factor[m.outcon] = 2.2 + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + + with self.assertRaises(NotImplementedError): + fs = egb_nlp.get_obj_scaling() + with self.assertRaises(NotImplementedError): + xs = egb_nlp.get_primals_scaling() + + cs = egb_nlp.get_constraints_scaling() + self.assertIsNone(cs) + + def test_scaling_greybox_only(self): + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputsScaleBoth() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + + comparison_c_order = [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + ] + c_order = egb_nlp.constraint_names() + + with self.assertRaises(NotImplementedError): + fs = egb_nlp.get_obj_scaling() + with self.assertRaises(NotImplementedError): + xs = egb_nlp.get_primals_scaling() + + cs = egb_nlp.get_constraints_scaling() + comparison_cs = np.asarray([3.1, 3.2, 4.1, 4.2], dtype=np.float64) + check_vectors_specific_order( + self, cs, c_order, comparison_cs, comparison_c_order + ) + + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputsScaleEqualities() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + cs = egb_nlp.get_constraints_scaling() + comparison_cs = np.asarray([3.1, 3.2, 1, 1], dtype=np.float64) + check_vectors_specific_order( + self, cs, c_order, comparison_cs, comparison_c_order + ) + + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputsScaleOutputs() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + egb_nlp = _ExternalGreyBoxAsNLP(m.egb) + cs = egb_nlp.get_constraints_scaling() + comparison_cs = np.asarray([1, 1, 4.1, 4.2], dtype=np.float64) + check_vectors_specific_order( + self, cs, c_order, comparison_cs, comparison_c_order + ) + + +class TestPyomoNLPWithGreyBoxModels(unittest.TestCase): + def test_error_no_variables(self): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.PressureDropSingleOutput()) + m.obj = pyo.Objective(expr=1) + with self.assertRaises(ValueError): + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + def test_error_fixed_inputs_outputs(self): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.PressureDropSingleOutput()) + m.egb.inputs['Pin'].fix(100) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + with self.assertRaises(NotImplementedError): + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.PressureDropTwoOutputs()) + m.egb.outputs['P2'].fix(50) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + with self.assertRaises(NotImplementedError): + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + def test_pressure_drop_single_output(self): + self._test_pressure_drop_single_output( + ex_models.PressureDropSingleOutput(), False + ) + self._test_pressure_drop_single_output( + ex_models.PressureDropSingleOutputWithHessian(), True + ) + + def _test_pressure_drop_single_output(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(0) + m.egb.outputs['Pout'].setub(100) + # m.dummy = pyo.Constraint(expr=sum(m.egb.inputs[i] for i in m.egb.inputs) + sum(m.egb.outputs[i] for i in m.egb.outputs) <= 1e6) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + self.assertEqual(4, pyomo_nlp.n_primals()) + self.assertEqual(1, pyomo_nlp.n_constraints()) + self.assertEqual(4, pyomo_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(4, pyomo_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.outputs[Pout]', + ] + x_order = pyomo_nlp.primals_names() + comparison_c_order = ['egb.output_constraints[Pout]'] + c_order = pyomo_nlp.constraint_names() + + xlb = pyomo_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = pyomo_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = pyomo_nlp.constraints_lb() + comparison_clb = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = pyomo_nlp.constraints_ub() + comparison_cub = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = pyomo_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = pyomo_nlp.init_duals() + comparison_duals_init = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(4, len(pyomo_nlp.create_new_vector('primals'))) + self.assertEqual(1, len(pyomo_nlp.create_new_vector('constraints'))) + self.assertEqual(1, len(pyomo_nlp.create_new_vector('duals'))) + + pyomo_nlp.set_primals(np.asarray([1, 2, 3, 4], dtype=np.float64)) + x = pyomo_nlp.get_primals() + self.assertTrue(np.array_equal(x, np.asarray([1, 2, 3, 4], dtype=np.float64))) + pyomo_nlp.set_primals(pyomo_nlp.init_primals()) + + pyomo_nlp.set_duals(np.asarray([42], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42], dtype=np.float64))) + pyomo_nlp.set_duals(np.asarray([21], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21], dtype=np.float64))) + + fac = pyomo_nlp.get_obj_factor() + self.assertEqual(fac, 1) + pyomo_nlp.set_obj_factor(42) + self.assertEqual(pyomo_nlp.get_obj_factor(), 42) + pyomo_nlp.set_obj_factor(1) + + f = pyomo_nlp.evaluate_objective() + self.assertEqual(f, 900) + + gradf = pyomo_nlp.evaluate_grad_objective() + comparison_gradf = np.asarray([0, 0, 0, 60], dtype=np.float64) + check_vectors_specific_order( + self, gradf, x_order, comparison_gradf, comparison_x_order + ) + c = pyomo_nlp.evaluate_constraints() + comparison_c = np.asarray([-22], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(1) + pyomo_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = pyomo_nlp.evaluate_jacobian() + comparison_j = np.asarray([[1, -36, -48, -1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + pyomo_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = pyomo_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (4, 4)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0], + [0, 0, -8 * 3 * 21, 0], + [0, -8 * 3 * 21, -8 * 2 * 21, 0], + [0, 0, 0, 2 * 1], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = pyomo_nlp.evaluate_hessian_lag() + + def test_pressure_drop_single_equality(self): + self._test_pressure_drop_single_equality( + ex_models.PressureDropSingleEquality(), False + ) + self._test_pressure_drop_single_equality( + ex_models.PressureDropSingleEqualityWithHessian(), True + ) + + def _test_pressure_drop_single_equality(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['Pout'].value = 50 + m.egb.inputs['Pout'].setlb(0) + m.egb.inputs['Pout'].setub(100) + m.obj = pyo.Objective(expr=(m.egb.inputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + self.assertEqual(4, pyomo_nlp.n_primals()) + self.assertEqual(1, pyomo_nlp.n_constraints()) + self.assertEqual(4, pyomo_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(4, pyomo_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[Pout]', + ] + x_order = pyomo_nlp.primals_names() + comparison_c_order = ['egb.eq_constraints[pdrop]'] + c_order = pyomo_nlp.constraint_names() + + xlb = pyomo_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = pyomo_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = pyomo_nlp.constraints_lb() + comparison_clb = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = pyomo_nlp.constraints_ub() + comparison_cub = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = pyomo_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = pyomo_nlp.init_duals() + comparison_duals_init = np.asarray([0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(4, len(pyomo_nlp.create_new_vector('primals'))) + self.assertEqual(1, len(pyomo_nlp.create_new_vector('constraints'))) + self.assertEqual(1, len(pyomo_nlp.create_new_vector('duals'))) + + pyomo_nlp.set_primals(np.asarray([1, 2, 3, 4], dtype=np.float64)) + x = pyomo_nlp.get_primals() + self.assertTrue(np.array_equal(x, np.asarray([1, 2, 3, 4], dtype=np.float64))) + pyomo_nlp.set_primals(pyomo_nlp.init_primals()) + + pyomo_nlp.set_duals(np.asarray([42], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42], dtype=np.float64))) + pyomo_nlp.set_duals(np.asarray([21], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21], dtype=np.float64))) + + fac = pyomo_nlp.get_obj_factor() + self.assertEqual(fac, 1) + pyomo_nlp.set_obj_factor(42) + self.assertEqual(pyomo_nlp.get_obj_factor(), 42) + pyomo_nlp.set_obj_factor(1) + + f = pyomo_nlp.evaluate_objective() + self.assertEqual(f, 900) + + gradf = pyomo_nlp.evaluate_grad_objective() + comparison_gradf = np.asarray([0, 0, 0, 60], dtype=np.float64) + check_vectors_specific_order( + self, gradf, x_order, comparison_gradf, comparison_x_order + ) + c = pyomo_nlp.evaluate_constraints() + comparison_c = np.asarray([22], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(1) + pyomo_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = pyomo_nlp.evaluate_jacobian() + comparison_j = np.asarray([[-1, 36, 48, 1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + pyomo_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = pyomo_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (4, 4)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0], + [0, 0, 8 * 3 * 21, 0], + [0, 8 * 3 * 21, 8 * 2 * 21, 0], + [0, 0, 0, 2 * 1], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = pyomo_nlp.evaluate_hessian_lag() + + def test_pressure_drop_two_outputs(self): + self._test_pressure_drop_two_outputs(ex_models.PressureDropTwoOutputs(), False) + self._test_pressure_drop_two_outputs( + ex_models.PressureDropTwoOutputsWithHessian(), True + ) + + def _test_pressure_drop_two_outputs(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.outputs['P2'].value = 80 + m.egb.outputs['P2'].setlb(10) + m.egb.outputs['P2'].setub(90) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(0) + m.egb.outputs['Pout'].setub(100) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + self.assertEqual(5, pyomo_nlp.n_primals()) + self.assertEqual(2, pyomo_nlp.n_constraints()) + self.assertEqual(8, pyomo_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(4, pyomo_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + ] + x_order = pyomo_nlp.primals_names() + comparison_c_order = [ + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + ] + c_order = pyomo_nlp.constraint_names() + + xlb = pyomo_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 10, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = pyomo_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 90, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = pyomo_nlp.constraints_lb() + comparison_clb = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = pyomo_nlp.constraints_ub() + comparison_cub = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = pyomo_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 80, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = pyomo_nlp.init_duals() + comparison_duals_init = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(5, len(pyomo_nlp.create_new_vector('primals'))) + self.assertEqual(2, len(pyomo_nlp.create_new_vector('constraints'))) + self.assertEqual(2, len(pyomo_nlp.create_new_vector('duals'))) + + pyomo_nlp.set_primals(np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + x = pyomo_nlp.get_primals() + self.assertTrue( + np.array_equal(x, np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + ) + pyomo_nlp.set_primals(pyomo_nlp.init_primals()) + + pyomo_nlp.set_duals(np.asarray([42, 10], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42, 10], dtype=np.float64))) + pyomo_nlp.set_duals(np.asarray([21, 5], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21, 5], dtype=np.float64))) + + fac = pyomo_nlp.get_obj_factor() + self.assertEqual(fac, 1) + pyomo_nlp.set_obj_factor(42) + self.assertEqual(pyomo_nlp.get_obj_factor(), 42) + pyomo_nlp.set_obj_factor(1) + + f = pyomo_nlp.evaluate_objective() + self.assertEqual(f, 900) + + gradf = pyomo_nlp.evaluate_grad_objective() + comparison_gradf = np.asarray([0, 0, 0, 0, 60], dtype=np.float64) + check_vectors_specific_order( + self, gradf, x_order, comparison_gradf, comparison_x_order + ) + c = pyomo_nlp.evaluate_constraints() + comparison_c = np.asarray([-16, -22], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(2) + pyomo_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = pyomo_nlp.evaluate_jacobian() + comparison_j = np.asarray([[1, -18, -24, -1, 0], [1, -36, -48, 0, -1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + pyomo_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = pyomo_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (5, 5)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0, 0], + [0, 0, (-4 * 3 * 21) + (-8 * 3 * 5), 0, 0], + [ + 0, + (-4 * 3 * 21) + (-8 * 3 * 5), + (-4 * 2 * 21) + (-8 * 2 * 5), + 0, + 0, + ], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 2 * 1], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = pyomo_nlp.evaluate_hessian_lag() + + def test_pressure_drop_two_equalities(self): + self._test_pressure_drop_two_equalities( + ex_models.PressureDropTwoEqualities(), False + ) + self._test_pressure_drop_two_equalities( + ex_models.PressureDropTwoEqualitiesWithHessian(), True + ) + + def _test_pressure_drop_two_equalities(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P2'].value = 80 + m.egb.inputs['P2'].setlb(10) + m.egb.inputs['P2'].setub(90) + m.egb.inputs['Pout'].value = 50 + m.egb.inputs['Pout'].setlb(0) + m.egb.inputs['Pout'].setub(100) + m.obj = pyo.Objective(expr=(m.egb.inputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + self.assertEqual(5, pyomo_nlp.n_primals()) + self.assertEqual(2, pyomo_nlp.n_constraints()) + self.assertEqual(8, pyomo_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(4, pyomo_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P2]', + 'egb.inputs[Pout]', + ] + x_order = pyomo_nlp.primals_names() + comparison_c_order = [ + 'egb.eq_constraints[pdrop2]', + 'egb.eq_constraints[pdropout]', + ] + c_order = pyomo_nlp.constraint_names() + + xlb = pyomo_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 10, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = pyomo_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 90, 100], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = pyomo_nlp.constraints_lb() + comparison_clb = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = pyomo_nlp.constraints_ub() + comparison_cub = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = pyomo_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 80, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = pyomo_nlp.init_duals() + comparison_duals_init = np.asarray([0, 0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(5, len(pyomo_nlp.create_new_vector('primals'))) + self.assertEqual(2, len(pyomo_nlp.create_new_vector('constraints'))) + self.assertEqual(2, len(pyomo_nlp.create_new_vector('duals'))) + + pyomo_nlp.set_primals(np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + x = pyomo_nlp.get_primals() + self.assertTrue( + np.array_equal(x, np.asarray([1, 2, 3, 4, 5], dtype=np.float64)) + ) + pyomo_nlp.set_primals(pyomo_nlp.init_primals()) + + pyomo_nlp.set_duals(np.asarray([42, 10], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([42, 10], dtype=np.float64))) + pyomo_nlp.set_duals(np.asarray([21, 5], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21, 5], dtype=np.float64))) + + fac = pyomo_nlp.get_obj_factor() + self.assertEqual(fac, 1) + pyomo_nlp.set_obj_factor(42) + self.assertEqual(pyomo_nlp.get_obj_factor(), 42) + pyomo_nlp.set_obj_factor(1) + + f = pyomo_nlp.evaluate_objective() + self.assertEqual(f, 900) + + gradf = pyomo_nlp.evaluate_grad_objective() + comparison_gradf = np.asarray([0, 0, 0, 0, 60], dtype=np.float64) + check_vectors_specific_order( + self, gradf, x_order, comparison_gradf, comparison_x_order + ) + c = pyomo_nlp.evaluate_constraints() + comparison_c = np.asarray([16, 6], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(2) + pyomo_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = pyomo_nlp.evaluate_jacobian() + comparison_j = np.asarray([[-1, 18, 24, 1, 0], [0, 18, 24, -1, 1]]) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + pyomo_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = pyomo_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (5, 5)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0, 0], + [0, 0, (4 * 3 * 21) + (4 * 3 * 5), 0, 0], + [0, (4 * 3 * 21) + (4 * 3 * 5), (4 * 2 * 21) + (4 * 2 * 5), 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 2 * 1], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = pyomo_nlp.evaluate_hessian_lag() + + def test_pressure_drop_two_equalities_two_outputs(self): + self._test_pressure_drop_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputs(), False + ) + self._test_pressure_drop_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputsWithHessian(), True + ) + + def _test_pressure_drop_two_equalities_two_outputs(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P1'].value = 80 + m.egb.inputs['P1'].setlb(10) + m.egb.inputs['P1'].setub(90) + m.egb.inputs['P3'].value = 70 + m.egb.inputs['P3'].setlb(20) + m.egb.inputs['P3'].setub(80) + m.egb.outputs['P2'].value = 75 + m.egb.outputs['P2'].setlb(15) + m.egb.outputs['P2'].setub(85) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(30) + m.egb.outputs['Pout'].setub(70) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + self.assertEqual(7, pyomo_nlp.n_primals()) + self.assertEqual(4, pyomo_nlp.n_constraints()) + self.assertEqual(16, pyomo_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(4, pyomo_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + ] + x_order = pyomo_nlp.primals_names() + comparison_c_order = [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + ] + c_order = pyomo_nlp.constraint_names() + + xlb = pyomo_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 10, 20, 15, 30], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = pyomo_nlp.primals_ub() + comparison_xub = np.asarray([150, 5, 5, 90, 80, 85, 70], dtype=np.float64) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = pyomo_nlp.constraints_lb() + comparison_clb = np.asarray([0, 0, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = pyomo_nlp.constraints_ub() + comparison_cub = np.asarray([0, 0, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = pyomo_nlp.init_primals() + comparison_xinit = np.asarray([100, 2, 3, 80, 70, 75, 50], dtype=np.float64) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = pyomo_nlp.init_duals() + comparison_duals_init = np.asarray([0, 0, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(7, len(pyomo_nlp.create_new_vector('primals'))) + self.assertEqual(4, len(pyomo_nlp.create_new_vector('constraints'))) + self.assertEqual(4, len(pyomo_nlp.create_new_vector('duals'))) + + pyomo_nlp.set_primals(np.asarray([1, 2, 3, 4, 5, 6, 7], dtype=np.float64)) + x = pyomo_nlp.get_primals() + self.assertTrue( + np.array_equal(x, np.asarray([1, 2, 3, 4, 5, 6, 7], dtype=np.float64)) + ) + pyomo_nlp.set_primals(pyomo_nlp.init_primals()) + + pyomo_nlp.set_duals(np.asarray([42, 10, 11, 12], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue( + np.array_equal(y, np.asarray([42, 10, 11, 12], dtype=np.float64)) + ) + pyomo_nlp.set_duals(np.asarray([21, 5, 6, 7], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue(np.array_equal(y, np.asarray([21, 5, 6, 7], dtype=np.float64))) + + fac = pyomo_nlp.get_obj_factor() + self.assertEqual(fac, 1) + pyomo_nlp.set_obj_factor(42) + self.assertEqual(pyomo_nlp.get_obj_factor(), 42) + pyomo_nlp.set_obj_factor(1) + + f = pyomo_nlp.evaluate_objective() + self.assertEqual(f, 900) + + gradf = pyomo_nlp.evaluate_grad_objective() + comparison_gradf = np.asarray([0, 0, 0, 0, 0, 0, 60], dtype=np.float64) + check_vectors_specific_order( + self, gradf, x_order, comparison_gradf, comparison_x_order + ) + c = pyomo_nlp.evaluate_constraints() + comparison_c = np.asarray([-2, 26, -13, -22], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(4) + pyomo_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = pyomo_nlp.evaluate_jacobian() + comparison_j = np.asarray( + [ + [-1, 9, 12, 1, 0, 0, 0], + [0, 18, 24, -1, 1, 0, 0], + [0, -9, -12, 1, 0, -1, 0], + [1, -36, -48, 0, 0, 0, -1], + ] + ) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + pyomo_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = pyomo_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (7, 7)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0], + [ + 0, + 0, + (2 * 3 * 21) + (4 * 3 * 5) + (-2 * 3 * 6) + (-8 * 3 * 7), + 0, + 0, + 0, + 0, + ], + [ + 0, + (2 * 3 * 21) + (4 * 3 * 5) + (-2 * 3 * 6) + (-8 * 3 * 7), + (2 * 2 * 21) + (4 * 2 * 5) + (-2 * 2 * 6) + (-8 * 2 * 7), + 0, + 0, + 0, + 0, + ], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 2 * 1], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = pyomo_nlp.evaluate_hessian_lag() + + def test_external_additional_constraints_vars(self): + self._test_external_additional_constraints_vars( + ex_models.PressureDropTwoEqualitiesTwoOutputs(), False + ) + self._test_external_additional_constraints_vars( + ex_models.PressureDropTwoEqualitiesTwoOutputsWithHessian(), True + ) + + def _test_external_additional_constraints_vars(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.hin = pyo.Var(bounds=(0, None), initialize=10) + m.hout = pyo.Var(bounds=(0, None)) + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.incon = pyo.Constraint(expr=0 <= m.egb.inputs['Pin'] - 10 * m.hin) + m.outcon = pyo.Constraint(expr=0 == m.egb.outputs['Pout'] - 10 * m.hout) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P1'].value = 80 + m.egb.inputs['P1'].setlb(10) + m.egb.inputs['P1'].setub(90) + m.egb.inputs['P3'].value = 70 + m.egb.inputs['P3'].setlb(20) + m.egb.inputs['P3'].setub(80) + m.egb.outputs['P2'].value = 75 + m.egb.outputs['P2'].setlb(15) + m.egb.outputs['P2'].setub(85) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(30) + m.egb.outputs['Pout'].setub(70) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + self.assertEqual(9, pyomo_nlp.n_primals()) + self.assertEqual(6, pyomo_nlp.n_constraints()) + self.assertEqual(20, pyomo_nlp.nnz_jacobian()) + if hessian_support: + self.assertEqual(4, pyomo_nlp.nnz_hessian_lag()) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + 'hin', + 'hout', + ] + x_order = pyomo_nlp.primals_names() + comparison_c_order = [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + 'incon', + 'outcon', + ] + c_order = pyomo_nlp.constraint_names() + + xlb = pyomo_nlp.primals_lb() + comparison_xlb = np.asarray([50, 1, 1, 10, 20, 15, 30, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, xlb, x_order, comparison_xlb, comparison_x_order + ) + xub = pyomo_nlp.primals_ub() + comparison_xub = np.asarray( + [150, 5, 5, 90, 80, 85, 70, np.inf, np.inf], dtype=np.float64 + ) + check_vectors_specific_order( + self, xub, x_order, comparison_xub, comparison_x_order + ) + clb = pyomo_nlp.constraints_lb() + comparison_clb = np.asarray([0, 0, 0, 0, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, clb, c_order, comparison_clb, comparison_c_order + ) + cub = pyomo_nlp.constraints_ub() + comparison_cub = np.asarray([0, 0, 0, 0, np.inf, 0], dtype=np.float64) + check_vectors_specific_order( + self, cub, c_order, comparison_cub, comparison_c_order + ) + + xinit = pyomo_nlp.init_primals() + comparison_xinit = np.asarray( + [100, 2, 3, 80, 70, 75, 50, 10, 0], dtype=np.float64 + ) + check_vectors_specific_order( + self, xinit, x_order, comparison_xinit, comparison_x_order + ) + duals_init = pyomo_nlp.init_duals() + comparison_duals_init = np.asarray([0, 0, 0, 0, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, duals_init, c_order, comparison_duals_init, comparison_c_order + ) + + self.assertEqual(9, len(pyomo_nlp.create_new_vector('primals'))) + self.assertEqual(6, len(pyomo_nlp.create_new_vector('constraints'))) + self.assertEqual(6, len(pyomo_nlp.create_new_vector('duals'))) + + pyomo_nlp.set_primals(np.asarray([1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=np.float64)) + x = pyomo_nlp.get_primals() + self.assertTrue( + np.array_equal(x, np.asarray([1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=np.float64)) + ) + pyomo_nlp.set_primals(pyomo_nlp.init_primals()) + + pyomo_nlp.set_duals(np.asarray([42, 10, 11, 12, 13, 14], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue( + np.array_equal(y, np.asarray([42, 10, 11, 12, 13, 14], dtype=np.float64)) + ) + pyomo_nlp.set_duals(np.asarray([0, 0, 21, 5, 6, 7], dtype=np.float64)) + y = pyomo_nlp.get_duals() + self.assertTrue( + np.array_equal(y, np.asarray([0, 0, 21, 5, 6, 7], dtype=np.float64)) + ) + + fac = pyomo_nlp.get_obj_factor() + self.assertEqual(fac, 1) + pyomo_nlp.set_obj_factor(42) + self.assertEqual(pyomo_nlp.get_obj_factor(), 42) + pyomo_nlp.set_obj_factor(1) + + f = pyomo_nlp.evaluate_objective() + self.assertEqual(f, 900) + + gradf = pyomo_nlp.evaluate_grad_objective() + comparison_gradf = np.asarray([0, 0, 0, 0, 0, 0, 60, 0, 0], dtype=np.float64) + check_vectors_specific_order( + self, gradf, x_order, comparison_gradf, comparison_x_order + ) + c = pyomo_nlp.evaluate_constraints() + comparison_c = np.asarray([-2, 26, -13, -22, 0, 50], dtype=np.float64) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + c = np.zeros(6) + pyomo_nlp.evaluate_constraints(out=c) + check_vectors_specific_order(self, c, c_order, comparison_c, comparison_c_order) + + j = pyomo_nlp.evaluate_jacobian() + comparison_j = np.asarray( + [ + [-1, 9, 12, 1, 0, 0, 0, 0, 0], + [0, 18, 24, -1, 1, 0, 0, 0, 0], + [0, -9, -12, 1, 0, -1, 0, 0, 0], + [1, -36, -48, 0, 0, 0, -1, 0, 0], + [1, 0, 0, 0, 0, 0, 0, -10, 0], + [0, 0, 0, 0, 0, 0, 1, 0, -10], + ] + ) + + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + j = 2.0 * j + pyomo_nlp.evaluate_jacobian(out=j) + check_sparse_matrix_specific_order( + self, + j, + c_order, + x_order, + comparison_j, + comparison_c_order, + comparison_x_order, + ) + + if hessian_support: + h = pyomo_nlp.evaluate_hessian_lag() + self.assertTrue(h.shape == (9, 9)) + comparison_h = np.asarray( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [ + 0, + 0, + (2 * 3 * 21) + (4 * 3 * 5) + (-2 * 3 * 6) + (-8 * 3 * 7), + 0, + 0, + 0, + 0, + 0, + 0, + ], + [ + 0, + (2 * 3 * 21) + (4 * 3 * 5) + (-2 * 3 * 6) + (-8 * 3 * 7), + (2 * 2 * 21) + (4 * 2 * 5) + (-2 * 2 * 6) + (-8 * 2 * 7), + 0, + 0, + 0, + 0, + 0, + 0, + ], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 2 * 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0], + ], + dtype=np.float64, + ) + check_sparse_matrix_specific_order( + self, + h, + x_order, + x_order, + comparison_h, + comparison_x_order, + comparison_x_order, + ) + else: + with self.assertRaises(NotImplementedError): + h = pyomo_nlp.evaluate_hessian_lag() + + @unittest.skipIf(not cyipopt_available, "CyIpopt needed to run tests with solve") + def test_external_greybox_solve(self): + self._test_external_greybox_solve( + ex_models.PressureDropTwoEqualitiesTwoOutputs(), False + ) + self._test_external_greybox_solve( + ex_models.PressureDropTwoEqualitiesTwoOutputsWithHessian(), True + ) + + def _test_external_greybox_solve(self, ex_model, hessian_support): + m = pyo.ConcreteModel() + m.mu = pyo.Var(bounds=(0, None), initialize=1) + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_model) + m.ccon = pyo.Constraint( + expr=m.egb.inputs['c'] == 128 / (3.14 * 1e-4) * m.mu * m.egb.inputs['F'] + ) + m.pcon = pyo.Constraint(expr=m.egb.inputs['Pin'] - m.egb.outputs['Pout'] <= 72) + m.pincon = pyo.Constraint(expr=m.egb.inputs['Pin'] == 100.0) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P1'].value = 80 + m.egb.inputs['P1'].setlb(10) + m.egb.inputs['P1'].setub(90) + m.egb.inputs['P3'].value = 70 + m.egb.inputs['P3'].setlb(20) + m.egb.inputs['P3'].setub(80) + m.egb.outputs['P2'].value = 75 + m.egb.outputs['P2'].setlb(15) + m.egb.outputs['P2'].setub(85) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(10) + m.egb.outputs['Pout'].setub(70) + m.obj = pyo.Objective( + expr=(m.egb.outputs['Pout'] - 20) ** 2 + (m.egb.inputs['F'] - 3) ** 2 + ) + + solver = pyo.SolverFactory('cyipopt') + if not hessian_support: + solver.config.options = {'hessian_approximation': 'limited-memory'} + status = solver.solve(m, tee=False) + + self.assertAlmostEqual(pyo.value(m.egb.inputs['F']), 3.0, places=3) + self.assertAlmostEqual(pyo.value(m.mu), 1.63542e-6, places=3) + self.assertAlmostEqual(pyo.value(m.egb.outputs['Pout']), 28.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.inputs['Pin']), 100.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.inputs['c']), 2.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.inputs['P1']), 82.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.inputs['P3']), 46.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.outputs['P2']), 64.0, places=3) + + def create_model_two_equalities_two_outputs(self, external_model): + m = pyo.ConcreteModel() + m.hin = pyo.Var(bounds=(0, None), initialize=10) + m.hout = pyo.Var(bounds=(0, None)) + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(external_model) + m.incon = pyo.Constraint(expr=0 <= m.egb.inputs['Pin'] - 10 * m.hin) + m.outcon = pyo.Constraint(expr=0 == m.egb.outputs['Pout'] - 10 * m.hout) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P1'].value = 80 + m.egb.inputs['P1'].setlb(10) + m.egb.inputs['P1'].setub(90) + m.egb.inputs['P3'].value = 70 + m.egb.inputs['P3'].setlb(20) + m.egb.inputs['P3'].setub(80) + m.egb.outputs['P2'].value = 75 + m.egb.outputs['P2'].setlb(15) + m.egb.outputs['P2'].setub(85) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(30) + m.egb.outputs['Pout'].setub(70) + return m + + def test_scaling_all_missing(self): + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputs() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + fs = pyomo_nlp.get_obj_scaling() + xs = pyomo_nlp.get_primals_scaling() + cs = pyomo_nlp.get_constraints_scaling() + self.assertIsNone(fs) + self.assertIsNone(xs) + self.assertIsNone(cs) + + def test_scaling_pyomo_model_only(self): + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputs() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + # m.scaling_factor[m.obj] = 0.1 # scale the objective + m.scaling_factor[m.egb.inputs['Pin']] = 1.1 # scale the variable + m.scaling_factor[m.egb.inputs['c']] = 1.2 # scale the variable + m.scaling_factor[m.egb.inputs['F']] = 1.3 # scale the variable + # m.scaling_factor[m.egb.inputs['P1']] = 1.4 # scale the variable + m.scaling_factor[m.egb.inputs['P3']] = 1.5 # scale the variable + m.scaling_factor[m.egb.outputs['P2']] = 1.6 # scale the variable + m.scaling_factor[m.egb.outputs['Pout']] = 1.7 # scale the variable + # m.scaling_factor[m.hin] = 1.8 + m.scaling_factor[m.hout] = 1.9 + # m.scaling_factor[m.incon] = 2.1 + m.scaling_factor[m.outcon] = 2.2 + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + 'hin', + 'hout', + ] + x_order = pyomo_nlp.primals_names() + comparison_c_order = [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + 'incon', + 'outcon', + ] + c_order = pyomo_nlp.constraint_names() + + fs = pyomo_nlp.get_obj_scaling() + self.assertEqual(fs, 1.0) + + xs = pyomo_nlp.get_primals_scaling() + comparison_xs = np.asarray( + [1.1, 1.2, 1.3, 1.0, 1.5, 1.6, 1.7, 1.0, 1.9], dtype=np.float64 + ) + check_vectors_specific_order( + self, xs, x_order, comparison_xs, comparison_x_order + ) + + cs = pyomo_nlp.get_constraints_scaling() + comparison_cs = np.asarray([1, 1, 1, 1, 1, 2.2], dtype=np.float64) + check_vectors_specific_order( + self, cs, c_order, comparison_cs, comparison_c_order + ) + + def test_scaling_greybox_only(self): + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputsScaleBoth() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + 'hin', + 'hout', + ] + x_order = pyomo_nlp.primals_names() + comparison_c_order = [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + 'incon', + 'outcon', + ] + c_order = pyomo_nlp.constraint_names() + + fs = pyomo_nlp.get_obj_scaling() + self.assertEqual(fs, 1.0) + + xs = pyomo_nlp.get_primals_scaling() + comparison_xs = np.asarray([1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=np.float64) + check_vectors_specific_order( + self, xs, x_order, comparison_xs, comparison_x_order + ) + + cs = pyomo_nlp.get_constraints_scaling() + comparison_cs = np.asarray([3.1, 3.2, 4.1, 4.2, 1, 1], dtype=np.float64) + check_vectors_specific_order( + self, cs, c_order, comparison_cs, comparison_c_order + ) + + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputsScaleEqualities() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + cs = pyomo_nlp.get_constraints_scaling() + comparison_cs = np.asarray([3.1, 3.2, 1, 1, 1, 1], dtype=np.float64) + check_vectors_specific_order( + self, cs, c_order, comparison_cs, comparison_c_order + ) + + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputsScaleOutputs() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + cs = pyomo_nlp.get_constraints_scaling() + comparison_cs = np.asarray([1, 1, 4.1, 4.2, 1, 1], dtype=np.float64) + check_vectors_specific_order( + self, cs, c_order, comparison_cs, comparison_c_order + ) + + def test_scaling_pyomo_model_and_greybox(self): + m = self.create_model_two_equalities_two_outputs( + ex_models.PressureDropTwoEqualitiesTwoOutputsScaleBoth() + ) + m.obj = pyo.Objective(expr=(m.egb.outputs['Pout'] - 20) ** 2) + m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + # m.scaling_factor[m.obj] = 0.1 # scale the objective + m.scaling_factor[m.egb.inputs['Pin']] = 1.1 # scale the variable + m.scaling_factor[m.egb.inputs['c']] = 1.2 # scale the variable + m.scaling_factor[m.egb.inputs['F']] = 1.3 # scale the variable + # m.scaling_factor[m.egb.inputs['P1']] = 1.4 # scale the variable + m.scaling_factor[m.egb.inputs['P3']] = 1.5 # scale the variable + m.scaling_factor[m.egb.outputs['P2']] = 1.6 # scale the variable + m.scaling_factor[m.egb.outputs['Pout']] = 1.7 # scale the variable + # m.scaling_factor[m.hin] = 1.8 + m.scaling_factor[m.hout] = 1.9 + # m.scaling_factor[m.incon] = 2.1 + m.scaling_factor[m.outcon] = 2.2 + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m) + + comparison_x_order = [ + 'egb.inputs[Pin]', + 'egb.inputs[c]', + 'egb.inputs[F]', + 'egb.inputs[P1]', + 'egb.inputs[P3]', + 'egb.outputs[P2]', + 'egb.outputs[Pout]', + 'hin', + 'hout', + ] + x_order = pyomo_nlp.primals_names() + comparison_c_order = [ + 'egb.eq_constraints[pdrop1]', + 'egb.eq_constraints[pdrop3]', + 'egb.output_constraints[P2]', + 'egb.output_constraints[Pout]', + 'incon', + 'outcon', + ] + c_order = pyomo_nlp.constraint_names() + + fs = pyomo_nlp.get_obj_scaling() + self.assertEqual(fs, 1.0) + + xs = pyomo_nlp.get_primals_scaling() + comparison_xs = np.asarray( + [1.1, 1.2, 1.3, 1.0, 1.5, 1.6, 1.7, 1.0, 1.9], dtype=np.float64 + ) + check_vectors_specific_order( + self, xs, x_order, comparison_xs, comparison_x_order + ) + + cs = pyomo_nlp.get_constraints_scaling() + comparison_cs = np.asarray([3.1, 3.2, 4.1, 4.2, 1, 2.2], dtype=np.float64) + check_vectors_specific_order( + self, cs, c_order, comparison_cs, comparison_c_order + ) + + @unittest.skipIf(not cyipopt_available, "CyIpopt needed to run tests with solve") + def test_external_greybox_solve_scaling(self): + m = pyo.ConcreteModel() + m.mu = pyo.Var(bounds=(0, None), initialize=1) + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model( + ex_models.PressureDropTwoEqualitiesTwoOutputsScaleBoth() + ) + m.ccon = pyo.Constraint( + expr=m.egb.inputs['c'] == 128 / (3.14 * 1e-4) * m.mu * m.egb.inputs['F'] + ) + m.pcon = pyo.Constraint(expr=m.egb.inputs['Pin'] - m.egb.outputs['Pout'] <= 72) + m.pincon = pyo.Constraint(expr=m.egb.inputs['Pin'] == 100.0) + m.egb.inputs['Pin'].value = 100 + m.egb.inputs['Pin'].setlb(50) + m.egb.inputs['Pin'].setub(150) + m.egb.inputs['c'].value = 2 + m.egb.inputs['c'].setlb(1) + m.egb.inputs['c'].setub(5) + m.egb.inputs['F'].value = 3 + m.egb.inputs['F'].setlb(1) + m.egb.inputs['F'].setub(5) + m.egb.inputs['P1'].value = 80 + m.egb.inputs['P1'].setlb(10) + m.egb.inputs['P1'].setub(90) + m.egb.inputs['P3'].value = 70 + m.egb.inputs['P3'].setlb(20) + m.egb.inputs['P3'].setub(80) + m.egb.outputs['P2'].value = 75 + m.egb.outputs['P2'].setlb(15) + m.egb.outputs['P2'].setub(85) + m.egb.outputs['Pout'].value = 50 + m.egb.outputs['Pout'].setlb(10) + m.egb.outputs['Pout'].setub(70) + m.obj = pyo.Objective( + expr=(m.egb.outputs['Pout'] - 20) ** 2 + (m.egb.inputs['F'] - 3) ** 2 + ) + + m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + m.scaling_factor[m.obj] = 0.1 # scale the objective + m.scaling_factor[m.egb.inputs['Pin']] = 1.1 # scale the variable + m.scaling_factor[m.egb.inputs['c']] = 1.2 # scale the variable + m.scaling_factor[m.egb.inputs['F']] = 1.3 # scale the variable + # m.scaling_factor[m.egb.inputs['P1']] = 1.4 # scale the variable + m.scaling_factor[m.egb.inputs['P3']] = 1.5 # scale the variable + m.scaling_factor[m.egb.outputs['P2']] = 1.6 # scale the variable + m.scaling_factor[m.egb.outputs['Pout']] = 1.7 # scale the variable + m.scaling_factor[m.mu] = 1.9 + m.scaling_factor[m.pincon] = 2.2 + + with TempfileManager.new_context() as temp: + logfile = temp.create_tempfile('_cyipopt-external-greybox-scaling.log') + solver = pyo.SolverFactory('cyipopt') + solver.config.options = { + 'hessian_approximation': 'limited-memory', + 'nlp_scaling_method': 'user-scaling', + 'output_file': logfile, + 'file_print_level': 10, + 'max_iter': 0, + } + status = solver.solve(m, tee=False) + + with open(logfile, 'r') as fd: + solver_trace = fd.read() + + self.assertIn('nlp_scaling_method = user-scaling', solver_trace) + self.assertIn(f'output_file = {logfile}', solver_trace) + self.assertIn('objective scaling factor = 0.1', solver_trace) + self.assertIn('x scaling provided', solver_trace) + self.assertIn('c scaling provided', solver_trace) + self.assertIn('d scaling provided', solver_trace) + # x_order: ['egb.inputs[F]', 'egb.inputs[P1]', 'egb.inputs[P3]', 'egb.inputs[Pin]', 'egb.inputs[c]', 'egb.outputs[P2]', 'egb.outputs[Pout]', 'mu'] + # c_order: ['ccon', 'pcon', 'pincon', 'egb.eq_constraints[pdrop1]', 'egb.eq_constraints[pdrop3]', 'egb.output_constraints[P2]', 'egb.output_constraints[Pout]'] + self.assertIn('DenseVector "x scaling vector" with 8 elements:', solver_trace) + self.assertIn( + 'x scaling vector[ 1]= 1.3000000000000000e+00', solver_trace + ) # F + self.assertIn( + 'x scaling vector[ 8]= 1.8999999999999999e+00', solver_trace + ) # mu + self.assertIn( + 'x scaling vector[ 7]= 1.7000000000000000e+00', solver_trace + ) # Pout + self.assertIn( + 'x scaling vector[ 4]= 1.1000000000000001e+00', solver_trace + ) # Pin + self.assertIn( + 'x scaling vector[ 5]= 1.2000000000000000e+00', solver_trace + ) # c + self.assertIn( + 'x scaling vector[ 2]= 1.0000000000000000e+00', solver_trace + ) # P1 + self.assertIn( + 'x scaling vector[ 3]= 1.5000000000000000e+00', solver_trace + ) # P3 + self.assertIn( + 'x scaling vector[ 6]= 1.6000000000000001e+00', solver_trace + ) # P2 + self.assertIn('DenseVector "c scaling vector" with 6 elements:', solver_trace) + self.assertIn( + 'c scaling vector[ 1]= 1.0000000000000000e+00', solver_trace + ) # ccon + self.assertIn( + 'c scaling vector[ 2]= 2.2000000000000002e+00', solver_trace + ) # pincon + self.assertIn( + 'c scaling vector[ 3]= 3.1000000000000001e+00', solver_trace + ) # eq_constraints[pdrop1] + self.assertIn( + 'c scaling vector[ 4]= 3.2000000000000002e+00', solver_trace + ) # eq_constraints[pdrop3] + self.assertIn( + 'c scaling vector[ 5]= 4.0999999999999996e+00', solver_trace + ) # P2_con + self.assertIn( + 'c scaling vector[ 6]= 4.2000000000000002e+00', solver_trace + ) # Pout_con + self.assertIn('DenseVector "d scaling vector" with 1 elements:', solver_trace) + self.assertIn( + 'd scaling vector[ 1]= 1.0000000000000000e+00', solver_trace + ) # pcon + + @unittest.skipIf(not cyipopt_available, "CyIpopt needed to run tests with solve") + def test_duals_after_solve(self): + m = pyo.ConcreteModel() + m.p = pyo.Var(initialize=1) + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.OneOutput()) + m.con = pyo.Constraint(expr=4 * m.p - 2 * m.egb.outputs['o'] == 0) + m.obj = pyo.Objective(expr=10 * m.p**2) + + # we want to check dual information so we need the suffixes + m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT) + m.ipopt_zL_out = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT) + m.ipopt_zU_out = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT) + + solver = pyo.SolverFactory('cyipopt') + status = solver.solve(m, tee=False) + + self.assertAlmostEqual(pyo.value(m.p), 10.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.inputs['u']), 4.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.outputs['o']), 20.0, places=3) + self.assertAlmostEqual(pyo.value(m.dual[m.con]), 50.0, places=3) + self.assertAlmostEqual( + m.dual[m.egb]['egb.output_constraints[o]'], -100.0, places=3 + ) + self.assertAlmostEqual( + pyo.value(m.ipopt_zL_out[m.egb.inputs['u']]), 500.0, places=3 + ) + self.assertAlmostEqual( + pyo.value(m.ipopt_zU_out[m.egb.inputs['u']]), 0.0, places=3 + ) + + del m.obj + m.obj = pyo.Objective(expr=-10 * m.p**2) + status = solver.solve(m, tee=False) + + self.assertAlmostEqual(pyo.value(m.p), 25.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.inputs['u']), 10.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.outputs['o']), 50.0, places=3) + self.assertAlmostEqual(pyo.value(m.dual[m.con]), -125.0, places=3) + self.assertAlmostEqual( + m.dual[m.egb]['egb.output_constraints[o]'], 250.0, places=3 + ) + self.assertAlmostEqual( + pyo.value(m.ipopt_zL_out[m.egb.inputs['u']]), 0.0, places=3 + ) + self.assertAlmostEqual( + pyo.value(m.ipopt_zU_out[m.egb.inputs['u']]), -1250.0, places=3 + ) + + m = pyo.ConcreteModel() + m.p = pyo.Var(initialize=1) + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model(ex_models.OneOutputOneEquality()) + m.con = pyo.Constraint(expr=4 * m.p - 2 * m.egb.outputs['o'] == 0) + m.obj = pyo.Objective(expr=10 * m.p**2) + + # we want to check dual information so we need the suffixes + m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT) + + solver = pyo.SolverFactory('cyipopt') + status = solver.solve(m, tee=False) + self.assertAlmostEqual(pyo.value(m.p), 2.5, places=3) + self.assertAlmostEqual(pyo.value(m.egb.inputs['u']), 1.0, places=3) + self.assertAlmostEqual(pyo.value(m.egb.outputs['o']), 5.0, places=3) + self.assertAlmostEqual(pyo.value(m.dual[m.con]), 12.5, places=3) + self.assertAlmostEqual( + m.dual[m.egb]['egb.output_constraints[o]'], -25.0, places=3 + ) + self.assertAlmostEqual( + m.dual[m.egb]['egb.eq_constraints[u2_con]'], 62.5, places=3 + ) + + +class TestGreyBoxObjectives(unittest.TestCase): + @unittest.skipIf(not cyipopt_available, "CyIpopt needed to run tests with solve") + def test_unconstrained(self): + solve_unconstrained() + + @unittest.skipIf(not cyipopt_available, "CyIpopt needed to run tests with solve") + def test_constrained(self): + solve_constrained() + + @unittest.skipIf(not cyipopt_available, "CyIpopt needed to run tests with solve") + def test_constrained_with_hessian(self): + solve_constrained_with_hessian() + + +# Regression tests to make sure PyomoNLPWithGreyBoxBlocks correctly handles variables that +# are external to the block when creating the NLP, but references in the constraints. +class TestPyomoNLPWithGreyBoxModelsExternalVars(unittest.TestCase): + def test_no_greybox_block(self): + m = pyo.ConcreteModel() + # Variable on the main model + m.x = pyo.Var(initialize=1) + + # One block contains constraints that reference the variable on the main model + m.b = pyo.Block() + m.b.y = pyo.Var(initialize=2) + + m.b.cons1 = pyo.Constraint(expr=m.x + 2 * m.b.y == 5) + m.b.cons2 = pyo.Constraint(expr=3 * m.x - 4 * m.b.y == -5) + + # Create NLP from m.b - should contain m.v even though it is external to the block + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m.b) + + self.assertEqual( + pyomo_nlp._pyomo_model_var_names_to_datas, {'x': m.x, 'b.y': m.b.y} + ) + + jac = pyomo_nlp.evaluate_jacobian().tocsr() + + # Due to external variable, the order is m.b.y, m.x + self.assertEqual(jac.shape, (2, 2)) + self.assertEqual(jac[0, 0], 2.0) + self.assertEqual(jac[0, 1], 1.0) + self.assertEqual(jac[1, 0], -4.0) + self.assertEqual(jac[1, 1], 3.0) + + def test_greybox_block_w_external_var(self): + m = pyo.ConcreteModel() + m.v = pyo.Var() + + m.b = pyo.Block() + m.b.egb = ExternalGreyBoxBlock() + m.b.egb.set_external_model(ex_models.PressureDropSingleOutput()) + + # Set egb variable values + m.b.egb.inputs['Pin'].value = 100 + m.b.egb.inputs['c'].value = 2 + m.b.egb.inputs['F'].value = 3 + m.b.egb.outputs['Pout'].value = 80 + + m.b.cons = pyo.Constraint(expr=m.v == m.b.egb.inputs['F']) + + # Create NLP from m.b - should contain m.v even though it is external to the block + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m.b) + + self.assertEqual( + pyomo_nlp._pyomo_model_var_names_to_datas, + { + 'v': m.v, + 'b.egb.inputs[Pin]': m.b.egb.inputs['Pin'], + 'b.egb.inputs[c]': m.b.egb.inputs['c'], + 'b.egb.inputs[F]': m.b.egb.inputs['F'], + 'b.egb.outputs[Pout]': m.b.egb.outputs['Pout'], + }, + ) + + jac = pyomo_nlp.evaluate_jacobian().tocsr() + + self.assertEqual(jac.shape, (2, 5)) + primals = pyomo_nlp.primals_names() + constraints = pyomo_nlp.constraint_names() + + expected = { + ('b.cons', 'v'): 1.0, + ('b.cons', 'b.egb.inputs[F]'): -1.0, + ('b.cons', 'b.egb.inputs[Pin]'): 0.0, + ('b.cons', 'b.egb.inputs[c]'): 0.0, + ('b.cons', 'b.egb.outputs[Pout]'): 0.0, + ('b.egb.output_constraints[Pout]', 'v'): 0.0, + ('b.egb.output_constraints[Pout]', 'b.egb.inputs[F]'): -48.0, # -4*c*2*F + ('b.egb.output_constraints[Pout]', 'b.egb.inputs[Pin]'): 1.0, + ('b.egb.output_constraints[Pout]', 'b.egb.inputs[c]'): -36.0, # -4*F**2 + ('b.egb.output_constraints[Pout]', 'b.egb.outputs[Pout]'): -1.0, + } + + for (c, v), val in expected.items(): + print(c, v) + self.assertAlmostEqual(jac[constraints.index(c), primals.index(v)], val) + + def test_greybox_block_w_constraints_w_external_var(self): + m = pyo.ConcreteModel() + m.v = pyo.Var() + + m.b = pyo.Block() + m.b.egb = ExternalGreyBoxBlock() + m.b.egb.set_external_model(ex_models.PressureDropSingleOutput()) + + # Set egb variable values + m.b.egb.inputs['Pin'].value = 100 + m.b.egb.inputs['c'].value = 2 + m.b.egb.inputs['F'].value = 3 + m.b.egb.outputs['Pout'].value = 80 + + m.b.cons = pyo.Constraint(expr=m.v == m.b.egb.inputs['F']) + + # Create NLP from m.b - should contain m.v even though it is external to the block + pyomo_nlp = PyomoNLPWithGreyBoxBlocks(m.b) + + self.assertEqual( + pyomo_nlp._pyomo_model_var_names_to_datas, + { + 'v': m.v, + 'b.egb.inputs[Pin]': m.b.egb.inputs['Pin'], + 'b.egb.inputs[c]': m.b.egb.inputs['c'], + 'b.egb.inputs[F]': m.b.egb.inputs['F'], + 'b.egb.outputs[Pout]': m.b.egb.outputs['Pout'], + }, + ) + + jac = pyomo_nlp.evaluate_jacobian().tocsr() + + self.assertEqual(jac.shape, (2, 5)) + primals = pyomo_nlp.primals_names() + constraints = pyomo_nlp.constraint_names() + + expected = { + ('b.cons', 'v'): 1.0, + ('b.cons', 'b.egb.inputs[F]'): -1.0, + ('b.cons', 'b.egb.inputs[Pin]'): 0.0, + ('b.cons', 'b.egb.inputs[c]'): 0.0, + ('b.cons', 'b.egb.outputs[Pout]'): 0.0, + ('b.egb.output_constraints[Pout]', 'v'): 0.0, + ('b.egb.output_constraints[Pout]', 'b.egb.inputs[F]'): -48.0, # -4*c*2*F + ('b.egb.output_constraints[Pout]', 'b.egb.inputs[Pin]'): 1.0, + ('b.egb.output_constraints[Pout]', 'b.egb.inputs[c]'): -36.0, # -4*F**2 + ('b.egb.output_constraints[Pout]', 'b.egb.outputs[Pout]'): -1.0, + } + + for (c, v), val in expected.items(): + self.assertAlmostEqual(jac[constraints.index(c), primals.index(v)], val) + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/pynumero/interfaces/tests/test_utils.py b/pyomo/contrib/pynumero/interfaces/tests/test_utils.py index 54a801edce8..0a6884b1d75 100644 --- a/pyomo/contrib/pynumero/interfaces/tests/test_utils.py +++ b/pyomo/contrib/pynumero/interfaces/tests/test_utils.py @@ -6,7 +6,11 @@ # 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.environ as pyo import pyomo.common.unittest as unittest from pyomo.contrib.pynumero.dependencies import ( numpy as np, @@ -19,6 +23,14 @@ raise unittest.SkipTest("Pynumero needs scipy and numpy to run NLP tests") import pyomo.contrib.pynumero.interfaces.utils as utils +from pyomo.contrib.pynumero.interfaces.pyomo_grey_box_nlp import ( + PyomoNLPWithGreyBoxBlocks, +) + +from pyomo.contrib.pynumero.asl import AmplInterface + +if not AmplInterface.available(): + raise unittest.SkipTest("Pynumero ASL interface is not available") class TestCondensedSparseSummation(unittest.TestCase): @@ -83,6 +95,46 @@ def test_repeated_row_col(self): self.assertTrue(np.array_equal(expected_row, C.row)) self.assertTrue(np.array_equal(expected_col, C.col)) + def test_empty_hessian_linear_model_does_not_crash(self): + """ + Regression test for CondensedSparseSummation bug where an empty + union of nonzeros caused unpack failure. + + Linear model -> Lagrangian Hessian has no nonzeros. + Previously this crashed inside CondensedSparseSummation._build_maps(). + """ + + m = pyo.ConcreteModel() + + # Linear model (same structure as user's example) + m.v1 = pyo.Var(initialize=1e-8) + m.v2 = pyo.Var(initialize=1) + m.v3 = pyo.Var(initialize=1) + + m.c1 = pyo.Constraint(expr=m.v1 == m.v2) + m.c2 = pyo.Constraint(expr=m.v1 == 1e-8 * m.v3) + m.c3 = pyo.Constraint(expr=1e8 * m.v1 + 1e10 * m.v2 == 1e-6 * m.v3) + + # PyNumero requires an objective + m.obj = pyo.Objective(expr=0) + + # This used to crash during initialization due to empty + # union of nonzeros in the Hessian of the Lagrangian + nlp = PyomoNLPWithGreyBoxBlocks(m) + + # Explicitly evaluate Hessian of Lagrangian + hess = nlp.evaluate_hessian_lag() + + # Shape should match number of primals + n = nlp.n_primals() + self.assertEqual(hess.shape, (n, n)) + + # For a purely linear model, Hessian must be structurally empty + self.assertEqual(hess.nnz, 0) + + # Also verify the sparse data vector is empty + self.assertEqual(len(hess.data), 0) + if __name__ == '__main__': TestCondensedSparseSummation().test_condensed_sparse_summation() diff --git a/pyomo/contrib/pynumero/interfaces/utils.py b/pyomo/contrib/pynumero/interfaces/utils.py index d18bd18378b..271e757b25f 100644 --- a/pyomo/contrib/pynumero/interfaces/utils.py +++ b/pyomo/contrib/pynumero/interfaces/utils.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. +# ___________________________________________________________________________________ import numpy as np from scipy.sparse import coo_matrix from pyomo.contrib.pynumero.sparse import BlockVector, BlockMatrix @@ -166,6 +169,27 @@ def _build_maps(self, list_of_matrices): nz_tuples.update(zip(m.row, m.col)) nz_tuples = sorted(nz_tuples) self._nz_tuples = nz_tuples + + # Check to ensure we have some non-zeroes. + # If not, we still need to maintain consistent internal state and maps so sum() returns a valid empty matrix. + if len(nz_tuples) == 0: + # All matrices have empty nonzero structure (e.g., linear models with empty Hessian). + # Maintain consistent internal state and maps so sum() returns a valid empty matrix. + self._row = np.array([], dtype=np.int64) + self._col = np.array([], dtype=np.int64) + self._shape = None + self._maps = [] + for m in list_of_matrices: + nnz = len(m.data) + # Map from each matrix's data vector (length nnz, usually 0) into the + # condensed data vector (length 0). This must exist for sum(). + self._maps.append(coo_matrix((0, nnz))) + if self._shape is None: + self._shape = m.shape + else: + assert self._shape == m.shape + return + self._row, self._col = list(zip(*nz_tuples)) row_col_to_nz_map = {t: i for i, t in enumerate(nz_tuples)}