Debugging a numeric singularity using block triangularization

We start with some imports. To debug a numeric singularity, we will need PyomoNLP from PyNumero to get the constraint Jacobian, and will need NumPy to compute condition numbers.

>>> import pyomo.environ as pyo
>>> from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP
>>> from pyomo.contrib.incidence_analysis import IncidenceGraphInterface
>>> import numpy as np

We now build the model we would like to debug. Compared to the model in Debugging a structural singularity with the Dulmage-Mendelsohn partition, we have converted the sum equation to use a sum over component flow rates rather than a sum over mass fractions.

>>> 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)
>>> # This equation is new!
>>> m.sum_flow_eqn = pyo.Constraint(
...     expr=sum(m.flow_comp[j] for j in m.components) == m.flow
... )
>>> 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
... })

We now construct the incidence graph and check unmatched variables and constraints to validate structural nonsingularity.

>>> igraph = IncidenceGraphInterface(m, include_inequality=False)
>>> var_dmp, con_dmp = igraph.dulmage_mendelsohn()
>>> print(len(var_dmp.unmatched))
0
>>> print(len(con_dmp.unmatched))
0

Our system is structurally nonsingular. Now we check whether we are numerically nonsingular (well-conditioned) by checking the condition number. Admittedly, deciding if a matrix is “singular” by looking at its condition number is somewhat of an art. We might define “numerically singular” as having a condition number greater than the inverse of machine precision (approximately 1e16), but poorly conditioned matrices can cause problems even if they don’t meet this definition. Here we use 1e10 as a somewhat arbitrary condition number threshold to indicate a problem in our system.

>>> # PyomoNLP requires exactly one objective function
>>> m._obj = pyo.Objective(expr=0.0)
>>> nlp = PyomoNLP(m)
>>> cond_threshold = 1e10
>>> cond = np.linalg.cond(nlp.evaluate_jacobian_eq().toarray())
>>> print(cond > cond_threshold)
True

The system is poorly conditioned. Now we can check diagonal blocks of a block triangularization to determine which blocks are causing the poor conditioning.

>>> var_blocks, con_blocks = igraph.block_triangularize()
>>> for i, (vblock, cblock) in enumerate(zip(var_blocks, con_blocks)):
...     submatrix = nlp.extract_submatrix_jacobian(vblock, cblock)
...     cond = np.linalg.cond(submatrix.toarray())
...     print(f"block {i}: {cond}")
...     if cond > cond_threshold:
...         for var in vblock:
...             print(f"  {var.name}")
...         for con in cblock:
...             print(f"  {con.name}")
block 0: 24.492504515710433
block 1: 1.2480741394486336e+17
  flow
  flow_comp[1]
  flow_comp[2]
  flow_comp[3]
  sum_flow_eqn
  flow_eqn[1]
  flow_eqn[2]
  flow_eqn[3]

We see that the second block is causing the singularity, and that this block contains the sum equation that we modified for this example. This suggests that converting this equation to sum over flow rates rather than mass fractions just converted a structural singularity to a numeric singularity, and didn’t really solve our problem. To see a fix that does resolve the singularity, see Debugging a structural singularity with the Dulmage-Mendelsohn partition.