from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("glpk")

T = ['A', 'B', 'C', 'D', 'E']

U = {('I','A'):2,
     ('I','B'):3,
     ('I','C'):9,
     ('A','D'):5,
     ('B','A'):7,
     ('B','E'):9,
     ('E','O'):9,
     ('C','B'):2,
     ('C','E'):1,
     ('D','O'):4,
     ('D','E'):1,
     ('A','E'):6}

A = list(U.keys())

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

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 == OutFlow
    
model.transshipment = Constraint(T, rule=flow_rule)

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

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

model.objective = Objective(expr = sum(model.f[i,j] for (i,j) in A if j=='O'), sense=maximize)

model.dual = Suffix(direction=Suffix.IMPORT)
results = opt.solve(model)
model.objective.expr()

8.0

for (i,j) in A:
    print ((i,j), 
           model.dual[model.capacity[i,j]], 
           model.capacity[i,j].uslack())

('I', 'A') 1.0 0.0
('I', 'B') 1.0 0.0
('I', 'C') 0.0 6.0
('A', 'D') 0.0 3.0
('B', 'A') 0.0 7.0
('B', 'E') 0.0 4.0
('E', 'O') 0.0 3.0
('C', 'B') 1.0 0.0
('C', 'E') 1.0 0.0
('D', 'O') 0.0 2.0
('D', 'E') 0.0 1.0
('A', 'E') 0.0 6.0