IPOPT restoration failed

[Bodhi Biswas]

Hello,

I apologize in advance for the following question(s).  I am not quite sure what the source of my problem is, so I will provide all information that may be necessary:

I am trying to model a cell signaling pathway.  It comprises of 25 ODE constraints (over 1500 time points) and ~25 rate coefficients which I am trying to estimate using Pyomo and certain experimental values.

The model is large with ~90,000 variables.  I am using Pyomo along with IPOPT and Pardiso (seems to be the fastest solver).

I already have known rate coefficient values that are accurate/feasible.  I can enter them into the model in a normal Python script and receive expected values.  In this sense, I know that my model is correct.  

Now, I have converted this model into a Pyomo script and initialized starting values using the known rate coefficients.

The Pyomo model does not work.  

First red flag: from the first iteration, the unscaled constraint violation is ~1e9.  I would expect this to be nearly 0, as I have initialized with feasible values. (I also have the "warm start" option on).

Second red flag: After many iterations (~300), I receive values that are completely incorrect.  Shouldn't Pyomo just use the values I have initialized with, since they are the closest to convergence feasibly?  (If so, the solver should also come to an answer in very few iterations, no?) If I narrow the bounds too much to basically force the answer I want, I receive "restoration failed".  Other times, the solver will return rate coefficient values that are correct (due to aggressive bounding), but all other values depending on the rate coefficients will be way off!

Third red flag:  The model values shouldn't exactly hit the experimental values.  In this case, a reasonable convergence tolerance is ~20.  I have set 'tol' to a generous 1e3 just to see if the model will terminate.  IPOPT is repeatedly iterating solutions with objective convergence values of ~1e2, but it STILL does not terminate.  The only way my model ever terminates is through "restoration failed", "converging to a point of local infeasibility,"  or by hitting a set max iter.  Am I not specifying the convergence tolerance correctly?

This is how I start the solver:

results = opt.solve(disc_instance,tee=True, options={'max_iter': 300, 'acceptable_tol':1e3,'tol': 1e03, 'dual_inf_tol': 1e10, 'compl_inf_tol': 1e-1, 'warm_start_init_point': 'yes', 'linear_solver': "pardiso" 'bound_frac': 1e-20, 'bound_push': 1e-20}, report_timing=True)

Conclusion:

Any suggestion as to what my problem could be?  Do I have too many constraints?  I don't think I'll be able to decrease the # of constraints much further...

Thank you, any suggestion is greatly appreciated!

Bodhi

PS: I am also open to sharing my code if needed.

[jose santiago rodriguez]

Hello,

Regarding the first "flag" you could use the following methods to check which constraints are the ones that are being violated.

def CheckInstanceFeasibility(instance, tolerance):

    print

    print "*** Feasibility check:"

    infeasibility_found = False

    for block in instance.block_data_objects():

        for con in block.component_data_objects(Constraint, active=True):

            resid = ComputeConstraintResid(con)

            if (resid > tolerance):

                infeasibility_found = True

                print con.cname(True), con, resid

    if infeasibility_found == False:

        print "   No infeasibilities found."

    print "***"

    print

def ComputeConstraintResid(con):

    bodyval = value(con.body)

    upper_resid = 0

    if con.upper is not None:

        upper_resid = max(0, bodyval - value(con.upper))

    lower_resid = 0

    if con.lower is not None:

        lower_resid = max(0, value(con.lower) - bodyval)

    return  max(upper_resid, lower_resid)

initialize your instance with your feasible solution and see which constraints are the ones that are giving you problems. That will give you a start for debugging your model. 

For the second red flag:

Yes IPOPT should take very few iterations if you are initializing near the optimal solution. However, if some of your variables are being initialize near the bounds IPOPT might be changing that initial value. To see if that is happening you can pass the option 'bound_push' to the solver as follows:

opt = SolverFactory('ipopt')

opt.options['bound_push'] = ... 

However the fact that IPOPT is going into restoration suggests that there is some problem in your model.

Could it be that your ODE system is stiff? how are you discretizing the differential equations? To check if the discretization is the problem i would suggest looking into the pyomo DAE toolbox. The discretization step is usually error prompting so with pyomo DAE you can avoid discretizing the ode on your own plus you can easily try different discretization techniques.

Check the documentation

https://software.sandia.gov/downloads/pub/pyomo/PyomoOnlineDocs.html#_dae_toolbox

These are some slides that are useful too

http://www.mcs.anl.gov/MACS/files/2014/09/Nicholson_SASSY_slides-1-2.pdf

Hope that helps,

-Santiago

[Bodhi Biswas]

Hi Santiago,

I used the code you provided to check which constraints are being violated.  Using this, I was able to decrease unscaled constraint violation somewhat.  There are many constraint violations that I don't see how I can fix.

I am already using a very low 'bound push' options.

I do believe it has something to do with employing a stiff ODE system.  From the DAE toolbox, I have tried employing both the 'finite difference' and 'collocation discretization' methods.  Still, I am eventually met with 'restoration failed'.

[Bethany Nicholson]

Bodhi,

I'd recommend first trying to simulate your system using pyomo.dae to ensure your system of ODEs is solvable using a particular discretization scheme. Do this by fixing all your input variables or degrees of freedom. The problem you send to ipopt should have the same number of variables as constraints. Are you supplying boundary conditions or just initial conditions? Have you checked the index of your system?

Bethany

On Mar 20, 2016 6:29 PM, "Bodhi Biswas" <bodhi...@gmail.com> wrote:

Hi Santiago,

I used the code you provided to check which constraints are being violated.  Using this, I was able to decrease unscaled constraint violation somewhat.  There are many constraint violations that I don't see how I can fix.

I am already using a very low 'bound push' options.

I do believe it has something to do with employing a stiff ODE system.  From the DAE toolbox, I have tried employing both the 'finite difference' and 'collocation discretization' methods.  Still, I am eventually met with 'restoration failure'.

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.