Add example of debugging a structural singularity

[Robbybp]

Adds an end-to-end diagnostics example, fixing a real problem that we encountered a few years ago.

Some to-dos on this example:

• Don't rely on writing an ipopt.opt file in the user's working directory
• Update old models (ported from IDAES 1.7) to use units so we don't get so many warnings about them (or just turn off logging, but this seems like poor form for an example)
• Add tests

---

Legal Acknowledgement

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

I agree my contributions are submitted under the license terms described in the LICENSE.txt file at the top level of this directory.
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.

---

📚 Documentation preview 📚: https://idaes-examples--114.org.readthedocs.build/en/114/

[andrewlee94]

Thank you very much @Robbybp , this look great! Edit: Updated the discussion to better reflect both the BFB and MB models

A couple of suggestions to help users understand what is going on:

1. A description of what "positive definite" and other more advanced terms mean might be helpful (or pointers to references). Many chemical engineers don;t get advanced mathematical modeling courses but end up doing modeling anyway (like me).
2. A physical description of the singularity might help readers link the maths to the physics. Here is a first pass from me:

Another way to think of this is in terms of the physical behaviour of the system. Consider a finite element in the system with solids flowing in and out of the element. The total flowrate of solids leaving nay element is defined by F = v*(1-epsilon)*rho where v is the velocity of the solids, epsilon is the voidage of the particle bed, and rho is the density of the solid material.

v and epsilon are defined elsewhere by a set of hydrodynamic calculations (either using an empirical correlation based on experimental observation for a fluidized bed or a constant bed velocity for a moving bed) which effectively enforce conservation of volume in the system. We also know that the solids are mostly inerts and that the mass flowrate of inerts, F*x_inerts = v*(1-epsilon)*rho*x_inerts, must be conserved. Given that v and epsilon are defined to enforce conservation of volume, this implies that the overall mass concentration of inerts, rho*x_inerts, must also be constant; this implies that either rho and x_inerts must either both be constant or vary inversely to each other in order to conserve the inert materials.

If we look at the equations, we see that rho is defined as being constant, however x_inerts is expected to vary as the solid material reacts with the gas phase - this is the source of our singularity. When the original model was written, the modeler made the seemingly reasonable assumption that the adsorbed gases were negligible in comparison to the mass of the inert phase and that the overall density of the solids was thus effectively constant. However, when taken in combination with the other assumptions in the model (conservation of volume), this results in a set of incompatible assumptions which manifested as a structural singularity.

[Robbybp]

@andrewlee94 Thanks for the feedback! I will incorporate it when I get around to the other improvements I mentioned above.

[Robbybp]

Update old models (ported from IDAES 1.7) to use units so we don't get so many warnings about them (or just turn off logging, but this seems like poor form for an example)

I've spent a couple of days now trying to update units in these property packages, but it may not be worth the effort. For one, we will not be able to get rid of all the unit inconsistency warnings, as many are due to accumulation equations (see Pyomo/pyomo#1790). Also, changing units of some variables (pressure and solid phase energies) makes the model numerically unstable, and significantly increases the amount of time we have to wait for the model to initialize. I seem to remember this being pretty annoying to fix the last time I had to update units for this model.

For now, I'm planning to just leave the model as-is, disable the logger when running the structural diagnostics, and leave the remark about ignoring unit inconsistencies in the notebook.

[Robbybp]

A description of what "positive definite" and other more advanced terms mean might be helpful (or pointers to references).

I've tried to remove this jargon. I still assume some mathematical knowledge (e.g. "Jacobian is singular"), but I think it is more manageable.

A physical description of the singularity might help readers link the maths to the physics.

I've attempted to add this as well. The problem with particle porosity is that, as written, particle density and skeletal density are both uniquely determined by the material holdups, so there is no reason for the dens_mass_particle == (1 - porosity) * dens_mass_skeletal constraint. However, we know (or think) that this constraint should hold. Hence the contradiction. I've added a more detailed explanation to the notebook.

[andrewlee94]

@MarcusHolly @bpaul4 @JavalVyas2000 Would you be able to find some time to review this and provide comments? This would probably also be useful for you to learn more about structural singularities and how they can occur.

[Robbybp]

Some typos that I will fix when responding to other review comments.

[Robbybp]

When I build the docs locally, I get some slightly odd looking output:

Note how the output is broken up into separate cells. This gets worse when displaying solver output, where each line of the log is in a different output cell:

That said, this is probably due to the underlying Jupyter notebook conversion tool (i.e. not our fault), and I don't think it needs to hold up this PR.

[bpaul4]

Can these warnings be hidden, or resolved?

[Robbybp]

I'll update the doc notebook next time I'm on a computer with a working jupyter environment.

[Robbybp]

This is an old model, and only runs for the purpose of this example. It's probably not worth it to try to get rid of all of the warnings.

[bpaul4]

Do the gas-solid contactor files need to be duplicated here, or could they just be imported from https://github.com/IDAES/idaes-pse/tree/main/idaes/models_extra/gas_solid_contactors?

[Robbybp]

The files need to be duplicated. They are ported from IDAES 1.7, and include a "broken" version of the model, so they should not be distributed with the current version of IDAES. I also don't think it's a good idea (a) enforce that this notebook be run with IDAES 1.7 or (b) somehow import these file directly from 1.7.

[bpaul4]

Looks great, very interesting and I enjoyed reading through it. Nice job, @Robbybp!

[MarcusHolly]

I will try to read through this today or tomorrow

[agarciadiego]

LGTM, just a couple of minor typos

[MarcusHolly]

LGTM

[JavalVyas2000]

Looks good to me, apart from the first cell in structural_singularity.ipynb file. @Robbybp, another minor comment would be that the notebook are not ran but the markdown writeup is written such that it would click well with the outputs being shown. I would suggest running the notebook before we merge it.

[Robbybp]

Thanks for the reviews, all. I'll fix the couple remaining typos and update the auxiliary notebooks soon. We can merge then.

[Robbybp]

@JavalVyas2000 I think I'll leave the base notebook with output cleared. The _doc notebook will have output displayed for easy viewing, so I don't think the base notebook needs this as well.

[JavalVyas2000]

@Robbybp yes its fine. Even if one of the notebooks is ran it should be fine.

[Robbybp]

As far as I'm concerned, this is ready to merge once tests pass.