Stylesheet style.css not found, please contact the developer of "arctic" template.
[[personal:blog:2017:0203_jump_for_gams_users]]
You are here: home » personal » blog » 2017 » 0203_jump_for_gams_users
Show page
Recent ChangesSitemap

Differences

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

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
personal:blog:2017:0203_jump_for_gams_users [2017/02/03 15:44] – antonellopersonal:blog:2017:0203_jump_for_gams_users [2025/05/02 09:41] (current) – external edit 127.0.0.1
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 [[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.   +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.   
  
  
Line 30: Line 30:
  
 Run, only once, the following code to install JuMP language and a couple of open source solvers: Run, only once, the following code to install JuMP language and a couple of open source solvers:
-<code lang="julia"> +<code julia> 
-Pkg.update()                        # To refresh the list of newest packages +using Pkg               # Load the package manager 
-Pkg.add("JuMP")                     # The mathematical optimisation library +Pkg.update()            # To refresh the list of newest packages 
-Pkg.add("GLPKMathProgInterface")    # A lineaqr and MIP solver +Pkg.add("CSV")          # A library to work with Comma Separated Values 
-Pkg.add("Ipopt")                    # A non-linear solver +Pkg.add("DataFrames")   # A library to deal with dataframes (R like tabular data) 
-Pkg.add("DataFrames")               # A library to deal with dataframes (R-like tabular data)+Pkg.add("JuMP")         # The mathematical optimisation library 
 +Pkg.add("GLPK")         # A linear and MIP solver 
 +Pkg.add("Ipopt")        # A non-linear solver (not needed in this example)
 </code> </code>
  
Line 42: Line 44:
 ==== Importing the libraries ==== ==== Importing the libraries ====
  
-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''+You will need to import as a minima the ''JuMP'' module and a suitable solver. In this case the problem is linear, so we can use ''GLPK'' (''HiGHS'' is another popular alternative). If the problem would have been non-linear, you could have used the ''Ipopt'' solver/package
  
-<code  lang="julia"> +<code  julia> 
-# 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  +# Import of the JuMP, GLPK, CSV and DataFrames modules (the latter twos just to import the data from a header based table, as in the original trasnport example in GAMS  
-using JuMP, DataFrames+using CSV, DataFrames, GLPK, JuMP
 </code> </code>
  
Line 56: Line 58:
    
 <code julia> <code julia>
-## Define sets ##+# Define sets #
 #  Sets #  Sets
 #       i   canning plants   / seattle, san-diego / #       i   canning plants   / seattle, san-diego /
Line 70: Line 72:
  
 <code julia> <code julia>
-## Define parameters ##+# Define parameters #
 #   Parameters #   Parameters
 #       a(i)  capacity of plant i in cases #       a(i)  capacity of plant i in cases
Line 94: Line 96:
 #      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  ;
-d_table = wsv"""+d_table = CSV.read(IOBuffer("""
 plants     new_york  chicago  topeka plants     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
-"""+"""), DataFrame, delim=" ", ignorerepeated=true,copycols=true)
 d = Dict( (r[:plants],m) => r[Symbol(m)] for r in eachrow(d_table), m in markets) 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 "(plant, market) => distance" dictionary # Here we are converting the table in a "(plant, market) => distance" dictionary
Line 130: Line 132:
  
 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.\\ 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.\\ +The solver engine to use is given as argument of the ''Model()'' call.\\ 
-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.: +We could pass solver-specific options with the ''set_optimizer_attribute'' function, e.g.: 
-''mymodel = Model(solver=IpoptSolver(print_level=0))''+''set_optimizer_attribute(trmodel, "msg_lev", GLPK.GLP_MSG_ON)''
  
 <code julia> <code julia>
-# Model declaration +# Model declaration (transport model) 
-trmodel = Model() # transport model+trmodel = Model(GLPK.Optimizer) 
 </code> </code>
  
Line 190: Line 192:
 <code julia> <code julia>
 print(trmodel) print(trmodel)
-</code+</code>
  
 ==== Resolution of the model ==== ==== Resolution of the model ====
  
-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)+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 (''MOI.OPTIMAL'' if all went fine)
  
 <code julia> <code julia>
-status = solve(trmodel)+optimize!(trmodel) 
 +status = termination_status(trmodel)
 </code> </code>
  
Line 205: Line 208:
  
 <code julia> <code julia>
-if status == :Optimal +if status == MOI.OPTIMAL 
-    println("Objective value: ", getobjectivevalue(trmodel)) +    println("Objective value: ", objective_value(trmodel)) 
-    println(getvalue(x))+    println("Shipped quantities: ") 
 +    println(value.(x))
     println("Shadow prices of supply:")     println("Shadow prices of supply:")
-    [println("$p = $(getdual(supply[p]))") for p in plants]+    [println("$p = $(dual(supply[p]))") for p in plants]
     println("Shadow prices of demand:")     println("Shadow prices of demand:")
-    [println("$m = $(getdual(demand[m]))") for m in markets]+    [println("$m = $(dual(demand[m]))") for m in markets] 
 + 
 else else
     println("Model didn't solved")     println("Model didn't solved")
Line 221: Line 226:
 ==== Editing and running the script ==== ==== Editing and running the script ====
 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.\\ 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.\\
-If you want advanced features and debugging capabilities you can use a dedicated Julia IDE, like e.g. [[http://junolab.org/|Juno]].+If you want advanced features and debugging capabilities you can use a dedicated Julia IDE, like the [[https://www.julia-vscode.org/|Julia extension for VSCode]].
  
-If you are using instead the Julia console,  you can run the script as ''julia transport.jl''.+If you are using instead the Julia terminal,  you can run the script as ''julia transport.jl''.
  
 ===== Further help ===== ===== Further help =====
-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]].+Documentation of JuMP is available from [[https://jump.dev/|this page]], and community-based support is available on [[https://discourse.julialang.org/c/domain/opt|the Discourse forum]].
  
 Happy modelling with JuMP ;-) Happy modelling with JuMP ;-)
Line 234: Line 239:
 Here is the complete script:  Here is the complete script: 
  
-<code julia> +<code Julia> 
-#= +# Transport example
-Transposition in JuMP of the basic transport model used in the GAMS tutorial +
- +
-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+
  
 +# Transposition in JuMP of the basic transport model used in the GAMS tutorial
 +# 
 +# 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 CSV, DataFrames, GLPK, JuMP
 + 
 # Sets # Sets
 plants  = ["seattle","san_diego"]          # canning plants plants  = ["seattle","san_diego"]          # canning plants
 markets = ["new_york","chicago","topeka"]  # markets markets = ["new_york","chicago","topeka"]  # markets
 + 
 # Parameters # Parameters
 a = Dict(              # capacity of plant i in cases a = Dict(              # capacity of plant i in cases
Line 265: Line 270:
   "topeka"    => 275,   "topeka"    => 275,
 ) )
 + 
 #  distance in thousands of miles #  distance in thousands of miles
-d_table = wsv"""+d_table = CSV.read(IOBuffer("""
 plants     new_york  chicago  topeka plants     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
-"""+"""), DataFrame, delim=" ", ignorerepeated=true,copycols=true)
 d = Dict( (r[:plants],m) => r[Symbol(m)] for r in eachrow(d_table), m in markets) d = Dict( (r[:plants],m) => r[Symbol(m)] for r in eachrow(d_table), m in markets)
 + 
 f = 90 # freight in dollars per case per thousand miles f = 90 # freight in dollars per case per thousand miles
 + 
 c = Dict() # transport cost in thousands of dollars per case ; 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] [ c[p,m] = f * d[p,m] / 1000 for p in plants, m in markets]
 + 
 # Model declaration # Model declaration
-trmodel = Model() # transport model +trmodel = Model(GLPK.Optimizer) # transport model 
 + 
 # Variables # Variables
 @variables trmodel begin @variables trmodel begin
     x[p in plants, m in markets] >= 0 # shipment quantities in cases     x[p in plants, m in markets] >= 0 # shipment quantities in cases
 end end
 + 
 # Constraints # Constraints
 @constraints trmodel begin @constraints trmodel begin
Line 294: Line 299:
         sum(x[p,m] for p in plants)  >=  b[m]         sum(x[p,m] for p in plants)  >=  b[m]
 end end
 + 
 # Objective # Objective
 @objective trmodel Min begin @objective trmodel Min begin
     sum(c[p,m]*x[p,m] for p in plants, m in markets)     sum(c[p,m]*x[p,m] for p in plants, m in markets)
 end end
 + 
 print(trmodel) print(trmodel)
 + 
 +optimize!(trmodel)
 +status = termination_status(trmodel)
  
-status = solve(trmodel) +if status == MOI.OPTIMAL 
- +    println("Objective value: ", objective_value(trmodel))
-if status == :Optimal +
-    println("Objective value: ", getobjectivevalue(trmodel))+
     println("Shipped quantities: ")     println("Shipped quantities: ")
-    println(getvalue(x))+    println(value.(x))
     println("Shadow prices of supply:")     println("Shadow prices of supply:")
-    [println("$p = $(getdual(supply[p]))") for p in plants]+    [println("$p = $(dual(supply[p]))") for p in plants]
     println("Shadow prices of demand:")     println("Shadow prices of demand:")
-    [println("$m = $(getdual(demand[m]))") for m in markets] +    [println("$m = $(dual(demand[m]))") for m in markets] 
 + 
 else else
     println("Model didn't solved")     println("Model didn't solved")
     println(status)     println(status)
 end end
- +
 # Expected result: # Expected result:
 # obj= 153.675 # obj= 153.675
personal/blog/2017/0203_jump_for_gams_users.1486136650.txt.gz · Last modified: (external edit)
Show pageOld revisions
Media ManagerBack to top
CC Attribution-Noncommercial-Share Alike 4.0 International
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0