To debug Implicit ATR Flowsheet

[sbugosen]

This is the first iteration of the ATR Implicit Flowsheet. Everything runs up until the solver, where my kernel automatically dies.

Potential issues:

1. In better_implicit_flowsheet.py I don't build arcs to the Gibbs Reactor when I first build the flowsheet. For this ATR flowsheet, I'm not exactly sure how to handle this when there are recycle streams.
2. Lagrange multipliers for N and Ar.

[Robbybp]

The error I get when running this flowsheet is function cbrt not available. Is this what you get as well? cbrt is an external cube root function, and this error seems to be coming from C (ASL). Somewhere we have an external function that is not being properly communicated to a .nl file... I'll look into this.

[sbugosen]

I did not get any helpful error description. I got this:

Cannot execute code, session has been disposed. Please try restarting the Kernel. The Kernel crashed while executing code in the the current cell or a previous cell. Please review the code in the cell

[Robbybp]

Yep, that's almost definitly a segfault caused by this error happening in ASL. We can solve this flowsheet with Ipopt via the nl interface, so we know we both have the cbrt function installed and callable by ASL under the right circumstances, I just have to figure out how the PyNumero/CyIpopt interface differs from the nl interface.

There might be some difference caused by us using IDAES's Ipopt? I hope not.

[sbugosen]

I will try to keep searching for possible modeling errors in this flowsheet.

[Robbybp]

I've narrowed this down and it seems to be a bug in ASL when reading two nl files and the second uses different external function libraries than the first. See https://github.com/ampl/asl/issues/13. Hopefully the ASL maintainers can fix this quickly, otherwise we'll have to implement some workaround in Pyomo (not that difficult, but a little ugly).

[Robbybp]

The workaround isn't horrible, just the following at the top of the file:

import os
os.environ["AMPLFUNC"] = "\n".join(
    [
        "/path/to/.idaes/bin/cubic_roots.so",
        "/path/to/.idaes/bin/functions.so",
    ]
)

Ideally, we get these library paths directly from IDAES. The question is how to identify all the external functions we need before instantiating the implicit function. identify_external_functions from pyomo.util.subsystems may be useful.

This can lead to multiple libraries defining the same external function (because ASL will see e.g. the cubic_roots.so library twice). But I think this will not horribly break anything as it doesn't matter which library ASL chooses to call the function from.

[Robbybp]

After I added the workaround for the external function libraries, I got a KeyError because of external variables being included in the constraints of the "explicit model". I remedied this by just using the Gibbs reactor outlet arc equations as the residual equations. Now we don't have to add dummy variables and replace the "original variables" with them in the rest of the model.

However, the model has trouble converging. I have tried with a couple different initializations and variations of the inner Newton solver, but I either get an evaluation error, or a residual that Ipopt is unable to converge. This could be a bug in how we are constructing the implicit function (or how I am evaluating the function and its derivatives under the hood), or it could be something difficult about the implicit function system. I will continue to debug.

[sbugosen]

Sounds good. This is what I plan to do next week:

1. Monday: I will present to you a first draft of the abstract + poster for the LANL Symposium.
2. Rest of the week: close that ALAMO ATR PR I have pending, write the FOCAPD abstract, take a look at the implicit flowsheet again to see if I can come up with something, evaluate using OMLT-NN for the ATR Flowsheet.

[Robbybp]

I found a bug in the setup of the implicit function. I had the recuperator inlet arc's pressure equality included in the residual equations, even though pressure was not set as an input. This script now works for me and solves to the same value as the full-space flowsheet. @sbugosen can you try the script and let me know if it works for you?

The solve only seems to converge with good initialization (i.e. after solving a square problem via the decomposition). Otherwise I get an evaluation error. I'll work on handling evaluation errors in PyNumero and CyIpopt, but until then, if we want to use this, we'll have to use the "good initialization". If we do so, we should use this initialization for the other formulations (full-space and surrogate) to be fair.

[sbugosen]

That's awesome! It solves for me, but it takes some time. It solves in 32 iterations. (I used conversion of 0.94, same results as in full space).

What I don't understand is why there are 842 variables in the implicit formulation. The full space formulation has 897 variables. It makes sense that the implicit formulation has less variables, but I thought the difference would be more significant. Am I looking at this correctly?

[Robbybp]

That iteration count sounds right to me. Yep, it does take longer, especially once the imprecision in the implicit function solve starts to effect Ipopt's ability to reduce infeasibility, around 1e-5 or so.

And yes, we expect the flowsheet with the implicit reactor to have the number of variables of the original flowsheet minus the number of external variables. I believe there are 897 variables in the original flowsheet and 55 external variables, so 842 looks like the right number of variables to me.

[sbugosen]

Got it, yes, I forgot that was the number of external vars.

So main modifications to make the flowsheet work:

1. Lagrange multipliers of N and Ar (inerts).
2. Take out pressure from the residual eqn set because it is fixed.
3. Solve the square problem and use this as initialization for the optimization.
4. Change reformer outlet temperature constraint to shell inlet temperature because reformer outlet temperature is treated implicitly. (Not use dummy variables for the implicit function). Use the variables that exist within the flowsheet.

5. Add external libraries for ASL.

   I want to document this so that I'm clear on what to expect when we do a more complex flowsheet.

[Robbybp]

My to-dos on this PR:

• [x] remove redundant code (e.g. the "build model" function in the implicit flowsheet file)
• [ ] incorporate @sbugosen's initialization routine instead of solving strongly connected components
• [ ] add the mixer into the implicit function formulation

@sbugosen Can you (a) confirm that the surrogate models replace the mixer as well as the reformer and (b) point me to the initialization routine you have been using?

[sbugosen]

@Robbybp (a) confirm that the surrogate models replace the mixer as well as the reformer

Yes

(b) point me to the initialization routine you have been using?

When I used this good initialization routine for Full Space and Implicit:

def make_simulation_model(T,P,initialize=True):
    m = pyo.ConcreteModel(name="ATR_Flowsheet")
    m.fs = FlowsheetBlock(dynamic=False)
    build_atr_flowsheet(m)
    set_atr_flowsheet_inputs(m,T,P)
    if initialize:
        initialize_atr_flowsheet(m)

    m.fs.reformer_bypass.split_fraction[0, "bypass_outlet"].fix(0.3)
    m.fs.reformer_mix.steam_inlet.flow_mol.fix(230)
    solver = get_solver() # max iterations = 200
    solver.solve(m, tee=True)
    return m

The Full Space and the Implicit converged all instances. The make_simulation found an optimal solution in every instance, and this helped a lot for the optimization to be successful. To get around this small issue, I read your Implicit Function paper again and found this:

"Compared with simultaneous dynamic optimization, using implicit functions to solve optimization problems in a reduced space has several advantages. While simultaneous dynamic optimization formulations and nonlinear optimization solvers are relatively mature, their convergence is sensitive to an initial guess and to scaling of variables and constraints. We note that a reduced space formulation leaves a user with fewer variables to initialize and scale."

So to prove that initialization shouldn't matter much for the implicit function, I used this routine instead:

def make_simulation_model(T,P,initialize=True):
    m = pyo.ConcreteModel(name="ATR_Flowsheet")
    m.fs = FlowsheetBlock(dynamic=False)
    build_atr_flowsheet(m)
    set_atr_flowsheet_inputs(m,T,P)
    if initialize:
        initialize_atr_flowsheet(m)

    m.fs.reformer_bypass.split_fraction[0, "bypass_outlet"].fix(0.3)
    m.fs.reformer_mix.steam_inlet.flow_mol.fix(466.7)
    solver = get_solver() # max iterations = 200
    solver.solve(m, tee=True)
    return m

In many cases, this routine did not find an optimal solution, but the implicit optimization did, regardless of this. This is good.

My results are that full space converged only 38 out of 64 optimization instances, while implicit converged 53 out of 64. I tested this routine for ALAMO, and again, even if the initialization routine maxxed out and could not find an optimal solution, the ALAMO formulation was able to converge every instance.

[sbugosen]

This paper supports the fact that initializing the optimization problem with a feasible solution really helps convergence.

https://www.sciencedirect.com/science/article/pii/S0140700713002909#bib5