Differences
This shows you the differences between two versions of the page.
Next revision | Previous revision Next revision Both sides next revision | ||
personal:blog:2017:0203_jump_for_gams_users [2017/02/03 14:03] antonello |
personal:blog:2017:0203_jump_for_gams_users [2018/06/18 15:11] 127.0.0.1 external edit |
||
---|---|---|---|
Line 1: | Line 1: | ||
- | ====== | + | ====== |
[[https:// | [[https:// | ||
Line 15: | Line 15: | ||
You have plenty of development environment to choose from (e.g. Jupiter, Juno), a clear modern language, the possibility to interface your model with third party libraries.. all of this basically for free.\\ | You have plenty of development environment to choose from (e.g. Jupiter, Juno), a clear modern language, the possibility to interface your model with third party libraries.. all of this basically for free.\\ | ||
It is also, at least for my user case, much faster than GAMS. Aside the preparation of the model to pass to the solver, where it is roughly equivalent, in the solver execution I can benefit of having on my system a version of IPOPT compiled with the much more performing ma27 linear solver, while for GAMS I would have to rely on the embedded version that is compiled with the MUMPS linear solver. That's part of the flexibility you gain in using JuMP in place of GAMS. | It is also, at least for my user case, much faster than GAMS. Aside the preparation of the model to pass to the solver, where it is roughly equivalent, in the solver execution I can benefit of having on my system a version of IPOPT compiled with the much more performing ma27 linear solver, while for GAMS I would have to rely on the embedded version that is compiled with the MUMPS linear solver. That's part of the flexibility you gain in using JuMP in place of GAMS. | ||
- | That's said, for people that don't need such flexibility, | + | That's said, for people that don't need such flexibility, |
- | 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 is named " | + | 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 is 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:// | 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 ===== | ===== Installation ===== | ||
- | **< | + | **Step 1:** |
+ | * Option a: Get an account on [[https:// | ||
+ | * Option b: Install Julia for your platform ([[http:// | ||
- | //This isn't true any more with Pyomo 4, where support for Python 3.x has been added..// | + | **Step 2:** |
- | ==== Ubuntu ==== | + | Run, only once, the following code to install JuMP language and a couple of open source solvers: |
- | //(tested in Ubuntu 14.04 LTS)// | + | <code julia> |
+ | Pkg.update() # To refresh the list of newest packages | ||
+ | Pkg.add(" | ||
+ | Pkg.add(" | ||
+ | Pkg.add(" | ||
+ | Pkg.add(" | ||
+ | </code> | ||
- | * **Install the python pre-requisites: | + | ===== Model components ===== |
- | * '' | + | |
- | * **Install pyomo:** | + | |
- | * '' | + | |
- | * '' | + | |
- | * **Install solvers: | + | |
- | * //linear and MIP solver (glpk)//: '' | + | |
- | * // | + | |
- | ==== Windows and Mac ==== | + | ==== Importing the libraries |
- | Please refer to the [[https:// | + | |
- | ===== Model components ===== | + | You will need to import as a minima the '' |
- | ==== Creation of the Model ==== | + | < |
- | In pyomo everything is an object. The various components | + | # Import |
- | There are two type of models in pyomo: A '' | + | using JuMP, DataFrames |
- | The first thing to do in the script is hence to load the pyomo library and to create | + | |
- | <code python> | + | |
- | # Import of the pyomo module | + | |
- | from pyomo.environ import * | + | |
- | + | ||
- | # Creation of a Concrete Model | + | |
- | model = ConcreteModel() | + | |
</ | </ | ||
- | ==== Set definition | + | ==== Defining the " |
- | Sets are created as attributes object | + | |
- | < | + | JuMP doesn' |
- | ## Define sets ## | + | While many works with position-based lists, I find more readable using dictionaries instead. So the " |
+ | One note: it seems that Julia/JuMP don't like much the " | ||
+ | |||
+ | < | ||
+ | # Define sets # | ||
# Sets | # Sets | ||
# | # | ||
# | # | ||
- | model.i = Set(initialize=['seattle',' | + | plants |
- | model.j = Set(initialize=[' | + | markets |
</ | </ | ||
- | ==== Parameters | + | |
- | Parameter objects are created specifying the sets over which they are defined and are initialised with either | + | ==== Definition of the " |
- | < | + | |
- | ## Define parameters | + | Capacity of plants and demand of markets |
+ | |||
+ | < | ||
+ | # Define parameters # | ||
# | # | ||
# | # | ||
# / | # / | ||
# san-diego | # san-diego | ||
+ | a = Dict( # capacity of plant i in cases | ||
+ | " | ||
+ | " | ||
+ | ) | ||
+ | |||
# | # | ||
# / | # / | ||
# chicago | # chicago | ||
# topeka | # topeka | ||
- | model.a | + | b = Dict( # demand at market j in cases |
- | model.b | + | " |
- | # Table d(i, | + | "chicago" |
+ | "topeka" | ||
+ | ) | ||
+ | |||
+ | # Table d(i, | ||
# new-york | # new-york | ||
# seattle | # seattle | ||
# san-diego | # san-diego | ||
- | dtab = { | + | d_table |
- | | + | plants |
- | (' | + | seattle |
- | (' | + | san_diego |
- | | + | """ |
- | (' | + | d = Dict( (r[:plants],m) => r[Symbol(m)] for r in eachrow(d_table), m in markets) |
- | (' | + | # Here we are converting the table in a " |
- | } | + | # r[: |
- | model.d = Param(model.i, model.j, initialize=dtab, doc=' | + | # m: the second key |
- | # Scalar f freight in dollars per case per thousand miles /90/ ; | + | # r[Symbol(m)]: |
- | model.f = Param(initialize=90, doc=' | + | |
- | </ | + | # Scalar f freight in dollars per case per thousand miles /90/ ; |
- | A third, powerful way to initialize a parameter is using a user-defined function.\\ | + | f = 90 # freight |
- | 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, |
- | # Parameter c(i, | + | |
# c(i,j) = f * d(i,j) / 1000 ; | # c(i,j) = f * d(i,j) / 1000 ; | ||
- | def c_init(model, i, j): | + | # We first declare an empty dictionary and then we fill it with the values |
- | | + | c = Dict() # transport cost in thousands of dollars per case ; |
- | model.c = Param(model.i, | + | [ c[p,m] = f * d[p,m] / 1000 for p in plants, m in markets] |
</ | </ | ||
+ | The above code take advantage of [[http:// | ||
+ | If we take the creation of the d dictionary as example, without List Comprehensions we would have had to write a nested for loop like: | ||
+ | <code julia> | ||
+ | d = Dict() | ||
+ | for r in eachrow(d_table) | ||
+ | for m in markets | ||
+ | d = (r[: | ||
+ | end | ||
+ | end | ||
+ | </ | ||
+ | Using List Comprehension is however quicker to code and more readable. | ||
- | ==== Variables | + | |
- | Similar | + | ==== Declaration of the model ==== |
+ | |||
+ | Here we declare a JuML optimisation model and we give it a name. This name will be then passed as first argument | ||
+ | We can, if we wish, works with several models at the same time.\\ | ||
+ | If we do not specify | ||
+ | '' | ||
+ | |||
+ | <code julia> | ||
+ | # Model declaration (transport model) | ||
+ | trmodel = Model() | ||
+ | |||
+ | </code> | ||
+ | |||
+ | ==== Declaration of the model variables ==== | ||
+ | |||
+ | Variables can have multiple-dimensions - that is, being indexed under several indexes -, and bounds | ||
Differently from GAMS, we don't need to define the variable that is on the left hand side of the objective function. | Differently from GAMS, we don't need to define the variable that is on the left hand side of the objective function. | ||
- | < | + | |
+ | < | ||
## Define variables ## | ## Define variables ## | ||
# Variables | # Variables | ||
Line 116: | Line 151: | ||
# | # | ||
# Positive Variable x ; | # Positive Variable x ; | ||
- | model.x = Var(model.i, model.j, bounds=(0.0,None), doc=' | + | @variables trmodel begin |
+ | | ||
+ | end | ||
</ | </ | ||
- | ==== Constrains | + | ==== Declaration of the model constraints |
- | 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. | + | |
- | < | + | As in GAMS, each constraint can actually |
+ | |||
+ | < | ||
## Define contrains ## | ## Define contrains ## | ||
# supply(i) | # supply(i) | ||
# supply(i) .. sum (j, x(i,j)) =l= a(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) | ||
# demand(j) .. sum(i, x(i,j)) =g= b(j); | # demand(j) .. sum(i, x(i,j)) =g= b(j); | ||
- | def demand_rule(model, j): | + | @constraints trmodel begin |
- | | + | supply[p in plants], # observe supply limit at plant p |
- | model.demand | + | sum(x[p,m] for m in markets) <= a[p] |
+ | demand[m in markets], # satisfy | ||
+ | sum(x[p,m] for p in plants) > | ||
+ | end | ||
</ | </ | ||
- | The above code take advantage | + | |
- | 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: | + | ==== Declaration |
- | < | + | |
- | def supply_rule(model, | + | Contrary to constraints and variables, the objective is always |
- | | + | |
- | for j in model.j: | + | < |
- | | + | # Objective |
- | | + | @objective trmodel Min begin |
+ | | ||
+ | end | ||
</ | </ | ||
- | Using List Comprehension is however quicker to code and more readable. | ||
- | ==== Objective & solving ==== | + | ==== Human-readable visualisation |
- | The definition | + | |
- | <code python> | + | If we wish we can get the optimisation model printed in a human-readable fashion, so we can expect all is like it should be |
- | ## Define Objective and solve ## | + | |
- | # cost define objective function | + | < |
- | # cost .. z =e= sum((i,j), c(i, | + | print(trmodel) |
- | # Model transport /all/ ; | + | |
- | # Solve transport using lp minimizing z ; | + | |
- | def objective_rule(model): | + | |
- | return sum(model.c[i, | + | |
- | model.objective | + | |
- | </ | + | |
- | 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: | + | |
- | < | + | |
- | def objective_rule(model): | + | |
- | obj = 0.0 | + | |
- | for ki in model.i: | + | |
- | for kj in model.j: | + | |
- | obj += model.c[ki, | + | |
- | return obj | + | |
</ | </ | ||
- | ==== Retrieving | + | ==== Resolution of the model ==== |
- | To retrieve | + | |
- | This function is called by pyomo after the solver has finished. | + | It is at this point that the solver is called |
- | < | + | |
- | ## Display of the output ## | + | < |
- | # Display x.l, x.m ; | + | status |
- | def pyomo_postprocess(options=None, instance=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 | + | ==== Visualisation of the results |
- | Differently from GAMS you can use whatever editor environment | + | While you can do any fancy output |
- | If you want advanced features and debugging capabilities | + | Notice that you can also easily retrieve the dual value associated to the constraint with '' |
- | You will normally run the script as '' | + | < |
- | If you want to run the script as '' | + | if status |
- | < | + | |
- | # This is an optional code path that allows the script to be run outside of | + | |
- | # pyomo command-line. | + | |
- | if __name__ | + | |
- | | + | |
- | | + | |
- | | + | else |
- | | + | |
- | | + | |
- | | + | end |
- | # | + | |
- | | + | |
- | | + | |
</ | </ | ||
- | Finally, if you are very lazy and want to run the script | + | |
- | <code python> | + | ==== Editing |
- | #!/usr/bin/env python | + | Differently from GAMS you can use whatever editor environment you wish to code a JuMP script. If you don't need debugging features, a simple text editor like Notepad++ (in windows), gedit or kate (in Linux) |
- | # -*- coding: utf-8 -*- | + | If you want advanced features and debugging capabilities you can use a dedicated Julia IDE, like e.g. [[http://junolab.org/|Juno]]. |
- | </ | + | |
+ | If you are using instead the Julia console, | ||
===== Further help ===== | ===== Further help ===== | ||
- | Documentation of pyomo is available from [[https://software.sandia.gov/trac/coopr/wiki/ | + | Documentation of JuMP is available from [[https://jump.readthedocs.io/en/latest/|this page]]. However if you want to do serious things with juMP, it is most likely that you will have to either look at the source code or consult the [[https://discourse.julialang.org/c/domain/opt|discussion |
+ | |||
+ | Happy modelling with JuMP ;-) | ||
- | Happy modelling with pyomo ;-) | ||
===== Complete script ===== | ===== Complete script ===== | ||
Here is the complete script: | Here is the complete script: | ||
- | < | + | < |
- | #!/ | + | # Transposition in JuMP of the basic transport model used in the GAMS tutorial |
- | # -*- coding: utf-8 -*- | + | # |
- | + | # This problem finds a least cost shipping schedule that meets | |
+ | # requirements at markets and supplies at factories. | ||
+ | # | ||
+ | # - Original formulation: | ||
+ | # Princeton University Press, Princeton, New Jersey, 1963. | ||
+ | # - Gams implementation: This formulation is described in detail in: | ||
+ | # Rosenthal, R E, Chapter 2: A GAMS Tutorial. In GAMS: A User's Guide. | ||
+ | # The Scientific Press, Redwood City, California, 1988. | ||
+ | # - JuMP implementation: | ||
+ | |||
+ | using JuMP, DataFrames | ||
+ | |||
+ | # Sets | ||
+ | plants | ||
+ | markets = [" | ||
+ | |||
+ | # Parameters | ||
+ | a = Dict( # capacity of plant i in cases | ||
+ | " | ||
+ | " | ||
+ | ) | ||
+ | b = Dict( # demand at market j in cases | ||
+ | " | ||
+ | " | ||
+ | " | ||
+ | ) | ||
+ | |||
+ | # distance in thousands of miles | ||
+ | d_table = wsv""" | ||
+ | plants | ||
+ | seattle | ||
+ | san_diego | ||
""" | """ | ||
- | Basic example of transport model from GAMS model library translated to Pyomo | + | d = Dict( (r[:plants],m) => r[Symbol(m)] for r in eachrow(d_table), m in markets) |
- | + | ||
- | To run this you need pyomo and a linear solver installed. | + | f = 90 # freight in dollars per case per thousand miles |
- | When these dependencies are installed you can solve this example in one of | + | |
- | this ways (glpk is the default solver): | + | c = Dict() # transport cost in thousands of dollars per case ; |
- | + | [ c[p,m] = f * d[p,m] / 1000 for p in plants, m in markets] | |
- | ./ | + | |
- | python transport.py | + | # Model declaration |
- | pyomo solve transport.py | + | trmodel |
- | pyomo solve --solver=glpk transport.py | + | |
- | + | # Variables | |
- | To display the results: | + | @variables trmodel begin |
- | + | x[p in plants, m in markets] >= 0 # shipment quantities | |
- | cat results.json | + | end |
- | 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 pyomo.environ import * | + | |
- | + | ||
- | # Creation of a Concrete Model | + | |
- | model = ConcreteModel() | + | |
- | + | ||
- | ## Define sets ## | + | |
- | # Sets | + | |
- | # | + | |
- | # | + | |
- | model.i = Set(initialize=[' | + | |
- | model.j = Set(initialize=[' | + | |
- | + | ||
- | ## Define parameters ## | + | |
- | # | + | |
- | # a(i) | + | |
- | # / | + | |
- | # san-diego | + | |
- | # b(j) | + | |
- | # / | + | |
- | # chicago | + | |
- | # topeka | + | |
- | model.a = Param(model.i, initialize={' | + | |
- | model.b = Param(model.j, | + | |
- | # Table d(i, | + | |
- | # new-york | + | |
- | # seattle | + | |
- | # san-diego | + | |
- | dtab = { | + | |
- | (' | + | |
- | (' | + | |
- | (' | + | |
- | (' | + | |
- | (' | + | |
- | (' | + | |
- | } | + | |
- | model.d = Param(model.i, | + | |
- | # | + | |
- | model.f | + | |
- | # | + | |
- | # | + | |
- | def c_init(model, | + | |
- | return model.f * model.d[i,j] / 1000 | + | |
- | model.c = Param(model.i, | + | |
- | + | ||
- | ## Define variables ## | + | |
- | # Variables | + | |
- | # x(i,j) shipment quantities | + | |
- | # | + | |
- | # | + | |
- | model.x | + | |
- | + | ||
- | ## Define contrains ## | + | |
- | # supply(i) | + | |
- | # supply(i) .. sum (j, x(i,j)) =l= a(i) | + | |
- | def supply_rule(model, | + | |
- | | + | |
- | model.supply | + | |
- | # 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 ## | + | # Constraints |
- | # | + | @constraints trmodel begin |
- | # cost .. z =e= | + | supply[p in plants], |
- | # Model transport /all/ ; | + | sum(x[p,m] for m in markets) < |
- | # Solve transport using lp minimizing z ; | + | |
- | def objective_rule(model): | + | sum(x[p,m] for p in plants) > |
- | return | + | end |
- | model.objective = Objective(rule=objective_rule, | + | |
- | + | # Objective | |
- | ## Display of the output ## | + | @objective trmodel Min begin |
- | # Display x.l, x.m ; | + | sum(c[p,m]*x[p,m] for p in plants, m in markets) |
- | def pyomo_postprocess(options=None, instance=None, results=None): | + | end |
- | | + | |
- | + | print(trmodel) | |
- | # This is an optional code path that allows the script to be run outside of | + | |
- | # pyomo command-line. | + | status = solve(trmodel) |
- | if __name__ | + | |
- | + | if status | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | else |
+ | | ||
+ | | ||
+ | end | ||
# Expected result: | # Expected result: | ||
Line 353: | Line 326: | ||
# | # | ||
</ | </ | ||
- | |||
- | |||
- | |||
- | ~~DISCUSSION~~ | ||
- | |||
- | |||
- | |||
~~DISCUSSION~~ | ~~DISCUSSION~~ | ||