pyomo icon indicating copy to clipboard operation
pyomo copied to clipboard

Support Gurobi nonlinear expressions in the direct/persistent solvers

Open torressa opened this issue 10 months ago • 2 comments

Summary

Gurobi v12 supports nonlinear constraints, it would be good if the direct/persistent solvers supported adding these directly to the gurobipy.Model object.

Rationale

Right now, to use this feature through Pyomo, one must use the gurobinl solver which uses the AMPL interface via an NL file.

Description

I would like to leverage the Pyomo expressions to either create an expression tree to pass this to gurobipy via Model.addGenConstrNLAdv or even create a NLExpr whatever is easiest. This requires a bit of knowledge of the Pyomo expressions and I haven't had time to dig into it yet.

Additional information

This can be used to translate nonlinear functions such as pyomo.kernel.sin as well as polynomials, hence it would solve #3249

torressa avatar Mar 12 '25 17:03 torressa

I am very interested in this functionality and would be happy to contribute in some way. I can certainly help with testing and maybe with some coding, if someone with more knowledge of the internals of Pyomo and the gurobipy API could sketch a design.

grahamsparrow-8451 avatar Mar 12 '25 18:03 grahamsparrow-8451

I'm working on a direct Gurobi interface that supports MINLP as a contrib package on this branch: https://github.com/emma58/pyomo/tree/gurobi-minlp

It's not quite ready for people to beat on yet, but I can update this issue when it is.

emma58 avatar Mar 12 '25 19:03 emma58

With Gurobi v12 supporting expression trees, I believe a simpler approach is possible. Since interest in MINLP using Pyomo+Gurobi is increasing over the past weeks, I've made a start in https://github.com/ronaldvdv-gurobi/pyomo (sin/cos/exp/sum/product) and added two unit tests as well. Will continue to add the missing expression types and functions and add more tests. Please feel free to comment on this approach!

ronaldvdv-gurobi avatar Sep 26 '25 11:09 ronaldvdv-gurobi

@ronaldvdv-gurobi, thanks for starting this off. I am very interested in testing this on our MINLP problems. I had previously tried using the AMPL interface, but I could not be sure that was always using the most efficient methods of interacting with Gurobi and since I have my problem coded in Pyomo I am keen to try this out. I believe that there will be significant improvements in MINLP support when Gurobi 13 is release in the next couple of months.

grahamsparrow-8451 avatar Sep 26 '25 12:09 grahamsparrow-8451

Yes, MINLP is very much a focus area for v13. What non-linear functions are you most interested in testing? I could definitely use a few (small) examples as testcases.

ronaldvdv-gurobi avatar Sep 26 '25 14:09 ronaldvdv-gurobi

Here is an example of a very small problem: https://github.com/coin-or/SHOT/issues/166 if it is of any use. Generally, we have objective and constraint functions that have pow() or log()/exp().

grahamsparrow-8451 avatar Sep 26 '25 15:09 grahamsparrow-8451

@grahamsparrow-8451 I've just added that issue as a unit test and fixed a small bug.

ronaldvdv-gurobi avatar Sep 29 '25 09:09 ronaldvdv-gurobi

@ronaldvdv-gurobi, thank you for looking into this! While your approach is indeed a far simpler implementation, mainly for performance reasons, this ultimately needs to be designed as a writer that walks the Pyomo expression tree and builds up the expression in gurobipy. In addition (and apologies that this isn't well advertised at the moment), the dev team has plans to eventually drop support for the solver interfaces currently in pyomo.solvers: We are working on a rewrite of all our interfaces in pyomo.contrib.solvers. The plan for the moment is to implement a direct interface to gurobi that supports MINLP in contrib.solvers, but keep it separate from the current gurobi direct rewrite (again because of performance). We will probably integrate them eventually, but are planning that as a separate step. I'm getting very close for a first draft of that interface on the branch I mentioned above. I'll merge in your tests from your branch--they're a great addition. Hoping to have a PR open this week!

emma58 avatar Sep 29 '25 16:09 emma58

Hey @emma58 , really glad to hear your branch is still WIP and could be available soon! Combining your Pyomo expertise with the testcases from my branch seems like the best of both worlds, especially if you believe a PR could be open (and merged) short-term. The reason for doing a simplified attempt was the request from several clients for a working MINLP interface to Gurobi - any approach that gets them started quickly would be great. If you want to discuss anything on this work, don't hesitate to reach out to me directly. As you've probably guessed, I work at Gurobi and would be happy to collaborate on this.

Just curious, what's the reason for implementing a new Gurobi direct interface in the contrib.solvers folder instead of extending the existing one with your walker-approach? I know from experience people are struggling with picking the right Pyomo-Gurobi interface and maintaining one more interface (until the full rewrite is completed and merged) can be tricky. How would the timeline for the full rewrites compare to the availability of your PR? Is it worth having a separate MINLP interface?

ronaldvdv-gurobi avatar Sep 29 '25 18:09 ronaldvdv-gurobi

@ronaldvdv-gurobi, the reason to do a new one is that the current gurobi_direct_v2 only supports linear problems and can handle very large linear problems since it is aware of templatized Pyomo expressions. So I'm hesitant to pollute it with catching the case where things go nonlinear at the moment since we want to do that in such a way that we don't lose these gains. We're also getting close to expanding all that to include quadratic, so I think it makes sense to get that in first. However, once we have, all the pieces will exist, so perhaps it's not such a big lift to integrate them all together? We'd be very happy to discuss that with you at some point--and we do love help! And in the long term, I agree that a separate MINLP interface is clunky, but I do think that's the path to getting something usable in the short (next few week) time horizon.

emma58 avatar Sep 29 '25 19:09 emma58

@ronaldvdv-gurobi, if you have a few minutes, could I have your help? I merged in your tests (on this branch), and 4 of them are failing with my current implementation: In particular, I dug into test_gurobi_minlp_sincosexp and I can't find what about the resulting model is incorrect. The Pyomo model is:

m = ConcreteModel(name="test")
m.x = Var(bounds=(-1, 4))
m.o = Objective(expr=sin(m.x) + cos(2 * m.x) + 1)
m.c = Constraint(expr=0.25 * exp(m.x) - m.x <= 0)

and if I write the lp file for the resulting Gurobi model (using this code):

from pyomo.contrib.solver.solvers.gurobi_direct_minlp import GurobiDirectMINLP
from pyomo.opt import WriterFactory
grb_model, var_map, obj, grb_cons, pyo_cons = WriterFactory(
    'gurobi_minlp'
).write(m, symbolic_solver_labels=True)
grb_model.write('debug.lp')

I get:

\ LP format - for model browsing. Use MPS format to capture full model detail.                                                                                                                                               
Minimize                                                                                                                                                                                                                     
  0 x + C1                                                                                                                                                                                                                   
Subject To                                                                                                                                                                                                                   
R0: C2 <= 0                                                                                                                                                                                                                 
Bounds                                                                                                                                                                                                                       
-1 <= x <= 4                                                                                                                                                                                                                
General Constraints                                                                                                                                                                                                          
\ C1 = (sin(x) + cos(2 * x)) + 1                                                                                                                                                                                             
GC0: C1 = NL : ( PLUS , -1 , -1 ) ( PLUS , -1 , 0 ) ( SIN , -1 , 1 )                                                                                                                                                        
  ( VARIABLE , x , 2 ) ( COS , -1 , 1 ) ( MULTIPLY , -1 , 4 )                                                                                                                                                                
  ( CONSTANT , 2 , 5 ) ( VARIABLE , x , 5 ) ( CONSTANT , 1 , 0 )                                                                                                                                                             
\ C2 = (0.25 * exp(x)) + (-1 * x)                                                                                                                                                                                            
GC1: C2 = NL : ( PLUS , -1 , -1 ) ( MULTIPLY , -1 , 0 )                                                                                                                                                                     
  ( CONSTANT , 0.25 , 1 ) ( EXP , -1 , 1 ) ( VARIABLE , x , 3 )                                                                                                                                                              
  ( MULTIPLY , -1 , 0 ) ( CONSTANT , -1 , 5 ) ( VARIABLE , x , 5 )                                                                                                                                                           
End             

And indeed if I call optimize on that gurobi model, I get the same wrong answer as the Pyomo interface is currently getting. But I'm not spotting the error?

(Also, if you'd like to take this offline, please email [email protected], and I can connect with you directly.)

emma58 avatar Sep 30 '25 18:09 emma58

Yep, I've had the same error: it's (likely) because the auxiliary variables generated by the solver interface, have their default lowerbound of 0 instead of -GRB.INFINITY . To confirm I made that change in my own interface and get the same (wrong) answer for the model above. And likewise, with infinite bounds, inspecting the solution shows that it requires the auxiliary (first) variable to have a negative value. Could you try if introducing infinite bounds solves the issue with your interface?

ronaldvdv-gurobi avatar Sep 30 '25 19:09 ronaldvdv-gurobi