Differences
This shows you the differences between two versions of the page.
personal:blog:2014:0904_pyomo_for_gams_users [2014/09/12 17:03] antonello [Creation of the Model] |
personal:blog:2014:0904_pyomo_for_gams_users [2018/06/18 15:11] |
||
---|---|---|---|
Line 1: | Line 1: | ||
- | ====== A Pyomo tutorial for GAMS users ====== | ||
- | |||
- | [[https:// | ||
- | It is developed by the Sandia National Laboratories and appeared in 2008 as an open source project. | ||
- | |||
- | [[http:// | ||
- | |||
- | This mini-tutorial is intended for gams users that want to try Pyomo. There may be two reasons for someone to with to use Pyomo instead of GAMS.\\ | ||
- | The most obvious one, even if often it isn't the key driver, is that GAMS is a commercial software while Pyomo being open-source is free both as freedom and as a free beer.\\ | ||
- | While for GAMS a licence for the underlying solver engine is often included with a particular version of GAMS, Pyomo would still require the user to buy a licence to use a specific commercial solvers. However Pyomo interfaces with both [[https:// | ||
- | The second reason (and, to me, the most important one) resides in the language features and in the availability of development environments. GAMS uses a VERY ODD syntax, somehow derived from the Cobol language, that is very distant from any programming language in use nowadays. For example a macro mechanism to provide an elementary way to structure the code in reusable components has been introduced only in GAMS 22.9. | ||
- | Its own editor is also very terrible, but as most text editors do not provide a gams syntax highlighting, | ||
- | |||
- | Pyomo, at the opposite, is both open source and.. it's python!\\ | ||
- | You have plenty of development environment to choose from, a clear modern language, the possibility to interface your model with third party libraries.. all of this basically for free.\\ | ||
- | While there are some reports of pyomo being somehow slower that GAMS it really depends. In my case it is actually much faster, as the IPOPT version that is embedded in GAMS uses the MUMPS linear solver, while on my system I have IPOPT compiled with the much more performing ma27 linear solver. That's part of the flexibility you gain in using pyomo in place of GAMS. | ||
- | |||
- | |||
- | So let's start. We will see how to code the trasnport.gms problem, the one that ship as default example in GAMS((yes, the default GAMS example in named " | ||
- | GAMS equivalent code is inserted as single-dash comments. The original GAMS code needs slightly different ordering of the commands and it's available at [[http:// | ||
- | |||
- | ===== Installation ===== | ||
- | |||
- | **Important: | ||
- | |||
- | ==== Ubuntu ==== | ||
- | //(tested in Ubuntu 14.04 LTS)// | ||
- | |||
- | * **Install the python pre-requisites: | ||
- | * '' | ||
- | * **Install pyomo:** | ||
- | * '' | ||
- | * '' | ||
- | * **Install solvers:** | ||
- | * //linear and MIP solver (glpk)//: '' | ||
- | * // | ||
- | |||
- | ==== Windows and Mac ==== | ||
- | Please refer to the [[https:// | ||
- | |||
- | ===== Model components ===== | ||
- | |||
- | ==== Creation of the Model ==== | ||
- | In pyomo everything is an object. The various components of the model (sets, parameters, variables, constrains, objective..) are all attributes of the main model object while being objects themselves.\\ | ||
- | There are two type of models in pyomo: A '' | ||
- | The first thing to do in the script is hence to load the pyomo library and to create a new ConcreteModel (we have little imagination here, and we call our model " | ||
- | <code python> | ||
- | # Import of the pyomo module | ||
- | from coopr.pyomo import * | ||
- | | ||
- | # Creation of a Concrete Model | ||
- | model = ConcreteModel() | ||
- | </ | ||
- | |||
- | ==== Set definition ==== | ||
- | Sets are created as attributes object of the main model objects and all the information is given as parameter in the constructor function. Specifically, | ||
- | <code python> | ||
- | ## Define sets ## | ||
- | # Sets | ||
- | # | ||
- | # | ||
- | model.i = Set(initialize=[' | ||
- | model.j = Set(initialize=[' | ||
- | </ | ||
- | |||
- | ==== Parameters ==== | ||
- | Parameter objects are created specifying the sets over which they are defined and are initialised with either a python dictionary or a scalar: | ||
- | <code python> | ||
- | ## Define parameters ## | ||
- | # | ||
- | # | ||
- | # / | ||
- | # san-diego | ||
- | # | ||
- | # / | ||
- | # chicago | ||
- | # topeka | ||
- | model.a = Param(model.i, | ||
- | model.b = Param(model.j, | ||
- | # Table d(i, | ||
- | # new-york | ||
- | # seattle | ||
- | # san-diego | ||
- | dtab = { | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | } | ||
- | model.d = Param(model.i, | ||
- | # Scalar f freight in dollars per case per thousand miles /90/ ; | ||
- | model.f = Param(initialize=90, | ||
- | </ | ||
- | A third, powerful way to initialize a parameter is using a user-defined function.\\ | ||
- | This function will be automatically called by pyomo with any possible (i,j) set. In this case pyomo will actually call c_init() six times in order to initialize the model.c parameter. | ||
- | <code python> | ||
- | # Parameter c(i, | ||
- | # c(i,j) = f * d(i,j) / 1000 ; | ||
- | def c_init(model, | ||
- | return model.f * model.d[i, | ||
- | model.c = Param(model.i, | ||
- | </ | ||
- | |||
- | ==== Variables ==== | ||
- | Similar to parameters, variables are created specifying their domain(s). For variables we can also specify the upper/lower bounds in the constructor.\\ | ||
- | Differently from GAMS, we don't need to define the variable that is on the left hand side of the objective function. | ||
- | <code python> | ||
- | ## Define variables ## | ||
- | # Variables | ||
- | # | ||
- | # | ||
- | # Positive Variable x ; | ||
- | model.x = Var(model.i, | ||
- | </ | ||
- | |||
- | ==== Constrains ==== | ||
- | At this point, it should not be a surprise that constrains are again defined as model objects with the required information passed as parameter in the constructor function. | ||
- | <code python> | ||
- | ## Define contrains ## | ||
- | # supply(i) | ||
- | # supply(i) .. sum (j, x(i,j)) =l= a(i) | ||
- | def supply_rule(model, | ||
- | return sum(model.x[i, | ||
- | model.supply = Constraint(model.i, | ||
- | # demand(j) | ||
- | # demand(j) .. sum(i, x(i,j)) =g= b(j); | ||
- | def demand_rule(model, | ||
- | return sum(model.x[i, | ||
- | model.demand = Constraint(model.j, | ||
- | </ | ||
- | The above code take advantage of [[https:// | ||
- | If we take the supply_rule as example, this is actually called two times by pyomo (once for each of the elements of i). Without List Comprehensions we would have had to write our function using a for loop, like: | ||
- | <code python> | ||
- | def supply_rule(model, | ||
- | supply = 0.0 | ||
- | for j in model.j: | ||
- | supply += model.x[i, | ||
- | return supply <= model.a[i] | ||
- | </ | ||
- | Using List Comprehension is however quicker to code and more readable. | ||
- | |||
- | ==== Objective & solving ==== | ||
- | The definition of the objective is similar to those of the constrains, except that most solvers require a scalar objective function, hence a unique function, and we can specify the sense (direction) of the optimisation. | ||
- | <code python> | ||
- | ## Define Objective and solve ## | ||
- | # cost define objective function | ||
- | # cost .. z =e= sum((i,j), c(i, | ||
- | # Model transport /all/ ; | ||
- | # Solve transport using lp minimizing z ; | ||
- | def objective_rule(model): | ||
- | return sum(model.c[i, | ||
- | model.objective = Objective(rule=objective_rule, | ||
- | </ | ||
- | As we are here looping over two distinct sets, we can see how List Comprehension really simplifies the code. The objective function could have being written without List Comprehension as: | ||
- | <code python> | ||
- | def objective_rule(model): | ||
- | obj = 0.0 | ||
- | for ki in model.i: | ||
- | for kj in model.j: | ||
- | obj += model.c[ki, | ||
- | return obj | ||
- | </ | ||
- | |||
- | ==== Retrieving the output ==== | ||
- | To retrieve the output and do something with it (either to just display it -like we do here-, to plot a graph with [[http:// | ||
- | This function is called by pyomo after the solver has finished. | ||
- | <code python> | ||
- | ## Display of the output ## | ||
- | # Display x.l, x.m ; | ||
- | def pyomo_postprocess(options=None, | ||
- | model.x.display() | ||
- | </ | ||
- | We can print model structure information with '' | ||
- | Results are also by default saved in a results.json file or, if PyYAML is installed in the system, in results.yml. | ||
- | |||
- | ==== Editing and running the script ==== | ||
- | Differently from GAMS you can use whatever editor environment you wish to code a pyomo script. If you don't need debugging features, a simple text editor like Notepad++ (in windows), gedit or kate (in Linux) will suffice. They already have syntax highlight for python.\\ | ||
- | If you want advanced features and debugging capabilities you can use a dedicated Python IDE, like e.g. [[https:// | ||
- | |||
- | You will normally run the script as '' | ||
- | If you want to run the script as '' | ||
- | <code python> | ||
- | # This is an optional code path that allows the script to be run outside of | ||
- | # pyomo command-line. | ||
- | if __name__ == ' | ||
- | #This replicates what the pyomo command-line tools does | ||
- | from coopr.opt import SolverFactory | ||
- | import coopr.environ | ||
- | opt = SolverFactory(" | ||
- | instance = model.create() | ||
- | results = opt.solve(instance) | ||
- | #sends results to stdout | ||
- | results.write() | ||
- | pyomo_postprocess(None, | ||
- | </ | ||
- | |||
- | Finally, if you are very lazy and want to run the script with just ./ | ||
- | <code python> | ||
- | # | ||
- | # -*- coding: utf-8 -*- | ||
- | </ | ||
- | |||
- | ===== Further help ===== | ||
- | Documentation of pyomo is available from [[https:// | ||
- | |||
- | Happy modelling with pyomo ;-) | ||
- | ===== Complete script ===== | ||
- | |||
- | Here is the complete script: | ||
- | |||
- | <code python> | ||
- | # | ||
- | # -*- coding: utf-8 -*- | ||
- | |||
- | """ | ||
- | Basic example of transport model from GAMS model library translated to Pyomo | ||
- | |||
- | To run this you need pyomo and a linear solver installed. | ||
- | When these dependencies are installed you can solve this example in one of | ||
- | this ways (glpk is the default solver): | ||
- | |||
- | ./ | ||
- | python transport.py | ||
- | pyomo transport.py | ||
- | pyomo --solver=glpk transport.py | ||
- | |||
- | To display the results: | ||
- | |||
- | cat results.json | ||
- | cat results.yml (if PyYAML is installed on your system) | ||
- | |||
- | GAMS equivalent code is inserted as single-dash comments. The original GAMS code | ||
- | needs slighly different ordering of the commands and it's available at | ||
- | http:// | ||
- | |||
- | Original problem formulation: | ||
- | Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions. | ||
- | Princeton University Press, Princeton, New Jersey, 1963. | ||
- | GAMS implementation: | ||
- | Rosenthal, R E, Chapter 2: A GAMS Tutorial. In GAMS: A User's Guide. | ||
- | The Scientific Press, Redwood City, California, 1988. | ||
- | Pyomo translation: | ||
- | Antonello Lobianco | ||
- | |||
- | This file is in the Public Domain | ||
- | """ | ||
- | |||
- | # Import | ||
- | from coopr.pyomo import * | ||
- | | ||
- | # Creation of a Concrete Model | ||
- | model = ConcreteModel() | ||
- | |||
- | ## Define sets ## | ||
- | # Sets | ||
- | # | ||
- | # | ||
- | model.i = Set(initialize=[' | ||
- | model.j = Set(initialize=[' | ||
- | |||
- | ## Define parameters ## | ||
- | # | ||
- | # | ||
- | # / | ||
- | # san-diego | ||
- | # | ||
- | # / | ||
- | # chicago | ||
- | # topeka | ||
- | model.a = Param(model.i, | ||
- | model.b = Param(model.j, | ||
- | # Table d(i, | ||
- | # new-york | ||
- | # seattle | ||
- | # san-diego | ||
- | dtab = { | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | } | ||
- | model.d = Param(model.i, | ||
- | # Scalar f freight in dollars per case per thousand miles /90/ ; | ||
- | model.f = Param(initialize=90, | ||
- | # Parameter c(i, | ||
- | # c(i,j) = f * d(i,j) / 1000 ; | ||
- | def c_init(model, | ||
- | return model.f * model.d[i, | ||
- | model.c = Param(model.i, | ||
- | |||
- | ## Define variables ## | ||
- | # Variables | ||
- | # | ||
- | # | ||
- | # Positive Variable x ; | ||
- | model.x = Var(model.i, | ||
- | |||
- | ## Define contrains ## | ||
- | # supply(i) | ||
- | # supply(i) .. sum (j, x(i,j)) =l= a(i) | ||
- | def supply_rule(model, | ||
- | return sum(model.x[i, | ||
- | model.supply = Constraint(model.i, | ||
- | # demand(j) | ||
- | # demand(j) .. sum(i, x(i,j)) =g= b(j); | ||
- | def demand_rule(model, | ||
- | return sum(model.x[i, | ||
- | model.demand = Constraint(model.j, | ||
- | |||
- | ## Define Objective and solve ## | ||
- | # cost define objective function | ||
- | # cost .. z =e= sum((i,j), c(i, | ||
- | # Model transport /all/ ; | ||
- | # Solve transport using lp minimizing z ; | ||
- | def objective_rule(model): | ||
- | return sum(model.c[i, | ||
- | model.objective = Objective(rule=objective_rule, | ||
- | |||
- | ## Display of the output ## | ||
- | # Display x.l, x.m ; | ||
- | def pyomo_postprocess(options=None, | ||
- | model.x.display() | ||
- | |||
- | # This is an optional code path that allows the script to be run outside of | ||
- | # pyomo command-line. | ||
- | if __name__ == ' | ||
- | |||
- | #This replicates what the pyomo command-line tools does | ||
- | from coopr.opt import SolverFactory | ||
- | import coopr.environ | ||
- | opt = SolverFactory(" | ||
- | instance = model.create() | ||
- | results = opt.solve(instance) | ||
- | #sends results to stdout | ||
- | results.write() | ||
- | pyomo_postprocess(None, | ||
- | |||
- | # Expected result: | ||
- | # obj= 153.675 | ||
- | # | ||
- | # | ||
- | # | ||
- | # | ||
- | # | ||
- | # | ||
- | </ | ||
- | |||
- | |||
- | |||
- | ~~DISCUSSION~~ | ||