Note

You can download this example as a Jupyter notebook or start it in interactive mode.

The Transport Problem#

Summary#

The goal of the Transport Problem is to select the quantities of a 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. The Pyomo version of the tutorial can be found at https://nbviewer.org/github/Pyomo/PyomoGallery/blob/master/transport/transport.ipynb.

For comparison, the variable names have been kept the same, though they are not pythonic. The equivalent GAMS code is included at the top of each cell as comments, while the Pyomo-equivalent is stated as text comments below.

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_{j \in J} x_{i,j} \leq a_{i}\), \(\forall i \in I\)

Satisfy demand at market j: \(\sum_{i \in I} 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\)

Linopy Formulation#

Creation of the Model#

In linopy, modeling is based on a central model object. While this object contains variables, constraints, expressions, and the objective function as attributes, it does not currently contain fixed data like sets and parameters. Thus, they have been defined separately. Since linopy is built heavily on xarray’s data structure, using xarray for that is convenient.

The first thing to do in the tutorial is to load the linopy and related libraries and create a new model object. We have little imagination here, and we call our model m. You can give it whatever name you want. However, if you give your model an other name, you also need to change all references to it throughout this tutorial.

[1]:

# Import of linopy and related modules
import linopy
import xarray as xr

# Creation of a Model
m = linopy.Model()

Equivalent in Pyomo

from pyomo.environ import *

model = ConcreteModel()

Set Definitions#

Sets can be created as xarray.DataArrays, where they require data in any form that numpy.ndarray could work with but can be provided additional data such as the name of the array, the names of its dimensions, or coordinates. Coordinates enable accessing data based on labels and aligning data. Since we want the sets i and j to function as coordinates later on, we only specify the name of their dimension here. Due to xarray’s flexibility, we can also define our coordinates as usual dictionaries, where the keys correspond to the names of the dimensions:

[2]:

## Define sets ##
#  Sets
#       i   canning plants   / seattle, san-diego /
#       j   markets          / new-york, chicago, topeka / ;

i = {"Canning Plants": ["seattle", "san-diego"]}
j = {"Markets": ["new-york", "chicago", "topeka"]}

Equivalent in Pyomo

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 as xarray.DataArrays and do specify over which sets they are created using coordinates. Since we set up i and j properly before, e.g. the supply a will automatically contain the dimension "Canning Plants" with the labels "seattle" and "san-diego" and can be accessed accordingly. This also works for multi-dimensional data such as the distance in thousands of miles d, where the dimension names are inferred again from the coordinates’ dimension names. Note, though, that xarray would prefer explicit statements of the form coords = [(dimension_name, data), ... ].

For the scalar parameter f, we simply define an integer variable.

[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  / ;

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

#  Scalar f  freight in dollars per case per thousand miles  /90/ ;

a = xr.DataArray([350, 600], coords=i, name="capacity of plant i in cases")
b = xr.DataArray([325, 300, 275], coords=j, name="demand at market j in cases")

d = xr.DataArray([[2.5, 1.7, 1.8], [2.5, 1.8, 1.4]], coords=i | j, name="distance in thousands of miles")

f = 90  # Freight in dollars per case per thousand miles

# Access data using e.g.:
# a.loc[{"Canning Plants":"seattle"}]
# d.loc[{"Canning Plants":"seattle", "Markets":"new-york"}]

Equivalent in Pyomo

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

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

model.f = Param(initialize=90, doc='Freight in dollars per case per thousand miles')

Working with xarrays enables us to get the resulting transport cost in thousands of dollars per case much simpler than in Pyomo:

[4]:

#  Parameter c(i,j)  transport cost in thousands of dollars per case ;
#            c(i,j) = f * d(i,j) / 1000 ;

c = d * f / 1000
c.name = "transport cost in thousands of dollars per case"

Equivalent in Pyomo

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#

Variables are created as model attributes. They get can get lower or upper bounds, a name, and the coordinates (and inferred dimensions) for which they are defined.

Differently from GAMS, we do not need to define the variable that is on the left-hand side of the objective function.

[5]:

## Define variables ##
#  Variables
#       x(i,j)  shipment quantities in cases
#       z       total transportation costs in thousands of dollars ;
#  Positive Variable x ;

x = m.add_variables(lower=0.0, coords=c.coords, name="Shipment quantities in cases")

# Inspect the variable by simply printing it:
# print(x)

Equivalent in Pyomo

model.x = Var(model.i, model.j, bounds=(0.0,None), doc='Shipment quantities in case')

Constraints#

Constraints are also defined as model attributes. Before assigning them to the model, though, it is very useful to check that they take the form we desire. Here, we make use of xarray’s .sum() functionality, for which we can specify the dimensions that should be summed over:

[6]:

x.sum(dims="Markets") <= a

[6]:

Constraint (unassigned) (Canning Plants: 2):
--------------------------------------------
[seattle]: +1 Shipment quantities in cases[seattle, new-york] + 1 Shipment quantities in cases[seattle, chicago] + 1 Shipment quantities in cases[seattle, topeka]         ≤ 350.0
[san-diego]: +1 Shipment quantities in cases[san-diego, new-york] + 1 Shipment quantities in cases[san-diego, chicago] + 1 Shipment quantities in cases[san-diego, topeka] ≤ 600.0

The output nicely confirms that this is indeed the constraint we want to implement, so we add it to the model:

[7]:

## Define contraints ##
# 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);

con = x.sum(dims="Markets") <= a
con1 = m.add_constraints(con, name="Observe supply limit at plant i")

con = x.sum(dims="Canning Plants") >= b
con2 = m.add_constraints(con, name="Satisfy demand at market j")

Equivalent in Pyomo

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

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

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
Linopy translation: Fridolin Glatter

The Transport Problem

Contents

The Transport Problem#

Summary#

Problem Statement#

Sets#

Parameters#

Variables#

Objective#

Constraints#

Linopy Formulation#

Creation of the Model#

Set Definitions#

Parameters#

Variables#

Constraints#

Objective and Solving#

References#