Differences

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

--- home_test_julia [2017/02/07 10:14]
antonello
+++ home_test_julia [2017/02/07 10:18]
antonello
@@ Line 1: / Line 1: @@
-===== 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 julia>
-Pkg.update()                        # To refresh the list of newest packages
-Pkg.add("JuMP")                     # The mathematical optimisation library
+Pkg.update()
-Pkg.add("GLPKMathProgInterface")    # A lineaqr and MIP solver
+Pkg.add("JuMP")           # asdasd
-Pkg.add("Ipopt")                    # A non-linear solver
+Pkg.add("GLPKMathProgInterface")
-Pkg.add("DataFrames")               # A library to deal with dataframes (R-like tabular data)
+Pkg.add("Ipopt")
+Pkg.add("DataFrames")
 </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 27: / 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 41: / 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>
-~~DISCUSSION~~