Efficiency improvements with Pyomo/Gurobi

[Goran V]

Hello,

My question is about running a scatter search algorithm with Python/Pyomo. The scatter search is done by fixing the value of a single variable, a threshold. I am running the 64-bit version of Pyomo and using the 64-bit version of Gurobi to solve the problem.

The example problem that I am trying to solve is an approximation of a multi-stage stochastic problem. The original problem is hourly in granularity and needs to have at least 24 periods, thus as a multi-stage stochastic problem it explodes really quickly. We generated a number of scenarios and are estimating the problem with a rolling-horizon approach within every scenario.

In the instance that I am now running, there are 5 scenarios, and within each scenario there are 24 IPs that are solved sequentially. So to evaluate a single value of the threshold, I need to solve 120 IPs. I will refer to each set of 120 IPs as one evaluation.

It takes about 23 seconds to solve one evaluation. I want to see if that time can be improved. The reason I think it can be improved is that when I look at the output of the solver, it finds the optimal solution of every IP at the root node (it doesn’t form a tree) in 0.01 seconds. Hence, for 120 IPs it should take about 1.2 seconds. I am guessing Pyomo or Gurobi or both spend a lot of time setting up the problem, feeding it to the solver, and getting the results back. Is there any way to cut on the setup time? From your experience, can you tell me which is more likely to cause a bottleneck here: Pyomo or Gurobi?

More details about my algorithm: I create only one abstract model as a global variable. Before I solve the instance, I update the values of some parameters (Python constructs used as parameters in the model), create the instance from the global abstract model, and activate/deactivate some of the constraints. Does the fact that I create a new instance every time explain the long setup time? And can it be accelerated?

Thank you very much!

Goran

[Watson, Jean-paul]

Goran:

I would suggest profiling your script (I'm assuming you are running a script), in order to get some solid numbers concerning what is taking the lion's share of the run-time. Please let us know if you need some pointers / code fragments on how to accomplish this.

jpw

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

[Siirola, John D]

JP is right: profiling the script would be useful.  That said, some possible options for speeding things up (depending on the details of your model) come to mind:

-          use mutable Params for the parameter values that you expect to change.  That way you only create the concrete model (from the abstract one) once, and then update the mutable params directly on the concrete model.

-          use the gurobi direct interface to Gurobi (that is the Python bindings and not the default LP file interface).

 

john

[Goran V]

Gentlemen:

Thank you for the quick replies and the advice. I was not familiar with profiling until now. I did some basic profiling of a modified version of my algorithm. This script is not the full scatter search, but it does make the same exact function calls when solving the 120 IPs per evaluation. Going off of what I read on the python.org documentation, when inspecting the output I believe I should be looking at the internal time and not the cumulative time (because the latter includes time spent on recursive function calls). Below is the list of the top 20 function calls by internal time. The one that stands out is obviously WaitForSingleObject. It is called 120 times, so every time I optimize an instance. I am not sure what this function does, can you help me understand? Is there anything else beyond this very basic analysis of the output I should be looking at?

John, I will be trying out your recommendations a little later and I will let you know how they work out.

Thank you very much for your time,

Goran

Fri Nov 29 14:41:17 2013    thrst

         24074128 function calls (23601523 primitive calls) in 38.625 seconds

   Ordered by: internal time

   List reduced from 723 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)

      120   11.791    0.098   11.791    0.098 {_subprocess.WaitForSingleObject}

     1704    5.558    0.003    5.558    0.003 {method 'acquire' of 'thread.lock' objects}

   364560    1.807    0.000    2.866    0.000 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\expr.py:935(generate_expression)

   200398    1.389    0.000    3.358    0.000 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:558(generate_canonical_repn)

   117958    1.367    0.000    1.833    0.000 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\io\cpxlp.py:158(_print_expr_canonical)

      720    1.156    0.002    1.156    0.002 {method 'close' of 'file' objects}

      120    1.009    0.008    4.379    0.036 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\io\cpxlp.py:365(_print_model_LP)

179669/167909    0.972    0.000    1.089    0.000 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:398(_collect_linear_sum)

     4560    0.799    0.000    4.586    0.001 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\preprocess\compute_canonical_repn.py:110(preprocess_constraint)

   105960    0.687    0.000    1.443    0.000 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\constraint.py:290(add)

   919794    0.406    0.000    0.406    0.000 {method 'write' of 'file' objects}

    88200    0.395    0.000    0.945    0.000 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\expr.py:1327(generate_relational_expression)

      120    0.338    0.003    0.561    0.005 d:\packages\coopr\src\coopr.plugins\coopr\plugins\solvers\GUROBI.py:300(process_soln_file)

316320/120    0.293    0.000    0.949    0.008 D:\Python 27\lib\copy.py:145(deepcopy)

  1075613    0.255    0.000    0.255    0.000 {isinstance}

179669/167909    0.255    0.000    0.338    0.000 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\expr.py:396(polynomial_degree)

  2742466    0.230    0.000    0.230    0.000 {id}

      120    0.223    0.002    0.375    0.003 d:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\PyomoModel.py:664(_load_solution)

      840    0.216    0.000    0.216    0.000 {nt.remove}

      480    0.212    0.000    0.212    0.000 {nt.open}

[Watson, Jean-paul]

The WaitForSingleObject must be the call to the external solver, i.e., the system call. That I understand. What I don't understand is the 6 seconds associated thread.lock. Are you multi-threading in your script, using the Python threading libraries?

The former time could be mitigated by using (at least in the case of Gurobi) the python bindings – our "gurobi-direct" version. The CPLEX python bindings appear to be no faster, at least for us, than the file-based interface. Which solver are you using? In any case, using the Gurobi direct interface would mitigate much of the print-lp function time. 

The other chunk of time relates to the preprocessing time within pyomo. While necessary, this time can be eliminated if you are incrementally modifying and solving a particular instance. Does your scatter search fall into this category, e.g., are you making minor modifications to your sub-problems as you go? If so, then we can discuss ways to mitigate that time.

jpw

 

From: Goran V <goranvo...@gmail.com>
Reply-To: "coopr...@googlegroups.com" <coopr...@googlegroups.com>
Date: Friday, November 29, 2013 1:19 PM
To: "coopr...@googlegroups.com" <coopr...@googlegroups.com>

--

[Goran V]

I am using Gurobi. I believe John already recommended I use the gurobi-direct version, and I do want to try it. I am not too familiar with this, can you point me in the right direction about how I would implement it?  Are you referring to the below Gurobi/Python interface?

http://www.gurobi.com/documentation/5.6/quick-start-guide/py_python_interface

To your second point, yes, I am making changes to my sub-problems and re-solving them a certain number of times. Particularly, I need to update the values of some parameters in the model (I do this through a numpy array and three global floats that I use in the model), and activate/deactivate some constraints. As explained in my original post, every time I update the parameters I also re-create an instance of the model. If there are ways to cut on the preprocessing time, I would like to try them out.

I have also just tried creating one instance of the model, solving it, updating the parameters and constraints, and then solving it again and repeating the process in a loop. Instead of creating an instance at every iteration of a loop, this just creates one instance for the whole loop, but the results of the optimizations are not consistent (when compared to creating an instance for every iteration of the loop). More precisely, the results of the first optimization are the same but the following results are not.

Please let me know if some code or pseudo-code would make it easier to understand, and I can send it by email.

[Gabriel Hackebeil]

Goran, 

I am using Gurobi. I believe John already recommended I use the gurobi-direct version, and I do want to try it. I am not too familiar with this, can you point me in the right direction about how I would implement it?  Are you referring to the below Gurobi/Python interface?

I think what John was referring to is Pyomo's solver interface to the Gurobi Python interface. If you have the Gurobi python module (gurobipy) installed, then Pyomo can send your model to Gurobi using this interface rather than by writing an LP file. So you don't actually have to implement anything new, just check out the gurobi docs for how to install this interface into your python packages. If the gurobipy module can be successfully imported then you can create an instance of the gurobi-python solver plugin by adding the solver_io='python' keyword declaration to SolverFactory, e.g.,

# Python interface

# Returns None if gurobipy module cannot be imported

opt = SolverFactory("gurobi", solver_io='python')

# LP file interface

# Requires gurobi.sh be in your PATH

opt = SolverFactory("gurobi") # ,solver_io='lp')

I have also just tried creating one instance of the model, solving it, updating the parameters and constraints, and then solving it again and repeating the process in a loop. Instead of creating an instance at every iteration of a loop, this just creates one instance for the whole loop, but the results of the optimizations are not consistent (when compared to creating an instance for every iteration of the loop). More precisely, the results of the first optimization are the same but the following results are not.

I think I have an idea of what's going wrong here. From what you describe, it sounds like you are not using a Param for the coefficients that you mean to update each time you solve your model. In order to update coefficients that appear in Pyomo Constraint and Objective expressions (after the model is constructed), one must place these values inside mutable Params. E.g.,

# A scalar Param

model.p = Param(mutable=True, ...)

# An indexed Param

model.s = Set()

model.q = Param(model.s, mutable=True, ...)

Once you update the value of a mutable parameter on your concrete instance, you must call instance.preprocess() before solving again. E.g.,

instance = model.create()

...

# update coefficients

instance.p = 5.0

instance.q[1] = 2.3

instance.preprocess()

# Solve again

results = opt.solve(instance)

instance.load(results)

Gabe

[Goran V]

Gabe, JP, and John:

Thank you for your help. I implemented the mutable Params as you recommended and it resulted in an improvement in speed of about 25%. The implementation of the gurobi direct interface improved the speed further resulting in a total improvement of about 60%. Given that I need to solve different instances of this problem many times, this is a huge help. Thank you once again for the quick replies!

Goran

On Wednesday, November 27, 2013 4:46:30 PM UTC-5, Goran V wrote:

[Watson, Jean-paul]

Great! Out of curiosity, would you mind re-posting your profile – following these modifications? Just wondering what was left in terms of performance hot-spots…

jpw

From: Goran V <goranvo...@gmail.com>
Reply-To: "coopr...@googlegroups.com" <coopr...@googlegroups.com>
Date: Friday, December 6, 2013 6:43 AM
To: "coopr...@googlegroups.com" <coopr...@googlegroups.com>
Subject: [EXTERNAL] Re: Efficiency improvements with Pyomo/Gurobi

--

[Goran V]

JP,

Attached are three profiles labeled accordingly. I had to rerun them all because the computer on which I ran the profile a week or so ago crashed. The machine on which I run the profile now takes much longer to go through the code. There is a big reduction in the number of function calls and solving time. Thanks again!

Original profile with creation of 120 instances per evaluation and LP file interface

Mon Dec 09 08:52:00 2013    thrst

         25119174 function calls (24752447 primitive calls) in 102.638 seconds

   Ordered by: internal time

   List reduced from 824 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)

      120   16.574    0.138   16.574    0.138 {_subprocess.WaitForSingleObject}

   364560    6.556    0.000   12.947    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\expr.py:935(generate_expression)

   200398    4.376    0.000   13.674    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:559(generate_canonical_repn)

   117958    4.341    0.000    7.376    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\plugins\io\cpxlp.py:158(_print_expr_canonical)

      120    4.129    0.034   23.514    0.196 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\plugins\io\cpxlp.py:371(_print_model_LP)

179669/167909    2.978    0.000    4.395    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:398(_collect_linear_sum)

  2698546    2.976    0.000    2.976    0.000 {id}

     4560    2.307    0.001   18.470    0.004 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\preprocess\compute_canonical_repn.py:110(preprocess_constraint)

   105960    2.180    0.000    5.951    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\constraint.py:332(add)

    88200    1.749    0.000    4.980    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\expr.py:1327(generate_relational_expression)

   919914    1.678    0.000    1.678    0.000 {method 'write' of 'file' objects}

      722    1.576    0.002    1.576    0.002 {method 'close' of 'file' objects}

267480/120    1.444    0.000    4.045    0.034 C:\Python27\lib\copy.py:145(deepcopy)

   818855    1.417    0.000    1.417    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\sparse_indexed_component.py:175(__getitem__)

   906557    1.302    0.000    1.302    0.000 {isinstance}

179669/167909    1.184    0.000    1.845    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\expr.py:396(polynomial_degree)

953050/939164    1.088    0.000    1.126    0.000 {len}

   740636    0.974    0.000    0.974    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:589(<genexpr>)

   834011    0.973    0.000    0.973    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\numvalue.py:306(is_indexed)

   764687    0.962    0.000    0.962    0.000 {method 'append' of 'list' objects}

Profile with creation of only 1 instance per evaluation instead of 120 (using mutable Params) and LP file interface

Mon Dec 09 09:08:10 2013    thrst

         12196593 function calls (12149317 primitive calls) in 60.038 seconds

   Ordered by: internal time

   List reduced from 803 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)

      120   16.824    0.140   16.824    0.140 {_subprocess.WaitForSingleObject}

   117958    4.298    0.000    7.315    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\plugins\io\cpxlp.py:158(_print_expr_canonical)

      120    4.122    0.034   23.515    0.196 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\plugins\io\cpxlp.py:371(_print_model_LP)

    95202    2.158    0.000    7.767    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:559(generate_canonical_repn)

  1760973    1.942    0.000    1.942    0.000 {id}

   919914    1.664    0.000    1.664    0.000 {method 'write' of 'file' objects}

94839/88910    1.614    0.000    3.038    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:398(_collect_linear_sum)

      603    1.519    0.003    1.519    0.003 {method 'close' of 'file' objects}

     2299    1.168    0.001   10.156    0.004 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\preprocess\compute_canonical_repn.py:110(preprocess_constraint)

   797448    1.113    0.000    1.113    0.000 {isinstance}

   170878    0.905    0.000    1.509    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\symbol_map.py:250(createSymbol)

      120    0.897    0.007    1.940    0.016 c:\packages\coopr\src\coopr.solvers\coopr\solvers\plugins\solvers\GUROBI.py:300(process_soln_file)

   266648    0.826    0.000    1.258    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\sparse_indexed_component.py:170(iteritems)

   248398    0.751    0.000    2.397    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\block.py:1021(active_subcomponents_data_generator)

      120    0.741    0.006    1.880    0.016 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\PyomoModel.py:689(_load_solution)

   437552    0.734    0.000    0.734    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\sparse_indexed_component.py:175(__getitem__)

   153120    0.702    0.000    1.128    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\var.py:178(ub)

   435120    0.684    0.000    0.684    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\var.py:109(domain)

94839/88910    0.664    0.000    1.144    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\expr.py:396(polynomial_degree)

   212276    0.659    0.000    1.026    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:633(canonical_degree)

Profile with creation of only 1 instance per evaluation and gurobi direct interface

Mon Dec 09 09:01:04 2013    thrst

         9614432 function calls (9568416 primitive calls) in 39.365 seconds

   Ordered by: internal time

   List reduced from 674 to 20 due to restriction <20>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)

      120   10.889    0.091   21.008    0.175 c:\packages\coopr\src\coopr.solvers\coopr\solvers\plugins\solvers\gurobi_direct.py:178(_populate_gurobi_instance)

    95202    2.120    0.000    7.592    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:559(generate_canonical_repn)

94839/88910    1.582    0.000    2.965    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:398(_collect_linear_sum)

  1260575    1.405    0.000    1.405    0.000 {id}

      120    1.327    0.011    1.991    0.017 c:\packages\coopr\src\coopr.solvers\coopr\solvers\plugins\solvers\gurobi_direct.py:568(_apply_solver)

     2299    1.140    0.000    9.933    0.004 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\preprocess\compute_canonical_repn.py:110(preprocess_constraint)

   170878    0.995    0.000    1.695    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\symbol_map.py:250(createSymbol)

      120    0.961    0.008    1.576    0.013 c:\packages\coopr\src\coopr.solvers\coopr\solvers\plugins\solvers\gurobi_direct.py:636(_postsolve)

   694259    0.828    0.000    0.828    0.000 {method 'append' of 'list' objects}

   264238    0.736    0.000    1.024    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\symbol_map.py:141(getSymbol)

   458640    0.735    0.000    0.735    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\var.py:109(domain)

      120    0.732    0.006    1.896    0.016 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\PyomoModel.py:689(_load_solution)

   153000    0.712    0.000    1.156    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\var.py:121(lb)

94839/88910    0.652    0.000    1.109    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\expr.py:396(polynomial_degree)

   170878    0.503    0.000    0.503    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\symbol_map.py:37(__call__)

   254792    0.488    0.000    0.488    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\sparse_indexed_component.py:175(__getitem__)

   361740    0.477    0.000    0.477    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\expr\canonical_repn.py:589(<genexpr>)

   132240    0.461    0.000    0.709    0.000 c:\packages\coopr\src\coopr.pyomo\coopr\pyomo\base\sparse_indexed_component.py:165(itervalues)

    88354    0.459    0.000    0.769    0.000 C:\Python27\lib\site-packages\six.py:437(iteritems)

   318613    0.454    0.000    0.454    0.000 {isinstance}

[Watson, Jean-paul]

Hi Goran,

I've stared at these profiles for a bit, and the general progression is logical – or at least consistent with other profiles I've examined. 

Looking at the last profile, the bottleneck is still in Pyomo – gurobi is spending only 2 seconds actually solving the instances. The rest is overhead in either Gurobi/Python (setting up the problem instance) or Pyomo (compiling, for lack of a better name, the instance into a flattened representation – which facilitates transformation into Gurobi/Python). 

To address this issue – it's popped up in some of my applications work and those of collaborators – I'm developing something called "persistent solver plugins". These are worth using in situations where only a small fraction of the instance changes between solves – specifically, a fraction of the constraints. Does the interaction pattern between your (mutable) parameters and constraints fall into this category? In other words, can you easily identify the subset of constraints to re-compile following a parameter update?

I expect to have a prototype of this functionality available soon after the winter break, assuming it would be applicable in your case.

jpw

From: Goran V <goranvo...@gmail.com>
Reply-To: "coopr...@googlegroups.com" <coopr...@googlegroups.com>
Date: Monday, December 9, 2013 7:16 AM
To: "coopr...@googlegroups.com" <coopr...@googlegroups.com>
Subject: [EXTERNAL] Re: Efficiency improvements with Pyomo/Gurobi

[Goran V]

Hi JP,

Yes, the interaction pattern in my case falls in that category. It would be great to have a chance to try out your prototype. I will check back with you in a month or so on its availability. Please let me know if you need a test case, I'll be more than happy to try implementing it.

Goran

[Watson, Jean-paul]

Sounds like an excellent plan. 

[Goran V]

Hi JP,

I'm very sorry about not reaching back sooner, I realize a ridiculous amount of time has passed, I got "distracted" by my dissertation. I started using Pyomo again recently and I just remembered this. Anyway, if the above is still relevant, I am running an implementation of the L-shaped method, and so I need to add a cut at each iteration, and therefore preprocess. If I could preprocess just the constraint I am adding, it would save some time. I got some ideas from the forum and I'll be replying to another thread about it, but I also wanted to check whether your persistent solver plugins could be applicable in this case.

Thanks!

[Watson, Jean-Paul]

We unfortunately don’t yet have support for per-constraint preprocessing in the persistent interfaces. It will be a summer project to re-design those, to support both incremental processing and other functionality…

jpw

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.