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.

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.

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 $x_{ij}$, the flow on the arc $(i,j)$.

model = ConcreteModel()
model.x = Var(A, within=NonNegativeReals)

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

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

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

def capacity_rule(model, i, j):
    return model.x[i,j] <= U[i,j]

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

The objective is to minimize the total cost of transportation.

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

Now we can solve and inspect the results:

model.dual = Suffix(direction=Suffix.IMPORT)

results = opt.solve(model)
model.x.get_values()

{('F1', 'F2'): 0.0,
 ('F1', 'DC'): 40.0,
 ('F1', 'W1'): 10.0,
 ('F2', 'DC'): 40.0,
 ('DC', 'W2'): 80.0,
 ('W1', 'W2'): 0.0,
 ('W2', 'W1'): 20.0}

Not ethe integer solution property: All $x_{ij}$ are integer. Even more, since all $b_i$ and $u_{ij}$ came in units of $10$, all $x_{ij}$ come in units of $10$ as well.

model.cost.expr()

49000.0

Let us also look at the shadow prices:

for i in N:
    print('Shadow price for', i, ':', model.dual[model.flow[i]])

Shadow price for F1 : 700.0
Shadow price for F2 : 600.0
Shadow price for DC : 300.0
Shadow price for W1 : -200.0
Shadow price for W2 : 0.0

So, for example: Suppose we could increase the supply by one unit in either F1 or F2, and the demand in either W1 or W2. Which combination should we choose to increase the cost the least? Increasing the b by one unit at either F1 or F2 will lead to an increase of 700 or 600. Increasing the demand at W1 or W2 means adding -1 to the respective b's. So this comes at a price of $(-1)*(-200) = 200$ for W1 and 0 for W2 (Careful here: We decrease the b, so we need to multiply the shadow prices by -1.)

Conclusion: We should choose F2 and W2, because this would increase the transportation costs the least, namely by 600.

for i in V:
    print('Shadow price for', i, ':', model.dual[model.capacity[i]])

Shadow price for ('F1', 'F2') : 0.0
Shadow price for ('DC', 'W2') : -200.0

This means: Increasing the capacity limit for arc (F1,F2) will not reduce the cost at all. But increasing the capacity limit at arc (DC,W2) would reduce the cost.