This file shows how to solve the minimum cost flow problem for the "Distribution Unlimited Co." prototype example from Hillier and Lieberman, Sections 3.4/9.6.

In [14]:
from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("glpk")

We have a distribution network with five nodes, two of which are sources (factories F1 and F2), one is a transshipment node (distribution center DC), and two are sinks (warehouses W1 and W2). The inflow/outflow at each node is given by $b_i$, the $C_{ij}$ denote the unit cost of transportation from node $i$ to node $j$ and the $U_{ij}$ are upper bounds on the capacity on some of the arcs.

In [15]:
b = {'F1':50,
     'F2':40,
     'DC':0,
     'W1':-30,
     'W2':-60}

C = {('F1','F2'):200,
     ('F1','DC'):400,
     ('F1','W1'):900,
     ('F2','DC'):300,
     ('DC','W2'):100,
     ('W1','W2'):300,
     ('W2','W1'):200}

U = {('F1','F2'):10,
     ('DC','W2'):80}

N = list(b.keys())
A = list(C.keys())
V = list(U.keys())

The decision variables are the $f_{ij}$, the flow on the arc $(i,j)$.

In [16]:
model = ConcreteModel()
model.f = Var(A, within=NonNegativeReals)

The first set of constraints reflects the requirement that inflow plus sources equals outflow at each node.

In [17]:
def flow_rule(model, n):
    InFlow  = sum(model.f[i,j] for (i,j) in A if j==n)
    OutFlow = sum(model.f[i,j] for (i,j) in A if i==n)
    return InFlow + b[n] == OutFlow
    
model.flow = Constraint(N, rule=flow_rule)

The second set of constraints implements the capacity restrictions on the arcs.

In [18]:
def capacity_rule(model, i, j):
    return U[i,j] >= model.f[i,j]

model.capacity = Constraint(V, rule=capacity_rule)

The objective is to minimize the total cost of transportation.

In [19]:
model.cost = Objective(expr = sum(model.f[a]*C[a] for a in A), sense=minimize)

Now we can solve and inspect the results:

In [20]:
results = opt.solve(model)
model.f.get_values()
Out[20]:
{('DC', 'W2'): 80.0,
 ('F1', 'DC'): 40.0,
 ('F1', 'F2'): 0.0,
 ('F1', 'W1'): 10.0,
 ('F2', 'DC'): 40.0,
 ('W1', 'W2'): 0.0,
 ('W2', 'W1'): 20.0}
In [21]:
model.cost.expr()
Out[21]:
49000.0