fe_factory_fix bug

[dthierry]

Hey @mchlshort I have tested both python 2 and 3 and it doesn't seem to be a problem. (Ex_5_sm_fe_by_fe.py).
Python 2

Python 3

Python 2

Python 3

Was there something in particular that didn't work out before?
~~Dav

[mchlshort]

Hi @dthierry, sorry for getting back to you late on this. Actually, I had the Michael's example working. The issue occurred on both of the aspirin sims. In the past I remember getting those profiles where the concentrations dropped below zero in an element, however they now seem to fail on the third element.

Can you please check those for me since mine are still failing? Additionally, while working on the FESimulator (the interface for Kipet and fe_factory) I noticed that these profiles have a strange jump at the second finite element. Everything seems ok from there onwards, but there may be an issue with the initialization or something?

I know you are away right now, so I will try to check it out. And enjoy your time away from the office

[dthierry]

Ah, I see. I think we weren't able to solve the Aspirin problem to begin with. There were issues like what you just reported.

[dthierry]

The problem is that the inputs don't match exactly. There is something making interpolations for us automatically. Turning off that element should help for now.

[mchlshort]

Hi @dthierry , I have been going over the PyomoSimulator class very closely today to try and figure out what the problem could be here. From what I can tell the problem is that, because the input data is not equally spaced, the interpolation is required in order to fix the input data in the correct places. Maybe just to give you the strategy that is being used here.

1. PyomoSimulator.apply_discretization set the places where all of the data needs to be solved for, i.e. model.time will consist of nfe*ncp +1 and partition the time space in relation to this (i.e. every ncp-th points is equally spaced)

2. PyomoSimulator._default_initialization fills the Z and X values with initializations (basically a constant at all model.time points that is the same as each components/states initial value) through using PyomoSimulator.initialize_from_trajectory.

3. PyomoSimulator.initialize_from_trajectory then in turn calls Simulator.interpolate_from_trajectory which just finds the nearest point to the data and does linear interpolation to ensure that the point has the value closest. For X and Z these are just constant values set at the initial value.

So then for the fix_from_trajectory this basically just uses the data file and fixes the relevant inputs to the values needed at the fe and ncp by using the interpolate_from_trajectory. Going through each dictionary, list and variable, I can't figure out what the problem in the set up is.

Maybe I am missing something here, but I am really struggling to figure out what the problem here could be and don't see how the interpolation could be turned off. Unless I am not quite getting which interpolations you are referring to.

[schenkch]

I also run into issues solving my problem without jumps increasing the number of finite elements which seems really strange to me. I can solve it for up to 30 finite elements but if I choose a higher number, the solution already fails in the zeroth element. It converges to a point of local infeasibility. I am not sure whether this is a pyomo problem or finite element factory problem though.

[mchlshort]

I am starting to think that it is a fe_factory issue. I have been carefully checking that all the parameters and variables are correctly initialized at each point in the pyomo model and it looks regular to me. When I input a low number of finite elements the model can solve, only if the PyomoSimulator solve is restricted to a single iteration, which makes me think that when it patches to the model something is off. As I increase the number of fe, it gets to a point where fe factory is able to work, however it fails when PyomoSimulator is called with an error evaluating one of the constraints. (again I can't seem to get the nl file for this since I am not sure how to get it to output when the solver fails).

Finally as the number of fes increase to a certain number (not very high) I get a failure at the first element. When this occurs the ipopt output looks quite strange:

"
***** Fe Factory: fe_initialize by DT @2018 *****
fe 0
Ipopt trunk:

---

This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt

---

This is Ipopt version trunk, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...: 331
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 54

Total number of variables............................: 82
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 82
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0

iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 1.80e+02 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
1 0.0000000e+00 1.77e+02 0.00e+00 -1.0 3.35e+02 - 1.00e+00 1.56e-02h 7
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
2 0.0000000e+00 1.75e+02 0.00e+00 -1.0 3.22e+02 - 1.00e+00 1.56e-02h 7

Warning: Cutting back alpha due to evaluation error
3 0.0000000e+00 1.75e+02 0.00e+00 -1.0 2.93e+02 - 1.00e+00 1.22e-04h 14
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
4 0.0000000e+00 1.75e+02 0.00e+00 -1.0 2.93e+02 - 1.00e+00 1.53e-05h 17
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
5r 0.0000000e+00 1.75e+02 9.99e+02 2.0 0.00e+00 - 0.00e+00 4.77e-07R 22
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Warning: Cutting back alpha due to evaluation error
6r 0.0000000e+00 1.75e+02 1.00e+03 2.0 9.21e+04 - 1.17e-02 2.55e-07f 13
Scaling factors are invalid - setting them all to 1.
Scaling factors are invalid - setting them all to 1.
WARNING: Problem in step computation; switching to emergency mode.
7r 0.0000000e+00 1.75e+02 1.00e+03 0.6 9.21e+04 20.0 0.00e+00 0.00e+00R 1
Scaling factors are invalid - setting them all to 1.
Scaling factors are invalid - setting them all to 1.
WARNING: Problem in step computation; switching to emergency mode.
Restoration phase is called at point that is almost feasible,
with constraint violation 0.000000e+00. Abort.
Restoration phase in the restoration phase failed.

Number of Iterations....: 7

                               (scaled)                 (unscaled)

Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 3.0096084412808239e-03 3.0096084412808239e-03
Constraint violation....: 9.6935181871682232e+01 1.7463551560456654e+02
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 9.6935181871682232e+01 1.7463551560456654e+02

Number of objective function evaluations = 14
Number of objective gradient evaluations = 6
Number of equality constraint evaluations = 85
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 8
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 8
Total CPU secs in IPOPT (w/o function evaluations) = 0.020
Total CPU secs in NLP function evaluations = 0.000

EXIT: Restoration Failed!
Traceback (most recent call last):

File "", line 1, in
runfile('/home/shortm/GrandKipetRepo/interp_off/kipet/kipet/examples/aspirin_sim_fe_factory_new_try.py', wdir='/home/shortm/GrandKipetRepo/interp_off/kipet/kipet/examples')
File "/home/shortm/anaconda2/lib/python2.7/site-packages/spyder/utils/site/sitecustomize.py", line 705, in runfile
execfile(filename, namespace)
File "/home/shortm/anaconda2/lib/python2.7/site-packages/spyder/utils/site/sitecustomize.py", line 94, in execfile
builtins.execfile(filename, *where)
File "/home/shortm/GrandKipetRepo/interp_off/kipet/kipet/examples/aspirin_sim_fe_factory_new_try.py", line 271, in
init.run()
File "/home/shortm/anaconda2/lib/python2.7/site-packages/kipet-1.0.0-py2.7.egg/kipet/library/fe_factory.py", line 441, in run
try:
File "/home/shortm/anaconda2/lib/python2.7/site-packages/kipet-1.0.0-py2.7.egg/kipet/library/fe_factory.py", line 333, in march_forward
"""
File "/home/shortm/anaconda2/lib/python2.7/site-packages/Pyomo-5.4.3-py2.7.egg/pyomo/opt/base/solvers.py", line 657, in solve
default_variable_value=self._default_variable_value)
File "/home/shortm/anaconda2/lib/python2.7/site-packages/Pyomo-5.4.3-py2.7.egg/pyomo/core/base/PyomoModel.py", line 249, in load_from
% str(results.solver.status))
ValueError: Cannot load a SolverResults object with bad status: error"

No idea how any of this is related at this stage. I think that when @dthierry is back on Tuesday we should sit down and try to figure out how all this fits together.

[dthierry]

We are still talking about the aspirin problem right?

[schenkch]

Yes but it seems so be the same problem for the aspirin problem and my problem.

[mchlshort]

Small update here:

Been testing with large numbers of elements and fe_factory will fail on the first element, even on the Michael's simulation, when you input a large number of elements(somewhere between 300 and 400). The same is true for the aspirin example. Something causes the first element to fail. I am not sure whether this is something to do with a derivative approaching a very large number or some other issue. When looking closely at the profiles for the Michael's sim, you can see a sharp change from the initial concentrations to another concentration, before a smooth profile is obtained. Is it related to this somehow? Perhaps fe_factory is not constructing the first element correctly?

[dthierry]

That is odd. I will check it out. Is this the fix_fe_fact branch?

[dthierry]

Hey Mike, I couldn't find a problem with the Ex_5; even with 2000 FEs the problem seems to be fine. Also, the jump seems to be consistent with the results.