Differences
This shows you the differences between two versions of the page.
Both sides previous revision Previous revision Next revision | Previous revision | ||
home_test_julia [2017/02/07 10:13] antonello created |
home_test_julia [2018/06/18 15:11] (current) |
||
---|---|---|---|
Line 1: | Line 1: | ||
- | ====== A JuMP tutorial for GAMS users ====== | ||
- | [[https:// | + | Run, only once, the following code to install JuMP language and a couple of open source |
- | It is developed by the MIT Operations Research Center and appeared in 2013 as an open source | + | <code julia> |
- | [[http:// | + | Pkg.update() |
- | + | Pkg.add(" | |
- | This mini-tutorial is intended for gams users that want to try JuMP. There may be two reasons for someone to with to use JuMP instead of GAMS.\\ | + | Pkg.add(" |
- | The most obvious one, even if often it isn't the key driver, is that GAMS is a commercial software while JuMP being open-source is free both as freedom and as a free beer.\\ | + | Pkg.add(" |
- | While for GAMS a licence for the underlying solver engine is often included with a particular version of GAMS, JuMP would still require the user to buy a licence to use a specific commercial solvers. However JuMP interfaces with both [[https:// | + | Pkg.add(" |
- | 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, | + | |
- | + | ||
- | JuMP, at the opposite, is both open source and it allow to write the model in a powerful general-purpose language like Julia\\ | + | |
- | 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. | + | |
- | 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 " | + | |
- | 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 ===== | + | |
- | + | ||
- | **Step 1:** | + | |
- | * Option a: Get an account on [[https:// | + | |
- | * Option b: Install Julia for your platform ([[http:// | + | |
- | + | ||
- | **Step 2:** | + | |
- | + | ||
- | Run, only once, the following code to install JuMP language and a couple of open source solvers: | + | |
- | <code lang=" | + | |
- | Pkg.update() | + | |
- | Pkg.add(" | + | |
- | Pkg.add(" | + | |
- | Pkg.add(" | + | |
- | Pkg.add(" | + | |
</ | </ | ||
- | ===== Model components ===== | ||
- | ==== Importing the libraries ==== | + | <code julia> |
- | + | ||
- | You will need to import as a minima the '' | + | |
- | + | ||
- | < | + | |
# 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 and DataFrames modules (the latter one just to import the data from a header-based table, as in the original trasnport example in GAMS | ||
using JuMP, DataFrames | using JuMP, DataFrames | ||
</ | </ | ||
- | ==== Defining the " | ||
- | |||
- | JuMP doesn' | ||
- | 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 " | ||
- | |||
<code julia> | <code julia> | ||
- | ## Define sets ## | + | # Define sets # |
# Sets | # Sets | ||
# | # | ||
Line 63: | Line 24: | ||
markets = [" | markets = [" | ||
</ | </ | ||
- | |||
- | |||
- | ==== Definition of the " | ||
- | |||
- | Capacity of plants and demand of markets are directly defined as dictionaries, | ||
- | |||
- | <code julia> | ||
- | ## Define parameters ## | ||
- | # | ||
- | # | ||
- | # / | ||
- | # san-diego | ||
- | a = Dict( # capacity of plant i in cases | ||
- | " | ||
- | " | ||
- | ) | ||
- | |||
- | # | ||
- | # / | ||
- | # chicago | ||
- | # topeka | ||
- | b = Dict( # demand at market j in cases | ||
- | " | ||
- | " | ||
- | " | ||
- | ) | ||
- | |||
- | # Table d(i, | ||
- | # new-york | ||
- | # seattle | ||
- | # san-diego | ||
- | d_table = wsv""" | ||
- | plants | ||
- | seattle | ||
- | san_diego | ||
- | """ | ||
- | d = Dict( (r[: | ||
- | # Here we are converting the table in a " | ||
- | # r[: | ||
- | # m: the second key | ||
- | # r[Symbol(m)]: | ||
- | |||
- | # Scalar f freight in dollars per case per thousand miles /90/ ; | ||
- | f = 90 # freight in dollars per case per thousand miles | ||
- | |||
- | # Parameter c(i, | ||
- | # c(i,j) = f * d(i,j) / 1000 ; | ||
- | # We first declare an empty dictionary and then we fill it with the values | ||
- | 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] | ||
- | </ | ||
- | 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. | ||
- | |||
- | |||
- | ==== 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.: | ||
- | '' | ||
- | |||
- | <code julia> | ||
- | # Model declaration | ||
- | trmodel = Model() # transport 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. | ||
- | |||
- | <code julia> | ||
- | ## Define variables ## | ||
- | # Variables | ||
- | # | ||
- | # | ||
- | # Positive Variable x ; | ||
- | @variables trmodel begin | ||
- | x[p in plants, m in markets] >= 0 # shipment quantities in cases | ||
- | end | ||
- | </ | ||
- | |||
- | ==== Declaration of the model constraints ==== | ||
- | |||
- | As in GAMS, each constraint can actually be a " | ||
- | |||
- | <code julia> | ||
- | ## Define contrains ## | ||
- | # supply(i) | ||
- | # supply(i) .. sum (j, x(i,j)) =l= a(i) | ||
- | # demand(j) | ||
- | # demand(j) .. sum(i, x(i,j)) =g= b(j); | ||
- | @constraints trmodel begin | ||
- | supply[p in plants], | ||
- | sum(x[p,m] for m in markets) | ||
- | demand[m in markets], | ||
- | sum(x[p,m] for p in plants) | ||
- | end | ||
- | </ | ||
- | |||
- | ==== Declaration of the model objective ==== | ||
- | |||
- | 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. | ||
- | |||
- | <code julia> | ||
- | # Objective | ||
- | @objective trmodel Min begin | ||
- | sum(c[p, | ||
- | end | ||
- | </ | ||
- | |||
- | ==== Human-readable visualisation of the model (optional) ==== | ||
- | |||
- | 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 | ||
- | |||
- | <code julia> | ||
- | print(trmodel) | ||
- | </ | ||
- | |||
- | ==== 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 (": | ||
- | |||
- | <code julia> | ||
- | status = solve(trmodel) | ||
- | </ | ||
- | |||
- | ==== Visualisation of the results ==== | ||
- | While you can do any fancy output you may wish after you retrieve the optimal value of the variables with '' | ||
- | Notice that you can also easily retrieve the dual value associated to the constraint with '' | ||
- | |||
- | <code julia> | ||
- | if status == :Optimal | ||
- | println(" | ||
- | println(getvalue(x)) | ||
- | println(" | ||
- | [println(" | ||
- | println(" | ||
- | [println(" | ||
- | else | ||
- | println(" | ||
- | println(status) | ||
- | end | ||
- | </ | ||
- | |||
- | |||
- | ==== 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.\\ | ||
- | If you want advanced features and debugging capabilities you can use a dedicated Julia IDE, like e.g. [[http:// | ||
- | |||
- | If you are using instead the Julia console, | ||
- | |||
- | ===== Further help ===== | ||
- | Documentation of JuMP is available from [[https:// | ||
- | |||
- | Happy modelling with JuMP ;-) | ||
- | |||
- | ===== Complete script ===== | ||
- | |||
- | Here is the complete script: | ||
- | |||
- | <code Julia> | ||
- | # 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: | ||
- | # 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. | ||
- | # - 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 | ||
- | """ | ||
- | d = Dict( (r[: | ||
- | |||
- | f = 90 # freight in dollars per case per thousand miles | ||
- | |||
- | 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] | ||
- | |||
- | # Model declaration | ||
- | trmodel = Model() # transport model | ||
- | |||
- | # Variables | ||
- | @variables trmodel begin | ||
- | x[p in plants, m in markets] >= 0 # shipment quantities in cases | ||
- | end | ||
- | |||
- | # Constraints | ||
- | @constraints trmodel begin | ||
- | supply[p in plants], | ||
- | sum(x[p,m] for m in markets) | ||
- | demand[m in markets], | ||
- | sum(x[p,m] for p in plants) | ||
- | end | ||
- | |||
- | # Objective | ||
- | @objective trmodel Min begin | ||
- | sum(c[p, | ||
- | end | ||
- | |||
- | print(trmodel) | ||
- | |||
- | status = solve(trmodel) | ||
- | |||
- | if status == :Optimal | ||
- | println(" | ||
- | println(" | ||
- | println(getvalue(x)) | ||
- | println(" | ||
- | [println(" | ||
- | println(" | ||
- | [println(" | ||
- | |||
- | else | ||
- | println(" | ||
- | println(status) | ||
- | end | ||
- | |||
- | # Expected result: | ||
- | # obj= 153.675 | ||
- | # | ||
- | # | ||
- | # | ||
- | # | ||
- | # | ||
- | # | ||
- | </ | ||
- | |||
- | ~~DISCUSSION~~ | ||
- |