Differences

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

--- home_test_julia [2017/02/07 09:27]
antonello
+++ home_test_julia [2018/06/18 15:11] (current)
@@ Line 1: / Line 1: @@
-====== 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.\\
+Run, only once, the following code to install JuMP language and a couple of open source solvers:
-It is developed by the MIT Operations Research Center and appeared in 2013 as an open source package of the relatively new Julia programming language.
+<code julia>
-[[http://www.gams.com/|GAMS]] (The General Algebraic Modeling System) does more or less the same things and appeared in the '70s as a project of the World Bank. GAMS is hence a very mature project (maybe //too// mature) with a lot of followers in the economic domain, where it is used mainly to solve equilibria problems.
+Pkg.update()
+Pkg.add("JuMP")           # asdasd
-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("GLPKMathProgInterface")
-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("Ipopt")
-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://www.gnu.org/software/glpk/|GLPK]] (for linear and mixed-integer programming) and [[https://projects.coin-or.org/Ipopt|IPOPT]] (for non-linear optimisation) open-source solvers, both of which are top on their classes, leaving the necessity to acquire a licence for a commercial solver to niche cases.\\
+Pkg.add("DataFrames")
-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, it's still the most common way to code in GAMS.\\
-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, 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 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]]
-===== Installation =====
-**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/]])
-**Step 2:**
-Run, only once, the following code to install JuMP language and a couple of open source solvers:
-<code lang="julia">
-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)
 </code>
-===== Model components =====
-==== 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''
 <code julia>
@@ Line 49: / Line 16: @@
 </code>
-==== Defining the "sets" ====
-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.\\
-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 ##
+# Define sets #
 #  Sets
 #       i   canning plants   / seattle, san-diego /
@@ Line 63: / Line 24: @@
 markets = ["new_york","chicago","topeka"]  # markets
 </code>
-==== Definition of the "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
-#       a(i)  capacity of plant i in cases
-#         /    seattle     350
-#              san-diego   600  /
-a = Dict(              # 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  / ;
-b = Dict(              # demand at market j in cases
-  "new_york"  => 325,
-  "chicago"   => 300,
-  "topeka"    => 275,
-)
-# 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  ;
-d_table = wsv"""
-plants     new_york  chicago  topeka
-seattle    2.5       1.7      1.8
-san_diego  2.5       1.8      1.4
-"""
-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
-# r[:plants]:   the first key, row field using a fixed header
-# m:            the second key
-# r[Symbol(m)]: the value, the row field with a dynamic header
-# Scalar f  freight in dollars per case per thousand miles  /90/ ;
-f = 90 # 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 ;
-# 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]
-</code>
-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:
-<code julia>
-d = Dict()
-for r in eachrow(d_table)
-  for m in markets
-    d = (r[:plants],m) => r[Symbol(m)]
-  end
-end
-</code>
-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.:
-''mymodel = Model(solver=IpoptSolver(print_level=0))''
-<code julia>
-# Model declaration
-trmodel = Model() # transport model
-</code>
-==== 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
-#       x(i,j)  shipment quantities in cases
-#       z       total transportation costs in thousands of dollars ;
-#  Positive Variable x ;
-@variables trmodel begin
-    x[p in plants, m in markets] >= 0 # shipment quantities in cases
-end
-</code>
-==== Declaration of the model constraints ====
-As in GAMS, each constraint can actually be a "family" of constraints:
-<code julia>
-## Define contrains ##
-# supply(i)   observe supply limit at plant i
-# supply(i) .. sum (j, x(i,j)) =l= a(i)
-# demand(j)   satisfy demand at market j ;
-# demand(j) .. sum(i, x(i,j)) =g= b(j);
-@constraints trmodel begin
-    supply[p in plants],   # observe supply limit at plant p
-        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
-</code>
-==== 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,m]*x[p,m] for p in plants, m in markets)
-end
-</code>
-==== 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)
-</code>
-==== 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)
-<code julia>
-status = solve(trmodel)
-</code>
-==== Visualisation of the results ====
-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.\\
-Notice that you can also easily retrieve the dual value associated to the constraint with ''getdual(constraint_name)''.
-<code julia>
-if status == :Optimal
-    println("Objective value: ", getobjectivevalue(trmodel))
-    println(getvalue(x))
-    println("Shadow prices of supply:")
-    [println("$p = $(getdual(supply[p]))") for p in plants]
-    println("Shadow prices of demand:")
-    [println("$m = $(getdual(demand[m]))") for m in markets]
-else
-    println("Model didn't solved")
-    println(status)
-end
-</code>
-==== 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://junolab.org/|Juno]].
-If you are using instead the Julia console,  you can run the script as ''julia transport.jl''.
-===== 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]].
-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: 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
-"""
-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
-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],   # observe supply limit at plant p
-        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
-# Objective
-@objective trmodel Min begin
-    sum(c[p,m]*x[p,m] for p in plants, m in markets)
-end
-print(trmodel)
-status = solve(trmodel)
-if status == :Optimal
-    println("Objective value: ", getobjectivevalue(trmodel))
-    println("Shipped quantities: ")
-    println(getvalue(x))
-    println("Shadow prices of supply:")
-    [println("$p = $(getdual(supply[p]))") for p in plants]
-    println("Shadow prices of demand:")
-    [println("$m = $(getdual(demand[m]))") for m in markets]
-else
-    println("Model didn't solved")
-    println(status)
-end
-# Expected result:
-# obj= 153.675
-#['seattle','new-york']   = 50
-#['seattle','chicago']    = 300
-#['seattle','topeka']     = 0
-#['san-diego','new-york'] = 275
-#['san-diego','chicago']  = 0
-#['san-diego','topeka']   = 275
-</code>