diff --git a/idaes/core/scaling/tests/test_get_jacobian_greybox.py b/idaes/core/scaling/tests/test_get_jacobian_greybox.py new file mode 100644 index 0000000000..0edc7cacdd --- /dev/null +++ b/idaes/core/scaling/tests/test_get_jacobian_greybox.py @@ -0,0 +1,504 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2026 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Tests for scaling utility functions. + +Author: Andrew Lee, Douglas Allan +""" + +import re + +import pytest + +from pyomo.environ import ( + assert_optimal_termination, + Constraint, + ConcreteModel, + SolverFactory, + Suffix, + Var, +) +from pyomo.contrib.pynumero.asl import AmplInterface +from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock +import pyomo.contrib.pynumero.interfaces.tests.external_grey_box_models as ex_models + +from idaes.core.scaling.util import ( + get_jacobian, + _get_jacobian_greybox_compatible, + jacobian_cond, + set_scaling_factor, +) +import idaes.logger as idaeslog + + +@pytest.mark.component +def test_get_jacobian_greybox_error_with_egb_constraints(): + m = ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model( + ex_models.PressureDropSingleOutput(), + build_implicit_constraint_objects=True, + ) + + # Add linking variables and constraints + m.Pin = Var(initialize=1e5) + m.c = Var(initialize=1e3) + m.F = Var(initialize=0.5) + m.Pout = Var(initialize=1e5) + + m.link_P_in = Constraint(expr=m.Pin == m.egb.inputs["Pin"]) + m.link_c = Constraint(expr=m.c == m.egb.inputs["c"]) + m.link_F = Constraint(expr=m.F == m.egb.inputs["F"]) + m.link_P_out = Constraint(expr=m.Pout == m.egb.outputs["Pout"]) + + m.Pin.fix() + m.c.fix() + m.F.fix() + + # Get Jacobian + with pytest.raises( + ValueError, + match="The model contains components that are not supported by the Pyomo NL writer. " + "This may be because the model contains a grey-box model. If you want to include " + "grey-box variables and constraints in the Jacobian, set " + "include_greybox=True when calling get_jacobian.", + ): + get_jacobian(m, include_greybox=False) + + +# Parity tests to ensure grey box compatible Jacobian matches legacy tools +# when no grey box is present +# Tests based on John Eslick's originals for the old scaling tools +@pytest.mark.skipif( + not AmplInterface.available(), reason="pynumero_ASL is not available" +) +class TestJacobianMethodsNoGreyBox: + @pytest.fixture + def model(self): + m = ConcreteModel() + x = m.x = Var(initialize=1e3) + y = m.y = Var(initialize=1e6) + z = m.z = Var(initialize=1e4) + m.c1 = Constraint(expr=0 == -x * y + z) + m.c2 = Constraint(expr=0 == 3 * x + 4 * y + 2 * z) + m.c3 = Constraint(expr=0 <= z**3) + return m + + @pytest.mark.unit + def test_jacobian(self, model): + """Make sure the Jacobian from Pynumero matches expectation. This is + mostly to ensure we understand the interface and catch if things change. + """ + m = model + jac, nlp = _get_jacobian_greybox_compatible(m) + assert not hasattr(m, "scaling_factor") + assert not hasattr(m, "scaling_hint") + + c1_row = nlp.constraint_names().index("c1") + c2_row = nlp.constraint_names().index("c2") + c3_row = nlp.constraint_names().index("c3") + x_col = nlp.primals_names().index("x") + y_col = nlp.primals_names().index("y") + z_col = nlp.primals_names().index("z") + + assert jac[c1_row, x_col] == pytest.approx(-1e6) + assert jac[c1_row, y_col] == pytest.approx(-1e3) + assert jac[c1_row, z_col] == pytest.approx(1) + + assert jac[c2_row, x_col] == pytest.approx(3) + assert jac[c2_row, y_col] == pytest.approx(4) + assert jac[c2_row, z_col] == pytest.approx(2) + + assert jac[c3_row, z_col] == pytest.approx(3e8) + + # Make sure scaling factors don't affect the result + set_scaling_factor(m.c1, 1e-6) + set_scaling_factor(m.x, 1e-3) + set_scaling_factor(m.y, 1e-6) + set_scaling_factor(m.z, 1e-4) + jac, _ = _get_jacobian_greybox_compatible(m, include_scaling_factors=False) + assert len(m.scaling_factor) == 4 + assert not hasattr(m, "scaling_hint") + assert jac[c1_row, x_col] == pytest.approx(-1e6) + + # Check the scaled jacobian calculation + jac_scaled, _ = _get_jacobian_greybox_compatible(m) + assert len(m.scaling_factor) == 4 + assert not hasattr(m, "scaling_hint") + assert jac_scaled[c1_row, x_col] == pytest.approx(-1000) + assert jac_scaled[c1_row, y_col] == pytest.approx(-1000) + assert jac_scaled[c1_row, z_col] == pytest.approx(0.01) + + @pytest.mark.unit + def test_scale_no_var_scale(self, model): + m = model + jac_scaled, nlp = _get_jacobian_greybox_compatible( + m, + include_scaling_factors=False, + include_ipopt_autoscaling=True, + min_scale=1e-6, + ) + # get_scaling_factor isn't called here so the suffix shouldn't exist + assert not hasattr(m, "scaling_factor") + assert not hasattr(m, "scaling_hint") + + c1_row = nlp.constraint_names().index("c1") + c2_row = nlp.constraint_names().index("c2") + c3_row = nlp.constraint_names().index("c3") + x_col = nlp.primals_names().index("x") + y_col = nlp.primals_names().index("y") + z_col = nlp.primals_names().index("z") + + assert jac_scaled[c1_row, x_col] == pytest.approx(-100) + assert jac_scaled[c1_row, y_col] == pytest.approx(-0.1) + assert jac_scaled[c1_row, z_col] == pytest.approx(1e-4) + + assert jac_scaled[c2_row, x_col] == pytest.approx(3) + assert jac_scaled[c2_row, y_col] == pytest.approx(4) + assert jac_scaled[c2_row, z_col] == pytest.approx(2) + + assert jac_scaled[c3_row, z_col] == pytest.approx(3e2) + + @pytest.mark.unit + def test_scale_with_var_scale(self, model): + m = model + m.scaling_factor = Suffix(direction=Suffix.EXPORT) + m.scaling_factor[m.x] = 1e-3 + m.scaling_factor[m.y] = 1e-6 + m.scaling_factor[m.z] = 1e-4 + + # In these tests derived from the old scaling tools, we use + # min_scale=1e-6 because that's what the old tools use as a + # default. However, we're using a new default of 1e-8 + # because that's what appears to be IPOPT's actual default. + jac_scaled, nlp = _get_jacobian_greybox_compatible( + m, + include_scaling_factors=True, + include_ipopt_autoscaling=True, + min_scale=1e-6, + ) + assert len(m.scaling_factor) == 3 + assert not hasattr(m, "scaling_hint") + + c1_row = nlp.constraint_names().index("c1") + c2_row = nlp.constraint_names().index("c2") + c3_row = nlp.constraint_names().index("c3") + x_col = nlp.primals_names().index("x") + y_col = nlp.primals_names().index("y") + z_col = nlp.primals_names().index("z") + + assert jac_scaled[c1_row, x_col] == pytest.approx(-1000) + assert jac_scaled[c1_row, y_col] == pytest.approx(-1000) + assert jac_scaled[c1_row, z_col] == pytest.approx(1e-2) + + assert jac_scaled[c2_row, x_col] == pytest.approx(0.075) + assert jac_scaled[c2_row, y_col] == pytest.approx(100) + assert jac_scaled[c2_row, z_col] == pytest.approx(0.5) + + assert jac_scaled[c3_row, z_col] == pytest.approx(3e6) + + @pytest.mark.unit + def test_exclude_scaling_factors_variables(self, model): + m = model + m.scaling_factor = Suffix(direction=Suffix.EXPORT) + m.scaling_factor[m.x] = 1e-3 + m.scaling_factor[m.y] = 1e-6 + m.scaling_factor[m.z] = 1e-4 + + jac_scaled, nlp = _get_jacobian_greybox_compatible( + m, + include_scaling_factors=False, + include_ipopt_autoscaling=True, + min_scale=1e-6, + ) + assert len(m.scaling_factor) == 3 + assert not hasattr(m, "scaling_hint") + + c1_row = nlp.constraint_names().index("c1") + c2_row = nlp.constraint_names().index("c2") + c3_row = nlp.constraint_names().index("c3") + x_col = nlp.primals_names().index("x") + y_col = nlp.primals_names().index("y") + z_col = nlp.primals_names().index("z") + + assert jac_scaled[c1_row, x_col] == pytest.approx(-100) + assert jac_scaled[c1_row, y_col] == pytest.approx(-0.1) + assert jac_scaled[c1_row, z_col] == pytest.approx(1e-4) + + assert jac_scaled[c2_row, x_col] == pytest.approx(3) + assert jac_scaled[c2_row, y_col] == pytest.approx(4) + assert jac_scaled[c2_row, z_col] == pytest.approx(2) + + assert jac_scaled[c3_row, z_col] == pytest.approx(3e2) + + @pytest.mark.unit + def test_condition_number(self, model, caplog): + """Calculate the condition number of the Jacobian""" + m = model + m.scaling_factor = Suffix(direction=Suffix.EXPORT) + m.scaling_factor[m.x] = 1e-3 + m.scaling_factor[m.y] = 1e-6 + m.scaling_factor[m.z] = 1e-4 + m.scaling_factor[m.c1] = 1e-6 + m.scaling_factor[m.c2] = 1e-6 + m.scaling_factor[m.c3] = 1e-12 + + jac, _ = _get_jacobian_greybox_compatible( + m, + include_scaling_factors=False, + include_ipopt_autoscaling=False, + ) + jac_scaled, _ = _get_jacobian_greybox_compatible( + m, + include_scaling_factors=True, + include_ipopt_autoscaling=False, + min_scale=1e-6, + ) + + n = jacobian_cond(m, jac=jac_scaled) + assert n == pytest.approx(687.47, rel=1e-3) + n = jacobian_cond(m, jac=jac) + assert n == pytest.approx(7.50567e7, rel=1e-3) + + # Nonsquare condition number + m.c3.deactivate() + jac, _ = _get_jacobian_greybox_compatible( + m, + include_scaling_factors=False, + include_ipopt_autoscaling=False, + ) + jac_scaled, _ = _get_jacobian_greybox_compatible( + m, + include_scaling_factors=True, + include_ipopt_autoscaling=False, + min_scale=1e-6, + ) + + # Scaled + with caplog.at_level(idaeslog.INFO): + n = jacobian_cond(m, jac=jac_scaled) + assert ( + "Nonsquare Jacobian. Using pseudoinverse to calculate Frobenius norm." + ) in caplog.text + assert n == pytest.approx(500.367, rel=1e-3) + + # Unscaled + with caplog.at_level(idaeslog.INFO): + n = jacobian_cond(m, jac=jac) + assert ( + "Nonsquare Jacobian. Using pseudoinverse to calculate Frobenius norm." + ) in caplog.text + assert n == pytest.approx(2.23741e5, rel=1e-3) + + +@pytest.mark.skipif( + not AmplInterface.available(), reason="pynumero_ASL is not available" +) +class TestJacobianMethodsWithGreyBox: + @pytest.fixture + def model(self): + m = ConcreteModel() + m.egb = ExternalGreyBoxBlock() + m.egb.set_external_model( + ex_models.PressureDropSingleOutput(), + build_implicit_constraint_objects=True, + ) + + # Add linking variables and constraints + m.Pin = Var(initialize=1e5) + m.c = Var(initialize=1e3) + m.F = Var(initialize=0.5) + m.Pout = Var(initialize=1e5) + + m.link_P_in = Constraint(expr=m.Pin == m.egb.inputs["Pin"]) + m.link_c = Constraint(expr=m.c == m.egb.inputs["c"]) + m.link_F = Constraint(expr=m.F == m.egb.inputs["F"]) + m.link_P_out = Constraint(expr=m.Pout == m.egb.outputs["Pout"]) + + m.Pin.fix() + m.c.fix() + m.F.fix() + + solver = SolverFactory("cyipopt") + results = solver.solve(m, tee=False) + assert_optimal_termination(results) + + return m + + @pytest.mark.unit + def test_jacobian(self, model): + """Make sure the Jacobian from Pynumero matches expectation. This is + mostly to ensure we understand the interface and catch if things change. + """ + m = model + jac, nlp = _get_jacobian_greybox_compatible(m, include_scaling_factors=False) + assert not hasattr(m, "scaling_factor") + assert not hasattr(m, "scaling_hint") + + expected_jac = { + ("link_P_in", "egb.inputs[Pin]"): -1.0, + ("link_c", "egb.inputs[c]"): -1.0, + ("link_F", "egb.inputs[F]"): -1.0, + ("link_P_out", "Pout"): 1.0, + ("link_P_out", "egb.outputs[Pout]"): -1.0, + ("egb.Pout_constraint", "egb.inputs[Pin]"): 1.0, + ("egb.Pout_constraint", "egb.inputs[c]"): -1.0, # -4*0.5**2 + ("egb.Pout_constraint", "egb.inputs[F]"): -4000.0, # -4*1e5*2*0.5 + ("egb.Pout_constraint", "egb.outputs[Pout]"): -1.0, + } + + jac_coo = jac.tocoo() + assert len(jac_coo.data) == 9 + for i, j, val in zip(jac_coo.row, jac_coo.col, jac_coo.data): + if (nlp.constraint_names()[i], nlp.primals_names()[j]) not in expected_jac: + assert val == 0 + else: + assert val == pytest.approx( + expected_jac[nlp.constraint_names()[i], nlp.primals_names()[j]] + ) + + # Make sure scaling factors don't affect the result + set_scaling_factor(m.link_P_in, 1e-5) + set_scaling_factor(m.link_c, 1e-3) + set_scaling_factor(m.link_F, 10) + set_scaling_factor(m.link_P_out, 1e-5) + set_scaling_factor(m.Pin, 1e-5) + set_scaling_factor(m.c, 1e-3) + set_scaling_factor(m.F, 10) + set_scaling_factor(m.Pout, 1e-5) + jac, _ = _get_jacobian_greybox_compatible(m, include_scaling_factors=False) + assert len(m.scaling_factor) == 8 + assert not hasattr(m, "scaling_hint") + + # Check the scaled jacobian calculation + jac_coo = jac.tocoo() + assert len(jac_coo.data) == 9 + for i, j, val in zip(jac_coo.row, jac_coo.col, jac_coo.data): + if (nlp.constraint_names()[i], nlp.primals_names()[j]) not in expected_jac: + assert val == 0 + else: + assert val == pytest.approx( + expected_jac[nlp.constraint_names()[i], nlp.primals_names()[j]] + ) + + @pytest.mark.unit + def test_jacobian_equality_only(self, model): + """Make sure the equality_only behaviour works as expected.""" + m = model + + # Add an inequality + m.ineq = Constraint(expr=m.Pout <= 1e6) + + jac, nlp = _get_jacobian_greybox_compatible(m, equality_constraints_only=False) + assert not hasattr(m, "scaling_factor") + assert not hasattr(m, "scaling_hint") + + expected_jac = { + ("link_P_in", "egb.inputs[Pin]"): -1.0, + ("link_c", "egb.inputs[c]"): -1.0, + ("link_F", "egb.inputs[F]"): -1.0, + ("link_P_out", "Pout"): 1.0, + ("link_P_out", "egb.outputs[Pout]"): -1.0, + ("egb.Pout_constraint", "egb.inputs[Pin]"): 1.0, + ("egb.Pout_constraint", "egb.inputs[c]"): -1.0, # -4*0.5**2 + ("egb.Pout_constraint", "egb.inputs[F]"): -4000.0, # -4*1e5*2*0.5 + ("egb.Pout_constraint", "egb.outputs[Pout]"): -1.0, + ("ineq", "Pout"): 1.0, + } + + jac_coo = jac.tocoo() + assert len(jac_coo.data) == 10 + for i, j, val in zip(jac_coo.row, jac_coo.col, jac_coo.data): + if (nlp.constraint_names()[i], nlp.primals_names()[j]) not in expected_jac: + assert val == 0 + else: + assert val == pytest.approx( + expected_jac[nlp.constraint_names()[i], nlp.primals_names()[j]] + ) + + # Test with equality_only=True to make sure the inequality is excluded + jac, nlp = _get_jacobian_greybox_compatible(m, equality_constraints_only=True) + jac_coo = jac.tocoo() + assert len(jac_coo.data) == 9 + + # Create a new mapping for constraint names that leaves out `ineq` + eq_constraints = [c for c in nlp.constraint_names() if not c.startswith("ineq")] + for i, j, val in zip(jac_coo.row, jac_coo.col, jac_coo.data): + if (eq_constraints[i], nlp.primals_names()[j]) not in expected_jac: + assert val == 0 + else: + assert val == pytest.approx( + expected_jac[eq_constraints[i], nlp.primals_names()[j]] + ) + + @pytest.mark.unit + def test_jacobian_w_ipopt_scaling(self, model): + """Make sure the IPOPT scaling behaviour works as expected.""" + m = model + jac, nlp = _get_jacobian_greybox_compatible(m, include_ipopt_autoscaling=True) + assert not hasattr(m, "scaling_factor") + assert not hasattr(m, "scaling_hint") + + expected_jac = { + ("link_P_in", "egb.inputs[Pin]"): -1.0, + ("link_c", "egb.inputs[c]"): -1.0, + ("link_F", "egb.inputs[F]"): -1.0, + ("link_P_out", "Pout"): 1.0, + ("link_P_out", "egb.outputs[Pout]"): -1.0, + # All EGB constraints get multiplied by 0.025 + ("egb.Pout_constraint", "egb.inputs[Pin]"): 0.025, + ("egb.Pout_constraint", "egb.inputs[c]"): -0.025, + ("egb.Pout_constraint", "egb.inputs[F]"): -100.0, + ("egb.Pout_constraint", "egb.outputs[Pout]"): -0.025, + } + + jac_coo = jac.tocoo() + assert len(jac_coo.data) == 9 + for i, j, val in zip(jac_coo.row, jac_coo.col, jac_coo.data): + if (nlp.constraint_names()[i], nlp.primals_names()[j]) not in expected_jac: + assert val == 0 + else: + assert val == pytest.approx( + expected_jac[nlp.constraint_names()[i], nlp.primals_names()[j]] + ) + + @pytest.mark.unit + def test_condition_number(self, model): + """Calculate the condition number of the Jacobian""" + m = model + m.scaling_factor = Suffix(direction=Suffix.EXPORT) + m.scaling_factor[m.link_P_in] = 1e-5 + m.scaling_factor[m.link_c] = 1e-3 + m.scaling_factor[m.link_F] = 10 + m.scaling_factor[m.link_P_out] = 1e-5 + m.scaling_factor[m.Pin] = 1e-5 + m.scaling_factor[m.c] = 1e-3 + m.scaling_factor[m.F] = 10 + m.scaling_factor[m.Pout] = 1e-5 + + jac, _ = _get_jacobian_greybox_compatible( + m, + include_scaling_factors=False, + include_ipopt_autoscaling=False, + ) + n = jacobian_cond(m, jac=jac) + assert n == pytest.approx(2.26274262e7, rel=1e-8) + + jac_scaled, _ = _get_jacobian_greybox_compatible( + m, + include_scaling_factors=True, + include_ipopt_autoscaling=False, + min_scale=1e-6, + ) + n = jacobian_cond(m, jac=jac_scaled) + assert n == pytest.approx(5.657178098e8, rel=1e-3) diff --git a/idaes/core/scaling/tests/test_util.py b/idaes/core/scaling/tests/test_util.py index 19f4df062c..b4c60c0a9e 100644 --- a/idaes/core/scaling/tests/test_util.py +++ b/idaes/core/scaling/tests/test_util.py @@ -37,6 +37,9 @@ from pyomo.common.tempfiles import TempfileManager from pyomo.dae import ContinuousSet, DerivativeVar from pyomo.contrib.pynumero.asl import AmplInterface +from pyomo.contrib.pynumero.interfaces.external_grey_box import ( + ExternalGreyBoxBlock, +) from idaes.core.scaling.util import ( get_jacobian, @@ -66,6 +69,10 @@ from idaes.core.util.scaling import get_constraint_transform_applied_scaling_factor import idaes.logger as idaeslog +from idaes.core.util.diagnostics_tools.tests.test_diagnostics_toolbox_greybox import ( + ResidualGrayBox, +) + currdir = this_file_dir() @@ -1779,6 +1786,19 @@ def test_get_scaling_factor_deactivated_suffix(self): m.scaling_factor.deactivate() assert get_scaling_factor(m.v) is None + @pytest.mark.unit + def test_get_scaling_factor_egb_constraint(self): + m = ConcreteModel() + + m.gb = ExternalGreyBoxBlock( + external_model=ResidualGrayBox(), build_implicit_constraint_objects=True + ) + + assert get_scaling_factor(m.gb.c1) == 1.0 + + # Test with default value + assert get_scaling_factor(m.gb.c1, default=42.0) == 42.0 + class TestSetScalingFactor: @pytest.mark.unit @@ -2398,6 +2418,7 @@ def test_find_unscaled_vars_and_constraints(): @pytest.mark.skipif( not AmplInterface.available(), reason="pynumero_ASL is not available" ) +@pytest.mark.parametrize("include_greybox", [True, False]) class TestJacobianMethods: @pytest.fixture def model(self): @@ -2411,7 +2432,7 @@ def model(self): return m @pytest.mark.unit - def test_jacobian(self, model): + def test_jacobian(self, model, include_greybox): """Make sure the Jacobian from Pynumero matches expectation. This is mostly to ensure we understand the interface and catch if things change. """ @@ -2444,13 +2465,15 @@ def test_jacobian(self, model): set_scaling_factor(m.x, 1e-3) set_scaling_factor(m.y, 1e-6) set_scaling_factor(m.z, 1e-4) - jac, _ = get_jacobian(m, include_scaling_factors=False) + jac, _ = get_jacobian( + m, include_scaling_factors=False, include_greybox=include_greybox + ) assert len(m.scaling_factor) == 4 assert not hasattr(m, "scaling_hint") assert jac[c1_row, x_col] == pytest.approx(-1e6) # Check the scaled jacobian calculation - jac_scaled, _ = get_jacobian(m) + jac_scaled, _ = get_jacobian(m, include_greybox=include_greybox) assert len(m.scaling_factor) == 4 assert not hasattr(m, "scaling_hint") assert jac_scaled[c1_row, x_col] == pytest.approx(-1000) @@ -2458,24 +2481,25 @@ def test_jacobian(self, model): assert jac_scaled[c1_row, z_col] == pytest.approx(0.01) @pytest.mark.unit - def test_scale_no_var_scale(self, model): + def test_scale_no_var_scale(self, model, include_greybox): m = model jac_scaled, nlp = get_jacobian( m, include_scaling_factors=False, include_ipopt_autoscaling=True, min_scale=1e-6, + include_greybox=include_greybox, ) # get_scaling_factor isn't called here so the suffix shouldn't exist assert not hasattr(m, "scaling_factor") assert not hasattr(m, "scaling_hint") - c1_row = nlp._condata_to_idx[m.c1] - c2_row = nlp._condata_to_idx[m.c2] - c3_row = nlp._condata_to_idx[m.c3] - x_col = nlp._vardata_to_idx[m.x] - y_col = nlp._vardata_to_idx[m.y] - z_col = nlp._vardata_to_idx[m.z] + c1_row = nlp.constraint_names().index(m.c1.name) + c2_row = nlp.constraint_names().index(m.c2.name) + c3_row = nlp.constraint_names().index(m.c3.name) + x_col = nlp.primals_names().index(m.x.name) + y_col = nlp.primals_names().index(m.y.name) + z_col = nlp.primals_names().index(m.z.name) assert jac_scaled[c1_row, x_col] == pytest.approx(-100) assert jac_scaled[c1_row, y_col] == pytest.approx(-0.1) @@ -2488,7 +2512,7 @@ def test_scale_no_var_scale(self, model): assert jac_scaled[c3_row, z_col] == pytest.approx(3e2) @pytest.mark.unit - def test_scale_with_var_scale(self, model): + def test_scale_with_var_scale(self, model, include_greybox): m = model m.scaling_factor = Suffix(direction=Suffix.EXPORT) m.scaling_factor[m.x] = 1e-3 @@ -2504,16 +2528,17 @@ def test_scale_with_var_scale(self, model): include_scaling_factors=True, include_ipopt_autoscaling=True, min_scale=1e-6, + include_greybox=include_greybox, ) assert len(m.scaling_factor) == 3 assert not hasattr(m, "scaling_hint") - c1_row = nlp._condata_to_idx[m.c1] - c2_row = nlp._condata_to_idx[m.c2] - c3_row = nlp._condata_to_idx[m.c3] - x_col = nlp._vardata_to_idx[m.x] - y_col = nlp._vardata_to_idx[m.y] - z_col = nlp._vardata_to_idx[m.z] + c1_row = nlp.constraint_names().index(m.c1.name) + c2_row = nlp.constraint_names().index(m.c2.name) + c3_row = nlp.constraint_names().index(m.c3.name) + x_col = nlp.primals_names().index(m.x.name) + y_col = nlp.primals_names().index(m.y.name) + z_col = nlp.primals_names().index(m.z.name) assert jac_scaled[c1_row, x_col] == pytest.approx(-1000) assert jac_scaled[c1_row, y_col] == pytest.approx(-1000) @@ -2526,7 +2551,7 @@ def test_scale_with_var_scale(self, model): assert jac_scaled[c3_row, z_col] == pytest.approx(3e6) @pytest.mark.unit - def test_exclude_scaling_factors_variables(self, model): + def test_exclude_scaling_factors_variables(self, model, include_greybox): m = model m.scaling_factor = Suffix(direction=Suffix.EXPORT) m.scaling_factor[m.x] = 1e-3 @@ -2537,17 +2562,18 @@ def test_exclude_scaling_factors_variables(self, model): m, include_scaling_factors=False, include_ipopt_autoscaling=True, + include_greybox=include_greybox, min_scale=1e-6, ) assert len(m.scaling_factor) == 3 assert not hasattr(m, "scaling_hint") - c1_row = nlp._condata_to_idx[m.c1] - c2_row = nlp._condata_to_idx[m.c2] - c3_row = nlp._condata_to_idx[m.c3] - x_col = nlp._vardata_to_idx[m.x] - y_col = nlp._vardata_to_idx[m.y] - z_col = nlp._vardata_to_idx[m.z] + c1_row = nlp.constraint_names().index(m.c1.name) + c2_row = nlp.constraint_names().index(m.c2.name) + c3_row = nlp.constraint_names().index(m.c3.name) + x_col = nlp.primals_names().index(m.x.name) + y_col = nlp.primals_names().index(m.y.name) + z_col = nlp.primals_names().index(m.z.name) assert jac_scaled[c1_row, x_col] == pytest.approx(-100) assert jac_scaled[c1_row, y_col] == pytest.approx(-0.1) @@ -2560,7 +2586,7 @@ def test_exclude_scaling_factors_variables(self, model): assert jac_scaled[c3_row, z_col] == pytest.approx(3e2) @pytest.mark.unit - def test_condition_number(self, model, caplog): + def test_condition_number(self, model, caplog, include_greybox): """Calculate the condition number of the Jacobian""" m = model m.scaling_factor = Suffix(direction=Suffix.EXPORT) @@ -2571,9 +2597,9 @@ def test_condition_number(self, model, caplog): m.scaling_factor[m.c2] = 1e-6 m.scaling_factor[m.c3] = 1e-12 - n = jacobian_cond(m, scaled=True) + n = jacobian_cond(m, scaled=True, include_greybox=include_greybox) assert n == pytest.approx(687.47, rel=1e-3) - n = jacobian_cond(m, scaled=False) + n = jacobian_cond(m, scaled=False, include_greybox=include_greybox) assert n == pytest.approx(7.50567e7, rel=1e-3) # Nonsquare condition number @@ -2581,7 +2607,7 @@ def test_condition_number(self, model, caplog): # Scaled with caplog.at_level(idaeslog.INFO): - n = jacobian_cond(m, scaled=True) + n = jacobian_cond(m, scaled=True, include_greybox=include_greybox) assert ( "Nonsquare Jacobian. Using pseudoinverse to calculate Frobenius norm." ) in caplog.text @@ -2589,14 +2615,14 @@ def test_condition_number(self, model, caplog): # Unscaled with caplog.at_level(idaeslog.INFO): - n = jacobian_cond(m, scaled=False) + n = jacobian_cond(m, scaled=False, include_greybox=include_greybox) assert ( "Nonsquare Jacobian. Using pseudoinverse to calculate Frobenius norm." ) in caplog.text assert n == pytest.approx(2.23741e5, rel=1e-3) @pytest.mark.unit - def test_condition_number_none(self, model): + def test_condition_number_none(self, model, include_greybox): with pytest.raises( RuntimeError, match=re.escape( @@ -2604,7 +2630,7 @@ def test_condition_number_none(self, model): "to calculate the condition number." ), ): - _ = jacobian_cond() + _ = jacobian_cond(include_greybox=include_greybox) def discretization_tester(transformation_method, scheme, t_skip, continuity_eqns=False): @@ -2741,3 +2767,176 @@ def diff_eqn(b, z, t): assert get_scaling_factor(m.xdot_disc_eq[0, i]) == approx(0.2) assert get_scaling_factor(m.xdot_disc_eq[1, i]) == approx(0.3) assert get_scaling_factor(m.xdot_disc_eq[2, i]) == approx(0.5) + + +class TestJacobianMethodsGreyBox: + @pytest.fixture + def model(self): + m = ConcreteModel() + + m.gb = ExternalGreyBoxBlock( + external_model=ResidualGrayBox(), build_implicit_constraint_objects=True + ) + + m.v1 = Var(initialize=5.0) + m.v2 = Var(initialize=5.0) + m.v3 = Var(initialize=-2.0) + m.v4 = Var(initialize=-4.0) + m.v5 = Var(initialize=2.0) + m.v6 = Var(initialize=0.0) + m.v7 = Var(initialize=1e-7) + + m.gb_link_v1 = Constraint(expr=m.gb.inputs["v1"] == m.v1) + m.gb_link_v2 = Constraint(expr=m.gb.inputs["v2"] == m.v2) + m.gb_link_v3 = Constraint(expr=m.gb.inputs["v3"] == m.v3) + m.gb_link_v4 = Constraint(expr=m.gb.inputs["v4"] == m.v4) + m.gb_link_v5 = Constraint(expr=m.gb.inputs["v5"] == m.v5) + m.gb_link_v6 = Constraint(expr=m.gb.inputs["v6"] == m.v6) + m.gb_link_v7 = Constraint(expr=m.gb.inputs["v7"] == m.v7) + + return m + + @pytest.mark.unit + def test_get_jacobian_include_greybox_false(self, model): + with pytest.raises( + ValueError, + match="The model contains components that are not supported by the Pyomo NL writer. " + "This may be because the model contains a grey-box model. If you want to include " + "grey-box variables and constraints in the Jacobian, set " + "include_greybox=True when calling get_jacobian.", + ): + get_jacobian(model, include_greybox=False) + + @pytest.mark.unit + def test_get_jacobian(self, model): + jac, nlp = get_jacobian(model, include_greybox=True) + + assert jac.shape == (11, 14) + + expected = { + ("gb_link_v1", "v1"): 1.0, + ("gb_link_v1", "gb.inputs[v1]"): -1.0, + ("gb_link_v2", "v2"): 1.0, + ("gb_link_v2", "gb.inputs[v2]"): -1.0, + ("gb_link_v3", "v3"): 1.0, + ("gb_link_v3", "gb.inputs[v3]"): -1.0, + ("gb_link_v4", "v4"): 1.0, + ("gb_link_v4", "gb.inputs[v4]"): -1.0, + ("gb_link_v5", "v5"): 1.0, + ("gb_link_v5", "gb.inputs[v5]"): -1.0, + ("gb_link_v6", "v6"): 1.0, + ("gb_link_v6", "gb.inputs[v6]"): -1.0, + ("gb_link_v7", "v7"): 1.0, + ("gb_link_v7", "gb.inputs[v7]"): -1.0, + ("gb.c1", "gb.inputs[v1]"): 1.0, + ("gb.c1", "gb.inputs[v2]"): 1.0, + ("gb.c2", "gb.inputs[v3]"): 1.0, + ("gb.c2", "gb.inputs[v4]"): -1.0, + ("gb.c2", "gb.inputs[v5]"): -1.0, + ("gb.c3", "gb.inputs[v3]"): 2.0, + ("gb.c3", "gb.inputs[v4]"): -3.0, + ("gb.c3", "gb.inputs[v5]"): -4.0, + ("gb.c3", "gb.inputs[v6]"): -1.0, + ("gb.c4", "gb.inputs[v1]"): -2e-08, + ("gb.c4", "gb.inputs[v7]"): 1.0, + } + + for i in range(11): + for j in range(14): + col = nlp.primals_names()[j] + row = nlp.constraint_names()[i] + assert jac[i, j] == pytest.approx( + expected.get((row, col), 0.0), rel=1e-6, abs=1e-12 + ) + + @pytest.mark.unit + def test_get_jacobian_w_scaling(self, model): + model.scaling_factor = Suffix(direction=Suffix.EXPORT) + model.scaling_factor[model.v7] = 1e7 + + jac, nlp = get_jacobian( + model, include_greybox=True, include_scaling_factors=True + ) + + assert jac.shape == (11, 14) + + expected = { + ("gb_link_v1", "v1"): 1.0, + ("gb_link_v1", "gb.inputs[v1]"): -1.0, + ("gb_link_v2", "v2"): 1.0, + ("gb_link_v2", "gb.inputs[v2]"): -1.0, + ("gb_link_v3", "v3"): 1.0, + ("gb_link_v3", "gb.inputs[v3]"): -1.0, + ("gb_link_v4", "v4"): 1.0, + ("gb_link_v4", "gb.inputs[v4]"): -1.0, + ("gb_link_v5", "v5"): 1.0, + ("gb_link_v5", "gb.inputs[v5]"): -1.0, + ("gb_link_v6", "v6"): 1.0, + ("gb_link_v6", "gb.inputs[v6]"): -1.0, + ("gb_link_v7", "v7"): 1e-07, + ("gb_link_v7", "gb.inputs[v7]"): -1.0, + ("gb.c1", "gb.inputs[v1]"): 1.0, + ("gb.c1", "gb.inputs[v2]"): 1.0, + ("gb.c2", "gb.inputs[v3]"): 1.0, + ("gb.c2", "gb.inputs[v4]"): -1.0, + ("gb.c2", "gb.inputs[v5]"): -1.0, + ("gb.c3", "gb.inputs[v3]"): 2.0, + ("gb.c3", "gb.inputs[v4]"): -3.0, + ("gb.c3", "gb.inputs[v5]"): -4.0, + ("gb.c3", "gb.inputs[v6]"): -1.0, + ("gb.c4", "gb.inputs[v1]"): -2e-08, + ("gb.c4", "gb.inputs[v7]"): 1.0, + } + + for i in range(11): + for j in range(14): + col = nlp.primals_names()[j] + row = nlp.constraint_names()[i] + print(row, col, jac[i, j]) + assert jac[i, j] == pytest.approx( + expected.get((row, col), 0.0), rel=1e-6, abs=1e-12 + ) + + @pytest.mark.unit + def test_get_jacobian_include_ipopt_autoscaling(self, model): + jac, nlp = get_jacobian( + model, include_greybox=True, include_ipopt_autoscaling=True + ) + + assert jac.shape == (11, 14) + + expected = { + ("gb_link_v1", "v1"): 1.0, + ("gb_link_v1", "gb.inputs[v1]"): -1.0, + ("gb_link_v2", "v2"): 1.0, + ("gb_link_v2", "gb.inputs[v2]"): -1.0, + ("gb_link_v3", "v3"): 1.0, + ("gb_link_v3", "gb.inputs[v3]"): -1.0, + ("gb_link_v4", "v4"): 1.0, + ("gb_link_v4", "gb.inputs[v4]"): -1.0, + ("gb_link_v5", "v5"): 1.0, + ("gb_link_v5", "gb.inputs[v5]"): -1.0, + ("gb_link_v6", "v6"): 1.0, + ("gb_link_v6", "gb.inputs[v6]"): -1.0, + ("gb_link_v7", "v7"): 1.0, + ("gb_link_v7", "gb.inputs[v7]"): -1.0, + ("gb.c1", "gb.inputs[v1]"): 1.0, + ("gb.c1", "gb.inputs[v2]"): 1.0, + ("gb.c2", "gb.inputs[v3]"): 1.0, + ("gb.c2", "gb.inputs[v4]"): -1.0, + ("gb.c2", "gb.inputs[v5]"): -1.0, + ("gb.c3", "gb.inputs[v3]"): 2.0, + ("gb.c3", "gb.inputs[v4]"): -3.0, + ("gb.c3", "gb.inputs[v5]"): -4.0, + ("gb.c3", "gb.inputs[v6]"): -1.0, + ("gb.c4", "gb.inputs[v1]"): -2e-08, + ("gb.c4", "gb.inputs[v7]"): 1.0, + } + + for i in range(11): + for j in range(14): + col = nlp.primals_names()[j] + row = nlp.constraint_names()[i] + assert jac[i, j] == pytest.approx( + expected.get((row, col), 0.0), rel=1e-6, abs=1e-12 + ) diff --git a/idaes/core/scaling/util.py b/idaes/core/scaling/util.py index e4aa98ebfd..57d36e71c5 100644 --- a/idaes/core/scaling/util.py +++ b/idaes/core/scaling/util.py @@ -64,6 +64,13 @@ from pyomo.dae.flatten import slice_component_along_sets from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP from pyomo.contrib.pynumero.asl import AmplInterface +from pyomo.contrib.pynumero.interfaces.pyomo_grey_box_nlp import ( + PyomoNLPWithGreyBoxBlocks, +) +from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlockData +from pyomo.contrib.pynumero.interfaces.external_grey_box_constraint import ( + ExternalGreyBoxConstraint, +) from idaes.core.util.exceptions import BurntToast import idaes.logger as idaeslog @@ -496,6 +503,10 @@ def get_scaling_factor(component, default: float = None, warning: bool = False): sfx_finder = SuffixFinder("scaling_factor") elif isinstance(component, ExpressionData): sfx_finder = SuffixFinder("scaling_hint") + elif isinstance(component, ExternalGreyBoxConstraint): + # External Grey Box Constraints cannot be given scaling factors, but we need them to behave as if + # they were Constraints, so we will return the default value here (or a value of 1.0 if no default) + return default if default is not None else 1.0 else: raise TypeError( f"Can get scaling factors for only VarData, ConstraintData, and (hints from) ExpressionData. " @@ -1137,6 +1148,7 @@ def get_jacobian( equality_constraints_only=False, max_grad=100, min_scale=1e-8, + include_greybox=False, ): """ Get the Jacobian matrix at the current model values. This function also @@ -1158,22 +1170,102 @@ def get_jacobian( min_scale: The smallest scaling factor that IPOPT gradient-based scaling is allowed to produce. This option corresponds to the IPOPT option nlp_scaling_min_value. Used only if include_ipopt_autoscaling is True. + include_greybox: if True, include grey-box variables and constraints in the Jacobian Returns: Jacobian matrix in Scipy CSR format, Pynumero nlp """ # Pynumero requires an objective, but I don't, so let's see if we have one - n_obj = 0 + has_obj = False for _ in m.component_data_objects(Objective, active=True): - n_obj += 1 + has_obj = True + break # Do not need to keep looking after we find one + # Add an objective if there isn't one - if n_obj == 0: + if not has_obj: dummy_objective_name = unique_component_name(m, "objective") setattr(m, dummy_objective_name, Objective(expr=0)) + # Create NLP and calculate the objective if not AmplInterface.available(): raise RuntimeError("Pynumero not available.") + + # Check to see if we should use legacy pathway or grey box pathway + if not include_greybox: + # Check to catch case where we have a grey box present but the + # user forgot to set include_greybox to True + try: + jac, nlp = _get_jacobian_legacy_pyomonlp( + m, + include_scaling_factors, + include_ipopt_autoscaling, + equality_constraints_only, + max_grad, + min_scale, + ) + except ValueError as err: + # If we get a ValueError check to see if it is because of unrecognized components + if "contains the following active components that the NL writer " in str( + err + ): + raise ValueError( + "The model contains components that are not supported by the Pyomo NL writer. " + "This may be because the model contains a grey-box model. If you want to include " + "grey-box variables and constraints in the Jacobian, set " + "include_greybox=True when calling get_jacobian." + ) from err + # Otherwise, raise the original error + raise err + else: + jac, nlp = _get_jacobian_greybox_compatible( + m, + include_scaling_factors, + include_ipopt_autoscaling, + equality_constraints_only, + max_grad, + min_scale, + ) + + # delete dummy objective + if not has_obj: + delattr(m, dummy_objective_name) + return jac, nlp + + +def _get_jacobian_legacy_pyomonlp( + m, + include_scaling_factors: bool = True, + include_ipopt_autoscaling: bool = False, + equality_constraints_only: bool = False, + max_grad: float = 100.0, + min_scale: float = 1e-8, +): + """ + Get the Jacobian matrix at the current model values. This function also + returns the Pynumero NLP which can be used to identify the constraints and + variables corresponding to the rows and columns. + + Args: + m: model to get Jacobian from + include_scaling_factors: if True scale the rows and columns of the Jacobian + using the user-defined scaling factors in the scaling_factor suffix. + include_ipopt_autoscaling: if True, include the gradient-based autoscaling + step that IPOPT uses to scale the rows of the Jacobian + equality_constraints_only: Option to include only equality constraints + in the calculated Jacobian. + max_grad: The value of the norm of the constraint gradient above which + gradient-based autoscaling will be performed by IPOPT. This option + corresponds to the IPOPT option nlp_scaling_max_gradient. Used + only if include_ipopt_autoscaling is True. + min_scale: The smallest scaling factor that IPOPT gradient-based scaling + is allowed to produce. This option corresponds to the IPOPT option + nlp_scaling_min_value. Used only if include_ipopt_autoscaling is True. + + Returns: + Jacobian matrix in Scipy CSR format, Pynumero nlp + """ nlp = PyomoNLP(m) + if equality_constraints_only: jac = nlp.evaluate_jacobian_eq().tocsr() else: @@ -1213,15 +1305,155 @@ def get_jacobian( jac = diags(grad_scaling_factors) @ jac - # delete dummy objective - if n_obj == 0: - delattr(m, dummy_objective_name) return jac, nlp +def _get_jacobian_greybox_compatible( + m, + include_scaling_factors: bool = True, + include_ipopt_autoscaling: bool = False, + equality_constraints_only: bool = False, + max_grad: float = 100.0, + min_scale: float = 1e-8, + eq_atol: float = 1e-12, + eq_rtol: float = 0.0, +): + """ + Grey-box compatible Jacobian utility using PyomoNLPWithGreyBoxBlocks only. + + Key behaviors vs legacy: + - Numeric Jacobian comes from PyomoNLPWithGreyBoxBlocks (supports ExternalGreyBoxBlock). + - Equality-only supported via per-row classification (prefer constraint equality property/type). + - Scaling: + * Do NOT apply scaling to ExternalGreyBoxConstraint rows. + * Do NOT apply scaling to variables inside ExternalGreyBoxBlock columns. + * Apply legacy suffix scaling to normal Pyomo vars/cons. + - Works for models without any grey-box blocks and should match legacy results + *except* for any differences in ordering introduced by PyomoNLPWithGreyBoxBlocks + (your tests will reveal whether ordering is identical). + + + Args: + m: model to get Jacobian from + include_scaling_factors: if True scale the rows and columns of the Jacobian + using the user-defined scaling factors in the scaling_factor suffix. + include_ipopt_autoscaling: if True, include the gradient-based autoscaling + step that IPOPT uses to scale the rows of the Jacobian + equality_constraints_only: Option to include only equality constraints + in the calculated Jacobian. Equality constraints are identified using a robust + classification that checks for an explicit equality property/type on the constraint + as well as checks the bounds with tight tolerances to account for numerical issues. + max_grad: The value of the norm of the constraint gradient above which + gradient-based autoscaling will be performed by IPOPT. This option + corresponds to the IPOPT option nlp_scaling_max_gradient. Used + only if include_ipopt_autoscaling is True. + min_scale: The smallest scaling factor that IPOPT gradient-based scaling + is allowed to produce. This option corresponds to the IPOPT option + nlp_scaling_min_value. Used only if include_ipopt_autoscaling is True. + eq_atol: Absolute tolerance for classifying equality constraints based on bounds. + eq_rtol: Relative tolerance for classifying equality constraints based on bounds. + """ + nlp = PyomoNLPWithGreyBoxBlocks(m) + + # Full Jacobian in NLP row/col order + jac_full = nlp.evaluate_jacobian().tocsr() + + # Lists in NLP order + vlist = nlp.vlist = list(nlp._pyomo_model_var_datas) + clist_all = list(nlp._constraint_datas) + + def _is_var_in_external_greybox_block(var): + """ + True if variable is declared under an ExternalGreyBoxBlock. + Conservative: check parent_block type name. + """ + try: + pb = var.parent_block() + except Exception: + return False + if pb is None: + return False + + if isinstance(pb, ExternalGreyBoxBlockData): + return True + return False + + def _is_equality_constraint(con, lb, ub): + """ + Classify equality constraints robustly. + - ExternalGreyBoxConstraint: always equality + - Otherwise, prefer con.equality if it exists and is boolean + - Otherwise use bounds test (finite & isclose) + """ + eq_attr = getattr(con, "equality", None) + if isinstance(eq_attr, (bool, np.bool_)): + return bool(eq_attr) + + finite = np.isfinite(lb) and np.isfinite(ub) + return bool(finite and np.isclose(lb, ub, atol=eq_atol, rtol=eq_rtol)) + + # --- Equality-only row selection --- + eq_rows = None + if equality_constraints_only: + lb_all = np.asarray(nlp.constraints_lb(), dtype=float) + ub_all = np.asarray(nlp.constraints_ub(), dtype=float) + + eq_mask = np.zeros(len(clist_all), dtype=bool) + for i, con in enumerate(clist_all): + eq_mask[i] = _is_equality_constraint(con, lb_all[i], ub_all[i]) + + eq_rows = np.where(eq_mask)[0] + jac = jac_full[eq_rows, :] + nlp.clist = [clist_all[i] for i in eq_rows] + else: + jac = jac_full + nlp.clist = clist_all + + # --- Scaling (suffix scaling for "normal" Pyomo parts; neutral for grey-box parts) --- + if include_scaling_factors: + var_sf = np.ones(len(vlist), dtype=float) + for j, var in enumerate(vlist): + if _is_var_in_external_greybox_block(var): + # Do NOT scale grey-box variables + var_sf[j] = 1.0 + else: + # Legacy behavior + var_sf[j] = get_scaling_factor(var, default=1, warning=False) + + con_sf_all = np.ones(len(clist_all), dtype=float) + for i, con in enumerate(clist_all): + if isinstance(con, tuple): + # Implicit Grey box constraints show up as tuples + # Do NOT scale grey-box constraints + con_sf_all[i] = 1.0 + else: + # ExternalGreyBoxConstraints will go down this path and be handled + # by get_scaling_factor + con_sf_all[i] = get_scaling_factor(con, default=1, warning=False) + + con_sf = con_sf_all if eq_rows is None else con_sf_all[eq_rows] + + # Guard against zeros / NaNs + var_sf = np.where((~np.isfinite(var_sf)) | (var_sf == 0.0), 1.0, var_sf) + con_sf = np.where((~np.isfinite(con_sf)) | (con_sf == 0.0), 1.0, con_sf) + + jac = diags(con_sf) @ (jac @ diags(1.0 / var_sf)) + + # --- IPOPT gradient-based autoscaling (row scaling) --- + if include_ipopt_autoscaling: + row_norms = spnorm(jac, ord=np.inf, axis=1) + grad_sf = np.ones(jac.shape[0], dtype=float) + for i in range(jac.shape[0]): + if row_norms[i] > max_grad: + grad_sf[i] = max(min_scale, max_grad / row_norms[i]) + jac = diags(grad_sf) @ jac + + return jac.tocsr(), nlp + + # TODO should we calculate the 2-norm condition number from the SVD # once #1566 is merged? -def jacobian_cond(m=None, scaled=True, jac=None): +def jacobian_cond(m=None, scaled=True, jac=None, include_greybox=False): """ Get the Frobenius condition number of the scaled or unscaled Jacobian matrix of a model. @@ -1230,6 +1462,7 @@ def jacobian_cond(m=None, scaled=True, jac=None): m: calculate the condition number of the Jacobian from this model. scaled: if True use scaled Jacobian, else use unscaled jac: (optional) previously calculated Jacobian + include_greybox: if True, include grey-box variables and constraints in the Jacobian and condition number calculation Returns: (float) Condition number @@ -1240,7 +1473,7 @@ def jacobian_cond(m=None, scaled=True, jac=None): "User must provide either a Pyomo model or a Jacobian " "to calculate the condition number." ) - jac, _ = get_jacobian(m, scaled) + jac, _ = get_jacobian(m, scaled, include_greybox=include_greybox) jac = jac.tocsc() if jac.shape[0] != jac.shape[1]: _log.info( diff --git a/idaes/core/util/diagnostics_tools/degeneracy_hunter.py b/idaes/core/util/diagnostics_tools/degeneracy_hunter.py index d255202fd3..71373a5f60 100644 --- a/idaes/core/util/diagnostics_tools/degeneracy_hunter.py +++ b/idaes/core/util/diagnostics_tools/degeneracy_hunter.py @@ -128,7 +128,7 @@ def __init__(self, model, **kwargs): ) if len(greybox_block_set(model)) != 0: raise NotImplementedError( - "Model contains Greybox models, which are not supported by Diagnostics toolbox at the moment" + "Model contains Greybox models, which are not supported by Degeneracy Hunter at the moment" ) self._model = model self.config = DHCONFIG(kwargs) diff --git a/idaes/core/util/diagnostics_tools/deprecated/tests/test_degeneracy_hunter_legacy.py b/idaes/core/util/diagnostics_tools/deprecated/tests/test_degeneracy_hunter_legacy.py index 9e4321e666..18035328a2 100644 --- a/idaes/core/util/diagnostics_tools/deprecated/tests/test_degeneracy_hunter_legacy.py +++ b/idaes/core/util/diagnostics_tools/deprecated/tests/test_degeneracy_hunter_legacy.py @@ -31,6 +31,8 @@ from idaes.core.util.diagnostics_tools.deprecated.degeneracy_hunter_legacy import ( DegeneracyHunter, ) + +# dummy_problem is used as a fixture the tests, but is not picked up by static analysis from idaes.core.util.diagnostics_tools.tests.utils import ( dummy_problem, ) diff --git a/idaes/core/util/diagnostics_tools/diagnostics_toolbox.py b/idaes/core/util/diagnostics_tools/diagnostics_toolbox.py index 18a1cd6121..d54dcc917a 100644 --- a/idaes/core/util/diagnostics_tools/diagnostics_toolbox.py +++ b/idaes/core/util/diagnostics_tools/diagnostics_toolbox.py @@ -44,13 +44,18 @@ from idaes.core.solvers.get_solver import get_solver from idaes.core.util.misc import compact_expression_to_string from idaes.core.util.model_statistics import ( - variables_in_activated_constraints_set, variables_not_in_activated_constraints_set, variables_with_none_value_in_activated_equalities_set, greybox_block_set, degrees_of_freedom, large_residuals_set, variables_near_bounds_set, + external_variables_set, + variables_fixed_to_zero_set, + variables_near_zero_set, + variables_with_none_value_set, + variables_violating_bounds_set, + variables_with_extreme_values_set, ) from idaes.core.scaling.util import ( get_jacobian, @@ -80,12 +85,6 @@ extreme_jacobian_columns, extreme_jacobian_rows, extreme_jacobian_entries, - var_in_block, - vars_fixed_to_zero, - vars_near_zero, - vars_violating_bounds, - vars_with_none_value, - vars_with_extreme_values, ) _log = idaeslog.getLogger(__name__) @@ -105,6 +104,7 @@ ), ) # TODO is a relative tolerance necessary if we're including scaling? +# AL: Yes, as you cannot guarantee a user will provide good scaling CONFIG.declare( "variable_bounds_relative_tolerance", ConfigValue( @@ -245,6 +245,26 @@ description="Feasibility tolerance for identifying infeasible constraints and bounds", ), ) +CONFIG.declare( + "include_grey_box_blocks", + ConfigValue( + default=True, + domain=bool, + description="Whether to include grey-box blocks in diagnostics checks. " + "Excluding Grey Box can reduce overhead for computing diagnostics, " + "but will result in apparent structural issues as the implicit " + "constraints will not be included in the analysis.", + ), +) +CONFIG.declare( + "apply_scaling", + ConfigValue( + default=True, + domain=bool, + description="Whether to apply scaling factors when evaluating diagnostics. " + " If True, any scaling factors present will be applied.", + ), +) @document_kwargs_from_configdict(CONFIG) @@ -300,13 +320,25 @@ def __init__(self, model: BlockData, **kwargs): "model argument must be an instance of a Pyomo BlockData object " "(either a scalar Block or an element of an indexed Block)." ) - if len(greybox_block_set(model)) != 0: - raise NotImplementedError( - "Model contains Greybox models, which are not supported by Diagnostics toolbox at the moment" - ) + self._model = model self.config = CONFIG(kwargs) + if len(greybox_block_set(model)) != 0: + if self.config.include_grey_box_blocks: + _log.warning( + "External Grey Box Models detected in model. If these were not constructed with " + "build_implicit_constraint_objects=True, then the implicit constraints " + "will not be detected by the Diagnostics Toolbox, which may lead to structural " + "issues being reported." + ) + else: + _log.warning( + "Grey Box Models detected in model and include_grey_box_blocks config option is set to False. " + "Grey Box blocks and their associated constraints will be ignored by the Diagnostics Toolbox, " + "which may lead to structural issues being reported." + ) + @property def model(self): """ @@ -326,17 +358,11 @@ def display_external_variables(self, stream=None): None """ - if stream is None: - stream = sys.stdout - - ext_vars = [] - for v in variables_in_activated_constraints_set(self._model): - if not var_in_block(v, self._model): - ext_vars.append(v.name) - write_report_section( stream=stream, - lines_list=ext_vars, + lines_list=external_variables_set( + self._model, include_greybox=self.config.include_grey_box_blocks + ), title="The following external variable(s) appear in constraints within the model:", header="=", footer="=", @@ -353,12 +379,11 @@ def display_unused_variables(self, stream=None): None """ - if stream is None: - stream = sys.stdout - write_report_section( stream=stream, - lines_list=variables_not_in_activated_constraints_set(self._model), + lines_list=variables_not_in_activated_constraints_set( + self._model, include_greybox=self.config.include_grey_box_blocks + ), title="The following variable(s) do not appear in any activated constraints within the model:", header="=", footer="=", @@ -375,12 +400,11 @@ def display_variables_fixed_to_zero(self, stream=None): None """ - if stream is None: - stream = sys.stdout - write_report_section( stream=stream, - lines_list=vars_fixed_to_zero(self._model), + lines_list=variables_fixed_to_zero_set( + self._model, include_greybox=self.config.include_grey_box_blocks + ), title="The following variable(s) are fixed to zero:", header="=", footer="=", @@ -398,16 +422,16 @@ def display_variables_at_or_outside_bounds(self, stream=None): None """ - if stream is None: - stream = sys.stdout - write_report_section( stream=stream, lines_list=[ f"{v.name} ({'fixed' if v.fixed else 'free'}): value={value(v)} bounds={v.bounds}" - for v in vars_violating_bounds( + for v in variables_violating_bounds_set( self._model, - tolerance=self.config.variable_bounds_violation_tolerance, + abs_tol=self.config.variable_bounds_violation_tolerance, + rel_tol=0.0, # For backward compatibility + include_greybox=self.config.include_grey_box_blocks, + apply_scaling=self.config.apply_scaling, ) ], title="The following variable(s) have values at or outside their bounds " @@ -427,12 +451,9 @@ def display_variables_with_none_value(self, stream=None): None """ - if stream is None: - stream = sys.stdout - write_report_section( stream=stream, - lines_list=vars_with_none_value(self._model), + lines_list=variables_with_none_value_set(self._model), title="The following variable(s) have a value of None:", header="=", footer="=", @@ -451,9 +472,6 @@ def display_variables_with_none_value_in_activated_constraints(self, stream=None None """ - if stream is None: - stream = sys.stdout - write_report_section( stream=stream, lines_list=[ @@ -467,7 +485,7 @@ def display_variables_with_none_value_in_activated_constraints(self, stream=None footer="=", ) - def _verify_active_variables_initialized(self, stream=None): + def _verify_active_variables_initialized(self): """ Validate that all variables are initialized (i.e., have values set to something other than None) before doing further numerical analysis. @@ -499,15 +517,15 @@ def display_variables_with_value_near_zero(self, stream=None): None """ - if stream is None: - stream = sys.stdout - write_report_section( stream=stream, lines_list=[ f"{v.name}: value={value(v)}" - for v in vars_near_zero( - self._model, self.config.variable_zero_value_tolerance + for v in variables_near_zero_set( + self._model, + tol=self.config.variable_zero_value_tolerance, + include_greybox=self.config.include_grey_box_blocks, + scale_variables=self.config.apply_scaling, ) ], title=f"The following variable(s) have a value close to zero " @@ -529,18 +547,17 @@ def display_variables_with_extreme_values(self, stream=None): None """ - if stream is None: - stream = sys.stdout - write_report_section( stream=stream, lines_list=[ f"{i.name}: {value(i)}" - for i in vars_with_extreme_values( - model=self._model, + for i in variables_with_extreme_values_set( + block=self._model, large=self.config.variable_large_value_tolerance, small=self.config.variable_small_value_tolerance, zero=self.config.variable_zero_value_tolerance, + include_greybox=self.config.include_grey_box_blocks, + apply_scaling=self.config.apply_scaling, ) ], title=f"The following variable(s) have extreme values " @@ -562,9 +579,6 @@ def display_variables_near_bounds(self, stream=None): None """ - if stream is None: - stream = sys.stdout - write_report_section( stream=stream, lines_list=[ @@ -573,6 +587,8 @@ def display_variables_near_bounds(self, stream=None): self._model, abs_tol=self.config.variable_bounds_absolute_tolerance, rel_tol=self.config.variable_bounds_relative_tolerance, + include_greybox=self.config.include_grey_box_blocks, + apply_scaling=self.config.apply_scaling, ) ], title=f"The following variable(s) have values close to their bounds " @@ -594,9 +610,6 @@ def display_components_with_inconsistent_units(self, stream=None): None """ - if stream is None: - stream = sys.stdout - write_report_section( stream=stream, lines_list=identify_inconsistent_units(self._model), @@ -619,13 +632,12 @@ def display_constraints_with_large_residuals(self, stream=None): None """ - if stream is None: - stream = sys.stdout - lrdict = large_residuals_set( self._model, tol=self.config.constraint_residual_tolerance, return_residual_values=True, + include_greybox=self.config.include_grey_box_blocks, + apply_scaling=self.config.apply_scaling, ) lrs = [] @@ -665,6 +677,12 @@ def compute_infeasibility_explanation(self, stream=None, solver=None, tee=False) None """ + if len(greybox_block_set(self._model)) != 0: + raise NotImplementedError( + "Model contains Greybox models, which are not supported by the " + "Infeasibility Explanation tools at the moment" + ) + if solver is None: solver = get_solver("ipopt_v2") if stream is None: @@ -797,12 +815,13 @@ def display_variables_with_extreme_jacobians(self, stream=None): None """ - self._verify_active_variables_initialized(stream=stream) - - if stream is None: - stream = sys.stdout + self._verify_active_variables_initialized() - jac, nlp = get_jacobian(self._model) + jac, nlp = get_jacobian( + self._model, + include_greybox=self.config.include_grey_box_blocks, + include_scaling_factors=self.config.apply_scaling, + ) xjc = extreme_jacobian_columns( jac=jac, @@ -836,12 +855,13 @@ def display_constraints_with_extreme_jacobians(self, stream=None): None """ - self._verify_active_variables_initialized(stream=stream) + self._verify_active_variables_initialized() - if stream is None: - stream = sys.stdout - - jac, nlp = get_jacobian(self._model) + jac, nlp = get_jacobian( + self._model, + include_greybox=self.config.include_grey_box_blocks, + include_scaling_factors=self.config.apply_scaling, + ) xjr = extreme_jacobian_rows( jac=jac, @@ -876,12 +896,13 @@ def display_extreme_jacobian_entries(self, stream=None): None """ - self._verify_active_variables_initialized(stream=stream) + self._verify_active_variables_initialized() - if stream is None: - stream = sys.stdout - - jac, nlp = get_jacobian(self._model, include_scaling_factors=True) + jac, nlp = get_jacobian( + self._model, + include_scaling_factors=self.config.apply_scaling, + include_greybox=self.config.include_grey_box_blocks, + ) xje = extreme_jacobian_entries( jac, nlp, @@ -912,10 +933,7 @@ def display_near_parallel_constraints(self, stream=None): None """ - self._verify_active_variables_initialized(stream=stream) - - if stream is None: - stream = sys.stdout + self._verify_active_variables_initialized() parallel = [ f"{i[0].name}, {i[1].name}" @@ -923,6 +941,7 @@ def display_near_parallel_constraints(self, stream=None): model=self._model, tolerance=self.config.parallel_component_tolerance, direction="row", + include_greybox=self.config.include_grey_box_blocks, ) ] @@ -946,10 +965,7 @@ def display_near_parallel_variables(self, stream=None): None """ - self._verify_active_variables_initialized(stream=stream) - - if stream is None: - stream = sys.stdout + self._verify_active_variables_initialized() parallel = [ f"{i[0].name}, {i[1].name}" @@ -957,6 +973,7 @@ def display_near_parallel_variables(self, stream=None): model=self._model, tolerance=self.config.parallel_component_tolerance, direction="column", + include_greybox=self.config.include_grey_box_blocks, ) ] @@ -986,6 +1003,8 @@ def _collect_constraint_mismatches(self, descend_into=True): List of strings summarising constraints with cancellations List of strings with constraint names where constraint contains no free variables """ + self._verify_active_variables_initialized() + walker = ConstraintTermAnalysisVisitor( term_mismatch_tolerance=self.config.constraint_term_mismatch_tolerance, term_cancellation_tolerance=self.config.constraint_term_cancellation_tolerance, @@ -1031,11 +1050,6 @@ def display_constraints_with_mismatched_terms(self, stream=None): None """ - self._verify_active_variables_initialized(stream=stream) - - if stream is None: - stream = sys.stdout - mismatch, _, _ = self._collect_constraint_mismatches() # Write the output @@ -1066,11 +1080,6 @@ def display_constraints_with_canceling_terms(self, stream=None): None """ - self._verify_active_variables_initialized(stream=stream) - - if stream is None: - stream = sys.stdout - _, canceling, _ = self._collect_constraint_mismatches() # Write the output @@ -1106,10 +1115,7 @@ def display_problematic_constraint_terms( None """ - self._verify_active_variables_initialized(stream=stream) - - if stream is None: - stream = sys.stdout + self._verify_active_variables_initialized() # Check that constraint is of correct type to give useful error message if not isinstance(constraint, ConstraintData): @@ -1125,7 +1131,11 @@ def display_problematic_constraint_terms( raise TypeError( f"{constraint.name} is not an instance of a Pyomo Constraint." ) - sf = get_scaling_factor(constraint, default=1, warning=False) + sf = ( + get_scaling_factor(constraint, default=1, warning=False) + if self.config.apply_scaling + else 1 + ) # Don't need to scale constraint_term_mismatch_tolerance and # constraint_term_cancellation_tolerance because they are both @@ -1205,11 +1215,6 @@ def display_constraints_with_no_free_variables(self, stream=None): # Although, in principle, this method doesn't require # all variables to be initialized, its current # implementation does. - self._verify_active_variables_initialized(stream=stream) - - if stream is None: - stream = sys.stdout - _, _, constant = self._collect_constraint_mismatches() # Write the output @@ -1295,7 +1300,9 @@ def _collect_structural_cautions(self): """ # Collect cautions cautions = [] - zero_vars = vars_fixed_to_zero(self._model) + zero_vars = variables_fixed_to_zero_set( + self._model, include_greybox=self.config.include_grey_box_blocks + ) if len(zero_vars) > 0: vstring = "variables" if len(zero_vars) == 1: @@ -1332,7 +1339,11 @@ def _collect_numerical_warnings( """ if jac is None or nlp is None: - jac, nlp = get_jacobian(self._model) + jac, nlp = get_jacobian( + self._model, + include_greybox=self.config.include_grey_box_blocks, + include_scaling_factors=self.config.apply_scaling, + ) warnings = [] next_steps = [] @@ -1355,8 +1366,12 @@ def _collect_numerical_warnings( next_steps.append(self.compute_infeasibility_explanation.__name__ + "()") # Variables outside bounds - violated_bounds = vars_violating_bounds( - self._model, tolerance=self.config.variable_bounds_violation_tolerance + violated_bounds = variables_violating_bounds_set( + self._model, + abs_tol=self.config.variable_bounds_violation_tolerance, + rel_tol=0.0, # For backward compatibility + include_greybox=self.config.include_grey_box_blocks, + apply_scaling=self.config.apply_scaling, ) if len(violated_bounds) > 0: cstring = "Variables" @@ -1446,7 +1461,11 @@ def _collect_numerical_cautions(self, jac=None, nlp=None): """ if jac is None or nlp is None: - jac, nlp = get_jacobian(self._model) + jac, nlp = get_jacobian( + self._model, + include_greybox=self.config.include_grey_box_blocks, + include_scaling_factors=self.config.apply_scaling, + ) cautions = [] @@ -1467,8 +1486,11 @@ def _collect_numerical_cautions(self, jac=None, nlp=None): ) # Variables near zero - near_zero = vars_near_zero( - self._model, self.config.variable_zero_value_tolerance + near_zero = variables_near_zero_set( + self._model, + tol=self.config.variable_zero_value_tolerance, + include_greybox=self.config.include_grey_box_blocks, + scale_variables=self.config.apply_scaling, ) if len(near_zero) > 0: cstring = "Variables" @@ -1480,11 +1502,13 @@ def _collect_numerical_cautions(self, jac=None, nlp=None): ) # Variables with extreme values - xval = vars_with_extreme_values( - model=self._model, + xval = variables_with_extreme_values_set( + block=self._model, large=self.config.variable_large_value_tolerance, small=self.config.variable_small_value_tolerance, zero=self.config.variable_zero_value_tolerance, + include_greybox=self.config.include_grey_box_blocks, + apply_scaling=self.config.apply_scaling, ) if len(xval) > 0: cstring = "Variables" @@ -1497,7 +1521,9 @@ def _collect_numerical_cautions(self, jac=None, nlp=None): ) # Variables with value None - none_value = vars_with_none_value(self._model) + none_value = variables_with_none_value_set( + self._model, include_greybox=self.config.include_grey_box_blocks + ) if len(none_value) > 0: cstring = "Variables" if len(none_value) == 1: @@ -1637,15 +1663,8 @@ def report_structural_issues(self, stream=None): None """ - if stream is None: - stream = sys.stdout - # Potential evaluation errors # TODO: High Index? - if len(greybox_block_set(self._model)) != 0: - raise NotImplementedError( - "Model contains Greybox models, which are not supported by Diagnostics toolbox at the moment" - ) stats = collect_model_statistics(self._model) warnings, next_steps = self._collect_structural_warnings() @@ -1689,11 +1708,13 @@ def report_numerical_issues(self, stream=None): None """ - self._verify_active_variables_initialized(stream=stream) + self._verify_active_variables_initialized() - if stream is None: - stream = sys.stdout - jac, nlp = get_jacobian(self._model) + jac, nlp = get_jacobian( + self._model, + include_greybox=self.config.include_grey_box_blocks, + include_scaling_factors=self.config.apply_scaling, + ) warnings, next_steps = self._collect_numerical_warnings(jac=jac, nlp=nlp) cautions = self._collect_numerical_cautions(jac=jac, nlp=nlp) @@ -1766,9 +1787,6 @@ def display_potential_evaluation_errors(self, stream=None): Returns: None """ - if stream is None: - stream = sys.stdout - warnings = self._collect_potential_eval_errors() write_report_section( stream=stream, diff --git a/idaes/core/util/diagnostics_tools/svd_toolbox.py b/idaes/core/util/diagnostics_tools/svd_toolbox.py index 792959b30d..f5c6b9d995 100644 --- a/idaes/core/util/diagnostics_tools/svd_toolbox.py +++ b/idaes/core/util/diagnostics_tools/svd_toolbox.py @@ -185,7 +185,7 @@ def __init__(self, model: BlockData, **kwargs): ) if len(greybox_block_set(model)) != 0: raise NotImplementedError( - "Model contains Greybox models, which are not supported by Diagnostics toolbox at the moment" + "Model contains Greybox models, which are not supported by the SVD Toolbox at the moment" ) self._model = model self.config = SVDCONFIG(kwargs) diff --git a/idaes/core/util/diagnostics_tools/tests/test_diagnostics_toolbox.py b/idaes/core/util/diagnostics_tools/tests/test_diagnostics_toolbox.py index eea4a72ad3..49812214ce 100644 --- a/idaes/core/util/diagnostics_tools/tests/test_diagnostics_toolbox.py +++ b/idaes/core/util/diagnostics_tools/tests/test_diagnostics_toolbox.py @@ -111,7 +111,7 @@ def model(self): return m @pytest.mark.component - def test_with_grey_box(self): + def test_with_grey_box(self, caplog): class BasicGrayBox(ExternalGreyBoxModel): def input_names(self): @@ -131,8 +131,22 @@ def evaluate_equality_constraints(self): m = ConcreteModel() m.gb = ExternalGreyBoxBlock(external_model=BasicGrayBox()) - with pytest.raises(NotImplementedError): - DiagnosticsToolbox(model=m) + + DiagnosticsToolbox(model=m) + assert ( + "External Grey Box Models detected in model. If these were not constructed with " + "build_implicit_constraint_objects=True, then the implicit constraints " + "will not be detected by the Diagnostics Toolbox, which may lead to structural " + "issues being reported." + ) in caplog.text + + # Test again with include_grey_box_blocks=False + DiagnosticsToolbox(model=m, include_grey_box_blocks=False) + assert ( + "Grey Box Models detected in model and include_grey_box_blocks config option is set to False. " + "Grey Box blocks and their associated constraints will be ignored by the Diagnostics Toolbox, " + "which may lead to structural issues being reported." + ) in caplog.text @pytest.mark.component def test_display_external_variables(self, model): @@ -497,7 +511,7 @@ def test_display_variables_with_extreme_jacobians(self): model.c2 = Constraint(expr=model.v1 == 1e-8 * model.v3) model.c3 = Constraint(expr=1e8 * model.v1 + 1e10 * model.v2 == 1e-6 * model.v3) - dt = DiagnosticsToolbox(model=model) + dt = DiagnosticsToolbox(model=model, include_grey_box_blocks=False) stream = StringIO() dt.display_variables_with_extreme_jacobians(stream) @@ -540,7 +554,7 @@ def test_display_constraints_with_extreme_jacobians(self): model.c2 = Constraint(expr=model.v1 == 1e-8 * model.v3) model.c3 = Constraint(expr=1e8 * model.v1 + 1e10 * model.v2 == 1e-6 * model.v3) - dt = DiagnosticsToolbox(model=model) + dt = DiagnosticsToolbox(model=model, include_grey_box_blocks=False) stream = StringIO() dt.display_constraints_with_extreme_jacobians(stream) @@ -580,7 +594,59 @@ def test_display_extreme_jacobian_entries(self): model.c2 = Constraint(expr=model.v1 == 1e-8 * model.v3) model.c3 = Constraint(expr=1e8 * model.v1 + 1e10 * model.v2 == 1e-6 * model.v3) - dt = DiagnosticsToolbox(model=model) + dt = DiagnosticsToolbox(model=model, include_grey_box_blocks=False) + + stream = StringIO() + dt.display_extreme_jacobian_entries(stream) + + expected = """==================================================================================== +The following constraint(s) and variable(s) are associated with extreme Jacobian +entries (<1.0E-04 or>1.0E+04): + + c3, v2: 1.000E+10 + c2, v3: 1.000E-08 + c3, v1: 1.000E+08 + c3, v3: 1.000E-06 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + # Test to make sure scaled Jacobian is used + set_scaling_factor(model.c3, 1e-8) + + stream = StringIO() + dt.display_extreme_jacobian_entries(stream) + + expected = """==================================================================================== +The following constraint(s) and variable(s) are associated with extreme Jacobian +entries (<1.0E-04 or>1.0E+04): + + c3, v3: 1.000E-14 + c2, v3: 1.000E-08 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + @pytest.mark.component + def test_display_extreme_jacobian_entries_w_include_grey_box(self): + # When adding grey box support, a bug was discovered in Pyomo that caused linear models + # to raise an Exception when using PyomoNLPWithGreyBoxBlocks. + # This test aims to confirm that remains fixed and that the include_grey_box_blocks option does + # not affect results for non-grey-box models. + model = ConcreteModel() + model.v1 = Var(initialize=1e-8) + model.v2 = Var(initialize=1) + model.v3 = Var(initialize=1) + + model.c1 = Constraint(expr=model.v1 == model.v2) + model.c2 = Constraint(expr=model.v1 == 1e-8 * model.v3) + model.c3 = Constraint(expr=1e8 * model.v1 + 1e10 * model.v2 == 1e-6 * model.v3) + + dt = DiagnosticsToolbox(model=model, include_grey_box_blocks=True) stream = StringIO() dt.display_extreme_jacobian_entries(stream) diff --git a/idaes/core/util/diagnostics_tools/tests/test_diagnostics_toolbox_greybox.py b/idaes/core/util/diagnostics_tools/tests/test_diagnostics_toolbox_greybox.py new file mode 100644 index 0000000000..49afeae404 --- /dev/null +++ b/idaes/core/util/diagnostics_tools/tests/test_diagnostics_toolbox_greybox.py @@ -0,0 +1,926 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2026 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +This module contains tests for the Diagnostics Toolbox using ExternalGreyBox models. +""" + +from io import StringIO +import logging +import re + +import pytest + +import numpy as np +import pytest +from scipy.sparse import coo_matrix + +from pyomo.contrib.pynumero.interfaces.external_grey_box import ( + ExternalGreyBoxBlock, + ExternalGreyBoxModel, +) +from pyomo.environ import ( + Block, + ConcreteModel, + Constraint, + SolverFactory, + units, + Var, +) + +from idaes.core.util.diagnostics_tools.diagnostics_toolbox import ( + DiagnosticsToolbox, +) + +logging.getLogger("cyipopt").setLevel(logging.WARNING) + + +# TODO: Pyomo NLP does not know how to handle Grey Box components +# TODO: Pyomo MIS does not handle Grey Box components + + +class OutputGrayBox(ExternalGreyBoxModel): + def input_names(self): + return ["v2", "v5", "v6"] + + def output_names(self): + return ["v1", "v3", "v4", "v7"] + + def set_input_values(self, input_values): + self._input_values = list(input_values) + + def finalize_block_construction(self, pyomo_block): + pyomo_block.outputs["v3"].setlb(0) # v3 + pyomo_block.outputs["v3"].setub(5) # v3 + pyomo_block.inputs["v5"].setlb(0) # v5 + pyomo_block.inputs["v5"].setub(5) # v5 + + def evaluate_outputs(self): + v1 = 10 - self._input_values[0] # c1 + # Linear combination of c2 and c3 to eliminate v3 and solve for v4 + v4 = -2 * self._input_values[1] - 2 * self._input_values[1] + v3 = v4 + self._input_values[1] # c2 + v7 = 2e-8 * v1 # c4 + return [v1, v3, v4, v7] + + def evaluate_jacobian_outputs(self): + row = np.zeros(6) # output index + col = np.zeros(6) # input index + data = np.zeros(6) # jacobian values + + # dv1/dv2 + row[0], col[0], data[0] = (0, 0, -1) + # dv4/dv5 + row[1], col[1], data[1] = (2, 1, -2) + # dv4/dv6 + row[2], col[2], data[2] = (2, 2, -1) + # dv3/dv5 + row[3], col[3], data[3] = (1, 1, -1) + # dv3/dv6 + row[4], col[4], data[4] = (1, 2, -1) + # dv7/dv2 + row[5], col[5], data[5] = (3, 0, -2e-8) + + return coo_matrix((data, (row, col)), shape=(4, 3)) + + +class ResidualGrayBox(ExternalGreyBoxModel): + def input_names(self): + return ["v1", "v2", "v3", "v4", "v5", "v6", "v7"] + + def output_names(self): + return [] + + def equality_constraint_names(self): + return ["c1", "c2", "c3", "c4"] + + def set_input_values(self, input_values): + self._input_values = list(input_values) + + def finalize_block_construction(self, pyomo_block): + pyomo_block.inputs["v3"].setlb(0) # v3 + pyomo_block.inputs["v3"].setub(5) # v3 + pyomo_block.inputs["v5"].setlb(0) # v5 + pyomo_block.inputs["v5"].setub(5) # v5 + + def evaluate_outputs(self): + raise NotImplementedError("This grey box only provides equality constraints") + + def evaluate_jacobian_outputs(self): + raise NotImplementedError("This grey box only provides equality constraints") + + def evaluate_equality_constraints(self): + # c1: m.v1 + m.b.v2 == 10 + # c2: m.b.v3 == m.b.v4 + m.b.v5 + # c3: 2 * m.b.v3 == 3 * m.b.v4 + 4 * m.b.v5 + m.b.v6 + # c4: m.b.v7 == 2e-8 * m.v1 + r = np.zeros(4) + r[0] = self._input_values[0] + self._input_values[1] - 10 + r[1] = self._input_values[2] - self._input_values[3] - self._input_values[4] + r[2] = ( + 2 * self._input_values[2] + - 3 * self._input_values[3] + - 4 * self._input_values[4] + - self._input_values[5] + ) + r[3] = self._input_values[6] - 2e-8 * self._input_values[0] + return r + + def evaluate_jacobian_equality_constraints(self): + row = np.zeros(11) # constraint index + col = np.zeros(11) # input index + data = np.zeros(11) # jacobian values + + # dc1/dv1 + row[0], col[0], data[0] = (0, 0, 1) + # dc1/dv2 + row[1], col[1], data[1] = (0, 1, 1) + # dc2/dv3 + row[2], col[2], data[2] = (1, 2, 1) + # dc2/dv4 + row[3], col[3], data[3] = (1, 3, -1) + # dc2/dv5 + row[4], col[4], data[4] = (1, 4, -1) + # dc3/dv3 + row[5], col[5], data[5] = (2, 2, 2) + # dc3/dv4 + row[6], col[6], data[6] = (2, 3, -3) + # dc3/dv5 + row[7], col[7], data[7] = (2, 4, -4) + # dc3/dv6 + row[8], col[8], data[8] = (2, 5, -1) + # dc4/dv1 + row[9], col[9], data[9] = (3, 0, -2e-8) + # dc4/dv7 + row[10], col[10], data[10] = (3, 6, 1) + + return coo_matrix((data, (row, col)), shape=(4, 7)) + + +@pytest.fixture() +def model(): + # Note: This model is deliberately infeasible and poorly scaled to + # test the diagnostics toolbox capabilities. + m = ConcreteModel() + m.b = Block() + + m.v1_1 = Var(units=units.m) # External variable + m.b.v2_1 = Var(units=units.m) + m.b.v3_1 = Var(bounds=(0, 5)) + m.b.v4_1 = Var() + m.b.v5_1 = Var(bounds=(0, 5)) + m.b.v6_1 = Var() + m.b.v7_1 = Var( + units=units.m, bounds=(0, 1) + ) # Poorly scaled variable with lower bound + + m.v1_2 = Var(units=units.m) # External variable + m.b.v2_2 = Var(units=units.m) + m.b.v3_2 = Var(bounds=(0, 5)) + m.b.v4_2 = Var() + m.b.v5_2 = Var(bounds=(0, 5)) + m.b.v6_2 = Var() + m.b.v7_2 = Var( + units=units.m, bounds=(0, 1) + ) # Poorly scaled variable with lower bound + + m.b.v8 = Var() # Unused variable + + m.b.v2_1.fix(5) + m.b.v5_1.fix(2) + m.b.v6_1.fix(0) + m.b.v2_2.fix(5) + m.b.v5_2.fix(2) + m.b.v6_2.fix(0) + + m.b.gb1 = ExternalGreyBoxBlock( + external_model=OutputGrayBox(), build_implicit_constraint_objects=True + ) + m.b.gb1_link_v1 = Constraint(expr=m.b.gb1.outputs["v1"] == m.v1_1) + m.b.gb1_link_v2 = Constraint(expr=m.b.gb1.inputs["v2"] == m.b.v2_1) + m.b.gb1_link_v3 = Constraint(expr=m.b.gb1.outputs["v3"] == m.b.v3_1) + m.b.gb1_link_v4 = Constraint(expr=m.b.gb1.outputs["v4"] == m.b.v4_1) + m.b.gb1_link_v5 = Constraint(expr=m.b.gb1.inputs["v5"] == m.b.v5_1) + m.b.gb1_link_v6 = Constraint(expr=m.b.gb1.inputs["v6"] == m.b.v6_1) + m.b.gb1_link_v7 = Constraint(expr=m.b.gb1.outputs["v7"] == m.b.v7_1) + + m.b.gb2 = ExternalGreyBoxBlock( + external_model=ResidualGrayBox(), + build_implicit_constraint_objects=True, + ) + m.b.gb2_link_v1 = Constraint(expr=m.b.gb2.inputs["v1"] == m.v1_2) + m.b.gb2_link_v2 = Constraint(expr=m.b.gb2.inputs["v2"] == m.b.v2_2) + m.b.gb2_link_v3 = Constraint(expr=m.b.gb2.inputs["v3"] == m.b.v3_2) + m.b.gb2_link_v4 = Constraint(expr=m.b.gb2.inputs["v4"] == m.b.v4_2) + m.b.gb2_link_v5 = Constraint(expr=m.b.gb2.inputs["v5"] == m.b.v5_2) + m.b.gb2_link_v6 = Constraint(expr=m.b.gb2.inputs["v6"] == m.b.v6_2) + m.b.gb2_link_v7 = Constraint(expr=m.b.gb2.inputs["v7"] == m.b.v7_2) + + # Solve model + solver = SolverFactory("cyipopt") + res = solver.solve(m) + + return m + + +@pytest.fixture() +def diagnostics_toolbox(model): + return DiagnosticsToolbox(model.b, include_grey_box_blocks=True) + + +@pytest.mark.unit +def test_display_external_variables(diagnostics_toolbox): + stream = StringIO() + + diagnostics_toolbox.display_external_variables(stream=stream) + + expected = """==================================================================================== +The following external variable(s) appear in constraints within the model: + + v1_1 + v1_2 + +==================================================================================== +""" + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_unused_variables(diagnostics_toolbox): + stream = StringIO() + + diagnostics_toolbox.display_unused_variables(stream=stream) + + expected = """==================================================================================== +The following variable(s) do not appear in any activated constraints within the model: + + b.v8 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_variables_fixed_to_zero(diagnostics_toolbox): + stream = StringIO() + + diagnostics_toolbox._model.gb1.outputs["v3"].fix(0) # Fix v3 to zero for testing + diagnostics_toolbox._model.gb2.inputs["v3"].fix(0) # Fix v3 to zero for testing + + diagnostics_toolbox.display_variables_fixed_to_zero(stream=stream) + + expected = """==================================================================================== +The following variable(s) are fixed to zero: + + b.v6_1 + b.v6_2 + b.gb1.outputs[v3] + b.gb2.inputs[v3] + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_variables_at_or_outside_bounds(diagnostics_toolbox): + stream = StringIO() + + diagnostics_toolbox.display_variables_at_or_outside_bounds(stream=stream) + + expected = """==================================================================================== +The following variable(s) have values at or outside their bounds (tol=0.0E+00): + + b.v3_1 (free): value=-9.997671315624657e-09 bounds=(0, 5) + b.v3_2 (free): value=-9.9943899712328e-09 bounds=(0, 5) + b.gb1.outputs[v3] (free): value=-9.998508737996996e-09 bounds=(0, 5) + b.gb2.inputs[v3] (free): value=-9.994692548467552e-09 bounds=(0, 5) + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_variables_with_none_value(diagnostics_toolbox): + stream = StringIO() + + # Set some variable values to None for testing + diagnostics_toolbox._model.gb1.outputs["v3"].set_value(None) + diagnostics_toolbox._model.gb2.inputs["v3"].set_value(None) + + diagnostics_toolbox.display_variables_with_none_value(stream=stream) + + expected = """==================================================================================== +The following variable(s) have a value of None: + + b.v8 + b.gb1.outputs[v3] + b.gb2.inputs[v3] + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_variables_with_none_value_in_activated_constraints( + diagnostics_toolbox, +): + stream = StringIO() + + # Set some variable values to None for testing + diagnostics_toolbox._model.v4_1.set_value(None) + diagnostics_toolbox._model.gb1.outputs["v3"].set_value(None) + diagnostics_toolbox._model.gb2.inputs["v3"].set_value(None) + + diagnostics_toolbox.display_variables_with_none_value_in_activated_constraints( + stream=stream + ) + + expected = """==================================================================================== +The following variable(s) have a value of None: + + b.gb1.outputs[v3] + b.v4_1 + b.gb2.inputs[v3] + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_variables_with_value_near_zero(diagnostics_toolbox): + stream = StringIO() + + diagnostics_toolbox.display_variables_with_value_near_zero(stream=stream) + + expected = """==================================================================================== +The following variable(s) have a value close to zero (tol=1.0E-08): + + b.v3_1: value=-9.997671315624657e-09 + b.v6_1: value=0 + b.v3_2: value=-9.9943899712328e-09 + b.v6_2: value=0 + b.gb1.outputs[v3]: value=-9.998508737996996e-09 + b.gb2.inputs[v3]: value=-9.994692548467552e-09 + b.gb2.inputs[v6]: value=-6.81818317855146e-13 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_variables_with_extreme_values(diagnostics_toolbox): + stream = StringIO() + + diagnostics_toolbox.display_variables_with_extreme_values(stream=stream) + + expected = """==================================================================================== +The following variable(s) have extreme values (<1.0E-04 or > 1.0E+04): + + b.v7_1: 1.0000000003005368e-07 + b.v7_2: 1.0000000003005368e-07 + b.gb1.outputs[v7]: 1.0000000001502684e-07 + b.gb2.inputs[v7]: 1.0000000001502684e-07 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_variables_near_bounds(diagnostics_toolbox): + stream = StringIO() + + diagnostics_toolbox.display_variables_near_bounds(stream=stream) + + expected = """==================================================================================== +The following variable(s) have values close to their bounds (abs=1.0E-04, rel=1.0E-04): + + b.v3_1: value=-9.997671315624657e-09 bounds=(0, 5) + b.v7_1: value=1.0000000003005368e-07 bounds=(0, 1) + b.v3_2: value=-9.9943899712328e-09 bounds=(0, 5) + b.v7_2: value=1.0000000003005368e-07 bounds=(0, 1) + b.gb1.outputs[v3]: value=-9.998508737996996e-09 bounds=(0, 5) + b.gb2.inputs[v3]: value=-9.994692548467552e-09 bounds=(0, 5) + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_components_with_inconsistent_units(diagnostics_toolbox): + stream = StringIO() + + diagnostics_toolbox.display_components_with_inconsistent_units(stream=stream) + + expected = """==================================================================================== +The following component(s) have unit consistency issues: + + b.gb1_link_v1 + b.gb1_link_v2 + b.gb1_link_v7 + b.gb2_link_v1 + b.gb2_link_v2 + b.gb2_link_v7 + +For more details on unit inconsistencies, import the assert_units_consistent method +from pyomo.util.check_units +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_constraints_with_large_residuals(diagnostics_toolbox): + stream = StringIO() + + diagnostics_toolbox.display_constraints_with_large_residuals(stream=stream) + + expected = """==================================================================================== +The following constraint(s) have large residuals (>1.0E-05): + + b.gb1_link_v5: 1.99380E+00 + b.gb1_link_v6: 1.11198E-02 + b.gb1.v3_constraint: 1.86107E-02 + b.gb2.c2: 6.66667E-01 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_underconstrained_set(diagnostics_toolbox): + stream = StringIO() + + # Create structural singularities + blk = diagnostics_toolbox._model + blk.v2_1.unfix() + blk.v2_2.unfix() + + diagnostics_toolbox.display_underconstrained_set(stream=stream) + + expected = """==================================================================================== +Dulmage-Mendelsohn Under-Constrained Set + + Independent Block 0: + + Variables: + + b.gb1.inputs[v2] + b.gb1.outputs[v1] + b.v2_1 + b.gb1.outputs[v7] + v1_1 + b.v7_1 + + Constraints: + + b.gb1.v1_constraint + b.gb1_link_v2 + b.gb1.v7_constraint + b.gb1_link_v1 + b.gb1_link_v7 + + Independent Block 1: + + Variables: + + b.gb2.inputs[v2] + b.gb2.inputs[v1] + b.v2_2 + v1_2 + b.gb2.inputs[v7] + b.v7_2 + + Constraints: + + b.gb2.c1 + b.gb2_link_v2 + b.gb2_link_v1 + b.gb2.c4 + b.gb2_link_v7 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.unit +def test_display_overconstrained_set(diagnostics_toolbox): + stream = StringIO() + + # Create structural singularities + blk = diagnostics_toolbox._model + blk.v4_1.fix() + blk.v4_2.fix() + + diagnostics_toolbox.display_overconstrained_set(stream=stream) + + expected = """==================================================================================== +Dulmage-Mendelsohn Over-Constrained Set + + Independent Block 0: + + Variables: + + b.gb1.outputs[v4] + b.gb1.inputs[v5] + b.gb1.inputs[v6] + + Constraints: + + b.gb1_link_v4 + b.gb1_link_v5 + b.gb1_link_v6 + b.gb1.v4_constraint + + Independent Block 1: + + Variables: + + b.gb2.inputs[v4] + b.gb2.inputs[v5] + b.gb2.inputs[v6] + b.gb2.inputs[v3] + + Constraints: + + b.gb2_link_v4 + b.gb2_link_v5 + b.gb2_link_v6 + b.gb2.c2 + b.gb2.c3 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.component +def test_display_near_parallel_constraints(diagnostics_toolbox): + m = diagnostics_toolbox.model + + # Add parallel constraints + m.parallel_con1 = Constraint(expr=m.gb1.outputs["v1"] == 10 - m.gb1.inputs["v2"]) + m.parallel_con2 = Constraint(expr=m.gb2.inputs["v1"] == 10 - m.gb2.inputs["v2"]) + + stream = StringIO() + + diagnostics_toolbox.display_near_parallel_constraints(stream=stream) + + expected = """==================================================================================== +The following pairs of constraints are nearly parallel: + + b.parallel_con1, b.gb1.v1_constraint + b.parallel_con2, b.gb2.c1 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.component +def test_display_near_parallel_variables(diagnostics_toolbox): + m = diagnostics_toolbox.model + + # Add parallel vars + m.pv1 = Var(initialize=m.gb1.outputs["v1"].value) # Parallel to gb1.outputs[v1] + m.pv2 = Var(initialize=m.gb1.inputs["v2"].value) # Parallel to gb1.inputs[v2] + m.parallel_con1 = Constraint(expr=m.pv1 == 10 - m.pv2) + + stream = StringIO() + + diagnostics_toolbox.display_near_parallel_variables(stream=stream) + + expected = """==================================================================================== +The following pairs of variables are nearly parallel: + + b.pv1, b.pv2 + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +# ==================================================================================== +# Reporting methods +@pytest.mark.component +def test_collect_structural_warnings(diagnostics_toolbox): + # Create structural singularities + blk = diagnostics_toolbox._model + blk.v4_1.fix() + blk.v4_2.fix() + + warnings, next_steps = diagnostics_toolbox._collect_structural_warnings() + + assert len(warnings) == 3 + assert "WARNING: -2 Degrees of Freedom" in warnings + assert "WARNING: 6 Components with inconsistent units" in warnings + assert """WARNING: Structural singularity found + Under-Constrained Set: 0 variables, 0 constraints + Over-Constrained Set: 7 variables, 9 constraints""" in warnings + + assert len(next_steps) == 2 + assert "display_components_with_inconsistent_units()" in next_steps + assert "display_overconstrained_set()" in next_steps + + +@pytest.mark.component +def test_collect_structural_cautions(diagnostics_toolbox): + # Create structural singularities + blk = diagnostics_toolbox._model + blk.v4_1.fix() + blk.v4_2.fix() + + cautions = diagnostics_toolbox._collect_structural_cautions() + + assert len(cautions) == 2 + assert "Caution: 2 variables fixed to 0" in cautions + assert "Caution: 1 unused variable (0 fixed)" in cautions + + +@pytest.mark.component +def test_collect_numerical_warnings(diagnostics_toolbox): + warnings, next_steps = diagnostics_toolbox._collect_numerical_warnings() + + assert len(warnings) == 2 + assert "WARNING: 4 Constraints with large residuals (>1.0E-05)" in warnings + assert "WARNING: 4 Variables at or outside bounds (tol=0.0E+00)" in warnings + + assert len(next_steps) == 3 + assert "display_constraints_with_large_residuals()" in next_steps + assert "compute_infeasibility_explanation()" in next_steps + assert "display_variables_at_or_outside_bounds()" in next_steps + + +@pytest.mark.component +def test_collect_numerical_cautions(diagnostics_toolbox): + cautions = diagnostics_toolbox._collect_numerical_cautions() + + for i in cautions: + print(i) + + assert len(cautions) == 5 + assert ( + "Caution: 6 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04)" + in cautions + ) + assert "Caution: 7 Variables with value close to zero (tol=1.0E-08)" in cautions + assert "Caution: 4 Variables with extreme value (<1.0E-04 or >1.0E+04)" in cautions + assert "Caution: 1 Variable with None value" in cautions + assert "Caution: 2 extreme Jacobian Entries (<1.0E-04 or >1.0E+04)" in cautions + + +@pytest.mark.component +def test_assert_no_structural_warnings(diagnostics_toolbox): + # Create structural singularities + blk = diagnostics_toolbox._model + blk.v4_1.fix() + blk.v4_2.fix() + + with pytest.raises(AssertionError, match=re.escape("Structural issues found (3).")): + diagnostics_toolbox.assert_no_structural_warnings() + + +@pytest.mark.component +def test_assert_no_numerical_warnings(diagnostics_toolbox): + with pytest.raises(AssertionError, match=re.escape("Numerical issues found (2).")): + diagnostics_toolbox.assert_no_numerical_warnings() + + +# ==================================================================================== +# Reporting methods +@pytest.mark.component +def test_report_structural_issues(diagnostics_toolbox): + # Create structural singularities + blk = diagnostics_toolbox._model + blk.v4_1.fix() + blk.v4_2.fix() + + stream = StringIO() + diagnostics_toolbox.report_structural_issues(stream=stream) + + expected = """==================================================================================== +Model Statistics + + Activated Blocks: 1 (Deactivated: 0) + Free Variables in Activated Constraints: 20 (External: 2) + Free Variables with only lower bounds: 0 + Free Variables with only upper bounds: 0 + Free Variables with upper and lower bounds: 8 + Fixed Variables in Activated Constraints: 8 (External: 0) + Activated Equality Constraints: 22 (Deactivated: 0) + Activated Inequality Constraints: 0 (Deactivated: 0) + Activated Objectives: 0 (Deactivated: 0) + GreyBox Statistics + Activated GreyBox models: 2 (Deactivated: 0) + Activated GreyBox Equalities: 8 (Deactivated: 0) + Free Variables in Activated GreyBox Equalities: 14 (Fixed: 0) + +------------------------------------------------------------------------------------ +3 WARNINGS + + WARNING: -2 Degrees of Freedom + WARNING: 6 Components with inconsistent units + WARNING: Structural singularity found + Under-Constrained Set: 0 variables, 0 constraints + Over-Constrained Set: 7 variables, 9 constraints + +------------------------------------------------------------------------------------ +2 Cautions + + Caution: 2 variables fixed to 0 + Caution: 1 unused variable (0 fixed) + +------------------------------------------------------------------------------------ +Suggested next steps: + + display_components_with_inconsistent_units() + display_overconstrained_set() + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.component +def test_report_numerical_issues(diagnostics_toolbox): + stream = StringIO() + diagnostics_toolbox.report_numerical_issues(stream=stream) + + print(stream.getvalue()) + expected = """==================================================================================== +Model Statistics + + Jacobian Condition Number: 7.849E+01 + +------------------------------------------------------------------------------------ +2 WARNINGS + + WARNING: 4 Constraints with large residuals (>1.0E-05) + WARNING: 4 Variables at or outside bounds (tol=0.0E+00) + +------------------------------------------------------------------------------------ +5 Cautions + + Caution: 6 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04) + Caution: 7 Variables with value close to zero (tol=1.0E-08) + Caution: 4 Variables with extreme value (<1.0E-04 or >1.0E+04) + Caution: 1 Variable with None value + Caution: 2 extreme Jacobian Entries (<1.0E-04 or >1.0E+04) + +------------------------------------------------------------------------------------ +Suggested next steps: + + display_constraints_with_large_residuals() + compute_infeasibility_explanation() + display_variables_at_or_outside_bounds() + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +# ==================================================================================== +# Functionality tests +# These tests are for methods that can/do not interact with Grey Box components, +# but we want to ensure work when Grey Box components are present in the model. +@pytest.mark.unit +def test_collect_constraint_mismatches(diagnostics_toolbox): + blk = diagnostics_toolbox._model + + blk.b3 = Block() + blk.b3.v1 = Var(initialize=2) + blk.b3.v2 = Var(initialize=3) + + # Constraint with no free variables + blk.b3.c1 = Constraint(expr=blk.b3.v1 == blk.b3.v2) + blk.b3.v1.fix() + blk.b3.v2.fix() + + # Constraint with mismatched terms + blk.b3.v3 = Var(initialize=10) + blk.b3.v4 = Var(initialize=10) + blk.b3.v5 = Var(initialize=1e-6) + blk.b3.c2 = Constraint(expr=blk.b3.v3 == blk.b3.v4 + blk.b3.v5) + + # Constraint with cancellation + blk.b3.v6 = Var(initialize=10) + blk.b3.c3 = Constraint(expr=blk.b3.v6 == 10 + blk.b3.v3 - blk.b3.v4) + + mismatch, canceling, constant = diagnostics_toolbox._collect_constraint_mismatches() + + print(mismatch) + print(canceling) + print(constant) + + assert mismatch == ["b.b3.c2: 1 mismatched term(s)"] + assert canceling == [ + "b.b3.c2: 1 potential canceling term(s)", + "b.b3.c3: 1 potential canceling term(s)", + ] + assert constant == ["b.b3.c1"] + + +@pytest.mark.component +def test_collect_potential_eval_errors(diagnostics_toolbox): + blk = diagnostics_toolbox._model + + blk.b3 = Block() + blk.b3.v1 = Var(initialize=2) + blk.b3.v2 = Var(initialize=0) + + # Constraint with potential division by zero + blk.b3.c1 = Constraint(expr=blk.b3.v1 == 10 / blk.b3.v2) + + errors = diagnostics_toolbox._collect_potential_eval_errors() + + assert len(errors) == 1 + assert "b.b3.c1: Potential division by 0 in 10/b.b3.v2" in errors[0] + + +@pytest.mark.component +def test_display_potential_evaluation_errors(diagnostics_toolbox): + blk = diagnostics_toolbox._model + + blk.b3 = Block() + blk.b3.v1 = Var(initialize=2) + blk.b3.v2 = Var(initialize=0) + + # Constraint with potential division by zero + blk.b3.c1 = Constraint(expr=blk.b3.v1 == 10 / blk.b3.v2) + + stream = StringIO() + diagnostics_toolbox.display_potential_evaluation_errors(stream=stream) + + expected = """==================================================================================== +1 WARNINGS + + b.b3.c1: Potential division by 0 in 10/b.b3.v2; Denominator bounds are (-inf, inf) + +==================================================================================== +""" + + assert stream.getvalue() == expected + + +@pytest.mark.component +def test_compute_infeasibility_explanation(diagnostics_toolbox): + with pytest.raises( + NotImplementedError, + match="Model contains Greybox models, which are not supported by the " + "Infeasibility Explanation tools at the moment", + ): + diagnostics_toolbox.compute_infeasibility_explanation() + + +@pytest.mark.component +def test_prepare_svd_toolbox(diagnostics_toolbox): + with pytest.raises( + NotImplementedError, + match="Model contains Greybox models, which are not supported by the " + "SVD Toolbox at the moment", + ): + diagnostics_toolbox.prepare_svd_toolbox() + + +@pytest.mark.component +def test_prepare_degeneracy_hunter(diagnostics_toolbox): + with pytest.raises( + NotImplementedError, + match="Model contains Greybox models, which are not supported by " + "Degeneracy Hunter at the moment", + ): + diagnostics_toolbox.prepare_degeneracy_hunter() diff --git a/idaes/core/util/diagnostics_tools/tests/test_svd_toolbox.py b/idaes/core/util/diagnostics_tools/tests/test_svd_toolbox.py index 00370b6c4b..a61b3d5ff9 100644 --- a/idaes/core/util/diagnostics_tools/tests/test_svd_toolbox.py +++ b/idaes/core/util/diagnostics_tools/tests/test_svd_toolbox.py @@ -15,10 +15,6 @@ """ from io import StringIO -from pyomo.contrib.pynumero.interfaces.external_grey_box import ( - ExternalGreyBoxBlock, - ExternalGreyBoxModel, -) import re import numpy as np @@ -32,12 +28,18 @@ ) from pyomo.contrib.pynumero.asl import AmplInterface from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP +from pyomo.contrib.pynumero.interfaces.external_grey_box import ( + ExternalGreyBoxBlock, + ExternalGreyBoxModel, +) from idaes.core.util.diagnostics_tools.svd_toolbox import ( SVDToolbox, svd_dense, svd_sparse, ) + +# dummy_problem is used as a fixture the tests, but is not picked up by static analysis from idaes.core.util.diagnostics_tools.tests.utils import ( dummy_problem, ) diff --git a/idaes/core/util/diagnostics_tools/tests/test_utils.py b/idaes/core/util/diagnostics_tools/tests/test_utils.py index f0672390ef..65b5f6f42a 100644 --- a/idaes/core/util/diagnostics_tools/tests/test_utils.py +++ b/idaes/core/util/diagnostics_tools/tests/test_utils.py @@ -25,15 +25,8 @@ Set, Var, ) -from pyomo.common.collections import ComponentSet from idaes.core.util.diagnostics_tools.utils import ( - var_in_block, - vars_fixed_to_zero, - vars_near_zero, - vars_violating_bounds, - vars_with_none_value, - vars_with_extreme_values, check_parallel_jacobian, extreme_jacobian_rows, extreme_jacobian_columns, @@ -43,213 +36,6 @@ from idaes.core.scaling.util import get_jacobian -@pytest.fixture -def model(): - m = ConcreteModel() - m.b = Block() - - m.v1 = Var() - m.v2 = Var() - m.v3 = Var() - m.v4 = Var() - - m.v1.fix(0) - m.v2.fix(3) - m.v3.set_value(0) - - return m - - -@pytest.mark.unit -def test_var_in_block(model): - assert var_in_block(model.v1, model) - assert not var_in_block(model.v1, model.b) - - -@pytest.mark.unit -def test_vars_fixed_to_zero(model): - zero_vars = vars_fixed_to_zero(model) - assert isinstance(zero_vars, ComponentSet) - assert len(zero_vars) == 1 - for i in zero_vars: - assert i is model.v1 - - -@pytest.mark.unit -def test_vars_near_zero(model): - model.v3.set_value(1e-5) - - near_zero_vars = vars_near_zero(model, variable_zero_value_tolerance=1e-5) - assert isinstance(near_zero_vars, ComponentSet) - assert len(near_zero_vars) == 2 - for i in near_zero_vars: - assert i.local_name in ["v1", "v3"] - - near_zero_vars = vars_near_zero(model, variable_zero_value_tolerance=1e-6) - assert isinstance(near_zero_vars, ComponentSet) - assert len(near_zero_vars) == 1 - for i in near_zero_vars: - assert i is model.v1 - - set_scaling_factor(model.v3, 1e5) - near_zero_vars = vars_near_zero(model, variable_zero_value_tolerance=1e-5) - assert isinstance(near_zero_vars, ComponentSet) - assert len(near_zero_vars) == 1 - for i in near_zero_vars: - assert i is model.v1 - - near_zero_vars = vars_near_zero(model, variable_zero_value_tolerance=1) - assert isinstance(near_zero_vars, ComponentSet) - assert len(near_zero_vars) == 2 - for i in near_zero_vars: - assert i.local_name in ["v1", "v3"] - - -@pytest.mark.unit -def test_vars_with_none_value(model): - none_value = vars_with_none_value(model) - - assert isinstance(none_value, ComponentSet) - assert len(none_value) == 1 - for i in none_value: - assert i is model.v4 - - -@pytest.mark.unit -def test_vars_with_bounds_issues(model): - model.v1.setlb(2) - model.v1.setub(6) - model.v2.setlb(0) - model.v2.setub(10) - model.v4.set_value(10) - model.v4.setlb(0) - model.v4.setub(1) - - bounds_issue = vars_violating_bounds(model, tolerance=0) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 2 - for i in bounds_issue: - assert i.local_name in ["v1", "v4"] - - m = ConcreteModel() - m.v = Var(initialize=-1e-4, bounds=(0, 1)) - - # Just outside lower bound - bounds_issue = vars_violating_bounds(m, tolerance=1e-6) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 1 - - bounds_issue = vars_violating_bounds(m, tolerance=1e3) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - set_scaling_factor(m.v, 1e-10, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e3) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - set_scaling_factor(m.v, 1e10, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e-6) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 1 - - # Just inside lower bound - m.v.set_value(1e-4) - - set_scaling_factor(m.v, 1, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e-6) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - bounds_issue = vars_violating_bounds(m, tolerance=1e3) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - set_scaling_factor(m.v, 1e-10, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e3) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - set_scaling_factor(m.v, 1e10, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e-6) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - # Just inside upper bound - m.v.set_value(1 - 1e-4) - - set_scaling_factor(m.v, 1, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e-6) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - bounds_issue = vars_violating_bounds(m, tolerance=1e3) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - set_scaling_factor(m.v, 1e-10, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e3) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - set_scaling_factor(m.v, 1e10, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e-6) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - # Just outside upper bound - m.v.set_value(1 + 1e-4) - - set_scaling_factor(m.v, 1, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e-6) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 1 - - bounds_issue = vars_violating_bounds(m, tolerance=1e3) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - set_scaling_factor(m.v, 1e-10, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e3) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 0 - - set_scaling_factor(m.v, 1e10, overwrite=True) - bounds_issue = vars_violating_bounds(m, tolerance=1e-6) - assert isinstance(bounds_issue, ComponentSet) - assert len(bounds_issue) == 1 - - -@pytest.mark.unit -def test_vars_with_extreme_values(): - m = ConcreteModel() - m.v1 = Var(initialize=1e-12) # below zero - m.v2 = Var(initialize=1e-8) # small - m.v3 = Var(initialize=1e-4) - m.v4 = Var(initialize=1e0) - m.v5 = Var(initialize=1e4) - m.v6 = Var(initialize=1e8) - m.v7 = Var(initialize=1e12) # large - - xvars = vars_with_extreme_values(m, large=1e9, small=1e-7, zero=1e-10) - - assert len(xvars) == 2 - for i in xvars: - assert i.name in ["v2", "v7"] - - set_scaling_factor(m.v1, 1e3) - set_scaling_factor(m.v2, 1e6) - set_scaling_factor(m.v4, 1e-12) - set_scaling_factor(m.v6, 1e3) - set_scaling_factor(m.v7, 1e-12) - - xvars = vars_with_extreme_values(m, large=1e9, small=1e-7, zero=1e-10) - - assert len(xvars) == 2 - for i in xvars: - assert i.name in ["v1", "v6"] - - class TestExtremeJacobianMethods: @pytest.fixture def model(self): diff --git a/idaes/core/util/diagnostics_tools/tests/test_writer_utils.py b/idaes/core/util/diagnostics_tools/tests/test_writer_utils.py index 5df78ff2bd..54644291a6 100644 --- a/idaes/core/util/diagnostics_tools/tests/test_writer_utils.py +++ b/idaes/core/util/diagnostics_tools/tests/test_writer_utils.py @@ -229,8 +229,12 @@ def evaluate_equality_constraints(self): m = ConcreteModel() - m.gb = ExternalGreyBoxBlock(external_model=BasicGrayBox()) - m.gb_inactive = ExternalGreyBoxBlock(external_model=BasicGrayBox()) + m.gb = ExternalGreyBoxBlock( + external_model=BasicGrayBox(), build_implicit_constraint_objects=True + ) + m.gb_inactive = ExternalGreyBoxBlock( + external_model=BasicGrayBox(), build_implicit_constraint_objects=True + ) m.gb_inactive.deactivate() m.a1 = Var(initialize=1) m.a1.fix() @@ -240,7 +244,7 @@ def evaluate_equality_constraints(self): m.o1 = Var(initialize=1) m.o1_eq = Constraint(expr=m.o1 == m.gb.outputs["o1"]) m.o1.fix() - stats = collect_model_statistics(m) + stats = collect_model_statistics(m, include_greybox=True) for k in stats: print(k) print(stats, len(stats)) diff --git a/idaes/core/util/diagnostics_tools/utils.py b/idaes/core/util/diagnostics_tools/utils.py index 4662530da5..3e59aa91a7 100644 --- a/idaes/core/util/diagnostics_tools/utils.py +++ b/idaes/core/util/diagnostics_tools/utils.py @@ -38,6 +38,7 @@ def check_parallel_jacobian( direction: str = "row", jac=None, nlp=None, + include_greybox: bool = False, ): """ Check for near-parallel rows or columns in the Jacobian. @@ -63,6 +64,7 @@ def check_parallel_jacobian( direction: 'row' (default, constraints) or 'column' (variables) jac: model Jacobian as a ``scipy.sparse.coo_matrix``, optional nlp: ``PyomoNLP`` of model, optional + include_greybox: whether to include grey-box components in the Jacobian Returns: list of 2-tuples containing parallel Pyomo components @@ -78,7 +80,7 @@ def check_parallel_jacobian( ) if jac is None or nlp is None: - jac, nlp = get_jacobian(model) + jac, nlp = get_jacobian(model, include_greybox=include_greybox) # Get vectors that we will check, and the Pyomo components # they correspond to. @@ -109,7 +111,7 @@ def check_parallel_jacobian( ) if nz in vectors_by_nz: # Store the index as well so we know what component this - # correrponds to. + # corresponds to. vectors_by_nz[nz].append((vec, vecidx)) else: vectors_by_nz[nz] = [(vec, vecidx)] diff --git a/idaes/core/util/diagnostics_tools/writer_utils.py b/idaes/core/util/diagnostics_tools/writer_utils.py index 13ab32dce9..208e7c2755 100644 --- a/idaes/core/util/diagnostics_tools/writer_utils.py +++ b/idaes/core/util/diagnostics_tools/writer_utils.py @@ -17,6 +17,8 @@ __author__ = "Alexander Dowling, Douglas Allan, Andrew Lee, Robby Parker, Ben Knueven" +from sys import stdout + from pyomo.common.collections import ComponentSet from idaes.core.util.model_statistics import ( @@ -37,19 +39,37 @@ unfixed_greybox_variables, greybox_variables, ) -from idaes.core.util.diagnostics_tools.utils import ( - var_in_block, -) MAX_STR_LENGTH = 84 TAB = " " * 4 -def collect_model_statistics(model): +def _var_in_block(var, block): + parent = var.parent_block() + while parent is not None: + if parent is block: + return True + parent = parent.parent_block() + return False + + +# TODO: This appears to duplicate much of the functionality in +# model_statistics.report_statistics +def collect_model_statistics(model, include_greybox=True): """ - Collects statistics about the model that are relevant for diagnostics tools. + Collect general statistics about the model, and format into a readable report. + + Args: + model: Pyomo model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + + Returns: + List of strings containing model statistics, formatted for display """ - vars_in_constraints = variables_in_activated_constraints_set(model) + vars_in_constraints = variables_in_activated_constraints_set( + model, include_greybox=include_greybox + ) fixed_vars_in_constraints = ComponentSet() free_vars_in_constraints = ComponentSet() free_vars_lb = ComponentSet() @@ -60,11 +80,11 @@ def collect_model_statistics(model): for v in vars_in_constraints: if v.fixed: fixed_vars_in_constraints.add(v) - if not var_in_block(v, model): + if not _var_in_block(v, model): ext_fixed_vars_in_constraints.add(v) else: free_vars_in_constraints.add(v) - if not var_in_block(v, model): + if not _var_in_block(v, model): ext_free_vars_in_constraints.add(v) if v.lb is not None: if v.ub is not None: @@ -98,8 +118,8 @@ def collect_model_statistics(model): f"(External: {len(ext_fixed_vars_in_constraints)})" ) stats.append( - f"{TAB}Activated Equality Constraints: {len(activated_equalities_set(model))+number_activated_greybox_equalities(model)} " - f"(Deactivated: {len(deactivated_equalities_set(model))+number_deactivated_greybox_equalities(model)})" + f"{TAB}Activated Equality Constraints: {len(activated_equalities_set(model, include_greybox=include_greybox))} " + f"(Deactivated: {len(deactivated_equalities_set(model, include_greybox=include_greybox))})" ) stats.append( f"{TAB}Activated Inequality Constraints: {len(activated_inequalities_set(model))} " @@ -153,6 +173,9 @@ def write_report_section( None """ + if stream is None: + stream = stdout + stream.write(f"{header * MAX_STR_LENGTH}\n") if title is not None: stream.write(f"{title}\n\n") diff --git a/idaes/core/util/model_statistics.py b/idaes/core/util/model_statistics.py index ccec7f7bc0..79bd2a25b5 100644 --- a/idaes/core/util/model_statistics.py +++ b/idaes/core/util/model_statistics.py @@ -26,6 +26,9 @@ from pyomo.common.collections import ComponentMap, ComponentSet from pyomo.common.deprecation import deprecation_warning from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock +from pyomo.contrib.pynumero.interfaces.external_grey_box_constraint import ( + ExternalGreyBoxConstraint, +) import idaes.logger as idaeslog from idaes.core.scaling import get_scaling_factor @@ -272,20 +275,30 @@ def number_deactivated_blocks(block): # ------------------------------------------------------------------------- # Basic Constraint methods -def total_constraints_set(block): +def total_constraints_set(block, include_greybox=True): """ Method to return a ComponentSet of all Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all Constraint components in block """ - return ComponentSet(activated_block_component_generator(block, ctype=Constraint)) + if include_greybox: + ctypes = (Constraint, ExternalGreyBoxConstraint) + else: + ctypes = (Constraint,) + return ComponentSet( + activated_block_component_generator( + block, ctype=ctypes, include_greybox=include_greybox + ) + ) -def number_total_constraints(block): +def number_total_constraints(block, include_greybox=True): """ Method to return the total number of Constraint components in a model. This will include the number of constraints provided by Greybox models using @@ -293,108 +306,147 @@ def number_total_constraints(block): Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of Constraint components in block """ - number_standard_constraints = sum( - 1 for _ in activated_block_component_generator(block, ctype=Constraint) - ) - number_greybox_constraints = number_activated_greybox_equalities(block) - return number_standard_constraints + number_greybox_constraints + return sum(1 for _ in total_constraints_set(block, include_greybox=include_greybox)) -def activated_constraints_generator(block): +def activated_constraints_generator(block, include_greybox=True): """ Generator which returns all activated Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A generator which returns all activated Constraint components block """ - for c in activated_block_component_generator(block, ctype=Constraint): + if include_greybox: + ctypes = (Constraint, ExternalGreyBoxConstraint) + else: + ctypes = (Constraint,) + for c in activated_block_component_generator( + block, ctype=ctypes, include_greybox=include_greybox + ): if c.active: yield c -def activated_constraints_set(block): +def activated_constraints_set(block, include_greybox=True): """ Method to return a ComponentSet of all activated Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all activated Constraint components in block """ - return ComponentSet(activated_constraints_generator(block)) + return ComponentSet( + activated_constraints_generator(block, include_greybox=include_greybox) + ) -def number_activated_constraints(block): +def number_activated_constraints(block, include_greybox=True): """ Method to return the number of activated Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of activated Constraint components in block """ - return sum(1 for _ in activated_constraints_generator(block)) + return sum( + 1 + for _ in activated_constraints_generator(block, include_greybox=include_greybox) + ) -def deactivated_constraints_generator(block): +def deactivated_constraints_generator(block, include_greybox=True): """ Generator which returns all deactivated Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A generator which returns all deactivated Constraint components block """ + # GreyBox constraints are a special case, as their activity depends on the ExternalGreyBoxBlock + # Using activated_block_component_generator will filter these out, so we need to handle them separately for c in activated_block_component_generator(block, ctype=Constraint): if not c.active: yield c + if include_greybox: + # we could potentially save time by looking only for deactivated greybox blocks here, + # but for robustness we will be thorough + for b in _iter_indexed_block_data_objects( + block, ctype=ExternalGreyBoxBlock, active=None, descend_into=True + ): + for c in b.component_data_objects( + ctype=ExternalGreyBoxConstraint, active=None, descend_into=False + ): + if not c.active: + yield c + -def deactivated_constraints_set(block): +def deactivated_constraints_set(block, include_greybox=True): """ Method to return a ComponentSet of all deactivated Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all deactivated Constraint components in block """ - return ComponentSet(deactivated_constraints_generator(block)) + return ComponentSet( + deactivated_constraints_generator(block, include_greybox=include_greybox) + ) -def number_deactivated_constraints(block): +def number_deactivated_constraints(block, include_greybox=True): """ Method to return the number of deactivated Constraint components in a - model. This will include number of deactivated equalities in a Greybox models - using number_deactivated_greybox_equalities function. + model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of deactivated Constraint components in block """ - standard_equalities = sum(1 for _ in deactivated_constraints_generator(block)) - greybox_equalities = number_deactivated_greybox_equalities(block) - return standard_equalities + greybox_equalities + return sum( + 1 + for _ in deactivated_constraints_generator( + block, include_greybox=include_greybox + ) + ) # ------------------------------------------------------------------------- # Equality Constraints -def total_equalities_generator(block): +def total_equalities_generator(block, include_greybox=True): """ Generator which returns all equality Constraint components in a model. @@ -408,44 +460,62 @@ def total_equalities_generator(block): if c.upper is not None and c.lower is not None and c.upper == c.lower: yield c + # GreyBox constraints are a special case, as they are always equalities + if include_greybox: + for b in _iter_indexed_block_data_objects( + block, ctype=ExternalGreyBoxBlock, active=None, descend_into=True + ): + for c in b.component_data_objects( + ctype=ExternalGreyBoxConstraint, active=None, descend_into=False + ): + yield c + -def total_equalities_set(block): +def total_equalities_set(block, include_greybox=True): """ Method to return a ComponentSet of all equality Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all equality Constraint components in block """ - return ComponentSet(total_equalities_generator(block)) + return ComponentSet( + total_equalities_generator(block, include_greybox=include_greybox) + ) -def number_total_equalities(block): +def number_total_equalities(block, include_greybox=True): """ Method to return the total number of equality Constraint components in a - model. This will include the number of activated equalities Greybox using the number_activated_greybox_equalities function. + model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of equality Constraint components in block """ - standard_equalities = sum(1 for _ in total_equalities_generator(block)) - greybox_equalities = number_activated_greybox_equalities(block) - return standard_equalities + greybox_equalities + return sum( + 1 for _ in total_equalities_generator(block, include_greybox=include_greybox) + ) -def activated_equalities_generator(block): +def activated_equalities_generator(block, include_greybox=True): """ Generator which returns all activated equality Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A generator which returns all activated equality Constraint components @@ -461,47 +531,59 @@ def activated_equalities_generator(block): ): yield c + # GreyBox constraints are a special case, as they are always equalities + if include_greybox: + for b in _iter_indexed_block_data_objects( + block, ctype=ExternalGreyBoxBlock, active=True, descend_into=True + ): + for c in b.component_data_objects( + ctype=ExternalGreyBoxConstraint, active=True, descend_into=False + ): + yield c -def activated_equalities_set(block): + +def activated_equalities_set(block, include_greybox=True): """ Method to return a ComponentSet of all activated equality Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all activated equality Constraint components in block """ - return ComponentSet(activated_equalities_generator(block)) + return ComponentSet( + activated_equalities_generator(block, include_greybox=include_greybox) + ) -def number_activated_equalities(block): +def number_activated_equalities(block, include_greybox=True): """ Method to return the number of activated equality Constraint components in - a model. This will include number of equalities in Greybox model using number_activated_greybox_equalities function. + a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of activated equality Constraint components in block """ return sum( - 1 for _ in activated_equalities_generator(block) - ) + number_activated_greybox_equalities(block) + 1 + for _ in activated_equalities_generator(block, include_greybox=include_greybox) + ) def number_activated_greybox_equalities(block) -> int: """ Function to compute total number of equality constraints for all GreyBox objects in this block. - A GreyBox model is always assumed to be 0DOFs where each output[i]==f(inputs) - where f is GreyBox model, this should be true regardless if - GreyBox model is doing internal optimization or not, as every output - is calculated through the GreyBox internal model using provided inputs. - Args: block : pyomo concrete model or pyomo block @@ -510,8 +592,10 @@ def number_activated_greybox_equalities(block) -> int: """ equalities = 0 for grey_box in activated_greybox_block_set(block): - equalities += len(grey_box.outputs) - equalities += grey_box.get_external_model().n_equality_constraints() + for _ in grey_box.component_data_objects( + ctype=ExternalGreyBoxConstraint, active=None, descend_into=False + ): + equalities += 1 return equalities @@ -519,11 +603,6 @@ def number_deactivated_greybox_equalities(block) -> int: """ Function to compute total number of equality constraints for all GreyBox objects in this block. - A GreyBox model is always assumed to be 0DOFs where each output[i]==f(inputs) - where f is GreyBox model, this should be true regardless if - GreyBox model is doing internal optimization or not, as every output - is calculated through a the GreyBox internal model using provided inputs. - Args: block : pyomo concrete model or pyomo block @@ -532,57 +611,70 @@ def number_deactivated_greybox_equalities(block) -> int: """ equalities = 0 for grey_box in deactivated_greybox_block_set(block): - equalities += len(grey_box.outputs) - equalities += grey_box.get_external_model().n_equality_constraints() + for _ in grey_box.component_data_objects( + ctype=ExternalGreyBoxConstraint, active=None, descend_into=False + ): + equalities += 1 return equalities -def deactivated_equalities_generator(block): +def deactivated_equalities_generator(block, include_greybox=True): """ Generator which returns all deactivated equality Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A generator which returns all deactivated equality Constraint components block """ - for c in total_equalities_generator(block): + for c in total_equalities_generator(block, include_greybox=include_greybox): if not c.active: yield c -def deactivated_equalities_set(block): +def deactivated_equalities_set(block, include_greybox=True): """ Method to return a ComponentSet of all deactivated equality Constraint components in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all deactivated equality Constraint components in block """ - return ComponentSet(deactivated_equalities_generator(block)) + return ComponentSet( + deactivated_equalities_generator(block, include_greybox=include_greybox) + ) -def number_deactivated_equalities(block): +def number_deactivated_equalities(block, include_greybox=True): """ Method to return the number of deactivated equality Constraint components in a model. This will include the number of deactivated equality constraints in Greybox models. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of deactivated equality Constraint components in block """ - standard_equalities = sum(1 for _ in deactivated_equalities_generator(block)) - greybox_equalities = number_deactivated_greybox_equalities(block) - return standard_equalities + greybox_equalities + return sum( + 1 + for _ in deactivated_equalities_generator( + block, include_greybox=include_greybox + ) + ) # ------------------------------------------------------------------------- @@ -729,127 +821,158 @@ def number_deactivated_inequalities(block): # Basic Variable Methods # Always use ComponentSets for Vars to avoid duplication of References # i.e. number methods should always use the ComponentSet, not a generator -def variables_set(block): +def variables_generator(block, include_greybox=True, active=True): + """ + Generator which returns all Var components in a model. + + Args: + block : model to be studied + include_greybox : whether to include greybox variables + active : whether to include only active sub-blocks (default = True) + + Returns: + A generator which returns all Var components in block + """ + for var in _iter_indexed_block_data_objects( + block, ctype=Var, active=active, descend_into=True + ): + yield var + + if include_greybox: + for egb in _iter_indexed_block_data_objects( + block, ctype=ExternalGreyBoxBlock, active=None, descend_into=True + ): + for v in egb.component_data_objects( + ctype=Var, active=None, descend_into=False + ): + yield v + + +def variables_set(block, include_greybox=True, active=True): """ Method to return a ComponentSet of all Var components in a model. Args: block : model to be studied + include_greybox : whether to include greybox variables + active : whether to include only active sub-blocks (default = True) Returns: A ComponentSet including all Var components in block """ var_set = ComponentSet() - for var in _iter_indexed_block_data_objects( - block, ctype=Var, active=True, descend_into=True + for var in variables_generator( + block, include_greybox=include_greybox, active=active ): var_set.add(var) - for var in greybox_variables(block): - var_set.add(var) return var_set -def number_variables(block): +def number_variables(block, include_greybox=True, active=True): """ Method to return the number of Var components in a model. Args: block : model to be studied + include_greybox : whether to include greybox variables + active : whether to include only active sub-blocks (default = True) Returns: Number of Var components in block """ - return len(variables_set(block)) + return len(variables_set(block, include_greybox=include_greybox, active=active)) -def fixed_variables_generator(block): +def fixed_variables_generator(block, include_greybox=True): """ Generator which returns all fixed Var components in a model. Args: block : model to be studied + include_greybox : whether to include greybox variables Returns: A generator which returns all fixed Var components block """ - for v in _iter_indexed_block_data_objects( - block, ctype=Var, active=True, descend_into=True - ): - if v.fixed: - yield v - # include greybox variables in set - for v in greybox_variables(block): + for v in variables_generator(block, include_greybox=include_greybox): if v.fixed: yield v -def fixed_variables_set(block): +def fixed_variables_set(block, include_greybox=True): """ Method to return a ComponentSet of all fixed Var components in a model. Args: block : model to be studied + include_greybox : whether to include greybox variables Returns: A ComponentSet including all fixed Var components in block """ - return ComponentSet(fixed_variables_generator(block)) + return ComponentSet( + fixed_variables_generator(block, include_greybox=include_greybox) + ) -def number_fixed_variables(block): +def number_fixed_variables(block, include_greybox=True): """ Method to return the number of fixed Var components in a model. Args: block : model to be studied + include_greybox : whether to include greybox variables Returns: Number of fixed Var components in block """ - return len(fixed_variables_set(block)) + return len(fixed_variables_set(block, include_greybox=include_greybox)) -def unfixed_variables_generator(block): +def unfixed_variables_generator(block, include_greybox=True): """ Generator which returns all unfixed Var components in a model. Args: block : model to be studied + include_greybox : whether to include greybox variables Returns: A generator which returns all unfixed Var components block """ - for v in _iter_indexed_block_data_objects( - block, ctype=Var, active=True, descend_into=True - ): + for v in variables_generator(block, include_greybox=include_greybox): if not v.fixed: yield v -def unfixed_variables_set(block): +def unfixed_variables_set(block, include_greybox=True): """ Method to return a ComponentSet of all unfixed Var components in a model. Args: block : model to be studied + include_greybox : whether to include greybox variables Returns: A ComponentSet including all unfixed Var components in block """ - return ComponentSet(unfixed_variables_generator(block)) + return ComponentSet( + unfixed_variables_generator(block, include_greybox=include_greybox) + ) -def number_unfixed_variables(block): +def number_unfixed_variables(block, include_greybox=True): """ Method to return the number of unfixed Var components in a model. Args: block : model to be studied + include_greybox : whether to include greybox variables Returns: Number of unfixed Var components in block """ - return len(unfixed_variables_set(block)) + return len(unfixed_variables_set(block, include_greybox=include_greybox)) def variables_near_bounds_generator( @@ -860,6 +983,8 @@ def variables_near_bounds_generator( skip_ub=False, abs_tol=1e-4, rel_tol=1e-4, + include_greybox=True, + apply_scaling=True, ): """ Generator which returns all Var components in a model which have a value @@ -871,6 +996,8 @@ def variables_near_bounds_generator( rel_tol : relative tolerance for inclusion in generator (default = 1e-4) skip_lb: Boolean to skip lower bound (default = False) skip_ub: Boolean to skip upper bound (default = False) + include_greybox : whether to include greybox variables (default = True) + apply_scaling: whether to apply scaling factors to the variable when determining if it is near a bound (default = True) Returns: A generator which returns all Var components block that are close to a @@ -893,13 +1020,11 @@ def variables_near_bounds_generator( abs_tol = tol rel_tol = tol - for v in _iter_indexed_block_data_objects( - block, ctype=Var, active=True, descend_into=True - ): + for v in variables_generator(block, include_greybox=include_greybox): # To avoid errors, check that v has a value if v.value is None: continue - sf = get_scaling_factor(v, default=1, warning=False) + sf = get_scaling_factor(v, default=1, warning=False) if apply_scaling else 1 # First, magnitude of variable if v.ub is not None and v.lb is not None: # Both upper and lower bounds, apply tol to (upper - lower) @@ -930,6 +1055,8 @@ def variables_near_bounds_set( skip_ub=False, abs_tol=1e-4, rel_tol=1e-4, + include_greybox=True, + apply_scaling=True, ): """ Method to return a ComponentSet of all Var components in a model which have @@ -941,6 +1068,9 @@ def variables_near_bounds_set( rel_tol : relative tolerance for inclusion in generator (default = 1e-4) skip_lb: Boolean to skip lower bound (default = False) skip_ub: Boolean to skip upper bound (default = False) + include_greybox : whether to include greybox variables (default = True) + apply_scaling: whether to apply scaling factors to the variable when determining + if it is near a bound (default = True) Returns: A ComponentSet including all Var components block that are close to a @@ -948,12 +1078,22 @@ def variables_near_bounds_set( """ return ComponentSet( variables_near_bounds_generator( - block, tol, relative, skip_lb, skip_ub, abs_tol, rel_tol + block, + tol, + relative, + skip_lb, + skip_ub, + abs_tol, + rel_tol, + include_greybox=include_greybox, + apply_scaling=apply_scaling, ) ) -def number_variables_near_bounds(block, tol=None, abs_tol=1e-4, rel_tol=1e-4): +def number_variables_near_bounds( + block, tol=None, abs_tol=1e-4, rel_tol=1e-4, include_greybox=True +): """ Method to return the number of all Var components in a model which have a value within tol (relative) of a bound. @@ -962,24 +1102,128 @@ def number_variables_near_bounds(block, tol=None, abs_tol=1e-4, rel_tol=1e-4): block : model to be studied abs_tol : absolute tolerance for inclusion in generator (default = 1e-4) rel_tol : relative tolerance for inclusion in generator (default = 1e-4) + include_greybox : whether to include greybox variables (default = True) Returns: Number of components block that are close to a bound """ return len( - variables_near_bounds_set(block, tol=tol, abs_tol=abs_tol, rel_tol=rel_tol) + variables_near_bounds_set( + block, + tol=tol, + abs_tol=abs_tol, + rel_tol=rel_tol, + include_greybox=include_greybox, + ) + ) + + +def variables_violating_bounds_generator( + block: Block, + abs_tol: float = 1e-4, + rel_tol: float = 1e-4, + include_greybox: bool = True, + apply_scaling: bool = True, +): + """ + Generator which returns all Var components in a model which have a value + that violates their bounds. + + Args: + block : model to be studied + abs_tol : absolute tolerance for violation of bounds + rel_tol : relative tolerance for violation of bounds + include_greybox : whether to include greybox variables (default = True) + apply_scaling: whether to apply scaling factors to the variable when determining if it is violating bounds (default = True) + + Returns: + A generator which returns all Var components block that are violating + their bounds + """ + for v in variables_generator(block, include_greybox=include_greybox): + if v.value is not None: + sf = get_scaling_factor(v, default=1, warning=False) if apply_scaling else 1 + tolerance = max(abs_tol, abs(value(v)) * rel_tol) + if v.lb is not None and sf * value(v) <= sf * v.lb - tolerance: + yield v + elif v.ub is not None and sf * value(v) >= sf * v.ub + tolerance: + yield v + + +def variables_violating_bounds_set( + block: Block, + abs_tol: float = 1e-4, + rel_tol: float = 1e-4, + include_greybox: bool = True, + apply_scaling: bool = True, +): + """ + Method to return a ComponentSet of all Var components in a model which have + a value that violates their bounds. + + Args: + block : model to be studied + abs_tol : absolute tolerance for violation of bounds + rel_tol : relative tolerance for violation of bounds + include_greybox : whether to include greybox variables (default = True) + apply_scaling: whether to apply scaling factors to the variable when determining if it is violating bounds (default = True) + Returns: + A ComponentSet including all Var components block that are violating + their bounds + """ + return ComponentSet( + variables_violating_bounds_generator( + block, + abs_tol=abs_tol, + rel_tol=rel_tol, + include_greybox=include_greybox, + apply_scaling=apply_scaling, + ) + ) + + +def number_variables_violating_bounds( + block: Block, + abs_tol: float = 1e-4, + rel_tol: float = 1e-4, + include_greybox: bool = True, + apply_scaling: bool = True, +): + """ + Method to return the number of all Var components in a model which have + a value that violates their bounds. + + Args: + block : model to be studied + abs_tol : absolute tolerance for violation of bounds + rel_tol : relative tolerance for violation of bounds + include_greybox : whether to include greybox variables (default = True) + apply_scaling: whether to apply scaling factors to the variable when determining if it is violating bounds (default = True) + Returns: + Number of components block that are violating their bounds + """ + return len( + variables_violating_bounds_set( + block, + abs_tol=abs_tol, + rel_tol=rel_tol, + include_greybox=include_greybox, + apply_scaling=apply_scaling, + ) ) # ------------------------------------------------------------------------- # Variables in Constraints -def variables_in_activated_constraints_set(block): +def variables_in_activated_constraints_set(block, include_greybox=True): """ Method to return a ComponentSet of all Var components which appear within a Constraint in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all Var components which appear within @@ -991,101 +1235,139 @@ def variables_in_activated_constraints_set(block): ): for v in identify_variables(c.body): var_set.add(v) - # include any vars in greyboxes - for v in greybox_variables(block): - var_set.add(v) + + if include_greybox: + for egb in _iter_indexed_block_data_objects( + block, ctype=ExternalGreyBoxBlock, active=True, descend_into=True + ): + # For simplicity and efficiency, we will assume all variables in greybox blocks are included + # For this to be False, the grey box itself would have to have a dangling variable + for v in egb.component_data_objects( + ctype=Var, active=None, descend_into=False + ): + var_set.add(v) return var_set -def number_variables_in_activated_constraints(block): +def number_variables_in_activated_constraints(block, include_greybox=True): """ Method to return the number of Var components that appear within active Constraints in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of Var components which appear within active Constraints in block """ - return len(variables_in_activated_constraints_set(block)) + return len( + variables_in_activated_constraints_set(block, include_greybox=include_greybox) + ) -def variables_not_in_activated_constraints_set(block): +def variables_not_in_activated_constraints_set(block, include_greybox=True): """ Method to return a ComponentSet of all Var components which do not appear within a Constraint in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all Var components which do not appear within activated Constraints in block """ - var_set = ComponentSet() - - active_vars = variables_in_activated_constraints_set(block) + all_vars = variables_set(block, include_greybox=include_greybox) + active_vars = variables_in_activated_constraints_set( + block, include_greybox=include_greybox + ) - for v in _iter_indexed_block_data_objects( - block, ctype=Var, active=True, descend_into=True - ): + var_set = ComponentSet() + for v in all_vars: if v not in active_vars: var_set.add(v) return var_set -def number_variables_not_in_activated_constraints(block): +def number_variables_not_in_activated_constraints(block, include_greybox=True): """ Method to return the number of Var components that do not appear within active Constraints in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of Var components which do not appear within active Constraints in block """ - return len(variables_not_in_activated_constraints_set(block)) + return len( + variables_not_in_activated_constraints_set( + block, include_greybox=include_greybox + ) + ) -def variables_in_activated_equalities_set(block): +def variables_in_activated_equalities_set(block, include_greybox=True): """ Method to return a ComponentSet of all Var components which appear within an equality Constraint in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all Var components which appear within activated equality Constraints in block """ var_set = ComponentSet() - for c in activated_equalities_generator(block): + # identify_variables does not work on ExternalGreyBoxConstraints, so we need to handle these separately + for c in activated_equalities_generator(block, include_greybox=False): for v in identify_variables(c.body): var_set.add(v) - # include any vars in greyboxes - for v in greybox_variables(block): - var_set.add(v) + + if include_greybox: + for egb in _iter_indexed_block_data_objects( + block, ctype=ExternalGreyBoxBlock, active=True, descend_into=True + ): + # For simplicity and efficiency, we will assume all variables in greybox blocks are included + # For this to be False, the grey box itself would have to have a dangling variable + # All grey box constraints are equalities, so we do not need to check for inequality constraints here + for v in egb.component_data_objects( + ctype=Var, active=None, descend_into=False + ): + var_set.add(v) return var_set -def number_variables_in_activated_equalities(block): +def number_variables_in_activated_equalities(block, include_greybox=True): """ Method to return the number of Var components which appear within activated equality Constraints in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + + block : model to be studied Returns: Number of Var components which appear within activated equality Constraints in block """ - return len(variables_in_activated_equalities_set(block)) + return len( + variables_in_activated_equalities_set(block, include_greybox=include_greybox) + ) def variables_in_activated_inequalities_set(block): @@ -1100,6 +1382,7 @@ def variables_in_activated_inequalities_set(block): A ComponentSet including all Var components which appear within activated inequality Constraints in block """ + # GreyBox constraints are all equalities, so we do not need to check for them here var_set = ComponentSet() for c in activated_inequalities_generator(block): for v in identify_variables(c.body): @@ -1156,54 +1439,68 @@ def number_variables_only_in_inequalities(block): # ------------------------------------------------------------------------- # Fixed Variables in Constraints -def fixed_variables_in_activated_equalities_set(block): +def fixed_variables_in_activated_equalities_set(block, include_greybox=True): """ Method to return a ComponentSet of all fixed Var components which appear within an equality Constraint in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all fixed Var components which appear within activated equality Constraints in block """ var_set = ComponentSet() - for v in variables_in_activated_equalities_set(block): + for v in variables_in_activated_equalities_set( + block, include_greybox=include_greybox + ): if v.fixed: var_set.add(v) return var_set -def number_fixed_variables_in_activated_equalities(block): +def number_fixed_variables_in_activated_equalities(block, include_greybox=True): """ Method to return the number of fixed Var components which appear within activated equality Constraints in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of fixed Var components which appear within activated equality Constraints in block """ - return len(fixed_variables_in_activated_equalities_set(block)) + return len( + fixed_variables_in_activated_equalities_set( + block, include_greybox=include_greybox + ) + ) -def unfixed_variables_in_activated_equalities_set(block): +def unfixed_variables_in_activated_equalities_set(block, include_greybox=True): """ Method to return a ComponentSet of all unfixed Var components which appear within an activated equality Constraint in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet of all unfixed Var components which appear within activated equality Constraints in block """ var_set = ComponentSet() - for v in variables_in_activated_equalities_set(block): + for v in variables_in_activated_equalities_set( + block, include_greybox=include_greybox + ): if not v.fixed: var_set.add(v) return var_set @@ -1219,6 +1516,7 @@ def unfixed_greybox_variables(block): Returns: A ComponentSet of all unfixed Var components which appear in Greybox models """ + # GreyBox variable should never be fixed - this function is useful to confirm that var_set = ComponentSet() for var in greybox_variables(block): if not var.fixed: @@ -1239,10 +1537,10 @@ def greybox_variables(block): """ var_set = ComponentSet() for grey_box in activated_greybox_block_set(block): - for in_var in grey_box.inputs: - var_set.add(grey_box.inputs[in_var]) - for out_var in grey_box.outputs: - var_set.add(grey_box.outputs[out_var]) + for v in grey_box.component_data_objects( + ctype=Var, active=None, descend_into=False + ): + var_set.add(v) return var_set @@ -1272,19 +1570,25 @@ def number_of_greybox_variables(block): return len(greybox_variables(block)) -def number_unfixed_variables_in_activated_equalities(block): +def number_unfixed_variables_in_activated_equalities(block, include_greybox=True): """ Method to return the number of unfixed Var components which appear within activated equality Constraints in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of unfixed Var components which appear within activated equality Constraints in block """ - return len(unfixed_variables_in_activated_equalities_set(block)) + return len( + unfixed_variables_in_activated_equalities_set( + block, include_greybox=include_greybox + ) + ) def fixed_variables_only_in_inequalities(block): @@ -1321,70 +1625,316 @@ def number_fixed_variables_only_in_inequalities(block): return len(fixed_variables_only_in_inequalities(block)) +def external_variables_set(block: Block, include_greybox: bool = True) -> ComponentSet: + """ + Method to return a ComponentSet of all Var components which appear within a + Constraint in a model, but do not appear in the variables_set of the model. + + Args: + block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + + Returns: + A ComponentSet including all Var components which appear within + activated Constraints in block, but do not appear in the variables_set of + the model + """ + ext_vars = ComponentSet() + local_var_set = variables_set(block, include_greybox=include_greybox, active=None) + conc_var_set = variables_in_activated_constraints_set( + block, include_greybox=include_greybox + ) + for v in conc_var_set: + if v not in local_var_set: + ext_vars.add(v) + return ext_vars + + +def number_external_variables(block: Block, include_greybox: bool = True) -> int: + """ + Method to return the number of Var components which appear within a + Constraint in a model, but do not appear in the variables_set of the model. + + Args: + block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + + Returns: + Number of Var components which appear within activated Constraints in block, + but do not appear in the variables_set of the model + """ + return len(external_variables_set(block, include_greybox=include_greybox)) + + +def variables_fixed_to_zero_set( + block: Block, include_greybox: bool = True +) -> ComponentSet: + """ + Method to return a ComponentSet of all Var components which are fixed to zero in a model. + + Args: + block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + + Returns: + A ComponentSet including all Var components which are fixed to zero in block + """ + var_set = ComponentSet() + for v in variables_set(block, include_greybox=include_greybox, active=None): + if v.fixed and value(v) == 0: + var_set.add(v) + return var_set + + +def number_variables_fixed_to_zero(block: Block, include_greybox: bool = True) -> int: + """ + Method to return the number of Var components which are fixed to zero in a model. + + Args: + block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + + Returns: + Number of Var components which are fixed to zero in block + """ + return len(variables_fixed_to_zero_set(block, include_greybox=include_greybox)) + + +def variables_near_zero_set( + block: Block, + tol: float = 1e-4, + include_greybox: bool = True, + scale_variables: bool = True, +) -> ComponentSet: + """ + Method to return a ComponentSet of all Var components which have a value within tol of zero in a model. + + Args: + block : model to be studied + tol: Tolerance for inclusion in generator (default = 1e-4) + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + scale_variables: Boolean to scale variables by their scaling factor when determining if they are near zero (default = True) + + Returns: + A ComponentSet including all Var components which have a value within tol of zero in block + """ + var_set = ComponentSet() + for v in variables_generator(block, include_greybox=include_greybox, active=None): + sf = get_scaling_factor(v, default=1, warning=False) if scale_variables else 1 + if v.value is not None and abs(value(v) * sf) <= tol: + var_set.add(v) + return var_set + + +def number_variables_near_zero( + block: Block, + tol: float = 1e-4, + include_greybox: bool = True, + scale_variables: bool = True, +) -> int: + """ + Method to return the number of Var components which have a value within tol of zero in a model. + + Args: + block : model to be studied + tol: Tolerance for inclusion in generator (default = 1e-4) + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + scale_variables: Boolean to scale variables by their scaling factor when determining if they are near zero (default = True) + + Returns: + Number of Var components which have a value within tol of zero in block + """ + return len( + variables_near_zero_set( + block, + tol=tol, + include_greybox=include_greybox, + scale_variables=scale_variables, + ) + ) + + +def variables_with_none_value_set(block, include_greybox=True): + """ + Method to return a ComponentSet of all Var components which have a value of None in a model. + + Args: + block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + + Returns: + A ComponentSet including all Var components which have a value of None in block + """ + var_set = ComponentSet() + for v in variables_generator(block, include_greybox=include_greybox, active=None): + if v.value is None: + var_set.add(v) + return var_set + + +def number_variables_with_none_value(block, include_greybox=True): + """ + Method to return the number of Var components which have a value of None in a model. + + Args: + block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + + Returns: + Number of Var components which have a value of None in block + """ + return len(variables_with_none_value_set(block, include_greybox=include_greybox)) + + +def variables_with_extreme_values_set( + block: Block, + large: float, + small: float, + zero: float, + include_greybox: bool = True, + apply_scaling: bool = True, +): + """Method to return a ComponentSet of all Var components which have a value with magnitude larger than large or smaller than small (but larger than zero) in a model. + + Args: + block : model to be studied + large: Threshold for large values + small: Threshold for small values + zero: Threshold for zero values (small values must be larger than this value to be included) + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + apply_scaling: Boolean to apply scaling factors to the variable when determining if it has an extreme value (default = True) + Returns: + A ComponentSet including all Var components which have a value with magnitude larger than large or smaller than small (but larger than zero) in block + """ + extreme_vars = ComponentSet() + for v in variables_generator(block, include_greybox=include_greybox, active=None): + sf = get_scaling_factor(v, default=1, warning=False) if apply_scaling else 1 + if v.value is not None: + mag = sf * abs(value(v)) + if mag > abs(large): + extreme_vars.add(v) + elif mag < abs(small) and mag > abs(zero): + extreme_vars.add(v) + + return extreme_vars + + +def number_variables_with_extreme_values( + block: Block, + large: float, + small: float, + zero: float, + include_greybox: bool = True, + apply_scaling: bool = True, +): + """Method to return the number of Var components which have a value with magnitude larger than large or smaller than small (but larger than zero) in a model. + + Args: + block : model to be studied + large: Threshold for large values + small: Threshold for small values + zero: Threshold for zero values (small values must be larger than this value to be included) + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) + apply_scaling: Boolean to apply scaling factors to the variable when determining if it has an extreme value (default = True) + Returns: + Number of Var components which have a value with magnitude larger than large or smaller than small (but larger than zero) in block + """ + return len( + variables_with_extreme_values_set( + block, + large, + small, + zero, + include_greybox=include_greybox, + apply_scaling=apply_scaling, + ) + ) + + # ------------------------------------------------------------------------- # Unused and un-Transformed Variables -def unused_variables_set(block): +def unused_variables_set(block, include_greybox=True): """ Method to return a ComponentSet of all Var components which do not appear within any activated Constraint in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all Var components which do not appear within any Constraints in block """ - return variables_set(block) - variables_in_activated_constraints_set(block) + return variables_set( + block, include_greybox=include_greybox + ) - variables_in_activated_constraints_set(block, include_greybox=include_greybox) -def number_unused_variables(block): +def number_unused_variables(block, include_greybox=True): """ Method to return the number of Var components which do not appear within any activated Constraint in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of Var components which do not appear within any activated Constraints in block """ - return len(unused_variables_set(block)) + return len(unused_variables_set(block, include_greybox=include_greybox)) -def fixed_unused_variables_set(block): +def fixed_unused_variables_set(block, include_greybox=True): """ Method to return a ComponentSet of all fixed Var components which do not appear within any activated Constraint in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all fixed Var components which do not appear within any Constraints in block """ var_set = ComponentSet() - for v in unused_variables_set(block): + for v in unused_variables_set(block, include_greybox=include_greybox): if v.fixed: var_set.add(v) return var_set -def number_fixed_unused_variables(block): +def number_fixed_unused_variables(block, include_greybox=True): """ Method to return the number of fixed Var components which do not appear within any activated Constraint in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of fixed Var components which do not appear within any activated Constraints in block """ - return len(fixed_unused_variables_set(block)) + return len(fixed_unused_variables_set(block, include_greybox=include_greybox)) def derivative_variables_set(block): @@ -1593,22 +2143,29 @@ def number_expressions(block): # ------------------------------------------------------------------------- # Other model statistics -def degrees_of_freedom(block): +def degrees_of_freedom(block, include_greybox=True): """ Method to return the degrees of freedom of a model. Args: block : model to be studied + include_greybox: Boolean to include GreyBox models (default = True) Returns: Number of degrees of freedom in block. """ return number_unfixed_variables_in_activated_equalities( - block - ) - number_activated_equalities(block) + block, include_greybox=include_greybox + ) - number_activated_equalities(block, include_greybox=include_greybox) -def large_residuals_set(block, tol=1e-5, return_residual_values=False): +def large_residuals_set( + block, + tol=1e-5, + return_residual_values=False, + include_greybox=True, + apply_scaling=True, +): """ Method to return a ComponentSet of all Constraint components with a residual greater than a given threshold which appear in a model. @@ -1618,6 +2175,9 @@ def large_residuals_set(block, tol=1e-5, return_residual_values=False): tol : residual threshold for inclusion in ComponentSet return_residual_values: boolean, if true return dictionary with residual values + include_greybox: boolean, if true include GreyBox constraints (default = True) + apply_scaling: boolean, if true apply scaling factors to the residuals when determining + if they are large (default = True) Returns: large_residual_set: A ComponentSet including all Constraint components @@ -1627,12 +2187,14 @@ def large_residuals_set(block, tol=1e-5, return_residual_values=False): return_residual_values is true) """ large_residuals_set = ComponentSet() - if return_residual_values: - residual_values = ComponentMap() - for c in _iter_indexed_block_data_objects( - block, ctype=Constraint, active=True, descend_into=True - ): - sf = get_scaling_factor(c, default=1, warning=False) + residual_values = ComponentMap() + + for c in activated_constraints_generator(block, include_greybox=include_greybox): + # Grey Box constraints cannot be scaled in the normal fashion + sf = 1 + if apply_scaling and not isinstance(c, ExternalGreyBoxConstraint): + sf = get_scaling_factor(c, default=1, warning=False) + try: val = value(c.body) except ValueError: @@ -1662,7 +2224,7 @@ def large_residuals_set(block, tol=1e-5, return_residual_values=False): return large_residuals_set -def number_large_residuals(block, tol=1e-5): +def number_large_residuals(block, tol=1e-5, include_greybox=True): """ Method to return the number Constraint components with a residual greater than a given threshold which appear in a model. @@ -1670,94 +2232,115 @@ def number_large_residuals(block, tol=1e-5): Args: block : model to be studied tol : residual threshold for inclusion in ComponentSet + include_greybox: boolean, if true include GreyBox constraints (default = True) Returns: Number of Constraint components with a residual greater than tol which appear in block """ - lr = 0 - for c in _iter_indexed_block_data_objects( - block, ctype=Constraint, active=True, descend_into=True - ): - if c.active and value(c.lower - c.body()) > tol: - lr += 1 - elif c.active and value(c.body() - c.upper) > tol: - lr += 1 - return lr + return len(large_residuals_set(block, tol=tol, include_greybox=include_greybox)) -def active_variables_in_deactivated_blocks_set(block): +def active_variables_in_deactivated_blocks_set(block, include_greybox=True): """ Method to return a ComponentSet of any Var components which appear within an active Constraint but belong to a deactivated Block in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including any Var components which belong to a deactivated Block but appear in an activate Constraint in block """ var_set = ComponentSet() + block_set = activated_blocks_set(block) - for v in variables_in_activated_constraints_set(block): + if include_greybox: + for c in activated_greybox_block_set(block): + block_set.add(c) + + for v in variables_in_activated_constraints_set( + block, include_greybox=include_greybox + ): if v.parent_block() not in block_set: var_set.add(v) return var_set -def number_active_variables_in_deactivated_blocks(block): +def number_active_variables_in_deactivated_blocks(block, include_greybox=True): """ Method to return the number of Var components which appear within an active Constraint but belong to a deactivated Block in a model. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of Var components which belong to a deactivated Block but appear in an activate Constraint in block """ - return len(active_variables_in_deactivated_blocks_set(block)) + return len( + active_variables_in_deactivated_blocks_set( + block, include_greybox=include_greybox + ) + ) -def variables_with_none_value_in_activated_equalities_set(block): +def variables_with_none_value_in_activated_equalities_set(block, include_greybox=True): """ Method to return a ComponentSet of all Var components which have a value of None in the set of activated constraints. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: A ComponentSet including all Var components which have a value of None in the set of activated constraints. """ var_set = ComponentSet() - for v in variables_in_activated_equalities_set(block): + for v in variables_in_activated_equalities_set( + block, include_greybox=include_greybox + ): if v.value is None: var_set.add(v) return var_set -def number_variables_with_none_value_in_activated_equalities(block): +def number_variables_with_none_value_in_activated_equalities( + block, include_greybox=True +): """ Method to return the number of Var components which have a value of None in the set of activated constraints. Args: block : model to be studied + include_greybox: Boolean to include implicit constraints from GreyBox + models (default = True) Returns: Number of Var components which have a value of None in the set of activated constraints. """ - return len(variables_with_none_value_in_activated_equalities_set(block)) + return len( + variables_with_none_value_in_activated_equalities_set( + block, include_greybox=include_greybox + ) + ) # ------------------------------------------------------------------------- # Reporting methods +# TODO: This appears to be duplicated in the diagnostics module - we should consider whether to merge these or keep them separate def report_statistics(block, ostream=None): """ Method to print a report of the model statistics for a Pyomo Block @@ -1851,7 +2434,7 @@ def report_statistics(block, ostream=None): # ------------------------------------------------------------------------- # Common sub-methods -def activated_block_component_generator(block, ctype): +def activated_block_component_generator(block, ctype, include_greybox=False): """ Generator which returns all the components of a given ctype which exist in activated Blocks within a model. @@ -1859,6 +2442,7 @@ def activated_block_component_generator(block, ctype): Args: block : model to be studied ctype : type of Pyomo component to be returned by generator. + include_greybox : whether to include greybox components Returns: A generator which returns all components of ctype which appear in @@ -1876,3 +2460,13 @@ def activated_block_component_generator(block, ctype): ): for c in b.component_data_objects(ctype=ctype, active=None, descend_into=False): yield c + + # If include_greybox is True, also yield components in active greybox sub-blocks + if include_greybox: + for b in _iter_indexed_block_data_objects( + block, ctype=ExternalGreyBoxBlock, active=True, descend_into=True + ): + for c in b.component_data_objects( + ctype=ctype, active=None, descend_into=False + ): + yield c diff --git a/idaes/core/util/tests/test_model_statistics.py b/idaes/core/util/tests/test_model_statistics.py index 10f7a6672b..63e3fc478a 100644 --- a/idaes/core/util/tests/test_model_statistics.py +++ b/idaes/core/util/tests/test_model_statistics.py @@ -999,6 +999,122 @@ def test_number_derivative_variables(): assert number_derivative_variables(m) == 0 +# ------------------------------------------------------------------------- +# Tests for new model_statistics functions +@pytest.mark.unit +def test_variables_fixed_to_zero_set(m): + m.v_zero = Var(initialize=0) + m.v_zero.fix(0) + var_set = variables_fixed_to_zero_set(m) + assert m.v_zero in var_set + assert len(var_set) == 1 + + +@pytest.mark.unit +def test_number_variables_fixed_to_zero(m): + m.v_zero = Var(initialize=0) + m.v_zero.fix(0) + assert number_variables_fixed_to_zero(m) == 1 + + +@pytest.mark.unit +def test_variables_near_zero_set(m): + m.v_small = Var(initialize=1e-5) + m.v_small.value = 1e-5 + var_set = variables_near_zero_set(m, tol=1e-4) + assert m.v_small in var_set + assert len(var_set) == 1 + + set_scaling_factor(m.v_small, 1e5) + var_set = variables_near_zero_set(m, tol=1e-4, scale_variables=True) + assert len(var_set) == 0 + var_set = variables_near_zero_set(m, tol=1e-4, scale_variables=False) + assert len(var_set) == 1 + assert m.v_small in var_set + + +@pytest.mark.unit +def test_number_variables_near_zero(m): + m.v_small = Var(initialize=1e-5) + m.v_small.value = 1e-5 + assert number_variables_near_zero(m, tol=1e-4) == 1 + set_scaling_factor(m.v_small, 1e5) + assert number_variables_near_zero(m, tol=1e-4, scale_variables=True) == 0 + assert number_variables_near_zero(m, tol=1e-4, scale_variables=False) == 1 + + +@pytest.mark.unit +def test_variables_with_none_value_set(m): + m.v_none = Var() + # Do not set value + var_set = variables_with_none_value_set(m) + assert m.v_none in var_set + for v in var_set: + if v is m.v_none: + assert v.value is None + else: + assert v.name.startswith("dv[") + assert len(var_set) == 12 + + +@pytest.mark.unit +def test_number_variables_with_none_value(m): + m.v_none = Var() + assert number_variables_with_none_value(m) == 12 + + +@pytest.mark.unit +def test_variables_with_extreme_values_set(m): + m.v_large = Var(initialize=1e6) + m.v_small = Var(initialize=1e-8) + m.v_zero = Var(initialize=0) + var_set = variables_with_extreme_values_set(m, large=1e5, small=1e-7, zero=1e-10) + assert m.v_large in var_set + assert m.v_small in var_set + assert m.v_zero not in var_set + assert len(var_set) == 2 + + set_scaling_factor(m.v_small, 1e8) + var_set = variables_with_extreme_values_set( + m, large=1e5, small=1e-7, zero=1e-10, apply_scaling=True + ) + assert m.v_large in var_set + assert m.v_small not in var_set + assert m.v_zero not in var_set + assert len(var_set) == 1 + var_set = variables_with_extreme_values_set( + m, large=1e5, small=1e-7, zero=1e-10, apply_scaling=False + ) + assert m.v_large in var_set + assert m.v_small in var_set + assert m.v_zero not in var_set + assert len(var_set) == 2 + + +@pytest.mark.unit +def test_number_variables_with_extreme_values(m): + m.v_large = Var(initialize=1e6) + m.v_small = Var(initialize=1e-8) + m.v_zero = Var(initialize=0) + assert ( + number_variables_with_extreme_values(m, large=1e5, small=1e-7, zero=1e-10) == 2 + ) + + set_scaling_factor(m.v_small, 1e8) + assert ( + number_variables_with_extreme_values( + m, large=1e5, small=1e-7, zero=1e-10, apply_scaling=True + ) + == 1 + ) + assert ( + number_variables_with_extreme_values( + m, large=1e5, small=1e-7, zero=1e-10, apply_scaling=False + ) + == 2 + ) + + @pytest.mark.unit def test_uninitialized_variables_in_activated_constraints(): m = ConcreteModel() @@ -1137,8 +1253,12 @@ def evaluate_equality_constraints(self): m = ConcreteModel() - m.gb = ExternalGreyBoxBlock(external_model=BasicGrayBox()) - m.gb_inactive = ExternalGreyBoxBlock(external_model=BasicGrayBox()) + m.gb = ExternalGreyBoxBlock( + external_model=BasicGrayBox(), build_implicit_constraint_objects=True + ) + m.gb_inactive = ExternalGreyBoxBlock( + external_model=BasicGrayBox(), build_implicit_constraint_objects=True + ) m.gb_inactive.deactivate() # test counting functions assert number_greybox_blocks(m) == 2 @@ -1170,7 +1290,9 @@ def evaluate_equality_constraints(self): assert degrees_of_freedom(m) == -1 assert number_variables_in_activated_constraints(m) == 7 assert number_total_constraints(m) == 5 - assert number_total_equalities(m) == 5 + # Total number of equalities = 8: 3 in each grey box plus 2 others + # Previous test value was incorrect - probably overlooked deactivated grey box + assert number_total_equalities(m) == 8 assert number_deactivated_equalities(m) == 3 assert number_deactivated_constraints(m) == 3 @@ -1367,3 +1489,239 @@ def test_number_active_variables_in_deactivated_blocks(m): @pytest.mark.unit def test_report_statistics(m): report_statistics(m) + + +@pytest.mark.unit +def test_external_variables_set(m): + m.b2["b"].cons_w_ext_var = Constraint(expr=m.b1.v1 == m.b2["b"].v1) + + assert len(external_variables_set(m)) == 0 + ext_vars = external_variables_set(m.b2) + assert len(ext_vars) == 1 + assert m.b1.v1 in ext_vars + + +@pytest.mark.unit +def test_number_external_variables(m): + m.b2["b"].cons_w_ext_var = Constraint(expr=m.b1.v1 == m.b2["b"].v1) + + assert number_external_variables(m) == 0 + assert number_external_variables(m.b2) == 1 + + +# ------------------------------------------------------------------------- +# Legacy tests for functions moved from model_diagnostics to model_statistics +class TestLegacyDiagnostics: + @pytest.fixture + def model(self): + m = ConcreteModel() + m.b = Block() + + m.v1 = Var() + m.v2 = Var() + m.v3 = Var() + m.v4 = Var() + + m.v1.fix(0) + m.v2.fix(3) + m.v3.set_value(0) + + return m + + @pytest.mark.unit + def test_vars_fixed_to_zero(self, model): + zero_vars = variables_fixed_to_zero_set(model) + assert isinstance(zero_vars, ComponentSet) + assert len(zero_vars) == 1 + for i in zero_vars: + assert i is model.v1 + + @pytest.mark.unit + def test_vars_near_zero(self, model): + model.v3.set_value(1e-5) + + near_zero_vars = variables_near_zero_set(model, tol=1e-5) + assert isinstance(near_zero_vars, ComponentSet) + assert len(near_zero_vars) == 2 + for i in near_zero_vars: + assert i.local_name in ["v1", "v3"] + + near_zero_vars = variables_near_zero_set(model, tol=1e-6) + assert isinstance(near_zero_vars, ComponentSet) + assert len(near_zero_vars) == 1 + for i in near_zero_vars: + assert i is model.v1 + + set_scaling_factor(model.v3, 1e5) + near_zero_vars = variables_near_zero_set(model, tol=1e-5) + assert isinstance(near_zero_vars, ComponentSet) + assert len(near_zero_vars) == 1 + for i in near_zero_vars: + assert i is model.v1 + + near_zero_vars = variables_near_zero_set(model, tol=1) + assert isinstance(near_zero_vars, ComponentSet) + assert len(near_zero_vars) == 2 + for i in near_zero_vars: + assert i.local_name in ["v1", "v3"] + + @pytest.mark.unit + def test_vars_near_zero_apply_scaling(self, model): + # v3 is 1e-4, set scaling so it appears near zero + model.v3.set_value(1e-4) + set_scaling_factor(model.v3, 1e2) + # With scaling, v3*sf = 1e-2 > tol, so not included + near_zero_vars = variables_near_zero_set(model, tol=1e-3, scale_variables=True) + assert model.v3 not in near_zero_vars + # Without scaling, v3 is included + near_zero_vars = variables_near_zero_set(model, tol=1e-3, scale_variables=False) + assert model.v3 in near_zero_vars + + @pytest.mark.unit + def test_vars_with_none_value(self, model): + none_value = variables_with_none_value_set(model) + + assert isinstance(none_value, ComponentSet) + assert len(none_value) == 1 + for i in none_value: + assert i is model.v4 + + @pytest.mark.unit + def test_vars_with_bounds_issues(self, model): + model.v1.setlb(2) + model.v1.setub(6) + model.v2.setlb(0) + model.v2.setub(10) + model.v4.set_value(10) + model.v4.setlb(0) + model.v4.setub(1) + + bounds_issue = variables_violating_bounds_set(model, abs_tol=0, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 2 + for i in bounds_issue: + assert i.local_name in ["v1", "v4"] + + m = ConcreteModel() + m.v = Var(initialize=-1e-4, bounds=(0, 1)) + + # Just outside lower bound + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e-6, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 1 + + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e3, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + set_scaling_factor(m.v, 1e-10, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e3, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + set_scaling_factor(m.v, 1e10, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e-6, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 1 + + # Just inside lower bound + m.v.set_value(1e-4) + + set_scaling_factor(m.v, 1, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e-6, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e3, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + set_scaling_factor(m.v, 1e-10, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e3, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + set_scaling_factor(m.v, 1e10, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e-6, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + # Just inside upper bound + m.v.set_value(1 - 1e-4) + + set_scaling_factor(m.v, 1, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e-6, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e3, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + set_scaling_factor(m.v, 1e-10, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e3, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + set_scaling_factor(m.v, 1e10, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e-6, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + # Just outside upper bound + m.v.set_value(1 + 1e-4) + + set_scaling_factor(m.v, 1, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e-6, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 1 + + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e3, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + set_scaling_factor(m.v, 1e-10, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e3, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 0 + + set_scaling_factor(m.v, 1e10, overwrite=True) + bounds_issue = variables_violating_bounds_set(m, abs_tol=1e-6, rel_tol=0) + assert isinstance(bounds_issue, ComponentSet) + assert len(bounds_issue) == 1 + + @pytest.mark.unit + def test_vars_with_extreme_values(self): + m = ConcreteModel() + m.v1 = Var(initialize=1e-12) # below zero + m.v2 = Var(initialize=1e-8) # small + m.v3 = Var(initialize=1e-4) + m.v4 = Var(initialize=1e0) + m.v5 = Var(initialize=1e4) + m.v6 = Var(initialize=1e8) + m.v7 = Var(initialize=1e12) # large + + xvars = variables_with_extreme_values_set(m, large=1e9, small=1e-7, zero=1e-10) + + assert len(xvars) == 2 + for i in xvars: + assert i.name in ["v2", "v7"] + + set_scaling_factor(m.v1, 1e3) + set_scaling_factor(m.v2, 1e6) + set_scaling_factor(m.v4, 1e-12) + set_scaling_factor(m.v6, 1e3) + set_scaling_factor(m.v7, 1e-12) + + xvars = variables_with_extreme_values_set(m, large=1e9, small=1e-7, zero=1e-10) + + assert len(xvars) == 2 + for i in xvars: + assert i.name in ["v1", "v6"] + + # Test apply_scaling argument + # With apply_scaling=False, only v2 and v7 should be included + xvars_noscale = variables_with_extreme_values_set( + m, large=1e9, small=1e-7, zero=1e-10, apply_scaling=False + ) + assert set([v.name for v in xvars_noscale]) == {"v2", "v7"} diff --git a/idaes/core/util/tests/test_model_statistics_grey_box.py b/idaes/core/util/tests/test_model_statistics_grey_box.py new file mode 100644 index 0000000000..cc8194604e --- /dev/null +++ b/idaes/core/util/tests/test_model_statistics_grey_box.py @@ -0,0 +1,2439 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2026 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +This module contains tests for the model statistics functions with the +presence of grey-box components. +""" + +import pytest + +import pyomo.environ as pyo +from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock +from pyomo.contrib.pynumero.interfaces.tests.external_grey_box_models import ( + PressureDropTwoEqualitiesTwoOutputsWithHessian, +) + +from idaes.core.util.model_statistics import * + + +@pytest.fixture +def model(): + m = pyo.ConcreteModel() + + # Add a Block containing a Grey Box model, with linking variables and constraints + m.b1 = pyo.Block() + m.b1.egb = ExternalGreyBoxBlock() + m.b1.egb.set_external_model( + PressureDropTwoEqualitiesTwoOutputsWithHessian(), + build_implicit_constraint_objects=True, + ) + + # Add Vars and linking constraints to m + m.b1.Pin = pyo.Var(initialize=101325, bounds=(0, 1e8)) # Not near bounds + m.b1.c = pyo.Var(initialize=1e-8, bounds=(0, 1)) # Near bounds + m.b1.F = pyo.Var() + m.b1.P1 = pyo.Var() + m.b1.P3 = pyo.Var() + m.b1.P2 = pyo.Var() + m.b1.Pout = pyo.Var() + + m.b1.link_Pin = pyo.Constraint(expr=m.b1.Pin == m.b1.egb.inputs["Pin"]) + m.b1.link_c = pyo.Constraint(expr=m.b1.c == m.b1.egb.inputs["c"]) + m.b1.link_F = pyo.Constraint(expr=m.b1.F == m.b1.egb.inputs["F"]) + m.b1.link_P1 = pyo.Constraint(expr=m.b1.P1 == m.b1.egb.inputs["P1"]) + m.b1.link_P3 = pyo.Constraint(expr=m.b1.P3 == m.b1.egb.inputs["P3"]) + m.b1.link_P2 = pyo.Constraint(expr=m.b1.P2 == m.b1.egb.outputs["P2"]) + m.b1.link_Pout = pyo.Constraint(expr=m.b1.Pout == m.b1.egb.outputs["Pout"]) + + # Add a second, unrelated Block to ensure that the model statistics + # functions are correctly identifying blocks and constraints in the presence of multiple blocks + m.b2 = pyo.Block() + m.b2.v1 = pyo.Var() + m.b2.c1 = pyo.Constraint(expr=m.b2.v1 == 1) + + # Add two inequalities and an objective to confirm behaviour + m.b1.ineq = pyo.Constraint(expr=m.b1.Pin >= 0) + m.b2.ineq = pyo.Constraint(expr=m.b2.v1 >= 0) + m.obj = pyo.Objective(expr=m.b1.Pout) + + # Add an Expression which uses variables from the Grey Box + m.expr = pyo.Expression(expr=m.b1.egb.inputs["Pin"] * m.b1.egb.inputs["c"]) + + # Set some values and bounds in the grey box to confirm that grey box + # variables are included in the variable counts + # Include both inputs and outputs to confirm that both are included + # in the variable counts + m.b1.egb.inputs["Pin"].set_value(101325) + m.b1.egb.inputs["Pin"].setlb(0) + m.b1.egb.inputs["Pin"].setub(1e8) + m.b1.egb.inputs["c"].set_value(1e-8) + m.b1.egb.inputs["c"].setlb(0) + m.b1.egb.inputs["c"].setub(1) + m.b1.egb.outputs["P2"].set_value(90000) + m.b1.egb.outputs["P2"].setlb(90000) + m.b1.egb.outputs["P2"].setub(1e8) + + return m + + +class TestBlockStatisticsGreyBox: + @pytest.mark.unit + def test_total_blocks_set_w_grey_box(self, model): + # Test that the total_blocks_set function correctly counts + # the number of blocks in the model + # Grey Box is not included as a normal block + assert len(total_blocks_set(model)) == 3 + for b in total_blocks_set(model): + assert b in [model, model.b1, model.b2] + + @pytest.mark.unit + def test_number_total_blocks_w_grey_box(self, model): + assert number_total_blocks(model) == 3 + + @pytest.mark.unit + def test_activated_blocks_set_w_grey_box(self, model): + # Test that the activated_blocks_set function correctly counts + # the number of activated blocks in the model + # Grey Box is not included as a normal block + assert len(activated_blocks_set(model)) == 3 + for b in activated_blocks_set(model): + assert b in [model, model.b1, model.b2] + + # Deactivate b1 and test again + model.b1.deactivate() + assert len(activated_blocks_set(model)) == 2 + for b in activated_blocks_set(model): + assert b in [model, model.b2] + + @pytest.mark.unit + def test_greybox_block_set_w_grey_box(self, model): + # Test that the grey_box_set function correctly identifies + # the Grey Box block in the model + gbs = greybox_block_set(model) + assert len(gbs) == 1 + assert model.b1.egb in gbs + + @pytest.mark.unit + def test_activated_greybox_block_set_w_grey_box(self, model): + # Test that the activated_greybox_block_set function correctly identifies + # the activated Grey Box block in the model + agbs = activated_greybox_block_set(model) + assert len(agbs) == 1 + assert model.b1.egb in agbs + + # Deactivate b1 and test again + model.b1.deactivate() + agbs = activated_greybox_block_set(model) + assert len(agbs) == 0 + + @pytest.mark.unit + def test_deactivated_greybox_block_set_w_grey_box(self, model): + # Test that the deactivated_greybox_block_set function correctly identifies + # the deactivated Grey Box block in the model + dgb = deactivated_greybox_block_set(model) + assert len(dgb) == 0 + + # Deactivate the grey box and test again + model.b1.egb.deactivate() + dgb = deactivated_greybox_block_set(model) + assert len(dgb) == 1 + assert model.b1.egb in dgb + + @pytest.mark.unit + def test_number_deactivated_greybox_block_w_grey_box(self, model): + # Test that the number_deactivated_greybox_block function correctly counts + # the number of deactivated Grey Box blocks in the model + assert number_deactivated_greybox_block(model) == 0 + + # Deactivate the grey box and test again + model.b1.egb.deactivate() + assert number_deactivated_greybox_block(model) == 1 + + @pytest.mark.unit + def test_number_greybox_blocks_w_grey_box(self, model): + # Test that the number_greybox_blocks function correctly counts + # the number of Grey Box blocks in the model + assert number_greybox_blocks(model) == 1 + + # Deactivate the grey box and test again (should not change the count) + model.b1.egb.deactivate() + assert number_greybox_blocks(model) == 1 + + @pytest.mark.unit + def test_number_activated_greybox_blocks_w_grey_box(self, model): + # Test that the number_activated_greybox_blocks function correctly counts + # the number of activated Grey Box blocks in the model + assert number_activated_greybox_blocks(model) == 1 + + # Deactivate the grey box and test again + model.b1.egb.deactivate() + assert number_activated_greybox_blocks(model) == 0 + + @pytest.mark.unit + def test_number_activated_blocks_w_grey_box(self, model): + # Test that the number_activated_blocks function correctly counts + # the number of activated blocks in the model + assert number_activated_blocks(model) == 3 + + # Deactivate b2 and test again + model.b2.deactivate() + assert number_activated_blocks(model) == 2 + + @pytest.mark.unit + def test_deactivated_blocks_set_w_grey_box(self, model): + # Test that the deactivated_blocks_set function correctly identifies + # the deactivated blocks in the model + dbs = deactivated_blocks_set(model) + assert len(dbs) == 0 + + # Deactivate b2 and test again + model.b2.deactivate() + dbs = deactivated_blocks_set(model) + assert len(dbs) == 1 + assert model.b2 in dbs + + @pytest.mark.unit + def test_number_deactivated_blocks_w_grey_box(self, model): + # Test that the number_deactivated_blocks function correctly counts the + # number of deactivated blocks in the model + assert number_deactivated_blocks(model) == 0 + + # Deactivate b2 and test again + model.b2.deactivate() + assert number_deactivated_blocks(model) == 1 + + +class TestConstraintStatisticsGreyBox: + @pytest.mark.unit + def test_total_constraints_set_w_grey_box(self, model): + # Test that the total_constraints_set function correctly counts the + # number of constraints in the model + # First, test with include_greybox = False + tcs = total_constraints_set(model, include_greybox=False) + assert len(tcs) == 10 + for c in tcs: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.ineq, + model.b2.ineq, + ] + + # Next, test with include_greybox = True + tcs = total_constraints_set(model, include_greybox=True) + assert len(tcs) == 14 + for c in tcs: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + model.b1.ineq, + model.b2.ineq, + ] + + @pytest.mark.unit + def test_number_total_constraints_w_grey_box(self, model): + # Test that the number_total_constraints function correctly counts the + # number of constraints in the model + # First, test with include_greybox = False + assert number_total_constraints(model, include_greybox=False) == 10 + + # Next, test with include_greybox = True + assert number_total_constraints(model, include_greybox=True) == 14 + + @pytest.mark.unit + def test_activated_constraints_generator_w_grey_box(self, model): + # Test that the activated_constraints_generator function correctly identifies + # the activated constraints in the model + # First, test with include_greybox = False + acg = activated_constraints_generator(model, include_greybox=False) + acg_list = list(acg) + assert len(acg_list) == 10 + for c in acg_list: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.ineq, + model.b2.ineq, + ] + + # Next, test with include_greybox = True + acg = activated_constraints_generator(model, include_greybox=True) + acg_list = list(acg) + assert len(acg_list) == 14 + for c in acg_list: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + model.b1.ineq, + model.b2.ineq, + ] + + # Now deactivate the grey box and test again with include_greybox = True + # (should not include grey box constraints) + model.b1.egb.deactivate() + acg = activated_constraints_generator(model, include_greybox=True) + acg_list = list(acg) + assert len(acg_list) == 10 + for c in acg_list: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.ineq, + model.b2.ineq, + ] + + @pytest.mark.unit + def test_activated_constraints_set_w_grey_box(self, model): + # Test that the activated_constraints_set function correctly identifies + # the activated constraints in the model + # First, test with include_greybox = False + acs = activated_constraints_set(model, include_greybox=False) + assert len(acs) == 10 + for c in acs: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.ineq, + model.b2.ineq, + ] + + # Next, test with include_greybox = True + acs = activated_constraints_set(model, include_greybox=True) + assert len(acs) == 14 + for c in acs: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + model.b1.ineq, + model.b2.ineq, + ] + + # Now deactivate the grey box and test again with include_greybox = True + # (should not include grey box constraints) + model.b1.egb.deactivate() + acs = activated_constraints_set(model, include_greybox=True) + assert len(acs) == 10 + for c in acs: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.ineq, + model.b2.ineq, + ] + + @pytest.mark.unit + def test_number_activated_constraints_w_grey_box(self, model): + # Test that the number_activated_constraints function correctly counts the number + # of activated constraints in the model + # First, test with include_greybox = False + assert number_activated_constraints(model, include_greybox=False) == 10 + + # Next, test with include_greybox = True + assert number_activated_constraints(model, include_greybox=True) == 14 + + # Now deactivate the grey box and test again with include_greybox = True + # (should not include grey box constraints) + model.b1.egb.deactivate() + assert number_activated_constraints(model, include_greybox=True) == 10 + + @pytest.mark.unit + def test_deactivated_constraints_generator_w_grey_box(self, model): + # Test that the deactivated_constraints_generator function correctly identifies + # the deactivated constraints in the model + # First, test with include_greybox = False + dcg = deactivated_constraints_generator(model, include_greybox=False) + dcg_list = list(dcg) + assert len(dcg_list) == 0 + + # Next, test with include_greybox = True + dcg = deactivated_constraints_generator(model, include_greybox=True) + dcg_list = list(dcg) + assert len(dcg_list) == 0 + + # Now deactivate the grey box and test again with include_greybox = True + # (should include grey box constraints) + model.b1.egb.deactivate() + dcg = deactivated_constraints_generator(model, include_greybox=True) + dcg_list = list(dcg) + assert len(dcg_list) == 4 + for c in dcg_list: + assert c in [ + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + ] + + @pytest.mark.unit + def test_deactivated_constraints_set_w_grey_box(self, model): + # Test that the deactivated_constraints_set function correctly identifies + # the deactivated constraints in the model + # First, test with include_greybox = False + dcs = deactivated_constraints_set(model, include_greybox=False) + assert len(dcs) == 0 + + # Next, test with include_greybox = True + dcs = deactivated_constraints_set(model, include_greybox=True) + assert len(dcs) == 0 + + # Now deactivate the grey box and test again with include_greybox = True + # (should include grey box constraints) + model.b1.egb.deactivate() + dcs = deactivated_constraints_set(model, include_greybox=True) + assert len(dcs) == 4 + for c in dcs: + assert c in [ + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + ] + + @pytest.mark.unit + def test_number_deactivated_constraints_w_grey_box(self, model): + # Test that the number_deactivated_constraints function correctly counts + # the number of deactivated constraints in the model + # First, test with include_greybox = False + assert number_deactivated_constraints(model, include_greybox=False) == 0 + + # Next, test with include_greybox = True + assert number_deactivated_constraints(model, include_greybox=True) == 0 + + # Now deactivate the grey box and test again with include_greybox = True + # (should include grey box constraints) + model.b1.egb.deactivate() + assert number_deactivated_constraints(model, include_greybox=True) == 4 + + @pytest.mark.unit + def test_total_equalities_generator_w_grey_box(self, model): + # Test that the total_equalities_generator function correctly identifies + # the equality constraints in the model + # First, test with include_greybox = False + teg = total_equalities_generator(model, include_greybox=False) + teg_list = list(teg) + assert len(teg_list) == 8 + for c in teg_list: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + ] + + # Next, test with include_greybox = True + teg = total_equalities_generator(model, include_greybox=True) + teg_list = list(teg) + assert len(teg_list) == 12 + for c in teg_list: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + model.b2.c1, + ] + + @pytest.mark.unit + def test_total_equalities_set_w_grey_box(self, model): + # Test that the total_equalities_set function correctly identifies the + # equality constraints in the model + # First, test with include_greybox = False + tes = total_equalities_set(model, include_greybox=False) + assert len(tes) == 8 + for c in tes: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + ] + + # Next, test with include_greybox = True + tes = total_equalities_set(model, include_greybox=True) + assert len(tes) == 12 + for c in tes: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + ] + + # Now deactivate the grey box and test again with include_greybox = True + # Greybox constraints should still be included in the count as we are not checking active + model.b1.egb.deactivate() + tes = total_equalities_set(model, include_greybox=True) + assert len(tes) == 12 + for c in tes: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + ] + + @pytest.mark.unit + def test_number_total_equalities_w_grey_box(self, model): + # Test that the number_total_equalities function correctly counts the number of + # equality constraints in the model + # First, test with include_greybox = False + assert number_total_equalities(model, include_greybox=False) == 8 + + # Next, test with include_greybox = True + assert number_total_equalities(model, include_greybox=True) == 12 + + # Now deactivate the grey box and test again with include_greybox = True + # Greybox constraints should still be included in the count as we are not checking active + model.b1.egb.deactivate() + assert number_total_equalities(model, include_greybox=True) == 12 + + @pytest.mark.unit + def test_activated_equalities_generator_w_grey_box(self, model): + # Test that the activated_equalities_generator function correctly identifies + # the activated equality constraints in the model + # First, test with include_greybox = False + aeg = activated_equalities_generator(model, include_greybox=False) + aeg_list = list(aeg) + assert len(aeg_list) == 8 + for c in aeg_list: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + ] + + # Next, test with include_greybox = True + aeg = activated_equalities_generator(model, include_greybox=True) + aeg_list = list(aeg) + assert len(aeg_list) == 12 + for c in aeg_list: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + ] + + # Now deactivate the grey box and test again with include_greybox = True + # (should not include grey box constraints) + model.b1.egb.deactivate() + aeg = activated_equalities_generator(model, include_greybox=True) + aeg_list = list(aeg) + assert len(aeg_list) == 8 + for c in aeg_list: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + ] + + @pytest.mark.unit + def test_activated_equalities_set_w_grey_box(self, model): + # Test that the activated_equalities_set function correctly identifies + # the activated equality constraints in the model + # First, test with include_greybox = False + aes = activated_equalities_set(model, include_greybox=False) + assert len(aes) == 8 + for c in aes: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + ] + + # Next, test with include_greybox = True + aes = activated_equalities_set(model, include_greybox=True) + assert len(aes) == 12 + for c in aes: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + ] + + # Now deactivate the grey box and test again with include_greybox = True + # (should not include grey box constraints) + model.b1.egb.deactivate() + aes = activated_equalities_set(model, include_greybox=True) + assert len(aes) == 8 + for c in aes: + assert c in [ + model.b1.link_Pin, + model.b1.link_c, + model.b1.link_F, + model.b1.link_P1, + model.b1.link_P3, + model.b1.link_P2, + model.b1.link_Pout, + model.b2.c1, + ] + + @pytest.mark.unit + def test_number_activated_equalities_w_grey_box(self, model): + # Test that the number_activated_equalities function correctly counts the + # number of activated equality constraints in the model + # First, test with include_greybox = False + assert number_activated_equalities(model, include_greybox=False) == 8 + + # Next, test with include_greybox = True + assert number_activated_equalities(model, include_greybox=True) == 12 + + # Now deactivate the grey box and test again with include_greybox = True + # (should not include grey box constraints) + model.b1.egb.deactivate() + assert number_activated_equalities(model, include_greybox=True) == 8 + + @pytest.mark.unit + def test_number_activated_greybox_equalities_w_grey_box(self, model): + # Test that the number_activated_greybox_equalities function correctly counts + # the number of activated equality constraints in the Grey Box blocks in the model + assert number_activated_greybox_equalities(model) == 4 + + # Now deactivate the grey box and test again (should be 0) + model.b1.egb.deactivate() + assert number_activated_greybox_equalities(model) == 0 + + @pytest.mark.unit + def test_number_deactivated_equalities_w_grey_box(self, model): + # Test that the number_deactivated_equalities function correctly counts + # the number of deactivated equality constraints in the Grey Box blocks in the model + assert number_deactivated_greybox_equalities(model) == 0 + + # Now deactivate the grey box and test again (should be 4) + model.b1.egb.deactivate() + assert number_deactivated_greybox_equalities(model) == 4 + + @pytest.mark.unit + def test_deactivated_equalities_generator_w_grey_box(self, model): + # Test that the deactivated_equalities_generator function correctly + # identifies the deactivated equality constraints in the Grey Box blocks in the model + deg = deactivated_equalities_generator(model, include_greybox=True) + deg_list = list(deg) + assert len(deg_list) == 0 + + # Now deactivate the grey box and test again (should include grey box constraints) + model.b1.egb.deactivate() + deg = deactivated_equalities_generator(model, include_greybox=True) + deg_list = list(deg) + assert len(deg_list) == 4 + for c in deg_list: + assert c in [ + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + ] + + @pytest.mark.unit + def test_deactivated_equalities_set_w_grey_box(self, model): + # Test that the deactivated_equalities_set function correctly identifies + # the deactivated equality constraints in the Grey Box blocks in the model + decs = deactivated_equalities_set(model, include_greybox=True) + assert len(decs) == 0 + + # Now deactivate the grey box and test again (should include grey box constraints) + model.b1.egb.deactivate() + decs = deactivated_equalities_set(model, include_greybox=True) + assert len(decs) == 4 + for c in decs: + assert c in [ + model.b1.egb.P2_constraint, + model.b1.egb.Pout_constraint, + model.b1.egb.pdrop1, + model.b1.egb.pdrop3, + ] + + @pytest.mark.unit + def test_number_deactivated_greybox_equalities_w_grey_box(self, model): + # Test that the number_deactivated_greybox_equalities function correctly + # counts the number of deactivated equality constraints in the Grey Box blocks in the model + assert number_deactivated_greybox_equalities(model) == 0 + + # Now deactivate the grey box and test again (should be 4) + model.b1.egb.deactivate() + assert number_deactivated_greybox_equalities(model) == 4 + + # Check inequality methods to ensure they work wit ha grey box present + @pytest.mark.unit + def test_total_inequalities_generator_w_grey_box(self, model): + tig = total_inequalities_generator(model) + tig_list = list(tig) + assert len(tig_list) == 2 + for c in tig_list: + assert c in [model.b1.ineq, model.b2.ineq] + + @pytest.mark.unit + def test_total_inequalities_set_w_grey_box(self, model): + tis = total_inequalities_set(model) + assert len(tis) == 2 + for c in tis: + assert c in [model.b1.ineq, model.b2.ineq] + + @pytest.mark.unit + def test_number_total_inequalities_w_grey_box(self, model): + assert number_total_inequalities(model) == 2 + + @pytest.mark.unit + def test_activated_inequalities_generator_w_grey_box(self, model): + aig = activated_inequalities_generator(model) + aig_list = list(aig) + assert len(aig_list) == 2 + for c in aig_list: + assert c in [model.b1.ineq, model.b2.ineq] + + @pytest.mark.unit + def test_activated_inequalities_set_w_grey_box(self, model): + aig = activated_inequalities_set(model) + assert len(aig) == 2 + for c in aig: + assert c in [model.b1.ineq, model.b2.ineq] + + @pytest.mark.unit + def test_number_activated_inequalities_w_grey_box(self, model): + assert number_activated_inequalities(model) == 2 + + @pytest.mark.unit + def test_deactivated_inequalities_generator_w_grey_box(self, model): + dig = deactivated_inequalities_generator(model) + dig_list = list(dig) + assert len(dig_list) == 0 + + # Deactivate one of the inequalities and test again + model.b1.ineq.deactivate() + dig = deactivated_inequalities_generator(model) + dig_list = list(dig) + assert len(dig_list) == 1 + for c in dig_list: + assert c in [model.b1.ineq] + + @pytest.mark.unit + def test_deactivated_inequalities_set_w_grey_box(self, model): + dis = deactivated_inequalities_set(model) + assert len(dis) == 0 + + # Deactivate one of the inequalities and test again + model.b1.ineq.deactivate() + dis = deactivated_inequalities_set(model) + assert len(dis) == 1 + for c in dis: + assert c in [model.b1.ineq] + + @pytest.mark.unit + def test_number_deactivated_inequalities_w_grey_box(self, model): + assert number_deactivated_inequalities(model) == 0 + + # Deactivate one of the inequalities and test again + model.b1.ineq.deactivate() + assert number_deactivated_inequalities(model) == 1 + + +class TestVariableStatisticsGreyBox: + @pytest.mark.unit + def test_variables_generator_w_grey_box(self, model): + # Start with include_greybox = False to confirm that we are not + # including grey box variables + vg = variables_generator(model, include_greybox=False) + vg_list = list(vg) + assert len(vg_list) == 8 + for v in vg_list: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + ] + + # Now test with include_greybox = True to confirm that we are + # including grey box variables + vg = variables_generator(model, include_greybox=True) + vg_list = list(vg) + assert len(vg_list) == 15 + for v in vg_list: + print(v.name) + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + @pytest.mark.unit + def test_variables_set_w_grey_box(self, model): + # Start with include_greybox = False to confirm that we are not + # including grey box variables + vs = variables_set(model, include_greybox=False) + assert len(vs) == 8 + for v in vs: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + ] + + # Now test with include_greybox = True to confirm that we are + # including grey box variables + vs = variables_set(model, include_greybox=True) + assert len(vs) == 15 + for v in vs: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + @pytest.mark.unit + def test_number_variables_w_grey_box(self, model): + # Start with include_greybox = False to confirm that we are not + # including grey box variables + assert number_variables(model, include_greybox=False) == 8 + + # Now test with include_greybox = True to confirm that we are + # including grey box variables + assert number_variables(model, include_greybox=True) == 15 + + @pytest.mark.unit + def test_fixed_variables_generator(self, model): + # Test that the fixed_variables_generator function correctly + # identifies the fixed variables in the model + fvg = fixed_variables_generator(model) + fvg_list = list(fvg) + assert len(fvg_list) == 0 + + # Fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + fvg = fixed_variables_generator(model) + fvg_list = list(fvg) + assert len(fvg_list) == 2 + for v in fvg_list: + assert v.name in ["b1.Pin", "b1.egb.inputs[Pin]"] + + # Test again with include_greybox = False to confirm that we are not including grey box variables + fvg = fixed_variables_generator(model, include_greybox=False) + fvg_list = list(fvg) + assert len(fvg_list) == 1 + for v in fvg_list: + assert v.name in ["b1.Pin"] + + @pytest.mark.unit + def test_fixed_variables_set(self, model): + # Test that the fixed_variables_set function correctly identifies + # the fixed variables in the model + fvs = fixed_variables_set(model) + assert len(fvs) == 0 + + # Fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + fvs = fixed_variables_set(model) + assert len(fvs) == 2 + for v in fvs: + assert v.name in ["b1.Pin", "b1.egb.inputs[Pin]"] + + # Test again with include_greybox = False to confirm that we are + # not including grey box variables + fvs = fixed_variables_set(model, include_greybox=False) + assert len(fvs) == 1 + for v in fvs: + assert v.name in ["b1.Pin"] + + @pytest.mark.unit + def test_number_fixed_variables(self, model): + # Test that the number_fixed_variables function correctly counts + # the number of fixed variables in the model + assert number_fixed_variables(model) == 0 + + # Fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + assert number_fixed_variables(model) == 2 + + # Test again with include_greybox = False to confirm that we are + # not including grey box variables + assert number_fixed_variables(model, include_greybox=False) == 1 + + @pytest.mark.unit + def test_unfixed_variables_generator(self, model): + # Test that the unfixed_variables_generator function correctly + # identifies the unfixed variables in the model + uv = unfixed_variables_generator(model) + uv_list = list(uv) + assert len(uv_list) == 15 + for v in uv_list: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Now fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + uv = unfixed_variables_generator(model) + uv_list = list(uv) + assert len(uv_list) == 13 + for v in uv_list: + assert v.name in [ + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Test again with include_greybox = False to confirm that we are + # not including grey box variables + uv = unfixed_variables_generator(model, include_greybox=False) + uv_list = list(uv) + assert len(uv_list) == 7 + for v in uv_list: + assert v.name in [ + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + ] + + @pytest.mark.unit + def test_unfixed_variables_set(self, model): + # Test that the unfixed_variables_set function correctly identifies + # the unfixed variables in the model + uvs = unfixed_variables_set(model) + assert len(uvs) == 15 + for v in uvs: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Now fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + uvs = unfixed_variables_set(model) + assert len(uvs) == 13 + for v in uvs: + assert v.name in [ + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Test again with include_greybox = False to confirm that we are + # not including grey box variables + uvs = unfixed_variables_set(model, include_greybox=False) + assert len(uvs) == 7 + for v in uvs: + assert v.name in [ + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + ] + + @pytest.mark.unit + def test_number_unfixed_variables(self, model): + # Test that the number_unfixed_variables function correctly counts the + # number of unfixed variables in the model + assert number_unfixed_variables(model) == 15 + + # Now fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + assert number_unfixed_variables(model) == 13 + + # Test again with include_greybox = False to confirm that we are not + # including grey box variables + assert number_unfixed_variables(model, include_greybox=False) == 7 + + @pytest.mark.unit + def test_variables_near_bounds_generator(self, model): + # Test that the variables_near_bounds_generator function correctly identifies + # the variables that are near their bounds in the model + vnbg = variables_near_bounds_generator(model, tol=1e-6) + vnbg_list = list(vnbg) + assert len(vnbg_list) == 3 + for v in vnbg_list: + assert v.name in ["b1.c", "b1.egb.inputs[c]", "b1.egb.outputs[P2]"] + + # Now test with include_greybox = False to confirm that we are not + # including grey box variables + vnbg = variables_near_bounds_generator(model, tol=1e-6, include_greybox=False) + vnbg_list = list(vnbg) + assert len(vnbg_list) == 1 + for v in vnbg_list: + assert v.name in ["b1.c"] + + @pytest.mark.unit + def test_variables_near_bounds_set(self, model): + # Test that the variables_near_bounds_set function correctly identifies the + # variables that are near their bounds in the model + vnbs = variables_near_bounds_set(model, tol=1e-6) + assert len(vnbs) == 3 + for v in vnbs: + assert v.name in ["b1.c", "b1.egb.inputs[c]", "b1.egb.outputs[P2]"] + + # Now test with include_greybox = False to confirm that we are not including grey box variables + vnbs = variables_near_bounds_set(model, tol=1e-6, include_greybox=False) + assert len(vnbs) == 1 + for v in vnbs: + assert v.name in ["b1.c"] + + @pytest.mark.unit + def test_number_variables_near_bounds(self, model): + # Test that the number_variables_near_bounds function correctly counts the + # number of variables that are near their bounds in the model + assert number_variables_near_bounds(model, tol=1e-6) == 3 + + # Now test with include_greybox = False to confirm that we are not including grey box variables + assert number_variables_near_bounds(model, tol=1e-6, include_greybox=False) == 1 + + @pytest.mark.unit + def test_variables_in_activated_constraints_set_w_grey_box(self, model): + # Test that the variables_in_activated_constraints_set function correctly + # identifies the variables that are in the activated constraints in the model + # First, test with include_greybox = False + # We should see all the variables, including those in the grey box, as they appear + # in the linking constraints. + vics = variables_in_activated_constraints_set(model, include_greybox=False) + assert len(vics) == 15 + for v in vics: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Next, deactivate the grey box - result should still be the same + model.b1.egb.deactivate() + vics = variables_in_activated_constraints_set(model, include_greybox=True) + assert len(vics) == 15 + for v in vics: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Now, turn off a linking constraint + model.b1.link_Pin.deactivate() + vics = variables_in_activated_constraints_set(model, include_greybox=True) + assert len(vics) == 14 + for v in vics: + assert v.name in [ + "b1.Pin", # This is in the inequality constraint + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + # "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Finally, turn the grey box back on + # All variables should be back, as the grey box variables are all considered + # to be in active constraints if the grey box is active + model.b1.egb.activate() + vics = variables_in_activated_constraints_set(model, include_greybox=True) + assert len(vics) == 15 + for v in vics: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + @pytest.mark.unit + def test_number_variables_in_activated_constraints_w_grey_box(self, model): + # Test that the number_variables_in_activated_constraints function correctly counts + # the number of variables that are in the activated constraints in the model + # First, test with include_greybox = False + assert ( + number_variables_in_activated_constraints(model, include_greybox=False) + == 15 + ) + + # Next, deactivate the grey box - result should still be the same + model.b1.egb.deactivate() + assert ( + number_variables_in_activated_constraints(model, include_greybox=True) == 15 + ) + + # Now, turn off a linking constraint + model.b1.link_Pin.deactivate() + assert ( + number_variables_in_activated_constraints(model, include_greybox=True) == 14 + ) + + # Finally, turn the grey box back on - should go back to 15 + model.b1.egb.activate() + assert ( + number_variables_in_activated_constraints(model, include_greybox=True) == 15 + ) + + @pytest.mark.unit + def test_variables_not_in_activated_constraints_set_w_grey_box(self, model): + # Test that the variables_not_in_activated_constraints_set function correctly identifies + # the variables that are not in the activated constraints in the model + # First, test with include_greybox = False + vniacs = variables_not_in_activated_constraints_set( + model, include_greybox=False + ) + assert len(vniacs) == 0 + + # Next, deactivate the grey box - result should still be the same as all vars appear + # in active linking constraints + model.b1.egb.deactivate() + vniacs = variables_not_in_activated_constraints_set(model, include_greybox=True) + assert len(vniacs) == 0 + + # Now, turn off a linking constraint - should see the Pin variable from the grey box as not + # being in an active constraint + model.b1.link_Pin.deactivate() + vniacs = variables_not_in_activated_constraints_set(model, include_greybox=True) + + assert len(vniacs) == 1 + for v in vniacs: + assert v.name in ["b1.egb.inputs[Pin]"] + + # Finally, turn the grey box back on - should go back to 0 + model.b1.egb.activate() + vniacs = variables_not_in_activated_constraints_set(model, include_greybox=True) + assert len(vniacs) == 0 + + @pytest.mark.unit + def test_number_variables_not_in_activated_constraints_w_grey_box(self, model): + # Test that the number_variables_not_in_activated_constraints function correctly counts + # the number of variables that are not in the activated constraints in the model + # First, test with include_greybox = False + assert ( + number_variables_not_in_activated_constraints(model, include_greybox=False) + == 0 + ) + + # Next, deactivate the grey box - result should still be the same + model.b1.egb.deactivate() + assert ( + number_variables_not_in_activated_constraints(model, include_greybox=True) + == 0 + ) + + # Now, turn off a linking constraint - should see the Pin variable from the grey box as not + # being in an active constraint + model.b1.link_Pin.deactivate() + assert ( + number_variables_not_in_activated_constraints(model, include_greybox=True) + == 1 + ) + + # Finally, turn the grey box back on - should go back to 0 + model.b1.egb.activate() + assert ( + number_variables_not_in_activated_constraints(model, include_greybox=True) + == 0 + ) + + @pytest.mark.unit + def test_variables_in_activated_equalities_set_w_grey_box(self, model): + # Test that the variables_in_activated_equalities_set function correctly identifies + # the variables that are in the activated equality constraints in the model + # First, test with include_greybox = False + # We should see all the variables, including those in the grey box, as they appear + # in the linking constraints. + vies = variables_in_activated_equalities_set(model, include_greybox=False) + assert len(vies) == 15 + for v in vies: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Next, deactivate the grey box - result should still be the same + model.b1.egb.deactivate() + vies = variables_in_activated_equalities_set(model, include_greybox=True) + assert len(vies) == 15 + for v in vies: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Now, turn off a linking constraint + model.b1.link_Pin.deactivate() + vies = variables_in_activated_equalities_set(model, include_greybox=True) + assert len(vies) == 13 + for v in vies: + assert v.name in [ + # "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + # "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Finally, turn the grey box back on + # GreyBox variables should be back, but b1.Pin should still be out as the linking constraint + # is still deactivated + model.b1.egb.activate() + vies = variables_in_activated_equalities_set(model, include_greybox=True) + assert len(vies) == 14 + for v in vies: + assert v.name in [ + # "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + @pytest.mark.unit + def test_number_variables_in_activated_equalities_w_grey_box(self, model): + # Test that the number_variables_in_activated_equalities function correctly counts + # the number of variables that are in the activated equality constraints in the model + # First, test with include_greybox = False + assert ( + number_variables_in_activated_equalities(model, include_greybox=False) == 15 + ) + + # Next, deactivate the grey box - result should still be the same + model.b1.egb.deactivate() + assert ( + number_variables_in_activated_equalities(model, include_greybox=True) == 15 + ) + + # Now, turn off a linking constraint + model.b1.link_Pin.deactivate() + assert ( + number_variables_in_activated_equalities(model, include_greybox=True) == 13 + ) + + # Finally, turn the grey box back on - should go back to 14 as only b1.Pin is out + model.b1.egb.activate() + assert ( + number_variables_in_activated_equalities(model, include_greybox=True) == 14 + ) + + @pytest.mark.unit + def test_variables_in_activated_inequalities_set_w_grey_box(self, model): + # Test that the variables_in_activated_inequalities_set function correctly identifies + # the variables that are in the activated inequality constraints in the model + vies = variables_in_activated_inequalities_set(model) + assert len(vies) == 2 + for v in vies: + assert v.name in ["b1.Pin", "b2.v1"] + + @pytest.mark.unit + def test_number_variables_in_activated_inequalities_w_grey_box(self, model): + # Test that the number_variables_in_activated_inequalities function correctly counts + # the number of variables that are in the activated inequality constraints in the model + assert number_variables_in_activated_inequalities(model) == 2 + + @pytest.mark.unit + def test_variables_only_in_inequalities_w_grey_box(self, model): + # Test that the variables_only_in_inequalities function correctly identifies + # the variables that are only in the inequality constraints in the model + vois = variables_only_in_inequalities(model) + assert len(vois) == 0 + + @pytest.mark.unit + def test_number_variables_only_in_inequalities_w_grey_box(self, model): + # Test that the number_variables_only_in_inequalities function correctly counts + # the number of variables that are only in the inequality constraints in the model + assert number_variables_only_in_inequalities(model) == 0 + + @pytest.mark.unit + def test_fixed_variables_in_activated_equalities_set_w_grey_box(self, model): + # Test that the fixed_variables_in_activated_equalities_set function correctly identifies + # the fixed variables that are in the activated equality constraints in the model + fvaes = fixed_variables_in_activated_equalities_set(model) + assert len(fvaes) == 0 + + # Now fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + fvaes = fixed_variables_in_activated_equalities_set(model) + assert len(fvaes) == 2 + for v in fvaes: + assert v.name in ["b1.Pin", "b1.egb.inputs[Pin]"] + + # Deactivate the GreyBox and linking constraint + model.b1.egb.deactivate() + model.b1.link_Pin.deactivate() + fvaes = fixed_variables_in_activated_equalities_set(model) + assert len(fvaes) == 0 + + @pytest.mark.unit + def test_number_fixed_variables_in_activated_equalities_w_grey_box(self, model): + # Test that the number_fixed_variables_in_activated_equalities function correctly counts + # the number of fixed variables that are in the activated equality constraints in the model + assert number_fixed_variables_in_activated_equalities(model) == 0 + + # Now fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + assert number_fixed_variables_in_activated_equalities(model) == 2 + + # Deactivate the GreyBox and linking constraint + model.b1.egb.deactivate() + model.b1.link_Pin.deactivate() + assert number_fixed_variables_in_activated_equalities(model) == 0 + + @pytest.mark.unit + def test_unfixed_variables_in_activated_equalities_set_w_grey_box(self, model): + # Test that the unfixed_variables_in_activated_equalities_set function correctly identifies + # the unfixed variables that are in the activated equality constraints in the model + uvaes = unfixed_variables_in_activated_equalities_set( + model, include_greybox=True + ) + assert len(uvaes) == 15 + for v in uvaes: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Now fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + uvaes = unfixed_variables_in_activated_equalities_set( + model, include_greybox=True + ) + assert len(uvaes) == 13 + for v in uvaes: + assert v.name in [ + "b1.c", + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b2.v1", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + @pytest.mark.unit + def test_unfixed_greybox_variables_w_grey_box(self, model): + # Test that the unfixed_greybox_variables function correctly identifies + # the unfixed grey box variables in the model + ugv = unfixed_greybox_variables(model) + ugv_list = list(ugv) + assert len(ugv_list) == 7 + for v in ugv_list: + assert v.name in [ + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + @pytest.mark.unit + def test_greybox_variables_w_grey_box(self, model): + # Test that the greybox_variables function correctly identifies + # the grey box variables in the model + gbv = greybox_variables(model) + gbv_list = list(gbv) + assert len(gbv_list) == 7 + for v in gbv_list: + assert v.name in [ + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + @pytest.mark.unit + def test_number_of_unfixed_greybox_variables_w_grey_box(self, model): + # Test that the number_of_unfixed_greybox_variables function correctly + # counts the number of unfixed grey box variables in the model + assert number_of_unfixed_greybox_variables(model) == 7 + # Now fix some of the grey box variables and test again + model.b1.egb.inputs["Pin"].fix(1) + assert number_of_unfixed_greybox_variables(model) == 6 + + @pytest.mark.unit + def test_number_of_greybox_variables_w_grey_box(self, model): + # Test that the number_of_greybox_variables function correctly counts + # the number of grey box variables in the model + assert number_of_greybox_variables(model) == 7 + + @pytest.mark.unit + def test_number_unfixed_variables_in_activated_equalities_w_grey_box(self, model): + # Test that the number_unfixed_variables_in_activated_equalities function correctly counts + # the number of unfixed variables that are in the activated equality constraints in the model + assert ( + number_unfixed_variables_in_activated_equalities( + model, include_greybox=True + ) + == 15 + ) + + # Now fix some of the variables and test again + model.b1.Pin.fix(1) + model.b1.egb.inputs["Pin"].fix(1) + assert ( + number_unfixed_variables_in_activated_equalities( + model, include_greybox=True + ) + == 13 + ) + + @pytest.mark.unit + def test_fixed_variables_only_in_inequalities_w_grey_box(self, model): + # Test that the fixed_variables_only_in_inequalities function correctly identifies + # the fixed variables that are only in the inequality constraints in the model + fvois = fixed_variables_only_in_inequalities(model) + assert len(fvois) == 0 + + # Fix Pin and deactivate the linking constraint so that Pin only appears + # in an inequality constraint + model.b1.Pin.fix(1) + model.b1.link_Pin.deactivate() + fvois = fixed_variables_only_in_inequalities(model) + assert len(fvois) == 1 + for v in fvois: + assert v.name in ["b1.Pin"] + + @pytest.mark.unit + def test_number_fixed_variables_only_in_inequalities_w_grey_box(self, model): + # Test that the number_fixed_variables_only_in_inequalities function correctly counts + # the number of fixed variables that are only in the inequality constraints in the model + assert number_fixed_variables_only_in_inequalities(model) == 0 + + # Fix Pin and deactivate the linking constraint so that Pin only appears + # in an inequality constraint + model.b1.Pin.fix(1) + model.b1.link_Pin.deactivate() + assert number_fixed_variables_only_in_inequalities(model) == 1 + + @pytest.mark.unit + def test_unused_variables_set(self, model): + # Test that the unused_variables_set function correctly identifies the + # unused variables in the model + uv = unused_variables_set(model, include_greybox=True) + assert len(uv) == 0 + + # Deactivate grey box and linking constraints + model.b1.egb.deactivate() + model.b1.link_Pin.deactivate() + uv = unused_variables_set(model, include_greybox=True) + assert len(uv) == 1 + for v in uv: + assert v.name in ["b1.egb.inputs[Pin]"] + + @pytest.mark.unit + def test_number_unused_variables(self, model): + # Test that the number_unused_variables function correctly counts the + # number of unused variables in the model + assert number_unused_variables(model, include_greybox=True) == 0 + + # Deactivate grey box and linking constraints + model.b1.egb.deactivate() + model.b1.link_Pin.deactivate() + assert number_unused_variables(model, include_greybox=True) == 1 + + @pytest.mark.unit + def test_fixed_unused_variables_set(self, model): + # Test that the fixed_unused_variables_set function correctly identifies the + # fixed unused variables in the model + fuv = fixed_unused_variables_set(model, include_greybox=True) + assert len(fuv) == 0 + + # Fix Pin and deactivate the linking constraint so that Pin is unused and fixed + model.b1.Pin.fix(1) + model.b1.link_Pin.deactivate() + model.b1.ineq.deactivate() # Deactivate the inequality constraint so that Pin + # is not in any active constraints + fuv = fixed_unused_variables_set(model, include_greybox=True) + assert len(fuv) == 1 + for v in fuv: + assert v.name in ["b1.Pin"] + + @pytest.mark.unit + def test_number_fixed_unused_variables(self, model): + # Test that the number_fixed_unused_variables function correctly counts + # the number of fixed unused variables in the model + assert number_fixed_unused_variables(model, include_greybox=True) == 0 + + # Fix Pin and deactivate the linking constraint so that Pin is unused and fixed + model.b1.Pin.fix(1) + model.b1.link_Pin.deactivate() + model.b1.ineq.deactivate() # Deactivate the inequality constraint so that + # Pin is not in any active constraints + assert number_fixed_unused_variables(model, include_greybox=True) == 1 + + @pytest.mark.unit + def test_derivative_variables_set(self, model): + # Test that the derivative_variables_set function correctly identifies + # the derivative variables in the model + dvs = derivative_variables_set(model) + assert len(dvs) == 0 + + @pytest.mark.unit + def test_number_derivative_variables(self, model): + # Test that the number_derivative_variables function correctly counts + # the number of derivative variables in the model + assert number_derivative_variables(model) == 0 + + +class TestObjectiveStatisticsGreyBox: + @pytest.mark.unit + def test_total_objectives_generator(self, model): + # Test that the total_objectives_generator function correctly identifies the + # objectives in the model + to = total_objectives_generator(model) + to_list = list(to) + assert len(to_list) == 1 + for o in to_list: + assert o.name in ["obj"] + + @pytest.mark.unit + def test_total_objectives_set(self, model): + # Test that the total_objectives_set function correctly identifies the + # objectives in the model + to = total_objectives_set(model) + assert len(to) == 1 + for o in to: + assert o.name in ["obj"] + + @pytest.mark.unit + def test_number_total_objectives(self, model): + # Test that the number_total_objectives function correctly counts the + # number of objectives in the model + assert number_total_objectives(model) == 1 + + @pytest.mark.unit + def test_activated_objectives_generator(self, model): + # Test that the activated_objectives_generator function correctly identifies the + # activated objectives in the model + ao = activated_objectives_generator(model) + ao_list = list(ao) + assert len(ao_list) == 1 + for o in ao_list: + assert o.name in ["obj"] + + @pytest.mark.unit + def test_activated_objectives_set(self, model): + # Test that the activated_objectives_set function correctly identifies the + # activated objectives in the model + ao = activated_objectives_set(model) + assert len(ao) == 1 + for o in ao: + assert o.name in ["obj"] + + @pytest.mark.unit + def test_number_activated_objectives(self, model): + # Test that the number_activated_objectives function correctly counts the + # number of activated objectives in the model + assert number_activated_objectives(model) == 1 + + @pytest.mark.unit + def test_deactivated_objectives_generator(self, model): + # Test that the deactivated_objectives_generator function correctly identifies the + # deactivated objectives in the model + model.obj.deactivate() + do = deactivated_objectives_generator(model) + do_list = list(do) + assert len(do_list) == 1 + for o in do_list: + assert o.name in ["obj"] + + @pytest.mark.unit + def test_deactivated_objectives_set(self, model): + # Test that the deactivated_objectives_set function correctly identifies the + # deactivated objectives in the model + model.obj.deactivate() + do = deactivated_objectives_set(model) + assert len(do) == 1 + for o in do: + assert o.name in ["obj"] + + @pytest.mark.unit + def test_number_deactivated_objectives(self, model): + # Test that the number_deactivated_objectives function correctly counts the + # number of deactivated objectives in the model + model.obj.deactivate() + assert number_deactivated_objectives(model) == 1 + + +class TestExpressionStatisticsGreyBox: + @pytest.mark.unit + def test_expressions_set(self, model): + # Test that the expressions_set function correctly identifies the + # expressions in the model + ex = expressions_set(model) + assert len(ex) == 1 + for e in ex: + assert e.name in ["expr"] + + @pytest.mark.unit + def test_number_expressions(self, model): + # Test that the number_expressions function correctly counts the + # number of expressions in the model + assert number_expressions(model) == 1 + + +# ------------------------------------------------------------------------- +# Tests for new model_statistics functions (grey box) +class TestNewStatisticsGreyBox: + @pytest.mark.unit + def test_external_variables_set(self, model): + # Add a cross-block constraint + model.b2.cross_cons = pyo.Constraint( + expr=model.b1.egb.inputs["c"] == model.b2.v1 + ) + + ext_vars = external_variables_set(model, include_greybox=True) + assert len(ext_vars) == 0 + + ext_vars = external_variables_set(model.b2, include_greybox=True) + assert len(ext_vars) == 1 + assert model.b1.egb.inputs["c"] in ext_vars + + ext_vars = external_variables_set(model.b2, include_greybox=False) + assert len(ext_vars) == 1 + assert model.b1.egb.inputs["c"] in ext_vars + + model.b2.cross_cons.deactivate() + ext_vars = external_variables_set(model.b2, include_greybox=True) + assert len(ext_vars) == 0 + + @pytest.mark.unit + def test_number_external_variables(self, model): + # Add a cross-block constraint + model.b2.cross_cons = pyo.Constraint( + expr=model.b1.egb.inputs["c"] == model.b2.v1 + ) + + assert number_external_variables(model, include_greybox=True) == 0 + assert number_external_variables(model.b2, include_greybox=True) == 1 + assert number_external_variables(model.b2, include_greybox=False) == 1 + model.b2.cross_cons.deactivate() + assert number_external_variables(model.b2, include_greybox=True) == 0 + + @pytest.mark.unit + def test_variables_fixed_to_zero_set(self, model): + model.b1.zero_var = pyo.Var(initialize=0) + model.b1.zero_var.fix(0) + # Fix a grey box var to zero (should never do this in practice) + model.b1.egb.inputs["F"].fix(0) + + var_set = variables_fixed_to_zero_set(model, include_greybox=True) + assert len(var_set) == 2 + for v in var_set: + assert v.name in ["b1.zero_var", "b1.egb.inputs[F]"] + + var_set = variables_fixed_to_zero_set(model, include_greybox=False) + assert len(var_set) == 1 + for v in var_set: + assert v.name in ["b1.zero_var"] + + @pytest.mark.unit + def test_number_variables_fixed_to_zero(self, model): + model.b1.zero_var = pyo.Var(initialize=0) + model.b1.zero_var.fix(0) + # Fix a grey box var to zero (should never do this in practice) + model.b1.egb.inputs["F"].fix(0) + + assert number_variables_fixed_to_zero(model, include_greybox=True) == 2 + assert number_variables_fixed_to_zero(model, include_greybox=False) == 1 + + @pytest.mark.unit + def test_variables_near_zero_set(self, model): + model.b1.near_zero = pyo.Var(initialize=1e-5) + model.b1.near_zero.value = 1e-5 + + var_set = variables_near_zero_set(model, tol=1e-4, include_greybox=True) + for v in var_set: + assert v.name in ["b1.near_zero", "b1.c", "b1.egb.inputs[c]"] + assert len(var_set) == 3 + var_set = variables_near_zero_set(model, tol=1e-4, include_greybox=False) + for v in var_set: + assert v.name in ["b1.near_zero", "b1.c"] + assert len(var_set) == 2 + + @pytest.mark.unit + def test_number_variables_near_zero(self, model): + model.b1.near_zero = pyo.Var(initialize=1e-5) + model.b1.near_zero.value = 1e-5 + assert number_variables_near_zero(model, tol=1e-4, include_greybox=True) == 3 + assert number_variables_near_zero(model, tol=1e-4, include_greybox=False) == 2 + + @pytest.mark.unit + def test_variables_with_none_value_set(self, model): + model.b1.none_var = pyo.Var() + var_set = variables_with_none_value_set(model, include_greybox=True) + for v in var_set: + assert v.name in [ + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b1.none_var", + "b2.v1", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + assert len(var_set) == 11 + var_set = variables_with_none_value_set(model, include_greybox=False) + for v in var_set: + assert v.name in [ + "b1.F", + "b1.P1", + "b1.P3", + "b1.P2", + "b1.Pout", + "b1.none_var", + "b2.v1", + ] + assert len(var_set) == 7 + + @pytest.mark.unit + def test_number_variables_with_none_value(self, model): + model.b1.none_var = pyo.Var() + assert number_variables_with_none_value(model, include_greybox=True) == 11 + assert number_variables_with_none_value(model, include_greybox=False) == 7 + + @pytest.mark.unit + def test_variables_with_extreme_values_set(self, model): + model.b1.large = pyo.Var(initialize=1e6) + model.b1.small = pyo.Var(initialize=1e-8) + model.b1.zero = pyo.Var(initialize=0) + var_set = variables_with_extreme_values_set( + model, large=1e5, small=1e-7, zero=1e-10, include_greybox=True + ) + for v in var_set: + assert v.name in [ + "b1.Pin", + "b1.c", + "b1.large", + "b1.small", + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + ] + assert len(var_set) == 6 + var_set = variables_with_extreme_values_set( + model, large=1e5, small=1e-7, zero=1e-10, include_greybox=False + ) + for v in var_set: + assert v.name in ["b1.Pin", "b1.c", "b1.large", "b1.small"] + assert len(var_set) == 4 + + @pytest.mark.unit + def test_number_variables_with_extreme_values(self, model): + model.b1.large = pyo.Var(initialize=1e6) + model.b1.small = pyo.Var(initialize=1e-8) + model.b1.zero = pyo.Var(initialize=0) + assert ( + number_variables_with_extreme_values( + model, large=1e5, small=1e-7, zero=1e-10, include_greybox=True + ) + == 6 + ) + assert ( + number_variables_with_extreme_values( + model, large=1e5, small=1e-7, zero=1e-10, include_greybox=False + ) + == 4 + ) + + +@pytest.mark.unit +def test_degrees_of_freedom(model): + # Test that the degree_of_freedom function correctly calculates the + # degree of freedom for the model + # Base degrees of freedom are 3: + # 15 vars - 7 linking constraints - 4 greybox constraints - 1 other equality = 3 + assert degrees_of_freedom(model, include_greybox=True) == 3 + + # Exclude grey box from check + assert degrees_of_freedom(model, include_greybox=False) == 3 + 4 + + # Deactivate grey box and check again - should be the same as excluding grey box + model.b1.egb.deactivate() + assert degrees_of_freedom(model, include_greybox=True) == 3 + 4 + + +@pytest.mark.unit +def test_large_residuals_set(model): + # Test that the large_residuals_set function correctly identifies the + # constraints with large residuals in the model + + # Set up the model so that we have some large residuals + model.b1.Pin.set_value(1e5) + model.b1.egb.inputs["Pin"].set_value(1.1e5) + + model.b1.c.set_value(0.9) + model.b1.egb.inputs["c"].set_value(0.9) # No residual + + model.b1.F.set_value(10) + model.b1.egb.inputs["F"].set_value(10.000009) + + # No value model.b1.P1 + model.b1.egb.inputs["P1"].set_value(2e5) + + model.b1.P2.set_value(1e5) + model.b1.egb.outputs["P2"].set_value(1e5) + model.b1.P3.set_value(1e5) + model.b1.egb.inputs["P3"].set_value(1e5) + model.b1.Pout.set_value(1e5) + model.b1.egb.outputs["Pout"].set_value(1e5) + + # Set b2.v1 to a value that should violate both equality and inequality constraints + model.b2.v1.set_value(-1) + + # Set input values inside the grey box + model.b1.egb.get_external_model().set_input_values( + [1.1e5, 0.9, 10.000009, 2e5, 1e5] + ) + + # Test with default tolerance - expect to see 8 constraints + lrs = large_residuals_set(model, include_greybox=True) + assert len(lrs) == 8 + for c in lrs: + assert c.name in [ + "b1.link_Pin", + "b1.link_P1", + "b2.c1", + "b2.ineq", + "b1.egb.pdrop1", + "b1.egb.pdrop3", + "b1.egb.P2_constraint", + "b1.egb.Pout_constraint", + ] + + # Lower tolerance - should see "b1.link_F" as well + lrs = large_residuals_set(model, tol=1e-8, include_greybox=True) + assert len(lrs) == 9 + for c in lrs: + assert c.name in [ + "b1.link_Pin", + "b1.link_P1", + "b1.link_F", + "b2.c1", + "b2.ineq", + "b1.egb.pdrop1", + "b1.egb.pdrop3", + "b1.egb.P2_constraint", + "b1.egb.Pout_constraint", + ] + + # Exclude grey box + lrs = large_residuals_set(model, tol=1e-8, include_greybox=False) + assert len(lrs) == 5 + for c in lrs: + assert c.name in [ + "b1.link_Pin", + "b1.link_P1", + "b1.link_F", + "b2.c1", + "b2.ineq", + ] + + # Deactivate grey box and check again - should be the same as excluding grey box + model.b1.egb.deactivate() + lrs = large_residuals_set(model, tol=1e-8, include_greybox=True) + assert len(lrs) == 5 + for c in lrs: + assert c.name in [ + "b1.link_Pin", + "b1.link_P1", + "b1.link_F", + "b2.c1", + "b2.ineq", + ] + + +@pytest.mark.unit +def test_large_residuals_set_w_value_w_grey_box(model): + # Test that the large_residuals_set function correctly identifies the + # constraints with large residuals in the model when we also check values + # Set up the model so that we have some large residuals + model.b1.Pin.set_value(1e5) + model.b1.egb.inputs["Pin"].set_value(1.1e5) + + model.b1.c.set_value(0.9) + model.b1.egb.inputs["c"].set_value(0.9) # No residual + + model.b1.F.set_value(10) + model.b1.egb.inputs["F"].set_value(10.000009) + + # No value model.b1.P1 + model.b1.egb.inputs["P1"].set_value(2e5) + + model.b1.P2.set_value(1e5) + model.b1.egb.outputs["P2"].set_value(1e5) + model.b1.P3.set_value(1e5) + model.b1.egb.inputs["P3"].set_value(1e5) + model.b1.Pout.set_value(1e5) + model.b1.egb.outputs["Pout"].set_value(1e5) + + # Set b2.v1 to a value that should violate both equality and inequality constraints + model.b2.v1.set_value(-1) + + # Set input values inside the grey box + model.b1.egb.get_external_model().set_input_values( + [1.1e5, 0.9, 10.000009, 2e5, 1e5] + ) + + # Test with default tolerance - expect to see 8 constraints + lrs = large_residuals_set(model, include_greybox=True, return_residual_values=True) + assert len(lrs) == 8 + + expected = { + "b1.link_Pin": 1e4, + "b1.link_P1": None, + "b2.c1": 2, + "b2.ineq": 1, + "b1.egb.pdrop1": 90090.00016200007, + "b1.egb.pdrop3": 99819.99967599986, + "b1.egb.P2_constraint": 99909.99983799993, + "b1.egb.Pout_constraint": 9639.999351999708, + } + + for c, r in lrs.items(): + print(f"Constraint: {c.name}, Residual: {r}") + if r is not None: + assert r == pytest.approx(expected[c.name], rel=1e-10) + else: + assert expected[c.name] is None + + +@pytest.mark.unit +def test_number_large_residuals_w_value_w_grey_box(model): + # Test that the number_large_residuals function correctly counts the + # number of constraints with large residuals in the model when we also check values + # Set up the model so that we have some large residuals + model.b1.Pin.set_value(1e5) + model.b1.egb.inputs["Pin"].set_value(1.1e5) + + model.b1.c.set_value(0.9) + model.b1.egb.inputs["c"].set_value(0.9) # No residual + + model.b1.F.set_value(10) + model.b1.egb.inputs["F"].set_value(10.000009) + + # No value model.b1.P1 + model.b1.egb.inputs["P1"].set_value(2e5) + + model.b1.P2.set_value(1e5) + model.b1.egb.outputs["P2"].set_value(1e5) + model.b1.P3.set_value(1e5) + model.b1.egb.inputs["P3"].set_value(1e5) + model.b1.Pout.set_value(1e5) + model.b1.egb.outputs["Pout"].set_value(1e5) + + # Set b2.v1 to a value that should violate both equality and inequality constraints + model.b2.v1.set_value(-1) + + # Set input values inside the grey box + model.b1.egb.get_external_model().set_input_values( + [1.1e5, 0.9, 10.000009, 2e5, 1e5] + ) + + # Test with default tolerance - expect to see 8 constraints + assert number_large_residuals(model, include_greybox=True) == 8 + + # Lower tolerance - should see "b1.link_F" as well + assert number_large_residuals(model, tol=1e-8, include_greybox=True) == 9 + + # Exclude grey box + assert number_large_residuals(model, tol=1e-8, include_greybox=False) == 5 + + # Deactivate grey box and check again - should be the same as excluding grey box + model.b1.egb.deactivate() + assert number_large_residuals(model, tol=1e-8, include_greybox=True) == 5 + + +@pytest.mark.unit +def test_active_variables_in_deactivated_blocks_set_w_grey_box(model): + # Test that the active_variables_in_deactivated_blocks_set function correctly identifies the + # active variables in the deactivated blocks in the model + avidbs = active_variables_in_deactivated_blocks_set(model) + assert len(avidbs) == 0 + + # Deactivate the grey box and test again - should see all the grey box variables as they are active + # but in a deactivated block + model.b1.egb.deactivate() + avidbs = active_variables_in_deactivated_blocks_set(model) + assert len(avidbs) == 7 + for v in avidbs: + assert v.name in [ + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + # Test with include_greybox = False - should be the same as deactivating the grey box + avidbs = active_variables_in_deactivated_blocks_set(model, include_greybox=False) + assert len(avidbs) == 7 + for v in avidbs: + assert v.name in [ + "b1.egb.inputs[Pin]", + "b1.egb.inputs[c]", + "b1.egb.inputs[F]", + "b1.egb.inputs[P1]", + "b1.egb.outputs[P2]", + "b1.egb.inputs[P3]", + "b1.egb.outputs[Pout]", + ] + + +@pytest.mark.unit +def test_number_active_variables_in_deactivated_blocks_w_grey_box(model): + # Test that the number_active_variables_in_deactivated_blocks function correctly counts the + # number of active variables in the deactivated blocks in the model + assert number_active_variables_in_deactivated_blocks(model) == 0 + + # Deactivate the grey box and test again - should see all the grey box variables as they are active + # but in a deactivated block + model.b1.egb.deactivate() + assert number_active_variables_in_deactivated_blocks(model) == 7 + + # Test with include_greybox = False - should be the same as deactivating the grey box + assert ( + number_active_variables_in_deactivated_blocks(model, include_greybox=False) == 7 + ) + + +@pytest.mark.unit +def test_variables_with_none_value_in_activated_equalities_set_w_grey_box(model): + # Test that the variables_with_none_value_in_activated_equalities_set function correctly identifies + # the variables with None value in the activated equality constraints in the model + vwnvieas = variables_with_none_value_in_activated_equalities_set( + model, include_greybox=True + ) + assert len(vwnvieas) == 10 + for v in vwnvieas: + assert v.name in [ + "b1.F", + "b1.egb.inputs[F]", + "b1.P1", + "b1.egb.inputs[P1]", + "b1.P3", + "b1.egb.inputs[P3]", + "b1.P2", + "b1.Pout", + "b1.egb.outputs[Pout]", + "b2.v1", + ] + + # Test with include_greybox = False + vwnvieas = variables_with_none_value_in_activated_equalities_set( + model, include_greybox=False + ) + assert len(vwnvieas) == 10 + for v in vwnvieas: + assert v.name in [ + "b1.F", + "b1.egb.inputs[F]", + "b1.P1", + "b1.egb.inputs[P1]", + "b1.P3", + "b1.egb.inputs[P3]", + "b1.P2", + "b1.Pout", + "b1.egb.outputs[Pout]", + "b2.v1", + ] + + # Deactivate the grey box and check again + model.b1.egb.deactivate() + vwnvieas = variables_with_none_value_in_activated_equalities_set( + model, include_greybox=False + ) + assert len(vwnvieas) == 10 + for v in vwnvieas: + assert v.name in [ + "b1.F", + "b1.egb.inputs[F]", + "b1.P1", + "b1.egb.inputs[P1]", + "b1.P3", + "b1.egb.inputs[P3]", + "b1.P2", + "b1.Pout", + "b1.egb.outputs[Pout]", + "b2.v1", + ] + + +@pytest.mark.unit +def test_number_variables_with_none_value_in_activated_equalities_w_grey_box(model): + # Test that the number_variables_with_none_value_in_activated_equalities function correctly counts + # the number of variables with None value in the activated equality constraints in the model + assert ( + number_variables_with_none_value_in_activated_equalities( + model, include_greybox=True + ) + == 10 + ) + + # Test with include_greybox = False + assert ( + number_variables_with_none_value_in_activated_equalities( + model, include_greybox=False + ) + == 10 + ) + + # Deactivate the grey box and check again + model.b1.egb.deactivate() + assert ( + number_variables_with_none_value_in_activated_equalities( + model, include_greybox=False + ) + == 10 + ) + + +@pytest.mark.unit +def test_external_variables_set_w_grey_box(model): + # Add a constraint on b2 that points to a Var in EGB + model.b2.cons_w_ext_var = Constraint(expr=model.b1.egb.inputs["Pin"] == model.b2.v1) + + assert len(external_variables_set(model)) == 0 + ext_vars = external_variables_set(model.b2) + assert len(ext_vars) == 1 + assert model.b1.egb.inputs["Pin"] in ext_vars + + +@pytest.mark.unit +def test_number_external_variables_w_grey_box(model): + model.b2.cons_w_ext_var = Constraint(expr=model.b1.egb.inputs["Pin"] == model.b2.v1) + + assert number_external_variables(model) == 0 + assert number_external_variables(model.b2) == 1