Debugging a Structural Singularity

###############################################################################
# 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-2023 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.
###############################################################################

Debugging a Structural Singularity#

Author: Robert Parker
Maintainer: Robert Parker
Updated: 2024-06-10

In this tutorial, we will use the IDAES Diagnostics Toolbox to diagnose and fix a structural singularity that is preventing us from solving an optimization problem.

Constructing the model#

Suppose a collaborator has given us a model to work with. They give us a square model and tell us what the degrees of freedom are. We construct an optimization problem and try to solve it. In this tutorial, we don’t want to worry too much about the details that go into constructing the model. This has been provided in the idaes_examples.mod.diagnostics.gas_solid_contactors.model module.

Model details (OKAY TO SKIP)#

The model we are trying to optimize is a dynamic model of a moving bed chemical looping combustion reactor. The model has been described by Okoli et al. and Parker and Biegler. This is a gas-solid reactor with counter-current flow. The degrees of freedom are gas and solid inlet flow rates, and we are trying to minimize the deviation from a desired operating point via a least-squares objective function.

Again, we don’t want to worry too much about the model. The make_model function will construct the optimization problem that we want to solve, and whenever we do something model-specific, we will explicitly make note of it.

Trying to solve the original model#

With that out of the way, let’s construct the model and try to solve it!

from idaes_examples.mod.diagnostics.gas_solid_contactors.model import make_model
import logging

# We'll turn off IDAES logging. This is not recommended in general, but this is an old model
# (from IDAES 1.7) that has been ported to work with the current version of IDAES. It generates
# a lot of warnings.
logging.getLogger("idaes").setLevel(logging.CRITICAL)
# We'll also turn off Pyomo logging. This will suppress unit inconsistency warnings later,
# which otherwise flood our console and slow down this notebook. We have unit inconsistencies
# as, in IDAES 1.7, we didn't rigorously enforce that models use units.
logging.getLogger("pyomo").setLevel(logging.CRITICAL)

# This constructs a dynamic model with degrees of freedom and an objective function.
model = make_model()

Before trying to solve the model, let’s make sure it conforms to our expectations, i.e. it (a) has degrees of freedom and (b) is well-initialized to a feasible point.

# Import some useful utilities from the model_statistics module.
# Degrees of freedom and constraint residuals are always good things to check before
# trying to solve a simulation or optimization problem.
from idaes.core.util.model_statistics import degrees_of_freedom, large_residuals_set

dof = degrees_of_freedom(model)
print(f"Degrees of freedom: {dof}")
has_large_residuals = bool(large_residuals_set(model, tol=1e-5))
print(f"Has large residuals: {has_large_residuals}")

In the above make_model function, the model has been “solved” to arrive at a feasible point, then degrees of freedom have been unfixed and an objective function has been added to give us an optimization problem. This looks good so far, so let’s try to solve the optimization problem.

# Import pyomo.environ for access to solvers
import pyomo.environ as pyo

solver = pyo.SolverFactory("ipopt")
solver.options["max_iter"] = 20
solver.options["print_user_options"] = "yes"
solver.options["OF_print_info_string"] = "yes"
res = solver.solve(model, tee=True)
print(f"Converged successfully: {pyo.check_optimal_termination(res)}")

IPOPT fails to solve the optimization problem… You can try increasing the iteration limit, but it is very unlikely that this model will ever solve. A telltale sign that something is wrong with our model is the persistence of regularization coefficients, that is, numbers in the lg(rg) column of the IPOPT log. These coefficients can have multiple causes. One is that the constraint Jacobian (partial derivative matrix) is singular, which indicates a problem with our model. We have set the print_info_string option in IPOPT to display “diagnostic tags” to help interpret these regularization coefficients. The “L” and “l” diagnostic tags, which appear repeatedly, indicate that the Jacobian is singular. For more information on IPOPT diagnostic tags, see the IPOPT documentation.

Debugging the original model#

Let’s run the diagnostics toolbox on the model and see what it has to say.

For good practice, we’ll first make sure the model we’re debugging is square. Remember that we’re assuming we already know how to toggle degrees of freedom in our model.

# Fix gas and solid flow rates at their respective inlets
model.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()
model.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()
# Part of our optimization problem was a set of constraints to enforce piecewise
# constant control inputs. We need to deactivate these as well.
model.piecewise_constant_constraints.deactivate()

Now we can run the diagnostics toolbox.

from idaes.core.util.model_diagnostics import DiagnosticsToolbox

dt = DiagnosticsToolbox(model)
dt.report_structural_issues()

Let’s look at the warnings we got:

  • Inconsistent units

  • Structural singularity

  • Potential evaluation errors

We’ll ignore the inconsistent units. The property package and unit model here were extracted from IDAES 1.7, before we rigorously enforced that all models use units. The potential evaluation errors we see here may be worth looking into, but looking at the failing IPOPT log above, we don’t notice any evaluation errors. (If evaluation errors occurred in IPOPT, we would see a message like “Error in AMPL evaluation” in the IPOPT iteration log, which we don’t see here.) The structural singularity looks like the most promising avenue to debug, especially as the IPOPT log displays persistent regularization coefficients that appear to be caused by a singular Jacobian.

Let’s follow the toolbox’s advice and display the under and over-constrained sets.

dt.display_underconstrained_set()
dt.display_overconstrained_set()

Over and under-constrained subsystems#

Structural singularities are characterized by the Dulmage-Mendelson decomposition, which partitions a system into minimal over and under-constrained subsystems. These subsystems contain the potentially unmatched constraints and variables, respectively. Here, “unmatched” effectively means “causing a singularity”. Pothen and Fan give a good overview of the Dulmage-Mendelsohn decomposition and Parker et al. give several examples.

The most straightforward way to fix a structural singularity is to fix variables that are in the under-constrained system and deactivate constraints in the over-constrained subsystem. However, this may not be applicable for every model. For example, we may need to add variables and constraints instead. What over and under-constrained subsystems are telling us is that something is wrong with our modeling assumptions. The particular fix that is appropriate will depend heavily on the model.

If the above output gives us any clues, we can go ahead and start trying to fix things. However, suppose it doesn’t. A good strategy is to try to break down the model into smaller, square subsystems that we think should be nonsingular. For a dynamic model like this one, a good candidate is the subsystem of variables and equations at each point in time.

# We've included a utility function to extract the subsystem of variables and equations
# at a specified point in time. If you are dealing with a process flowsheet, here you
# may want to extract each unit model individually.
from idaes_examples.mod.diagnostics.util import get_subsystem_at_time

# TemporarySubsystemManager is used to temporarily fix some variables to make sure
# we're debugging a square subsystem.
from pyomo.util.subsystems import TemporarySubsystemManager

# Let's start with t=0. Really, we'd probably want to do this in a loop and try all time points.
t0 = model.fs.time.first()
t_block, inputs = get_subsystem_at_time(model, model.fs.time, t0)
# We'll temporarily fix the "inputs" to make sure we have a square system while debugging
with TemporarySubsystemManager(to_fix=inputs):
    dt = DiagnosticsToolbox(t_block)
    dt.report_structural_issues()
    dt.display_underconstrained_set()
    dt.display_overconstrained_set()

These over and under-constrained subsystems aren’t much smaller, but now the over-constrained system decomposes into 10 small, independent blocks. These should be easier to debug.

Debugging the over-constrained subsystem#

To debug the over-constrained subsystem, we look for a constraint that is not calculating any of the variables in the subsystem. The “odd constraint out” here seems to be the mass fraction sum, sum_component_eqn. This must “solve for” one of the mass fractions, which means one of the material_holdup_calculation equations must “solve for” particle density rather than mass fraction. If we want to see what variables are contained in one of these constraints, we can always pprint it:

model.fs.MB.solid_phase.properties[0, 0.9].sum_component_eqn.pprint()
model.fs.MB.solid_phase.material_holdup_calculation[0, 0.9, "Sol", "Fe3O4"].pprint()

If one of these material_holdup_calculation equations is solving for particle density, then that means that density_particle_constraint is not actually solving for density. Maybe density_particle_constraint is over-determining our system?

model.fs.MB.solid_phase.properties[0, 0.9].density_particle_constraint.pprint()

But this looks like a very reasonable constraint. After some thought, which admittedly requires some knowledge of the process we are modeling, we decide that the right approach is to make particle porosity a variable. We have assumed that porosity is constant, but this overconstrained subsystem is telling us that this assumption is not valid.

How did we figure this out? (OKAY TO SKIP)#

Adding a variable (including by unfixing a parameter) to an over-constraining constraint will often remove that constraint from the over-constrained subsystem. But how did we know that this was the right thing to do? If you just care about using the diagnostics toolbox to extract as much information about a singularity as possible, you can skip this section. But if you are curious how we determined that particle porosity should not be constant, read on.

dens_mass_skeletal is determined purely by the composition of solid, which is made up of Fe2O3, Fe3O4, and inert Ti2O3. We can view the density_skeletal_constraint as follows:

model.fs.MB.solid_phase.properties[0, 0.9].density_skeletal_constraint.pprint()

If we assume a constant particle porosity, this gives us a particle porosity that is also uniquely determined by the solid composition by the above density_particle_constraint:

dens_mass_particle = (1 - porosity) * dens_mass_skeletal

But the composition of the solid is determined by the (somewhat misnamed) material_holdup_calculation constraints. While the name of these constraints implies they “calculate holdups,” material holdups at \(t=0\) are fixed as initial conditions (because holdups are the differential variables with respect to time in this model). At other time points, we assume that holdups are specified by differential and discretization equations of the model. This means that the material_holdup_calculation constraints actually calculate the solid phase mass fractions from the holdups. But as we hinted at above, the 4-by-4 system of holdup calculation constraints, sum_component_eqn (which simply constrains the sum of mass fractions to be one), mass fractions, and dens_mass_particle, uniquely solve for dens_mass_particle as well as the mass fractions. But if the holdup variables can be used to solve for the mass fractions, they also solve for dens_mass_skeletal. So both sides of density_particle_constraint are already uniquely determined! This implies that we don’t need this constraint at all, but we also know that this constraint has to hold. Something has to give. With this in mind, we actually have several options for how to resolve this overspecification:

  1. Remove density_particle_constraint. Then we would have dens_mass_particle and dens_mass_skeletal, with no relationship between them. This would leave us with a mathematically sound model, but with densities that contradict constant particle porosity that we have assumed (which is used elsewhere in the reaction rate calculation equations).

  2. Remove the constraints that calculate skeletal density from composition.

  3. Relax particle porosity from a parameter to a variable.

Options 2 and 3 are equally valid. We’ve chosen option 3, meaning we assume that the particle “evolves” with a density that is well determined from its constituent species, rather than changing density to accommodate whatever mass it accumulates via reaction without altering its volume. This exercise should remind us that all mathematical modeling is somewhat of an art. In the process of choosing the “least bad” model, it is fairly easy to over or under-specify something by making the wrong combination of assumptions, and the Dulmage-Mendelsohn decomposition is a great tool for detecting when this has happened.

Debugging the under-constrained subsystem#

The under-constrained system does not decompose into independent subsystems, making it more difficult to debug. However, by inspection, we notice that the same constraints and variables seem to be repeated at each point in the length domain. For each point in space, the “odd variable out” seems to be the total flow rate flow_mass. Using some intuition about this particular process model, we may conclude that this variable should be calculated from the solid phase velocity, which is constant. We expect an equation that looks like

flow_mass == velocity * area * density

But this equation isn’t here… so we need to add it.

Fixing the model#

We’ll start by creating a fresh copy of the model, so we don’t accidentally rely on IPOPT’s point of termination.

model2 = make_model()
# Make the model square while we try to fix the structural singularity
model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()
model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()
model2.piecewise_constant_constraints.deactivate()

Adding a new particle porosity variable#

model2.fs.MB.particle_porosity = pyo.Var(
    model2.fs.time,
    model2.fs.MB.length_domain,
    initialize=model2.fs.solid_properties.particle_porosity.value,
)

Now we need to replace the old particle porosity parameter with this new variable. Luckily, the old parameter is actually implemented as a fixed variable, so we can easily identify all the constraints it participates in with IncidenceGraphInterface:

from pyomo.contrib.incidence_analysis import IncidenceGraphInterface

igraph = IncidenceGraphInterface(model2, include_fixed=True)
porosity_param = model2.fs.solid_properties.particle_porosity
print(f"Constraints containing {porosity_param.name}:")
for con in igraph.get_adjacent_to(porosity_param):
    print(f"  {con.name}")

Particle porosity only appears in two constraints: the density constraint we saw above, and the reaction rate equation. We can replace particle porosity in these constraints using Pyomo’s replace_expressions function:

from pyomo.core.expr import replace_expressions

for t, x in model2.fs.time * model2.fs.MB.length_domain:
    substitution_map = {id(porosity_param): model2.fs.MB.particle_porosity[t, x]}
    sp = model2.fs.MB.solid_phase
    cons = [
        sp.properties[t, x].density_particle_constraint,
        sp.reactions[t, x].gen_rate_expression["R1"],
    ]
    for con in cons:
        con.set_value(
            replace_expressions(
                con.expr, substitution_map, descend_into_named_expressions=True
            )
        )

We have added a new particle_porosity variable, and are using it in the relevant locations. Now we can move on to adding the missing constraint.

Adding a new density-flow rate constraint#

@model2.fs.MB.Constraint(model2.fs.time, model2.fs.MB.length_domain)
def density_flowrate_constraint(mb, t, x):
    return (
        mb.velocity_superficial_solid[t]
        * mb.bed_area
        * mb.solid_phase.properties[t, x].dens_mass_particle
        == mb.solid_phase.properties[t, x].flow_mass
    )

Testing the new model#

Let’s see if these changes have fixed our model.

# Construct a new diagnostics toolbox
dt = DiagnosticsToolbox(model2)
dt.report_structural_issues()

The structural singularity seems to be gone! Let’s unfix our degrees of freedom and see if we can solve.

model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()
model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()
model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()
model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()
model2.piecewise_constant_constraints.activate()

res = solver.solve(model2, tee=True)
print(f"Converged successfully: {pyo.check_optimal_termination(res)}")

This doesn’t look much better. What’s going on? I thought we just fixed the issue?

Debugging the model, take two#

Let’s check the diagnostics toolbox for numerical issues.

model2.fs.MB.gas_phase.properties[:, 0].flow_mol.fix()
model2.fs.MB.solid_phase.properties[:, 1].flow_mass.fix()
model2.piecewise_constant_constraints.deactivate()
dt.report_numerical_issues()

Looks like we have “parallel constraints”, which are another form of singularity. Let’s follow the toolbox’s advice to see what they are.

dt.display_near_parallel_constraints()

density_flowrate_constraint is the constraint that we added. What is solid_super_vel?

model2.fs.MB.solid_super_vel[0].pprint()

This is the same as the constraint we just added! Looks like that constraint already existed at the solid inlet. We can easily deactivate the new constraints at this point in the length domain:

model2.fs.MB.density_flowrate_constraint[:, 1.0].deactivate();

But now we have removed constraints from a square model, and expect to have degrees of freedom. Let’s see what the diagnostics toolbox has to say.

dt = DiagnosticsToolbox(model2)
dt.report_structural_issues()

But this doesn’t help us very much. We have some extraneous degrees of freedom, but with 8881 variables in the under-constrained subsystem, it will be difficult to tell what they are. After some thought (and model-specific intuition), we land on the conclusion that maybe we need to fix particle porosity at the solid inlet. Here, total flow rate is specified, and the solid_super_vel equation is using it to compute velocity. So we need dens_mass_particle to be known, which means we need particle_porosity to be fixed.

model2.fs.MB.particle_porosity[:, 1.0].fix();

Let’s run the diagnostics toolbox as a sanity check.

dt = DiagnosticsToolbox(model2)
dt.report_structural_issues()
dt.report_numerical_issues()

Looks good! Now we can release our degrees of freedom and try to solve again.

model2.fs.MB.gas_phase.properties[:, 0].flow_mol.unfix()
model2.fs.MB.gas_phase.properties[0, 0].flow_mol.fix()
model2.fs.MB.solid_phase.properties[:, 1].flow_mass.unfix()
model2.fs.MB.solid_phase.properties[0, 1].flow_mass.fix()
model2.piecewise_constant_constraints.activate()

res = solver.solve(model2, tee=True)
print(f"Converged successfully: {pyo.check_optimal_termination(res)}")

It worked! For the simple optimization problem we have set up, this solve looks a lot more like what we expect.

Takeaways from this tutorial#

What have we learned?

  1. IPOPT using non-zero regularization coefficients hints at a singular Jacobian (especially when “L”/”l” diagnostic tags are present).

  2. When this happens, start by calling report_structural_issues to check for a structural singularity. If this looks good, call report_numerical_issues to check for a numerical singularity.

  3. When debugging a structural singularity, decomposing a problem into subsystems that each should be nonsingular (e.g. unit models or points in time) is very useful.

  4. The solution to a structural singularity is often to relax a fixed parameter, add a constraint that was forgotten, remove a constraint that was redundant, or fix an extraneous degree of freedom.

  5. Model-specific intuition is usually necessary to diagnose and fix modeling issues. (If you’re an algorithm developer, learn about the models you’re using! If you don’t understand your models, you don’t understand your algorithms!)

  6. A modeling issue doesn’t necessarily have a unique solution. This is especially true when the issue involves invalid assumptions.

  7. Debugging is an iterative process — fixing one issue can introduce another.

References#

[1] Okoli et al., “A framework for the optimization of chemical looping combustion processes”. Powder Tech, 2020.

[2] Parker and Biegler, “Dynamic modeling and nonlinear model predictive control of a moving bed chemical looping combustion reactor”. IFAC PapersOnline, 2022.

[3] Dulmage and Mendelsohn, “Coverings of bipartite graphs”. Can J. Math., 1958.

[4] Pothen and Fan, “Computing the block triangular form of a sparse matrix”. ACM Trans. Math. Softw., 1990.

[5] Parker et al., “Applications of the Dulmage-Mendelsohn decomposition for debugging nonlinear optimization problems”. Comp. Chem. Eng., 2023.

Contents