# Differences

This shows you the differences between two versions of the page.

 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] (current) 2017/02/07 10:26 antonello [Declaration of the model] 2017/02/07 10:24 antonello [Definition of the parameters] 2017/02/07 10:21 antonello [Defining the sets] 2017/02/07 10:20 antonello [Importing the libraries] 2017/02/07 10:20 antonello 2017/02/07 10:12 antonello 2017/02/07 09:05 antonello [Complete script] 2017/02/06 11:09 antonello [Human-readable visualisation of the model (optional)] 2017/02/03 16:44 antonello 2017/02/03 15:32 antonello 2017/02/03 14:45 antonello 2017/02/03 14:29 antonello 2017/02/03 14:09 antonello 2017/02/03 14:03 antonello Next revision Previous revision 2017/02/07 10:26 antonello [Declaration of the model] 2017/02/07 10:24 antonello [Definition of the parameters] 2017/02/07 10:21 antonello [Defining the sets] 2017/02/07 10:20 antonello [Importing the libraries] 2017/02/07 10:20 antonello 2017/02/07 10:12 antonello 2017/02/07 09:05 antonello [Complete script] 2017/02/06 11:09 antonello [Human-readable visualisation of the model (optional)] 2017/02/03 16:44 antonello 2017/02/03 15:32 antonello 2017/02/03 14:45 antonello 2017/02/03 14:29 antonello 2017/02/03 14:09 antonello 2017/02/03 14:03 antonello Line 1: Line 1: - ====== ​a JuMP tutorial for GAMS users ====== + ====== ​A JuMP tutorial for GAMS users ====== [[https://​github.com/​JuliaOpt/​JuMP.jl|JuMP]] (Julia for Mathematical Optimization) is an Algebraic Modelling Language (AML) that allows to write optimisation problems using a concise mathematical formulation,​ acting as interface to the specific solver engine API. For non-linear optimisation problems it allows to keep a high-level approach that doesn'​t require the modeller to compute the Jacobian or the Hessian.\\ [[https://​github.com/​JuliaOpt/​JuMP.jl|JuMP]] (Julia for Mathematical Optimization) is an Algebraic Modelling Language (AML) that allows to write optimisation problems using a concise mathematical formulation,​ acting as interface to the specific solver engine API. For non-linear optimisation problems it allows to keep a high-level approach that doesn'​t require the modeller to compute the Jacobian or the Hessian.\\ 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,​ the package automatically install a local pre-compiled version of the solver, so just adding the package relative to the solver is enough to start writing the model. Even more, for people that doesn'​t care too much about performances,​ there is a service on JuliaBox.com that allows to run Julia/JuMP scripts for free in the browser, without anything to install on the local computer. ​ + That's said, for people that don't need such flexibility,​ the package automatically install a local pre-compiled version of the solver, so just adding the package relative to the solver is enough to start writing the model. Even more, for people that doesn'​t care too much about performances,​ there is a service on [[https://​juliabox.com|JuliaBox.com]] that allows to run Julia/JuMP scripts for free in the browser, without anything to install on the local computer. ​ - 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 "​tra//​sn//​port"​ )), using Pyomo. For a fictions product, there are three canning plants and three markets and the objective of the model is to find the optimal allocation of products between plants and markets that minimises the (transport) costs.\\ + 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 "​tra//​sn//​port"​ )), using JuMP. For a fictions product, there are three canning plants and three markets and the objective of the model is to find the optimal allocation of products between plants and markets that minimises the (transport) costs.\\ 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://​www.gams.com/​mccarl/​trnsport.gms]] 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://​www.gams.com/​mccarl/​trnsport.gms]] ===== Installation ===== ===== Installation ===== - **<​del>​Important: Pyomo requires python 2.x. While python 3.x support is work in progress, at the moment only python 2.x is supported.** + **Step 1:** + * Option a: Get an account on [[https://​juliabox.com|JuliaBox.com]] to run julia/JuMP script without installing anything on the local computer + * Option b: Install Julia for your platform ([[http://​julialang.org/​downloads/​|http://​julialang.org/downloads/​]]) - //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)// + + Pkg.update()                        # To refresh the list of newest packages + Pkg.add("​JuMP"​)                     # The mathematical optimisation library + Pkg.add("​GLPKMathProgInterface"​) ​   # A lineaqr and MIP solver + Pkg.add("​Ipopt"​) ​                   # A non-linear solver + Pkg.add("​DataFrames"​) ​              # A library to deal with dataframes (R like tabular data) + - * **Install the python pre-requisites:​** + ===== Model components ===== - * ''​sudo apt-get install python-yaml,​ python-pip''​ + - * **Install pyomo:** + - * ''​sudo pip install pyomo''​ + - * ''​sudo pip install pyomo.extras''​ + - * **Install solvers:​** + - * //linear and MIP solver (glpk)//: ''​sudo apt-get install glpk36 glpk-utils''​ + - * //​non-linaer solver (ipopt)//: ''​sudo apt-get install coinor-libipopt1''​ + - ==== Windows and Mac ==== + ==== Importing the libraries ​==== - Please refer to the [[https://​software.sandia.gov/​downloads/​pub/​coopr/​CooprInstallGuide.html|Coopr installation guide]] + - ===== Model components ===== + You will need to import as a minima the ''​JuMP''​ module. If you wish to specify a solver engine rather than letting JuMP select a suitable one, you will need to import also the module relative to the solver, e.g. ''​Ipopt''​ or  ''​GLPKMathProgInterface''​ - ==== Creation of the Model ==== + <​code ​ julia> - 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.\\ + # Import ​of the JuMP and DataFrames modules ​(the latter ​one just to import ​the data from a header based table, as in the original trasnport example in GAMS - There are two type of models in pyomo: A ''​ConcreteModel''​ is one where all the data is defined at the model creation. We are going to use this type of model in this tutorial. Pyomo however supports also an ''​AbstractModel'',​ where the model structure is firstly generated and then particular instances of the model are generated with a particular set of data.\\ + using JuMP, DataFrames - 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 "​model"​. You can give it whatever name you want((However,​ if you give your model an other name, you also need to add a ''​pyomo_create_model(options=None,​ model_options=None)''​ function that returns your model))): + - ​ + - # Import of the pyomo module + - from pyomo.environ import * + - + - # Creation of a Concrete Model + - model = ConcreteModel() + ​ - ==== Set definition ​==== + ==== Defining the "​sets" ​==== - Sets are created as attributes object ​of the main model objects and all the information is given as parameter ​in the constructor function. Specifically, we are passing to the constructor ​the initial ​elements of the set and a documentation string to keep track on what our set represents: + - <​code ​python> + JuMP doesn'​t really have a concept ​of sets, but it uses the native containers available ​in the core Julia language\\Variables,​ parameters and constraints can be indexed using these containers.\\ - ## Define sets ## + While many works with position-based lists, I find more readable using dictionaries instead. So the "​sets" ​are represented as lists, but then everything else is a dictionary with the elements of the list as keys.\\ + One note: it seems that Julia/JuMP don't like much the "​-"​ symbol, so I replaced it to "​_"​.\\ + + <​code ​julia> + # Define sets # #  Sets #  Sets #       ​i ​  ​canning plants ​  / seattle, san-diego / #       ​i ​  ​canning plants ​  / seattle, san-diego / #       ​j ​  ​markets ​         / new-york, chicago, topeka / ; #       ​j ​  ​markets ​         / new-york, chicago, topeka / ; - model.i = Set(initialize=['seattle','​san-diego'​], doc='​Canning plans'​) + plants  ​= ["seattle","​san_diego"​]          # canning plants - model.j = Set(initialize=['​new-york'​,'chicago', 'topeka'], doc='​Markets'​) + markets ​= ["​new_york"​,"chicago","topeka"]  # markets ​ - ==== Parameters ​==== + - Parameter objects are created specifying the sets over which they are defined and are initialised with either ​a python ​dictionary ​or a scalar: + ==== Definition of the "​parameters" ​==== - <​code ​python> + - ## Define parameters ​## + Capacity of plants and demand of markets ​are directly ​defined ​as dictionaries,​ while the distance is first read as a DataFrame from a white-space separated table and then it is converted in a "​(plant,​ market) => value" ​dictionary. + + <​code ​julia> + # Define parameters # #   ​Parameters #   ​Parameters #       ​a(i) ​ capacity of plant i in cases #       ​a(i) ​ capacity of plant i in cases #         / ​   seattle ​    350 #         / ​   seattle ​    350 #              san-diego ​  ​600 ​ / #              san-diego ​  ​600 ​ / + a = Dict(              # capacity of plant i in cases + "​seattle" ​  => 350, + "​san_diego"​ => 600, + ) + #       ​b(j) ​ demand at market j in cases #       ​b(j) ​ demand at market j in cases #         / ​   new-york ​   325 #         / ​   new-york ​   325 #              chicago ​    300 #              chicago ​    300 #              topeka ​     275  / ; #              topeka ​     275  / ; - model.a ​= Param(model.i, initialize={'​seattle':​350,'​san-diego':​600},​ doc='​Capacity of plant i in cases') + b = Dict(              # demand at market j in cases - model.b ​= Param(model.j,​ initialize={'​new-york':​325,'chicago':300,'topeka':275}, doc='​Demand at market j in cases') + "​new_york"  ​=> 325, - #  Table d(i,​j) ​ distance in thousands of miles + "chicago" ​  ​=> ​300, + "topeka" ​   => 275, + ) + + # Table d(i,​j) ​ distance in thousands of miles #                    new-york ​      ​chicago ​     topeka #                    new-york ​      ​chicago ​     topeka #      seattle ​         2.5           ​1.7 ​         1.8 #      seattle ​         2.5           ​1.7 ​         1.8 #      san-diego ​       2.5           ​1.8 ​         1.4  ; #      san-diego ​       2.5           ​1.8 ​         1.4  ; - dtab = { + d_table ​= wsv"""​ - ​('​seattle', ​ '​new-york'​) : 2.5, + plants ​    ​new_york ​ chicago ​ topeka - ('​seattle', ​ '​chicago'​) ​ : 1.7, + seattle ​   ​2.5 ​      ​1.7      1.8 - ('​seattle', ​ '​topeka'​) ​  : ​1.8, + san_diego  ​2.5       ​1.8      1.4 - ​('​san-diego','​new-york'​) : 2.5, + """​ - ('​san-diego','​chicago'​) ​ : 1.8, + d = Dict( (r[:plants],m) => r[Symbol(m)] for r in eachrow(d_table), m in markets) - ('​san-diego','​topeka'​) ​  : ​1.4, + # Here we are converting the table in a "​(plant,​ market) => distance"​ dictionary - } + # r[:​plants]: ​  the first key, row field using a fixed header - model.d = Param(model.i, model.j, initialize=dtab, doc='​Distance ​in thousands of miles') + # m:            the second key - #  Scalar f  freight in dollars per case per thousand miles  /90/ ; + # r[Symbol(m)]:​ the value, the row field with a dynamic header - model.f = Param(initialize=90, doc='​Freight ​in dollars per case per thousand miles') + - ​ + # 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 ​in dollars per case per thousand miles - 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. + - + # Parameter c(i,​j) ​ transport cost in thousands of dollars per case ; - #  Parameter c(i,​j) ​ transport cost in thousands of dollars per case ; + #            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 - ​return model.f * model.d[i,j] / 1000 + c = Dict() # transport cost in thousands of dollars per case ; - model.c = Param(model.i,​ model.j, initialize=c_init,​ doc='​Transport cost in thousands of dollar per case') + [ c[p,m] = f * d[p,m] / 1000 for p in plants, m in markets] ​ ​ + The above code take advantage of [[http://​docs.julialang.org/​en/​stable/​manual/​arrays/#​comprehensions|List Comprehensions]],​ a powerful feature of the Julia language that provides a concise way to loop over a list. + 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: + + d = Dict() + for r in eachrow(d_table) + for m in markets + d = (r[:​plants],​m) => r[Symbol(m)] + end + end + ​ + Using List Comprehension is however quicker to code and more readable. - ==== Variables ​==== + - Similar ​to parameters, variables ​are created specifying their domain(s). For variables we can also specify the upper/lower bounds ​in the constructor.\\ + ==== 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 ​to all the subsequent operations, like creation of variables, constraints and objective function.\\ + We can, if we wish, works with several models at the same time.\\ + If we do not specify ​a solver, we let JuML use a suitable solver for the type of problem. Aside to specify the solver, we can also pass it solver-level options, e.g.: + ''​mymodel = Model(solver=IpoptSolver(print_level=0))''​ + + + # Model declaration (transport model) + trmodel = Model() + + + + ==== Declaration of the model variables ==== + + Variables can have multiple-dimensions - that is, being indexed under several indexes -, and bounds ​are given at the same time as their declaration.\\ 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. - <​code ​python> + + <​code ​julia> ## Define variables ## ## Define variables ## #  Variables #  Variables Line 116: Line 151: #       ​z ​      total transportation costs in thousands of dollars ; #       ​z ​      total transportation costs in thousands of dollars ; #  Positive Variable x ; #  Positive Variable x ; - model.x = Var(model.i, model.j, bounds=(0.0,None), doc='​Shipment ​quantities in case') + @variables trmodel begin + ​x[p in plants, m in markets] >= 0 # shipment ​quantities in cases + 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.  ​ + - <​code ​python> + As in GAMS, each constraint can actually ​be a "​family"​ of constraints:​ + + <​code ​julia> ## Define contrains ## ## Define contrains ## # supply(i) ​  ​observe supply limit at plant i # supply(i) ​  ​observe supply limit at plant 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,​ i): - return sum(model.x[i,​j] for j in model.j) <= model.a[i] - model.supply = Constraint(model.i,​ rule=supply_rule,​ doc='​Observe supply limit at plant i') # demand(j) ​  ​satisfy demand at market j ;  ​ # demand(j) ​  ​satisfy demand at market 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 - ​return ​sum(model.x[i,j] for i in model.i) >= model.b[j] + supply[p in plants],   # observe supply limit at plant p - model.demand ​= Constraint(model.j, rule=demand_rule,​ doc='​Satisfy ​demand at market ​j') + sum(x[p,m] for m in markets)  <=  a[p] + demand[m in markets],  # satisfy ​demand at market ​m + sum(x[p,m] for p in plants)  >​= ​ b[m] + end ​ - The above code take advantage ​of [[https://​docs.python.org/​2/​tutorial/​datastructures.html#​list-comprehensions|List Comprehensions]], a powerful feature of the python language that provides ​a concise way to loop over a list. + - 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 ​of the model objective ==== - <​code ​python> + - def supply_rule(model,​ i): + Contrary to constraints and variables, the objective is always ​a unique function. Note that it is at this point that we specify ​the direction ​of the optimisation. - ​supply = 0.0 + - for j in model.j: + <​code ​julia> - ​supply += model.x[i,j] + # Objective - ​return supply <= model.a[i] + @objective trmodel Min begin + ​sum(c[p,m]*x[p,m] for p in plants, m in markets) + end ​ - Using List Comprehension is however quicker to code and more readable. - ==== Objective & solving ==== + ==== Human-readable visualisation ​of the model (optional) ==== - 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. + - ​ + 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 + <​code ​julia> - #  cost ..        z  =e=  sum((i,j), c(i,​j)*x(i,​j)) ; + print(trmodel) - #  Model transport /all/ ; + - #  Solve transport using lp minimizing z ; + - def objective_rule(model): + - return sum(model.c[i,​j]*model.x[i,​j] for i in model.i for j in model.j) + - model.objective ​= Objective(rule=objective_rule,​ sense=minimize, doc='​Define objective function'​) + - ​ + - 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,​kj]*model.x[ki,​kj] + - return obj + ​ - ==== Retrieving ​the output ​==== + ==== Resolution of the model ==== - To retrieve ​the output ​and do something with it (either ​to just display it -like we do here-, to plot a graph with [[http://​matplotlib.org|matplotlib]] or to save it in a csv file) we use the ''​pyomo_postprocess()''​ function.\\ + - This function is called by pyomo after the solver has finished. ​ + It is at this point that the solver is called ​and the model is passed ​to the solver engine for its solution. The return value is the status of the optimisation ​(":​Optimal"​ if all went fine) - <​code ​python> + - ## Display of the output ## + <​code ​julia> - # Display x.l, x.m ; + status ​= solve(trmodel) - def pyomo_postprocess(options=None, instance=None,​ results=None):​ + - model.x.display() + ​ - We can print model structure information with ''​model.pprint()''​ ("​pprint"​ stand for "​pretty print"​).\\ - 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 ​==== + ==== Visualisation of the results ​==== - 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.\\ + While you can do any fancy output ​you may wish after you retrieve the optimal value of the variables with ''​getvalue(var_name)''​, you can just ''​println(getvalue(x))''​ to get a basic output.\\ - If you want advanced features and debugging capabilities ​you can use a dedicated Python IDE, like e.g. [[https://​code.google.com/​p/​spyderlib/​|Spyder]]. + Notice that you can also easily retrieve the dual value associated to the constraint with ''​getdual(constraint_name)''​. - You will normally run the script as ''​pyomo solve --solver=glpk transport.py''​. You can output solver specific output adding the option ''​--stream-output''​.\\ + <​code ​julia> - If you want to run the script as ''​python transport.py''​ add the following lines at the end:\\ + if status ​== :Optimal - <​code ​python> + ​println("​Objective value: ", getobjectivevalue(trmodel)) - # This is an optional code path that allows the script to be run outside of + ​println(getvalue(x)) - # pyomo command-line. ​ For example: ​ python transport.py + ​println("​Shadow prices of supply:"​) - if __name__ ​== '​__main__'​: + ​[println("​$p ​=$(getdual(supply[p]))"​) ​for p in plants] - ​#This replicates what the pyomo command-line tools does + ​println("​Shadow prices of demand:"​) - ​from pyomo.opt import SolverFactory + ​[println("​$m ​=$(getdual(demand[m]))") for m in markets] - ​import pyomo.environ + else - ​opt = SolverFactory("glpk") + ​println("Model didn't solved"​) - ​instance = model.create() + ​println(status) - ​results ​= opt.solve(instance) + end - #​sends results to stdout + - ​results.write() + - ​pyomo_postprocess(None, instance, results) + ​ - Finally, if you are very lazy and want to run the script ​with just ./transport.py (and you are in Linux) ​add the following lines at the top: + - + ==== Editing ​and running ​the script ​==== - #!/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) ​will suffice. They already have syntax highlight for Julia.\\ - # -*- 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, ​ you can run the script as ''​julia transport.jl''​. ===== Further help ===== ===== Further help ===== - Documentation of pyomo is available from [[https://software.sandia.gov/trac/coopr/wiki/​Documentation|this page]]. However if you want to do serious things with pyomo, it is most likely that you will have to either look at the source code or consult the [[https://groups.google.com/forum/#!forum/coopr-forum|mailing list]]. + 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 ​forum]]. + + Happy modelling with JuMP ;-) - Happy modelling with pyomo ;-) ===== Complete script ===== ===== Complete script ===== Here is the complete script: ​ Here is the complete script: ​ - <​code ​python> + <​code ​Julia> - #!/​usr/​bin/​env python + # 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:​ Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions. + # 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:​ Antonello Lobianco + + using JuMP, DataFrames + + # Sets + plants ​ = ["​seattle","​san_diego"​] ​         # canning plants + markets = ["​new_york","​chicago","​topeka"​] ​ # markets + + # Parameters + a = Dict(              # capacity of plant i in cases + "​seattle" ​  => 350, + "​san_diego"​ => 600, + ) + b = Dict(              # demand at market j in cases + "​new_york" ​ => 325, + "​chicago" ​  => 300, + "​topeka" ​   => 275, + ) + + #  distance in thousands of miles + d_table = wsv"""​ + plants ​    ​new_york ​ chicago ​ topeka + seattle ​   2.5       ​1.7 ​     1.8 + san_diego ​ 2.5       ​1.8 ​     1.4 """​ """​ - 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] - ./​transport.py ​(Linux only) + - python transport.py + # Model declaration - pyomo solve transport.py + trmodel ​= Model() # transport model - pyomo solve --solver=glpk transport.py + - + # Variables - To display the results: + @variables trmodel begin - + x[p in plants, m in markets] >= 0 # shipment quantities ​in cases - 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://​www.gams.com/​mccarl/​trnsport.gms + - + - 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 + - #       ​i ​  ​canning ​plants ​  / seattle, san-diego / + - #       ​j ​  ​markets ​         / new-york, chicago, topeka / ; + - model.i = Set(initialize=['​seattle','​san-diego'​], doc='​Canning plans') + - model.j = Set(initialize=['​new-york','​chicago',​ '​topeka'​],​ doc='​Markets'​) + - + - ## Define parameters ## + - #   ​Parameters + - #       a(i)  ​capacity of plant i in cases + - #         / ​   seattle ​    350 + - #              san-diego ​  ​600 ​ / + - #       b(j)  ​demand at market j in cases + - #         / ​   new-york ​   325 + - #              chicago ​    300 + - #              topeka ​     275  / ; + - model.a = Param(model.i, initialize={'​seattle':​350,'​san-diego':​600},​ doc='​Capacity of plant i in cases') + - model.b = Param(model.j,​ initialize={'​new-york':​325,'​chicago':​300,'​topeka':​275},​ doc='​Demand at market j in cases'​) + - #  Table d(i,​j) ​ distance in thousands of miles + - #                    new-york ​      ​chicago ​     topeka + - #      seattle ​         2.5           ​1.7 ​         1.8 + - #      san-diego ​       2.5           ​1.8 ​         1.4  ; + - dtab = { + - ('​seattle', ​ '​new-york'​) : 2.5, + - ('​seattle', ​ '​chicago'​) ​ : 1.7, + - ('​seattle', ​ '​topeka'​) ​  : 1.8, + - ('​san-diego','​new-york'​) : 2.5, + - ('​san-diego','​chicago'​) ​ : 1.8, + - ('​san-diego','​topeka'​) ​  : 1.4, + - } + - model.d = Param(model.i,​ model.j, initialize=dtab,​ doc='​Distance in thousands of miles'​) + - #  ​Scalar f  ​freight in dollars per case per thousand miles  /90/ ; + - model.f ​= Param(initialize=90,​ doc='​Freight in dollars per case per thousand miles') + - #  ​Parameter c(i,​j)  ​transport cost in thousands of dollars per case ; + - #            ​c(i,j) = f * d(i,j) / 1000 ; + - def c_init(model,​ i, j): + - return model.f * model.d[i,j] / 1000 + - model.c = Param(model.i,​ model.j, initialize=c_init,​ doc='​Transport cost in thousands of dollar per case'​) + - + - ## Define variables ## + - #  Variables + - #       x(i,j)  shipment quantities ​in cases + - #       ​z ​      total transportation costs in thousands of dollars ; + - #  ​Positive Variable x ; + - model.x ​= Var(model.i, model.j, bounds=(0.0,​None), doc='​Shipment quantities in case'​) + - + - ## Define contrains ## + - # supply(i) ​  ​observe supply limit at plant i + - # supply(i) .. sum (j, x(i,j)) =l= a(i) + - def supply_rule(model,​ i): + - ​return sum(model.x[i,j] for j in model.j) <= model.a[i] + - model.supply ​= Constraint(model.i,​ rule=supply_rule,​ doc='​Observe supply limit at plant i') + - # demand(j) ​  ​satisfy demand at market j ; + - # demand(j) .. sum(i, x(i,j)) =g= b(j); + - def demand_rule(model,​ j): + - return sum(model.x[i,​j] for i in model.i) >= model.b[j]  ​ + - model.demand = Constraint(model.j,​ rule=demand_rule,​ doc='​Satisfy demand at market j') + - ## Define Objective and solve ## + # Constraints - #  ​cost ​       define objective function + @constraints trmodel begin - #  cost ..        z  =e=  ​sum((i,j), c(i,j)*x(i,j)) ; + supply[p in plants], ​  # observe supply limit at plant p - #  Model transport /all/ ; + sum(x[p,m] for m in markets)  <​= ​ a[p] - #  Solve transport using lp minimizing z ; + ​demand[m in markets], ​ # ​satisfy demand at market m - def objective_rule(model):​ + sum(x[p,m] for p in plants)  >​= ​ b[m] - return ​sum(model.c[i,​j]*model.x[i,j] for i in model.i for j in model.j) + end - model.objective = Objective(rule=objective_rule,​ sense=minimize,​ doc='​Define objective function'​) + - + # 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 - ​model.x.display() + - + print(trmodel) - # This is an optional code path that allows the script to be run outside of + - # pyomo command-line. ​ For example: ​ python transport.py + status = solve(trmodel) - if __name__ ​== '​__main__'​: + - + if status ​== :Optimal - ​#This replicates what the pyomo command-line tools does + ​println("​Objective value: ", getobjectivevalue(trmodel)) - ​from pyomo.opt import SolverFactory + ​println("​Shipped quantities: ") - ​import pyomo.environ + ​println(getvalue(x)) - ​opt = SolverFactory("glpk") + ​println("Shadow prices of supply:") - ​instance ​= model.create() + ​[println("​$p ​=$(getdual(supply[p]))") for p in plants] - ​results = opt.solve(instance) + ​println("​Shadow prices of demand:"​) - ​#sends results to stdout + ​[println("​$m =$(getdual(demand[m]))"​) for m in markets] - ​results.write() + - ​pyomo_postprocess(None, instance, results) + else + ​println("Model didn't solved"​) + ​println(status) + end # Expected result: # Expected result: Line 353: Line 326: #​['​san-diego','​topeka'​] ​  = 275 #​['​san-diego','​topeka'​] ​  = 275 ​ - - - - ~~DISCUSSION~~ - - - ~~DISCUSSION~~ ~~DISCUSSION~~
personal/blog/2017/0203_jump_for_gams_users.1486126996.txt.gz · Last modified: 2018/06/18 15:10 (external edit)