# A Pyomo tutorial for GAMS users

• 2015.01.12: Updated to run under Pyomo 4. See previous revisions if you still use Pyomo 3.

Pyomo (Python Optimisation Modeling Object) 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 Sandia National Laboratories and appeared in 2008 as an open source project.

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 Pyomo. There may be two reasons for someone to with to use Pyomo instead of GAMS.
The most obvious one, even if often it isn't the key driver, is that GAMS is a commercial software while Pyomo 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, Pyomo would still require the user to buy a licence to use a specific commercial solvers. However Pyomo 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.

Pyomo, at the opposite, is both open source and.. it's python!
You have plenty of development environment to choose from, a clear modern language, the possibility to interface your model with third party libraries.. all of this basically for free.
While there are some reports of pyomo being somehow slower that GAMS it really depends. In my case it is actually much faster, as the IPOPT version that is embedded in GAMS uses the MUMPS linear solver, while on my system I have IPOPT compiled with the much more performing ma27 linear solver. That's part of the flexibility you gain in using pyomo in place of GAMS.

So let's start. We will see how to code the trasnport.gms problem, the one that ship as default example in GAMS1), using Pyomo. 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

Important: Pyomo requires python 2.x. While python 3.x support is work in progress, at the moment only python 2.x is supported.

This isn't true any more with Pyomo 4, where support for Python 3.x has been added..

### Ubuntu

(tested in Ubuntu 14.04 LTS)

• Install the python pre-requisites:
• sudo apt-get install python-yaml, python-pip
• Install pyomo:
• sudo pip install pyomo
• sudo pip install pyomo.extras
• Install solvers:
• linear and MIP solver (glpk): sudo apt-get install glpk36 glpk-utils
• non-linaer solver (ipopt): sudo apt-get install coinor-libipopt1

### Windows and Mac

Please refer to the Coopr installation guide

## Model components

### Creation of the Model

In pyomo everything is an object. The various components of the model (sets, parameters, variables, constrains, objective..) are all attributes of the main model object while being objects themselves.
There are two type of models in pyomo: A ConcreteModel is one where all the data is defined at the model creation. We are going to use this type of model in this tutorial. Pyomo however supports also an AbstractModel, where the model structure is firstly generated and then particular instances of the model are generated with a particular set of data.
The first thing to do in the script is hence to load the pyomo library and to create a new ConcreteModel (we have little imagination here, and we call our model “model”. You can give it whatever name you want2)):

# Import of the pyomo module
from pyomo.environ import *

# Creation of a Concrete Model
model = ConcreteModel()

### 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)
yes, the default GAMS example is named “trasnport”
2)
However, if you give your model an other name, you also need to add a pyomo_create_model(options=None, model_options=None) function that returns your model

## Discussion

, 2014/09/06 20:21

Just a few corrections: Gams allows to use online comments (See the manual), it also allows to define Marcos (functions) that can be reused and gams code can be written in any editor (I use Emacs and the great Gams-mode> Cheers Renger

, 2014/09/12 16:46, 2014/09/12 16:55

Thank you. I were not aware of the macro functionality. I updated the text.

I still defined macros as an “elementary” method to provide a way to structure the code in reusable components as, compared with proper “functions”, macros don't allow to define default parameters, to have multiple outputs passing the arguments by reference, to not say in OO programming the possibility to have function overriding depending on the type/number of the parameters..
That's said, I have also read of the specific features of the GAMS macros, like controlling the expansion of the arguments using the ampersands (&) in the macro body…

Concerning the editors, emacs is one of the few ones that has a syntax highlight mode for GAMS, and I suspect that due to the pretty distinctive GAMS language, the ones that do actually provide syntax highlight work with a small subset of the GAMS syntax.

, 2014/09/09 03:21

Looks useful to lookut for GAMS alternative, can you also include the freely available solvers (good quality) which can be used with python for LP, NLP, MIP, MINLP problems.

, 2014/09/12 16:54, 2014/09/12 16:54

Hello, I have experience with these two solvers:

• GLPK: LP/MIP problems
• IPOPT: LNP/MINLP problems

They both are excellent for their class of problems (IPOPT is comparable with CONOPT). The only real problem of them (especially IPOPT, GLPK is simpler) is that it is not very easy to compile them in Windows, and I didn't try to use a windows-compiled IPOPT with Pyomo (altought I have compiled and used its C++ API directly in windows).

, 2014/09/30 10:59

Thanks for the tutorial.

, 2015/04/28 20:42

Thank you! This helps a lot. Do you have a suggestion on how to build/define more complex tables as for example in Gams's gas transmission line example, as well?

, 2015/04/29 08:23

I think that more complex tables should be handled putting data in a separate file. If you point me to the gams example I'll see if I can get it in Pyomo..

, 2015/04/29 13:55

This would be great. You can find the Gams example I'm talking about here: http://www.gams.com/modlib/libhtml/gastrans.htm

, 2015/04/29 17:18, 2015/04/29 17:30

Hello, I think you have 2 ways.

If you want to keep the specific data in the same python model code you can use the ConcretModel() and define your sets as:

model = ConcreteModel()
m     = model # shorter alias

## Define sets ##
#Sets i Towns / Anderlues, Antwerpen, Arlon, Berneau, Blaregnies, Brugge, Dudzele,
#               Gent, Liege, Loenhout, Mons, Namur, Petange, Peronnes, Sinsin,
#               Voeren, Wanze, Warnand, Zeebrugge, Zomergem /
#     a Arcs  / 1*24 /
#     as(a) active arcs
#     ap(a) passive arcs
#     aij(a,i,i) arc description
#
#
#set nc node data column headers /
#       slo supply lower bound (mill M3 per day)
#       sup supply upper bound (mill M3 per day)
#       plo pressure lower bound (bar)
#       pup pressure upper bound (bar)
#       c   cost ($per MBTU) / # #set ac Arc data column headers / # D diameter (mm) # L length (km) # act type indicator ( 1 = active) / #alias (i,j) m.i = Set(initialize=['Anderlues', 'Antwerpen', 'Arlon', 'Berneau', 'Blaregnies', 'Brugge', 'Dudzele', 'Gent', 'Liege', 'Loenhout', 'Mons', 'Namur', 'Petange', 'Peronnes', 'Sinsin', 'Voeren', 'Wanze', 'Warnand', 'Zeebrugge', 'Zomergem'], name='Towns') m.a = RangeSet(24, name='Arcs') m.aa = Set(within=m.a, name='Active arcs') # 'as' is a reserved keyword m.ap = Set(within=m.a, name='Passive arcs') m.aij = Set(within=m.a*m.i*m.i, name='Arc description') m.nc = Set(initialize=['slo', # supply lower bound (mill M3 per day) 'sup', # supply upper bound (mill M3 per day) 'plo', # pressure lower bound (bar) 'pup', # pressure upper bound (bar) 'c'], # cost ($ per MBTU)
name='node data column headers')
m.ac  = Set(initialize=['D',   # diameter (mm)
'L',   # length (km)
'act'],# type indicator ( 1 = active)
name='Arc data column headers')
m.j   = SetOf(model.i)  # an alias

The problem is that then to create tables you would have to use something like:

NdataTab = {
('Anderlues', 'slo'): 0,
('Anderlues', 'sup'): 1.2,
...
}
m.Ndata = Param(m.i,m.nc, initialize=NdataTab, name='Node data')

(I would suggest the decroise libreoffice extension if you choose this way to unpivot tables that you can then simply copy/paste)

Alternatively you can declare an Abstract model and load data from a file using table AMPL syntax.

I am not aware of a concise method to load a tabled-formatted parameter directly in code. Any how, separate the logic from the data is a cleaner approach.

, 2015/04/29 18:37

Thank you so much!

, 2015/09/14 05:48

Thank you for posting this! I would like to try to compile IPOPT with a linear solver like MA27, but I can't find any good guide to do so.

Do you know of any good guide on how to do so?

Thanks again!

, 2015/09/14 07:30

Hi Pau, the best guide is the one in the IPOPT website (PDF).

I also wrote some instruction for our project, that use IPOPT,here.

In general there is a old way to obtain and compile MA27 and a newer one, depending on the version of the MA27 libraries that you have.

I also suggest you to read the readme that ships with the version of IPOPT you download and, in last istance, subscribe the very useful IPOPT mailing list.

, 2016/01/22 01:16

Hello antonello,

You can tell me please how would this

Abstract model of Pyomo.

thanks

, 2016/10/07 23:13

Hello Antonello,

I'm struggling to solve the transport problem by getting the data from .dat or .xls file, and honestly, I don't know if there is a better way to define the set, param, and var other than typing it one by one especially when dealing with a lot of data. I'm still new to pyomo, I used to do optimization by excel but it limited in terms of variables.

your help is highly appreciated!

Enter your comment. Wiki syntax is allowed:

personal/blog/2014/0904_pyomo_for_gams_users.txt · Last modified: 2018/06/18 15:11 (external edit)