*sn*port”

**This is an old revision of the document!**

### Table of Contents

# A JuMP tutorial for GAMS users

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.

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.

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.

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.

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.

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 GLPK (for linear and mixed-integer programming) and 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.

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 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^{1)}, 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 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/)

**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 <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> # Sets plants = ["seattle","san_diego"] # canning plants markets = ["new_york","chicago","topeka"] # markets

### Defining the "parameters"

Capacity of plants and demand of markets are directly defined as dictionaries, while the distance is first read from a white-space separated table and then it is converted in a “(plant, market) ⇒ value” dictionary.

# 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) # 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 f = 90 # freight in dollars per case per thousand miles # 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]

# 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(getvalue(x))

else

println("Model didn't solved") println(status)

end </code>

### Set definition

Sets are created as attributes object of the main model objects and all the information is given as parameter in the constructor function. Specifically, we are passing to the constructor the initial elements of the set and a documentation string to keep track on what our set represents:

## Define sets ## # Sets # i canning plants / seattle, san-diego / # j markets / new-york, chicago, topeka / ; model.i = Set(initialize=['seattle','san-diego'], doc='Canning plans') model.j = Set(initialize=['new-york','chicago', 'topeka'], doc='Markets')

### Parameters

Parameter objects are created specifying the sets over which they are defined and are initialised with either a python dictionary or a scalar:

## Define parameters ## # Parameters # a(i) 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 / ; model.a = Param(model.i, initialize={'seattle':350,'san-diego':600}, doc='Capacity of plant i in cases') model.b = Param(model.j, initialize={'new-york':325,'chicago':300,'topeka':275}, doc='Demand at market j in cases') # 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 ; dtab = { ('seattle', 'new-york') : 2.5, ('seattle', 'chicago') : 1.7, ('seattle', 'topeka') : 1.8, ('san-diego','new-york') : 2.5, ('san-diego','chicago') : 1.8, ('san-diego','topeka') : 1.4, } model.d = Param(model.i, model.j, initialize=dtab, doc='Distance in thousands of miles') # Scalar f freight in dollars per case per thousand miles /90/ ; model.f = Param(initialize=90, doc='Freight in dollars per case per thousand miles')

A third, powerful way to initialize a parameter is using a user-defined function.

This function will be automatically called by pyomo with any possible (i,j) set. In this case pyomo will actually call c_init() six times in order to initialize the model.c parameter.

# Parameter c(i,j) transport cost in thousands of dollars per case ; # c(i,j) = f * d(i,j) / 1000 ; def c_init(model, i, j): return model.f * model.d[i,j] / 1000 model.c = Param(model.i, model.j, initialize=c_init, doc='Transport cost in thousands of dollar per case')

### Variables

Similar to parameters, variables are created specifying their domain(s). For variables we can also specify the upper/lower bounds in the constructor.

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 ; model.x = Var(model.i, model.j, bounds=(0.0,None), doc='Shipment quantities in case')

### Constrains

At this point, it should not be a surprise that constrains are again defined as model objects with the required information passed as parameter in the constructor function.

## Define contrains ## # supply(i) observe supply limit at plant i # supply(i) .. sum (j, x(i,j)) =l= a(i) def supply_rule(model, i): return sum(model.x[i,j] for j in model.j) <= model.a[i] model.supply = Constraint(model.i, rule=supply_rule, doc='Observe supply limit at plant i') # demand(j) satisfy demand at market j ; # demand(j) .. sum(i, x(i,j)) =g= b(j); def demand_rule(model, j): return sum(model.x[i,j] for i in model.i) >= model.b[j] model.demand = Constraint(model.j, rule=demand_rule, doc='Satisfy demand at market j')

The above code take advantage of List Comprehensions, a powerful feature of the python language that provides a concise way to loop over a list. If we take the supply_rule as example, this is actually called two times by pyomo (once for each of the elements of i). Without List Comprehensions we would have had to write our function using a for loop, like:

def supply_rule(model, i): supply = 0.0 for j in model.j: supply += model.x[i,j] return supply <= model.a[i]

Using List Comprehension is however quicker to code and more readable.

### Objective & solving

The definition of the objective is similar to those of the constrains, except that most solvers require a scalar objective function, hence a unique function, and we can specify the sense (direction) of the optimisation.

## Define Objective and solve ## # cost define objective function # cost .. z =e= sum((i,j), c(i,j)*x(i,j)) ; # Model transport /all/ ; # Solve transport using lp minimizing z ; def objective_rule(model): return sum(model.c[i,j]*model.x[i,j] for i in model.i for j in model.j) model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')

As we are here looping over two distinct sets, we can see how List Comprehension really simplifies the code. The objective function could have being written without List Comprehension as:

def objective_rule(model): obj = 0.0 for ki in model.i: for kj in model.j: obj += model.c[ki,kj]*model.x[ki,kj] return obj

### Retrieving the output

To retrieve the output and do something with it (either to just display it -like we do here-, to plot a graph with matplotlib or to save it in a csv file) we use the `pyomo_postprocess()`

function.

This function is called by pyomo after the solver has finished.

## Display of the output ## # Display x.l, x.m ; def pyomo_postprocess(options=None, instance=None, results=None): model.x.display()

We can print model structure information with `model.pprint()`

(“pprint” stand for “pretty print”).

Results are also by default saved in a results.json file or, if PyYAML is installed in the system, in results.yml.

### Editing and running the script

Differently from GAMS you can use whatever editor environment you wish to code a pyomo 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 python.

If you want advanced features and debugging capabilities you can use a dedicated Python IDE, like e.g. Spyder.

You will normally run the script as `pyomo solve –solver=glpk transport.py`

. You can output solver specific output adding the option `–stream-output`

.

If you want to run the script as `python transport.py`

add the following lines at the end:

# This is an optional code path that allows the script to be run outside of # pyomo command-line. For example: python transport.py if __name__ == '__main__': #This replicates what the pyomo command-line tools does from pyomo.opt import SolverFactory import pyomo.environ opt = SolverFactory("glpk") instance = model.create() results = opt.solve(instance) #sends results to stdout results.write() pyomo_postprocess(None, instance, results)

Finally, if you are very lazy and want to run the script with just ./transport.py (and you are in Linux) add the following lines at the top:

#!/usr/bin/env python # -*- coding: utf-8 -*-

## Further help

Documentation of pyomo is available from this page. However if you want to do serious things with pyomo, it is most likely that you will have to either look at the source code or consult the mailing list.

Happy modelling with pyomo

## Complete script

Here is the complete script:

#!/usr/bin/env python # -*- coding: utf-8 -*- """ Basic example of transport model from GAMS model library translated to Pyomo To run this you need pyomo and a linear solver installed. When these dependencies are installed you can solve this example in one of this ways (glpk is the default solver): ./transport.py (Linux only) python transport.py pyomo solve transport.py pyomo solve --solver=glpk transport.py To display the results: cat results.json cat results.yml (if PyYAML is installed on your system) GAMS equivalent code is inserted as single-dash comments. The original GAMS code needs slighly different ordering of the commands and it's available at http://www.gams.com/mccarl/trnsport.gms Original problem formulation: Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions. 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. Pyomo translation: Antonello Lobianco This file is in the Public Domain """ # Import from pyomo.environ import * # Creation of a Concrete Model model = ConcreteModel() ## Define sets ## # Sets # i canning plants / seattle, san-diego / # j markets / new-york, chicago, topeka / ; model.i = Set(initialize=['seattle','san-diego'], doc='Canning plans') model.j = Set(initialize=['new-york','chicago', 'topeka'], doc='Markets') ## Define parameters ## # Parameters # a(i) 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 / ; model.a = Param(model.i, initialize={'seattle':350,'san-diego':600}, doc='Capacity of plant i in cases') model.b = Param(model.j, initialize={'new-york':325,'chicago':300,'topeka':275}, doc='Demand at market j in cases') # 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 ; dtab = { ('seattle', 'new-york') : 2.5, ('seattle', 'chicago') : 1.7, ('seattle', 'topeka') : 1.8, ('san-diego','new-york') : 2.5, ('san-diego','chicago') : 1.8, ('san-diego','topeka') : 1.4, } model.d = Param(model.i, model.j, initialize=dtab, doc='Distance in thousands of miles') # Scalar f freight in dollars per case per thousand miles /90/ ; model.f = Param(initialize=90, doc='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 ; def c_init(model, i, j): return model.f * model.d[i,j] / 1000 model.c = Param(model.i, model.j, initialize=c_init, doc='Transport cost in thousands of dollar per case') ## Define variables ## # Variables # x(i,j) shipment quantities in cases # z total transportation costs in thousands of dollars ; # Positive Variable x ; model.x = Var(model.i, model.j, bounds=(0.0,None), doc='Shipment quantities in case') ## Define contrains ## # supply(i) observe supply limit at plant i # supply(i) .. sum (j, x(i,j)) =l= a(i) def supply_rule(model, i): return sum(model.x[i,j] for j in model.j) <= model.a[i] model.supply = Constraint(model.i, rule=supply_rule, doc='Observe supply limit at plant i') # demand(j) satisfy demand at market j ; # demand(j) .. sum(i, x(i,j)) =g= b(j); def demand_rule(model, j): return sum(model.x[i,j] for i in model.i) >= model.b[j] model.demand = Constraint(model.j, rule=demand_rule, doc='Satisfy demand at market j') ## Define Objective and solve ## # cost define objective function # cost .. z =e= sum((i,j), c(i,j)*x(i,j)) ; # Model transport /all/ ; # Solve transport using lp minimizing z ; def objective_rule(model): return sum(model.c[i,j]*model.x[i,j] for i in model.i for j in model.j) model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function') ## Display of the output ## # Display x.l, x.m ; def pyomo_postprocess(options=None, instance=None, results=None): model.x.display() # This is an optional code path that allows the script to be run outside of # pyomo command-line. For example: python transport.py if __name__ == '__main__': #This replicates what the pyomo command-line tools does from pyomo.opt import SolverFactory import pyomo.environ opt = SolverFactory("glpk") instance = model.create() results = opt.solve(instance) #sends results to stdout results.write() pyomo_postprocess(None, instance, results) # 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

^{1)}

## Discussion

You can use a normal `if`/`else` statement or use the more concise ? ternary operator:

aVariable = CONDITION ? ExprIfTrue : ExprIfFalse

Obviously if the condition is whitin a model to optimise, the model will need to be declared nonlinear.