13

I am interested in the performance of Pyomo to generate an OR model with a huge number of constraints and variables (about 10e6). I am currently using GAMS to launch the optimizations but I would like to use the different python features and therefore use Pyomo to generate the model.

I made some tests and apparently when I write a model, the python methods used to define the constraints are called each time the constraint is instanciated. Before going further in my implementation, I would like to know if there exists a way to create directly a block of constraints based on numpy array data ? From my point of view, constructing constraints by block may be more efficient for large models.

Do you think it is possible to obtain performance comparable to GAMS or other AML languages with pyomo or other python modelling library ?

Thanks in advance for your help !

Dimitri Tomanos
  • 153
  • 2
  • 8
  • 1
    Usually GAMS is faster than Pyomo (and on my models also usually faster than AMPL especially when we have a lot of data manipulation steps but I know of models where AMPL is significantly faster). – Erwin Kalvelagen Apr 16 '17 at 18:31
  • Pyomo is very slow to generate models. In my problem, the program takes around 3 hours in total out of which solver (Gurobi) takes only 5 seconds. – Vinod Kumar Chauhan Jun 16 '20 at 10:09

3 Answers3

8

While you can use NumPy data when creating Pyomo constraints, you cannot currently create blocks of constraints in a single NumPy-style command with Pyomo. Fow what it's worth, I don't believe that you can in languages like AMPL or GAMS, either. While Pyomo may eventually support users defining constraints using matrix and vector operations, it is not likely that that interface would avoid generating the individual constraints, as the solver interfaces (e.g., NL, LP, MPS files) are all "flat" representations that explicit represent individual constraints. This is because Pyomo needs to explicitly generate representations of the algebra (i.e., the expressions) to send out to the solvers. In contrast, NumPy only has to calculate the result: it gets its efficiency by creating the data in a C/C++ backend (i.e., not in Python), relying on low-level BLAS operations to compute the results efficiently, and only bringing the result back to Python.

As far as performance and scalability goes, I have generated raw models with over 13e6 variables and 21e6 constraints. That said, Pyomo was designed for flexibility and extensibility over speed. Runtimes in Pyomo can be an order of magnitude slower than AMPL when using cPython (although that can shrink to within a factor of 4 or 5 using pypy). At least historically, AMPL has been faster than GAMS, so the gap between Pyomo and GAMS should be smaller.

jsiirola
  • 2,259
  • 2
  • 9
  • 12
  • Thanks a lot for your answer, I always had this doubt about Pyomo. But can you clarify the below part further? > [...] solver interfaces [...] are all "flat" representations that explicit[ly] represent individual constraints. [...] Pyomo needs to explicitly generate representations of the _algebra_ (i.e., the expressions) to send out to the solvers. In contrast, NumPy only has to calculate the result: it gets its efficiency by creating the data in a C/C++ backend [...] relying on low-level BLAS operations to compute the results efficiently, and only bringing the result back to Python. – Milind R Jun 21 '17 at 10:12
  • How is this different from the case with matrix representation vs. writing out simultaneous linear equations in scalar form? Surely the solver input files would have some matrix-ish way to denote large number of constraints? Or are you saying that's not the case because of sparsity? – Milind R Jun 21 '17 at 10:18
  • To my knowledge the standard solver input formats do not support a "matrix-ish" way to specify constraints. For linear solvers, the standard formats are LP (row based) and MPS (column based). Both explicitly list all the nonzeros in the "A" matrix. For general nonlinear, there really isn't a standard format. AMPL's NL format is the closest, and that is row based. Some convex solvers support matrix style input, but those are not formats that general nonlinear solvers typically take. – jsiirola Jun 22 '17 at 05:53
  • @jsiirola Just curious: how long does it take to generate your mentioned "raw models with over 13e6 variables and 21e6 constraints" in Pyomo? And how this compares to generation time AMPL? – Remis Jul 30 '19 at 15:15
6

I was also wondering the same when I came across this piece of code from Jonas Hörsch and Tom Brown and it was very useful to me:

https://github.com/FRESNA/PyPSA/blob/master/pypsa/opt.py

They define classes to define constraints more efficiently than the original Pyomo parser do. I did some tests on a large model that I have and it reduced the generation time considerably.

Jon Cardoso-Silva
  • 991
  • 11
  • 25
  • Thanks so much for the tip - do you know how much of a speed up it gave roughly for larger models and whether this was still the best appraoch for improving Pyomo's performance? Thanks! – Jack Simpson Sep 26 '22 at 11:54
3

You can build big linear (LP) and mix-integer (MILP) optimization problems in Python with the open-source tool Linopy. Linopy promises a speedup of times 4-6 and a memory reduction of roughly 50% reaching roughly Julia JUMP performance. See the benchmark:

enter image description here

The tool is part of the PyPSA ecosystem. This tool is the next-level version of the PyPSA 'opt.py' developments that Jon Cardodo mentioned. It has roughly the same speed, same performance but better usability -- reported by developers.

Max Parzen
  • 103
  • 1
  • 6