Implicit constraint

[andrewlee94]

Fixes None.

Summary/Motivation:

Pyomo currently has no way of explicitly representing the constraints within an ExternalGreyBoxBlock, which means that they cannot be seen by diagnostics tool (such as Incidence Analysis). This PR adds a new optional ExternalGreyBoxConstraint component which duck-types Constraint which can be added to ExternalGreyBoxBlocks to provide a representation of the grey box constraints. The incidence analysis and PyomoNLPWithGreyBoxBlocks code has been updated to recognise and understand these components.

Changes proposed in this PR:

• Adds ExternalGreyBoxConstraint objects
• Updates incidence analysis tools to look for ExternalGreyBoxConstraints and include these in the incidence analysis
• Updates PyomoNLPWithGreyBoxBlocks to include ExternalGreyBoxConstraints in constraint mappings

Legal Acknowledgement

By contributing to this software project, I have read the contribution guide and agree to the following terms and conditions for my contribution:

1. I agree my contributions are submitted under the BSD license.
2. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

[Robbybp]

Thanks for the PR! At first glance, this looks pretty good. It will take a little while to get through this.

1. Can you post a simple example of calling an Incidence Analysis method on a model with a gray box block
2. I probably would have tried to do this only with something like your EGBConstraintBody object, i.e., I would not have created ExternalGreyBoxConstraint. Did you consider this? What was the motivation to add an "official" Component object?

[andrewlee94]

@Robbybp

1. There are some tests which show how the Incidence Analysis tools work with ExternalGreyBoxBlocks, which are probably the best examples I have: https://github.com/Pyomo/pyomo/pull/3870/changes#diff-97563c1c5b8bf2ed2f37b88015b2fdb3004548d7c4b640dc7a68906a4f48b8d7. The big thing to note in the examples is the use of build_implicit_constraint_objects=True so that the ExternalGreyBoxConstraint objects are constructed - from there usage is identical and the results are what I expected them to look like.

2. The reason for a full component was to make it easy to integrate with things like the IDAES DiagnosticsToolbox. Having a component lets us use component_data_objects(ExternalGreyBoxConstraint) so it was easy to update the diagnostics tools to include these. It also lets us duck-type the implicit constraint objects so they look and behave much like constraints to these tools; the main difference being that they do not have an expr which raises an exception.

[Robbybp]

Some initial comments. I haven't had time to think through the implications of ExternalGreyBoxConstraint yet.

This is reminding me that our grey-box infrastructure could use some better documentation.

[Robbybp]

Any reason this is optional?

[andrewlee94]

Only to ensure there is backward compatibility. I decided to make sure there was still a way to guarantee the model would behave the same as before.

[Robbybp]

We can't assume that inputs are indexed by input_names(). If inputs are provided to set_external_model, they will be a reference-to-list of VarData. Here's an example that fails with the current implementation:

inputs-as-reference.py

import numpy as np
import pyomo.environ as pyo
from scipy.sparse import coo_matrix

from pyomo.contrib.pynumero.interfaces.external_grey_box import (
    ExternalGreyBoxBlock,
    ExternalGreyBoxModel,
)
from pyomo.contrib.incidence_analysis import IncidenceGraphInterface

class MyGreyBox(ExternalGreyBoxModel):
    def n_inputs(self):
        return 4

    def input_names(self):
        return [f"x[{i}]" for i in range(1, self.n_inputs() + 1)]

    def n_equality_constraints(self):
        return 3

    def equality_constraint_names(self):
        return [f"eq[{i}]" for i in range(1, self.n_equality_constraints() + 1)]

    def set_input_values(self, input_values):
        if len(input_values) != self.n_inputs():
            raise ValueError(
                f"Expected {self.n_inputs()} inputs, got {len(input_values)}."
            )
        self._inputs = np.asarray(input_values, dtype=float)

    def evaluate_equality_constraints(self):
        x1, x2, x3, x4 = self._inputs
        return np.array(
            [
                x2 * x3**1.5 * x4 - 5.0,
                x1 - x2 + x4 - 4.0,
                x1 * np.exp(-x3) - 1.0,
            ],
            dtype=float,
        )

    def evaluate_jacobian_equality_constraints(self):
        x1, x2, x3, x4 = self._inputs

        # Rows correspond to eq[1], eq[2], eq[3]
        # Cols correspond to x[1], x[2], x[3], x[4]
        rows = np.array([0, 0, 0, 1, 1, 1, 2, 2], dtype=int)
        cols = np.array([1, 2, 3, 0, 1, 3, 0, 2], dtype=int)
        data = np.array(
            [
                x3**1.5 * x4,  # d(eq1)/d(x2)
                1.5 * x2 * x3**0.5 * x4,  # d(eq1)/d(x3)
                x2 * x3**1.5,  # d(eq1)/d(x4)
                1.0,  # d(eq2)/d(x1)
                -1.0,  # d(eq2)/d(x2)
                1.0,  # d(eq2)/d(x4)
                np.exp(-x3),  # d(eq3)/d(x1)
                -x1 * np.exp(-x3),  # d(eq3)/d(x3)
            ],
            dtype=float,
        )
        return coo_matrix(
            (data, (rows, cols)),
            shape=(self.n_equality_constraints(), self.n_inputs()),
        )

m = pyo.ConcreteModel()
m.x = pyo.Var([1, 2, 3, 4], bounds=(0.0, None), initialize=1.0)
m.objective = pyo.Objective(
    expr=m.x[1] ** 2 + 2 * m.x[2] ** 2 + 3 * m.x[3] ** 2 + 4 * m.x[4] ** 2
)
m.grey_box = ExternalGreyBoxBlock()
m.grey_box.set_external_model(
    MyGreyBox(),
    inputs=[m.x[i] for i in range(1, 5)],
    build_implicit_constraint_objects=True,
)

igraph = IncidenceGraphInterface(m)
matching = igraph.maximum_matching()

[andrewlee94]

I was not even aware you could do that (hence your comment about documentation I guess). What is the canonical way to defining inputs and names then?

[Robbybp]

I'm not sure if there is a "canonical" way. Here's what I remember:

• Carl's original implementation created new variables on the ExternalGreyBoxBlock. These variables are indexed by the names defined on the ExternalGreyBoxModel. The names are also used under the hood in PyomoNLPWithGreyBoxBlocks (or somewhere) as keys (or something).
• I did not want to duplicate variable for every input, so I added the option to provide input (and output) variables if they already exist on the model. In this case, we create a Reference to the user-provided variables. Names are still used under the hood somewhere.

You can always assume the names exist, although personally I would prefer if we did not rely on them. inputs is a Var whose values are the VarData to be treated as inputs to the grey-box, in order.

[andrewlee94]

So, inputs.keys() or [v.local_name for v in inputs] then.

[Robbybp]

inputs.keys():

• If inputs were provided, these are ints
• If inputs were not provided, these are names

[v.local_name for v in inputs.values()]

• If inputs were provided these are the original var names
• If inputs were not provided, these are, e.g., inputs[input_name]

[Robbybp]

What is the canonical way to defining inputs and names then?

Names can always be accessed by ExternalGreyBoxModel.input_names. Input variables can always be accessed by ExternalGreyBoxBlockData.inputs.values. These are the canonical ways to access both of these data.

[Robbybp]

As above, we can't assume that inputs are indexed by names.

[Robbybp]

My preference would be to always use the Jacobian and to determine incidence using only the Jacobian's structure. I don't like using numeric values of the Jacobian as that implicitly depends on the input values we're evaluating.

In general, I'm open to other methods of determining variable-constraint incidence. Can you provide an example use-case for considering the variable-constraint incidence to be dense even if the Jacobian is sparse? Does this have to do with decomposability when solving strongly connected components?

[andrewlee94]

I think there are a few things going on here, and I am not sure there is a single correct answer. Right now, I think there are three possible paths (up from the two I documented in the code), and all three likely have cases where they would be "more correct".

Firstly, a concrete case I can think of where we would want to assume all variables are incident on all constraints: the block-triangularization solver. As the ExternalGreyBoxModel can only be called as a monolithic block, from the perspective of the block-triangularization solver it cannot be decomposed which effectively means all variables are incident on all constraints.

Beyond this however, things get a bit murkier and I'll start with an example.

The main use case I have for ExternalGreyBoxModels is to wrap a procedural thermodynamic flash algorithm which handles four phases (aqueous, organic, vapor and solids). Incidentally, the procedural code returns hard-zero values for non-physical phases (which results in numerical issues I need to deal with, but that is a separate issue). For the ExternalGreyBoxModel however, this raises an issue that numerically, the structure of the model is effectively changing as phases appear and disappear. When I initially implemented this I tried to filter out all those zero values in the wrapper for efficiency, but that mean I ran into issues in the ExternalGreyBoxModel because the shape of the Jacobian was changing from iteration to iteration. My solution was to make the Jacobian fully specified, and to keep all the zero values. However, this means that we cannot deduce the numerical structure of the problem from the Jacobian shape, and hence the idea to use the value of the Jacobian to get the local strucutre.

This would also be needed/useful for detecting local singularities, such as cases where a multiplier variable goes to zero and effectively makes a constraint under/overspecified.

TLDR: I see three possible approaches, each with their own value:

1. Assume full incidence: useful for the block-triangularization solver, or users of a grey box who don't want to know about what is inside it.
2. Incidence based on Jacobian structure: cheap and unambiguous, but requires that the Jacobian shape reflect the numerical structure of the problem. This is feasible for many problems but starts to fail in cases like mine, and requires the developer of the ExternalGreyBoxModel to implement the Jacobian correctly.
3. Incidence based on Jacobian value: requires some inference of what Jacobian value to consider a variable incident on a constraint (and thus scaling dependent), but can handle more complex cases where the numerical structure of the problem is state-dependent.

[Robbybp]

In that case, the "gray-box way" to do this would be to have a Jacobian that captures every variable-constraint pair that can possibly be nonzero. But you're saying that that is not so useful for debugging? I'll think about what I'd like to do here.

Incidence analysis was not designed for cases where the incidence is state-dependent, so if we support this here, I'd like to think about supporting this for regular constraints as well. We could generate the entire incidence graph then filter edges (maybe with a new method) if they correspond to a low value?

[andrewlee94]

I think there are two parts to that.

The first part if whether we can (and should) expect the grey box Jacobian to properly capture "every variable-constraint pair that can possibly be nonzero", as this depends on the developer of the grey box and how they implement the Jacobian. The best answer there might just be to make sure things are documented clearly.

The second part is whether incidence analysis should look for state-dependent structure, which as you note it was not designed to do (although I think it will at least exclude any expression branches with a fixed zero). Your suggestion about having a separate method for that (if we want to do it) is probably the cleanest approach for that.

However, I am not sure that addresses the point about the block-triangularization solver and the distinction between numerically decomposable and functionally decomposable.

[Robbybp]

For your use case, do you need the block-triangularization solver to support GreyBox?

[andrewlee94]

I would need to think about it more, but you have a better understanding of how the everything works and what might be possible. I agree that is it definitely possible in some cases, but I can see edge cases that we would need to take care of.

For example, say we had a grey box that was fundamentally two separate, block-triangularizable sub-systems, A and B (obviously poor design by the developer, but possible). If we ran it through the BT tools and solver, it is possible that it would identify that sub-block A could be solved for before the inputs for sub-block B were known (numerically at least), meaning that if we were to call the grey box then we could have undefined behaviour when the grey box solves sub-block B (e.g. what happens if the inputs have value of None? What if the current values are infeasible for sub-block B?).

However, as I am not expecting the use the block-triangularization solver myself, I would be fine with deferring this if we think it needs more thought or effort.

[Robbybp]

In this case, B would need to be initialized with well-defined inputs. (I don't think this is a restriction. I think A's inputs need to be valid values at the start of the solve as well.) It doesn't matter if the initial values are infeasible, we just ignore B's residual while solving A. If A's outputs end up making B infeasible, that's a general problem that can happen in any BT algorithm (and in my opinion is a useful feature of the algorithm).

But this doesn't need to happen in this PR. We can always raise NotImplementedError for now. I'm just looking for ways to break this apart and limit the functionality I promise to support moving forward.

[andrewlee94]

Agreed - should I add a check for ExternalGreyBoxBlocks in the BT solver and raise an exception then?

[Robbybp]

Leave it for now. First I want to (a) see what is required to support user-provided inputs and outputs and (b) decide what to do about state-dependent Jacobian structure in the GB.

[andrewlee94]

I think state-dependent structure is probably something for a separate PR too - as you noted the current tools do not consider state dependence on a pure Pyomo model, so grey box models should behave the same way. If we go with that, then all we need for incidence analysis is to look at the shape of the Jacobian that is returned which should be easy.

User inputs should also just be a matter of finding the right way to get the index-name-object mappings.

[andrewlee94]

@Robbybp I implemented an approach to use the Jacobian shape to determine incidence which turned out to be fairly simple, and also addressed the issue with input names. However, using your test case revealed two issues:

1. First and most importantly, to be able to use the Jacobian shape, it means we need to be able to evaluate the Jacobian. This means the input variables have to be initialized at least (but not converged). However, the incidence analysis tools are structural in nature, and thus are meant to be applied before a model is initialized and/or solved. Thus, there is a chance (depending on how the Grey Box is implemented), that we do not have values for the inputs and we will run into failures when we try to get the Jacobian. I ran into this with your test case, which did not initialize the variable values.

2. Secondly, and simpler, I assume output names can be assigned the same way, meaning we need a similar fix there as well?

[andrewlee94]

Digging a little deeper, it is actually a bit more complicated. The main error is coming from the fact that the method to get the Jacobain requires the ExternalGreyBoxModel to have the _inputs attribute. However, this does not get created until the set_input_values method has been called.

This means the issue really lies in the code for the grey box, however all the examples are written this way.

[andrewlee94]

@Robbybp I just pushed my latest version with the fixes for using Jacobian shape and allowing external input and output variables.

[blnicho]

@andrewlee94 @Robbybp is this ready for final reviews?

[Robbybp]

is this ready for final reviews?

I still want to go through the ImplicitConstraint implementation in more detail, but I don't think that should stop anybody else from reviewing.

[blnicho]

is this ready for final reviews?

I still want to go through the ImplicitConstraint implementation in more detail, but I don't think that should stop anybody else from reviewing.

@Robbybp will you have time to look at this PR this week. We would like to get it into the Pyomo May release (which has already been pushed back by a week). The dev team is starting to review it also.

[jsiirola]

Overall, this looks reasonable. Several questions, though...

[jsiirola]

Why does this need a special override (couldn't this be handled in get_incident_variables)?

[Robbybp]

I thought it was more natural to have a function that would handle "any" constraint rather than adding more logic to get_incident_variables, which is documented to operate on expressions.

[jsiirola]

My thought process was to make the EGBBody look more like an expression so that the special case handling was unnecessary.

For example, if EGBBody had an args property that returned the input / output variables, then the existing code should be able to handle it, right?

[jsiirola]

Should this use either NumericValue (or more likely NumericExpression) as a base class? This is effectively representing an anonymous named expression....

[andrewlee94]

What would this look like in practice? I have not looked into NumericExpression before so I am not familiar with what it looks like and what would need to be implemented.

[jsiirola]

NumericExpression is the base class for all Pyomo expression objects (operators, named expressions, etc). Among other things, this would let the "expression" that defines the row in the EGB to be treated like any other expression.

[Robbybp]

My first thought is that this sounds a bit dangerous. The only benefit I see is access to identify_variables, and there is some potentially confusing incompatibility that would come up. What would happen if you try to use this in another constraint or expression?

[jsiirola]

It depends on the writer. For all current writers, you would end up with an error from an unsupported operator in the constraint expression, which I think is what we would want.

[andrewlee94]

What benefits would we get from doing this then (besides identify_variables), and what would the "expression" actually look like - we can get a value and partial derivatives of the expression, but the form of numerical expression is not defined.

[jsiirola]

Why is this part of the container and not the Data?

[andrewlee94]

This is the equivalent of the rule for EGBConstraints, where we provide the list of IDs we need to create. I changed the name to implicit_constraint_ids to be slightly more specific, and updated the doc string to explain what it is and why it is needed.

[Robbybp]

Can you revert this just to clean up the diff slightly.

[Robbybp]

I assume these are added for compatibility with IncidenceGraphInterface? Or is there another reason beyond this?

[andrewlee94]

It may also have been needed for the IDAES diagnostics tools as they use PyomoNLP interface, and assumed the presence of some of these methods.

[Robbybp]

I'm not sure I understand the purpose of any of the tests in this file. The name suggests that we are testing PyomoNLPWithGreyBoxBlocks in the case where the underlying ExternalGreyBoxBlock has constructed ExternalGreyBoxConstraints.

1. This PR makes only minor changes to PyomoNLPWithGreyBoxBlocks, so I don't see the need for extensive testing of its variables, constraints, Jacobian, etc.
2. Since we now always construct the constraints, the existing tests for PyomoNLPWithGreyBoxBlocks test this already, right?

If anything, I think this file should just contain unit tests for the methods added to ducktype PyomoNLP (or these tests can go in the existing PyomoNLPWithGreyBoxNLP test file).

[andrewlee94]

Did you address this in the PR you opened? I believe this was originally added when EGBConstraint were optional - the old tests covered the previous case where EGBConstraints were absent, and these tests covered when they were present.

[Robbybp]

No. I only touched test_external_grey_box_model_with_constraints.py.

I recommend removing this file.

If they don't already exist, can you add tests for the new PyomoNLPWithGreyBoxBlocks methods to the existing test_pyomo_grey_box_nlp.py file? (If they do exist, can you point them out?)

[Robbybp]

Can some of these properties be removed (for now)? I think most of them are not necessary for the functionality this PR is adding. My understanding is that subclassing Component (instead of ducktyping) has the following advantages:

• Access to .pprint
• .name behaves like a normal Pyomo name
• Compatibility with component_data_objects

Everything else, to me, has marginal value for right now. My suggestion is to remove as many properties and methods as possible from this file. My rationale is that it will be relatively easy to add these later if we have a good reason to, but relatively annoying to change if we add them now.

[andrewlee94]

I am not sure if these are all strictly necessary, but they were added so that EGBConstraints would behave like normal Constraints in the IDAES diagnostics and statistics tools.

[Robbybp]

Can you show me the usage patterns where you would like this to behave like a regular constraint?

In general, I'm okay with branching on con.ctype, e.g.

[con.lower if con.ctype is Constraint else 0 for con in constraints]

instead of making this consistently ducktyped. (This is not a great example, as lower is one of the attributes I might actually keep.)

I actually prefer to shift the burden onto the user like this instead of adding more functionality that we need to support, especially when I'm not 100% sure the functionality is well-motivated.

I think at least the following attributes/methods are not necessary: strict_lower/upper, has_lb/ub, get/set_value, (l/u)slack, active, (de)activate, and one of ub/lb or upper/lower.

[andrewlee94]

I was going in the opposite direction here, especially seeing that as these are equalities these values are well defined. The problem with branching on constraint type is that it does require the user (or middle developer) to be aware of the difference, and I am not that confident in users of IDAES. Specifically, in IDAES there are a decent number of statistics and diagnostics functions that rely on bounds for making determinations (including determining if a constraint is an equality, as they were written before the equality property was added/documented). My goal was to make EGBConstraints fit seamlessly with that so that i) I did not have to rewrite all those functions, and ii) to prevent issues in future from devleopers overlooking EGBConstraints as a possibility.

[Robbybp]

How do these tests differ from those in the Incidence Analysis tests? I would lean towards putting all of the Incidence Analysis-Grey box integration tests in Incidence Analysis.

[Robbybp]

I just opened a PR into this branch, andrewlee94#2, where I remove many of the tests in test_external_grey_box_model_with_constraints.py.

[Robbybp]

What is the purpose of implicit_constraint_id?

[andrewlee94]

This is the "rule" for the constraint - it tells us which output or equality constraint key this constraint related to. This gets passed on to the EGBConstraintBody object so that it can look up the correct values from the EGBModel.

[Robbybp]

Okay, I see. What happens if an equality constraint and an output have the same name? Do we need additional data to tells us which "type" of implicit constraint we have?

[andrewlee94]

I had not considered that - my first question there is whether we should support that in general as that ambiguity might extend further than just this (i.e., you can't have two Python objects with the same name, so should we allow implicit constraints to have the same name)?

It would be possible to fully specify the constraint id by including both type and name/index, possibly just breaking the current argument into two (constraint_id and output_id) and allowing only one to be set.

[jsiirola]

Why not just use an integer to indicate the position in the list of constraints / output variables (where output variables are offset by the number of constraints)?

[Robbybp]

Why not just use an integer to indicate the position in the list of constraints / output variables

I would either go with this approach or use a tuple of (enum, int). I don't love that names are a mandatory part of ExternalGreyBoxModel and would like to rely on them as little as possible.