Debugging a structural singularity with the Dulmage-Mendelsohn partition

We start with some imports and by creating a Pyomo model we would like to debug. Usually the model is much larger and more complicated than this. This particular system appeared when debugging a dynamic 1-D partial differential-algebraic equation (PDAE) model representing a chemical looping combustion reactor.

>>> import pyomo.environ as pyo
>>> from pyomo.contrib.incidence_analysis import IncidenceGraphInterface

>>> m = pyo.ConcreteModel()
>>> m.components = pyo.Set(initialize=[1, 2, 3])
>>> m.x = pyo.Var(m.components, initialize=1.0/3.0)
>>> m.flow_comp = pyo.Var(m.components, initialize=10.0)
>>> m.flow = pyo.Var(initialize=30.0)
>>> m.density = pyo.Var(initialize=1.0)
>>> m.sum_eqn = pyo.Constraint(
...     expr=sum(m.x[j] for j in m.components) - 1 == 0
... )
>>> m.holdup_eqn = pyo.Constraint(m.components, expr={
...     j: m.x[j]*m.density - 1 == 0 for j in m.components
... })
>>> m.density_eqn = pyo.Constraint(
...     expr=1/m.density - sum(1/m.x[j] for j in m.components) == 0
... )
>>> m.flow_eqn = pyo.Constraint(m.components, expr={
...     j: m.x[j]*m.flow - m.flow_comp[j] == 0 for j in m.components
... })

To check this model for structural singularity, we apply the Dulmage-Mendelsohn partition. var_dm_partition and con_dm_partition are named tuples with fields for each of the four subsets defined by the partition: unmatched, overconstrained, square, and underconstrained.

>>> igraph = IncidenceGraphInterface(m)
>>> # Make sure we have a square system
>>> print(len(igraph.variables))
8
>>> print(len(igraph.constraints))
8
>>> var_dm_partition, con_dm_partition = igraph.dulmage_mendelsohn()

If any variables or constraints are unmatched, the (Jacobian of the) model is structurally singular.

>>> # Note that the unmatched variables/constraints are not mathematically
>>> # unique and could change with implementation!
>>> for var in var_dm_partition.unmatched:
...     print(var.name)
flow_comp[1]
>>> for con in con_dm_partition.unmatched:
...     print(con.name)
density_eqn

This model has one unmatched constraint and one unmatched variable, so it is structurally singular. However, the unmatched variable and constraint are not unique. For example, flow_comp[2] could have been unmatched instead of flow_comp[1]. The exact variables and constraints that are unmatched depends on both the order in which variables are identified in Pyomo expressions and the implementation of the matching algorithm. For a given implementation, however, these variables and constraints should be deterministic.

Unique subsets of variables and constraints that are useful when debugging a structural singularity are the underconstrained and overconstrained subsystems. The variables in the underconstrained subsystem are contained in the unmatched and underconstrained fields of the var_dm_partition named tuple, while the constraints are contained in the underconstrained field of the con_dm_partition named tuple. The variables in the overconstrained subsystem are contained in the overconstrained field of the var_dm_partition named tuple, while the constraints are contained in the overconstrained and unmatched fields of the con_dm_partition named tuple.

We now construct the underconstrained and overconstrained subsystems:

>>> uc_var = var_dm_partition.unmatched + var_dm_partition.underconstrained
>>> uc_con = con_dm_partition.underconstrained
>>> oc_var = var_dm_partition.overconstrained
>>> oc_con = con_dm_partition.overconstrained + con_dm_partition.unmatched

And display the variables and constraints contained in each:

>>> # Note that while these variables/constraints are uniquely determined,
>>> # their order is not!

>>> # Overconstrained subsystem
>>> for var in oc_var:
>>>     print(var.name)
x[1]
density
x[2]
x[3]
>>> for con in oc_con:
>>>     print(con.name)
sum_eqn
holdup_eqn[1]
holdup_eqn[2]
holdup_eqn[3]
density_eqn

>>> # Underconstrained subsystem
>>> for var in uc_var:
>>>     print(var.name)
flow_comp[1]
flow
flow_comp[2]
flow_comp[3]
>>> for con in uc_con:
>>>     print(con.name)
flow_eqn[1]
flow_eqn[2]
flow_eqn[3]

At this point we must use our intuition about the system being modeled to identify “what is causing” the singularity. Looking at the under and over- constrained systems, it appears that we are missing an equation to calculate flow, the total flow rate, and that density is over-specified as it is computed by both the bulk density equation and one of the component density equations.

With this knowledge, we can eventually figure out (a) that we need an equation to calculate flow from density and (b) that our “bulk density equation” is actually a skeletal density equation. Admittedly, this is difficult to figure out without the full context behind this particular system.

The following code constructs a new version of the model and verifies that it is structurally nonsingular:

>>> import pyomo.environ as pyo
>>> from pyomo.contrib.incidence_analysis import IncidenceGraphInterface

>>> m = pyo.ConcreteModel()
>>> m.components = pyo.Set(initialize=[1, 2, 3])
>>> m.x = pyo.Var(m.components, initialize=1.0/3.0)
>>> m.flow_comp = pyo.Var(m.components, initialize=10.0)
>>> m.flow = pyo.Var(initialize=30.0)
>>> m.dens_bulk = pyo.Var(initialize=1.0)
>>> m.dens_skel = pyo.Var(initialize=1.0)
>>> m.porosity = pyo.Var(initialize=0.25)
>>> m.velocity = pyo.Param(initialize=1.0)
>>> m.sum_eqn = pyo.Constraint(
...     expr=sum(m.x[j] for j in m.components) - 1 == 0
... )
>>> m.holdup_eqn = pyo.Constraint(m.components, expr={
...     j: m.x[j]*m.dens_bulk - 1 == 0 for j in m.components
... })
>>> m.dens_skel_eqn = pyo.Constraint(
...     expr=1/m.dens_skel - sum(1/m.x[j] for j in m.components) == 0
... )
>>> m.dens_bulk_eqn = pyo.Constraint(
...     expr=m.dens_bulk == (1 - m.porosity)*m.dens_skel
... )
>>> m.flow_eqn = pyo.Constraint(m.components, expr={
...     j: m.x[j]*m.flow - m.flow_comp[j] == 0 for j in m.components
... })
>>> m.flow_dens_eqn = pyo.Constraint(
...     expr=m.flow == m.velocity*m.dens_bulk
... )

>>> igraph = IncidenceGraphInterface(m, include_inequality=False)
>>> print(len(igraph.variables))
10
>>> print(len(igraph.constraints))
10
>>> var_dm_partition, con_dm_partition = igraph.dulmage_mendelsohn()

>>> # There are now no unmatched variables or equations
>>> print(len(var_dm_partition.unmatched))
0
>>> print(len(con_dm_partition.unmatched))
0