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.