The Transport Problem

Note: Adapted from https://github.com/Pyomo/PyomoGallery, see LICENSE.BSD

Summary

The goal of the Transport Problem is to select the quantities of an homogeneous good that has several production plants and several punctiform markets as to minimise the transportation costs.

It is the default tutorial for the GAMS language, and 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.

Problem Statement

The Transport Problem can be formulated mathematically as a linear programming problem using the following model.

Sets

$I$ = set of canning plants
$J$ = set of markets

Parameters

$a_i$ = capacity of plant $i$ in cases, $\forall i \in I$
$b_j$ = demand at market $j$ in cases, $\forall j \in J$
$d_{i,j}$ = distance in thousands of miles, $\forall i \in I, \forall j \in J$
$f$ = freight in dollars per case per thousand miles
$c_{i,j}$ = transport cost in thousands of dollars per case

$c_{i,j}$ is obtained exougenously to the optimisation problem as $c_{i,j} = f \cdot d_{i,j}$, $\forall i \in I, \forall j \in J$

Variables

$x_{i,j}$ = shipment quantities in cases
z = total transportation costs in thousands of dollars

Objective

Minimize the total cost of the shipments:
$\min_{x} z = \sum_{i \in I} \sum_{j \in J} c_{i,j} x_{i,j}$

Constraints

Observe supply limit at plant i:
$\sum_{i \in I} x_{i,j} \leq a_{i}$, $\forall i \in I$

Satisfy demand at market j:
$\sum_{j \in J} x_{i,j} \geq b_{j}$, $\forall j \in J$

Non-negative transportation quantities
$x_{i,j} \geq 0$, $\forall i \in I, \forall j \in J$

Pyomo Formulation

Creation of the Model

In pyomo everything is an object. The various components of the model (sets, parameters, variables, constraints, 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 to load the pyomo library and create a new ConcreteModel object. We have little imagination here, and we call our model "model". You can give it whatever name you want. However, if you give your model an other name, you also need to create a model object at the end of your script:


In [1]:
# Import of the pyomo module
from pyomo.environ import *
 
# Creation of a Concrete Model
model = ConcreteModel()

Set Definitions

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:


In [2]:
## 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:


In [3]:
## 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.


In [4]:
#  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.


In [5]:
## 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.


In [6]:
## 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:


In [7]:
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 and 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.


In [8]:
## 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 ;
#
# itertools.product() returns the Cartesian product of two or more iterables

import itertools
def objective_rule(model):
    return sum(model.c[i,j]*model.x[i,j] for i, j in itertools.product(model.i, 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:


In [9]:
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

We use the pyomo_postprocess() function to retrieve the output and do something with it. For example, we could display solution values (see below), plot a graph with matplotlib or save it in a csv file.

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


In [10]:
## 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:


In [11]:
# This emulates what the pyomo command-line tools does
from pyomo.opt import SolverFactory
import pyomo.environ

opt = SolverFactory("glpk")
results = opt.solve(model)

# sends results to stdout
results.write()

print("\nDisplaying Solution\n" + '-'*60)
pyomo_postprocess(None, None, results)


# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 153.675
  Upper bound: 153.675
  Number of objectives: 1
  Number of constraints: 6
  Number of variables: 7
  Number of nonzeros: 13
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
  Error rc: 0
  Time: 0.011461734771728516
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

Displaying Solution
------------------------------------------------------------
x : Shipment quantities in case
    Size=6, Index=x_index, Domain=Reals
    Key                       : Lower : Value : Upper : Fixed : Stale
     ('san-diego', 'chicago') :   0.0 :   0.0 :  None : False : False
    ('san-diego', 'new-york') :   0.0 : 325.0 :  None : False : False
      ('san-diego', 'topeka') :   0.0 : 275.0 :  None : False : False
       ('seattle', 'chicago') :   0.0 : 300.0 :  None : False : False
      ('seattle', 'new-york') :   0.0 :   0.0 :  None : False : False
        ('seattle', 'topeka') :   0.0 :   0.0 :  None : False : False

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:


In [12]:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

Complete script

Here is the complete script:


In [13]:
!cat transport.py


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

# 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 emulates what the pyomo command-line tools does
    from pyomo.opt import SolverFactory
    import pyomo.environ
    opt = SolverFactory("glpk")
    results = opt.solve(model)
    #sends results to stdout
    results.write()
    print("\nDisplaying Solution\n" + '-'*60)
    pyomo_postprocess(None, None, results)

Solutions

Running the model lead to the following output:


In [14]:
!pyomo solve --solver=glpk transport.py


[    0.00] Setting up Pyomo environment
[    0.00] Applying Pyomo preprocessing actions
[    0.00] Creating model
[    0.00] Applying solver
[    0.02] Processing results
    Number of solutions: 1
    Solution Information
      Gap: 0.0
      Status: feasible
      Function Value: 153.67499999999998
    Solver results file: results.json
[    0.02] Applying Pyomo postprocessing actions
x : Shipment quantities in case
    Size=6, Index=x_index, Domain=Reals
    Key                       : Lower : Value : Upper : Fixed : Stale
     ('san-diego', 'chicago') :   0.0 :   0.0 :  None : False : False
    ('san-diego', 'new-york') :   0.0 : 325.0 :  None : False : False
      ('san-diego', 'topeka') :   0.0 : 275.0 :  None : False : False
       ('seattle', 'chicago') :   0.0 : 300.0 :  None : False : False
      ('seattle', 'new-york') :   0.0 :   0.0 :  None : False : False
        ('seattle', 'topeka') :   0.0 :   0.0 :  None : False : False
[    0.02] Pyomo Finished

By default, the optimization results are stored in the file results.json:


In [15]:
!cat results.json


{
    "Problem": [
        {
            "Lower bound": 153.675,
            "Name": "unknown",
            "Number of constraints": 6,
            "Number of nonzeros": 13,
            "Number of objectives": 1,
            "Number of variables": 7,
            "Sense": "minimize",
            "Upper bound": 153.675
        }
    ],
    "Solution": [
        {
            "number of solutions": 1,
            "number of solutions displayed": 1
        },
        {
            "Constraint": "No values",
            "Gap": 0.0,
            "Message": null,
            "Objective": {
                "objective": {
                    "Value": 153.67499999999998
                }
            },
            "Problem": {},
            "Status": "feasible",
            "Variable": {
                "x[san-diego,new-york]": {
                    "Value": 325.0
                },
                "x[san-diego,topeka]": {
                    "Value": 275.0
                },
                "x[seattle,chicago]": {
                    "Value": 300.0
                }
            }
        }
    ],
    "Solver": [
        {
            "Error rc": 0,
            "Statistics": {
                "Branch and bound": {
                    "Number of bounded subproblems": 0,
                    "Number of created subproblems": 0
                }
            },
            "Status": "ok",
            "Termination condition": "optimal",
            "Time": 0.004419803619384766
        }
    ]
}

This solution shows that the minimum transport costs is attained when 300 cases are sent from the Seattle plant to the Chicago market, 50 cases from Seattle to New-York and 275 cases each are sent from San-Diego plant to New-York and Topeka markets.

The total transport costs will be $153,675.

References

  • 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