This is an old revision of the document!


Installation

Step 1:

Step 2:

Run, only once, the following code to install JuMP language and a couple of open source solvers:

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)

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

# 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

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 “_”.

## Define sets ##
#  Sets
#       i   canning plants   / seattle, san-diego /
#       j   markets          / new-york, chicago, topeka / ;
plants  = ["seattle","san_diego"]          # canning plants
markets = ["new_york","chicago","topeka"]  # markets

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.

## 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] 

The above code take advantage of 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.

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
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.

## 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

Declaration of the model constraints

As in GAMS, each constraint can actually be a “family” of constraints:

## 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

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.

# Objective
@objective trmodel Min begin
    sum(c[p,m]*x[p,m] for p in plants, m in markets)
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

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 (“:Optimal” if all went fine)

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 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).

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

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. 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 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 discussion forum.

Happy modelling with JuMP ;-)

Complete script

Here is the complete script:

# 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
home_test_julia.1486458864.txt.gz · Last modified: 2018/06/18 15:10 (external edit)
CC Attribution-Noncommercial-Share Alike 4.0 International
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0