Here is a variation of Problem 2 from Homework Sheet 6.

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

We use the numbers from the problem to set up the following LP transportation problem.

P = [1, 2, 3] # List of plants

DC = [1, 2, 3, 4] # List of distribution centers (DC)
d = {1:10, 2:10, 3:10, 4:10} # demand at the DCs
s = {1:12, 2:17, 3:11} # supply at the plants
# Here: total supply = total demand (both are 40)

# Distance from plant i to DC j
c = {(1,1):800, (1,2):1300, (1,3):400, (1,4):700,
(2,1):1100, (2,2):1400, (2,3):600, (2,4):1000,
(3,1):600, (3,2):1200, (3,3):800, (3,4):900}

model = ConcreteModel()

model.x = Var(P, DC, within=NonNegativeReals)

#Transportation cost given in the problem description
model.z = Objective(expr = sum((0.5*c[i,j]+100)*model.x[i,j] for i in P for j in DC), sense=minimize)

# Usual supply rule
def supply_rule (model, i):
    return sum(model.x[i,j] for j in DC) <= s[i]
model.supply = Constraint(P, rule=supply_rule)

# Usual demand rule
def demand_rule (model, j):
    return sum(model.x[i,j] for i in P) >= d[j]
model.demand = Constraint(DC, rule=demand_rule)

# Import dual variables and solve model
model.dual = Suffix(direction=Suffix.IMPORT)
results = opt.solve(model)

model.x.get_values()

{(1, 1): 0.0,
 (1, 2): 0.0,
 (1, 3): 2.0,
 (1, 4): 10.0,
 (2, 1): 0.0,
 (2, 2): 9.0,
 (2, 3): 8.0,
 (2, 4): 0.0,
 (3, 1): 10.0,
 (3, 2): 1.0,
 (3, 3): 0.0,
 (3, 4): 0.0}

model.z.expr()

20200.0

for j in DC:
    print(j, model.dual[model.demand[j]])

1 500.0
2 800.0
3 400.0
4 550.0

for i in P:
    print(i, model.dual[model.supply[i]])

1 -100.0
2 0.0
3 -100.0

Now suppose we can increase the supply at one of the plants by one unit. Thus, we can also increase the demand at one of the DCs by one unit.

Question: Which plant and DC should we choose in order to increase the cost the least?

Answer: The shadow price at DC3 is the smallest (400), and the shadow prices at either plant 1 or 3 are the most negative (-100). So if we choose to increase the supply by one unit at plants 1 or 3 and ship to DC3 we would increase the costs by 400-100 = 300 to a total of 20500.

Change the code above to see this result by changing the respective supplies and demands.

---

Next, let us solve part (b) from the homework sheet. Here, we introduce an extra dummy supply (plant 4), which has 0 cost to ship to any of the DCs.

# New parameters
P = [1, 2, 3, 4]
DC = [1, 2, 3, 4]
d = {1:15, 2:10, 3:10, 4:10}
s = {1:12, 2:17, 3:11, 4:5}
c = {(1,1):800, (1,2):1300, (1,3):400, (1,4):700,
(2,1):1100, (2,2):1400, (2,3):600, (2,4):1000,
(3,1):600, (3,2):1200, (3,3):800, (3,4):900,
(4,1):0, (4,2):0, (4,3):0, (4,4):0}

# Same code as above
model = ConcreteModel()
model.x = Var(P, DC, within=NonNegativeReals)
model.z = Objective(expr = sum((0.5*c[i,j]+100)*model.x[i,j] for i in P for j in DC), sense=minimize)
def supply_rule (model, i):
    return sum(model.x[i,j] for j in DC) <= s[i]
model.supply = Constraint(P, rule=supply_rule)
def demand_rule (model, j):
    return sum(model.x[i,j] for i in P) >= d[j]
model.demand = Constraint(DC, rule=demand_rule)
model.dual = Suffix(direction=Suffix.IMPORT)
results = opt.solve(model)

model.x.get_values()

{(1, 1): 4.0,
 (1, 2): 0.0,
 (1, 3): 0.0,
 (1, 4): 8.0,
 (2, 1): 0.0,
 (2, 2): 5.0,
 (2, 3): 10.0,
 (2, 4): 2.0,
 (3, 1): 11.0,
 (3, 2): 0.0,
 (3, 3): 0.0,
 (3, 4): 0.0,
 (4, 1): 0.0,
 (4, 2): 5.0,
 (4, 3): 0.0,
 (4, 4): 0.0}

All the Virtual "dummy" supply from the "dummy plant" 4 is shipped to DC2, so DC2 is undersupplied by 5 units.

model.z.expr()

19700.0