Search code examples
pythonoptimizationpyomo

Pyomo: Weird results for a simple 'Transportation theory' problem


I'm trying to solve the classical 'Transportation theory' to optimize the production and transportation of goods from 4 factories to 4 distribution centers, through 2 weeks.

There are:

  • 4 max capacities in each factory (items/week), the same ones for the 2 weeks
  • 4x4 transportation costs from each factory to each distribution center ($/item), the same ones for the 2 weeks
  • 4 weekly demand volume in each distribution center (items/week), the same ones for the 2 weeks

Code:

from pyomo.environ import *

model = ConcreteModel()

# DATA
costs = {\
    (1, 'Factory A1', 'DC 1'): 1000,(1, 'Factory A2', 'DC 1'): 1001, (1, 'Factory A3', 'DC 1'): 1002, (1, 'Factory B', 'DC 1'): 5000, \
    (1, 'Factory A1', 'DC 2'): 2000,(1, 'Factory A2', 'DC 2'): 2001, (1, 'Factory A3', 'DC 2'): 2002, (1, 'Factory B', 'DC 2'): 6000, \
    (1, 'Factory A1', 'DC 3'): 3000,(1, 'Factory A2', 'DC 3'): 3001, (1, 'Factory A3', 'DC 3'): 3002, (1, 'Factory B', 'DC 3'): 7000, \
    (1, 'Factory A1', 'DC 4'): 4000,(1, 'Factory A2', 'DC 4'): 4001, (1, 'Factory A3', 'DC 4'): 4002, (1, 'Factory B', 'DC 4'): 8000, \
    (2, 'Factory A1', 'DC 1'): 1000,(2, 'Factory A2', 'DC 1'): 1001, (2, 'Factory A3', 'DC 1'): 1002, (2, 'Factory B', 'DC 1'): 5000, \
    (2, 'Factory A1', 'DC 2'): 2000,(2, 'Factory A2', 'DC 2'): 2001, (2, 'Factory A3', 'DC 2'): 2002, (2, 'Factory B', 'DC 2'): 6000, \
    (2, 'Factory A1', 'DC 3'): 3000,(2, 'Factory A2', 'DC 3'): 3001, (2, 'Factory A3', 'DC 3'): 3002, (2, 'Factory B', 'DC 3'): 7000, \
    (2, 'Factory A1', 'DC 4'): 4000,(2, 'Factory A2', 'DC 4'): 4001, (2, 'Factory A3', 'DC 4'): 4002, (2, 'Factory B', 'DC 4'): 8000 \
    }
capacities = {\
    (1, 'Factory A1'): 1000, (1, 'Factory A2'): 1000, (1, 'Factory A3'): 1000, (1, 'Factory B'): 1100, \
    (2, 'Factory A1'): 1000, (2, 'Factory A2'): 1000, (2, 'Factory A3'): 1000, (2, 'Factory B'): 1100}
demand= {(1, 'DC 1'): 500, (1, 'DC 2'): 501, (1, 'DC 3'): 502, (1, 'DC 4'):503, \
    (2, 'DC 1'): 500, (2, 'DC 2'): 501, (2, 'DC 3'): 502, (2, 'DC 4'): 503}

M_prod = max(capacities.values())  # big M value for max production capacity for constraint use

# SETS
model.T = Set(initialize=[1, 2])  # Weeks
model.F = Set(initialize=['Factory A1','Factory A2', 'Factory A3', 'Factory B'])  # Factory
model.D = Set(initialize=["DC 1", "DC 2", "DC 3", "DC 4"])         # Destination

# PARAMS
model.cap = Param(model.T, model.F, initialize=capacities, default=0)
model.demand = Param(model.T, model.D, initialize=demand, default=0)
model.cost = Param(model.T, model.F, model.D, initialize=costs, default=10)

# VARS
model.x = Var(model.T, model.F, model.D, bounds=(0.0, None))  # Production
model.running =   Var(model.F, model.T, within=Binary)    # factory f running in period t
model.startup =   Var(model.F, model.T, within=Binary)    # startup cost induced in period t by factory f

# CONSTRAINTS
def f_supply(model, t, f):
    return sum(model.x[t, f, d] for d in model.D) <= model.cap[t, f]
model.supply = Constraint(model.T, model.F, rule=f_supply)

def f_demand(model, t, d):
    return sum(model.x[t, f, d] for f in model.F) >= model.demand[t, d]
model.meet_demand = Constraint(model.T, model.D, rule=f_demand)

# OBJETIVE
def f_objetive(model):
    return sum(
        model.cost[t, f, d] * model.x[t, f, d]
        for t in model.T for f in model.F for d in model.D
    )
model.objetive = Objective(rule=f_objetive, sense=minimize) 

opt = SolverFactory("glpk")
results = opt.solve(model)
print('\nTotal cost = ', model.objetive())
model.x.display()

Results:

Total cost =  10042024.0
x : Size=32, Index=x_index
    Key                       : Lower : Value : Upper : Fixed : Stale : Domain
    (1, 'Factory A1', 'DC 1') :   0.0 : 500.0 :  None : False : False :  Reals
    (1, 'Factory A1', 'DC 2') :   0.0 :   0.0 :  None : False : False :  Reals
    (1, 'Factory A1', 'DC 3') :   0.0 : 500.0 :  None : False : False :  Reals
    (1, 'Factory A1', 'DC 4') :   0.0 :   0.0 :  None : False : False :  Reals
    (1, 'Factory A2', 'DC 1') :   0.0 :   0.0 :  None : False : False :  Reals
    (1, 'Factory A2', 'DC 2') :   0.0 : 497.0 :  None : False : False :  Reals
    (1, 'Factory A2', 'DC 3') :   0.0 :   0.0 :  None : False : False :  Reals
    (1, 'Factory A2', 'DC 4') :   0.0 : 503.0 :  None : False : False :  Reals
    (1, 'Factory A3', 'DC 1') :   0.0 :   0.0 :  None : False : False :  Reals
    (1, 'Factory A3', 'DC 2') :   0.0 :   4.0 :  None : False : False :  Reals
    (1, 'Factory A3', 'DC 3') :   0.0 :   2.0 :  None : False : False :  Reals
    (1, 'Factory A3', 'DC 4') :   0.0 :   0.0 :  None : False : False :  Reals
     (1, 'Factory B', 'DC 1') :   0.0 :   0.0 :  None : False : False :  Reals
     (1, 'Factory B', 'DC 2') :   0.0 :   0.0 :  None : False : False :  Reals
     (1, 'Factory B', 'DC 3') :   0.0 :   0.0 :  None : False : False :  Reals
     (1, 'Factory B', 'DC 4') :   0.0 :   0.0 :  None : False : False :  Reals
    (2, 'Factory A1', 'DC 1') :   0.0 : 500.0 :  None : False : False :  Reals
    (2, 'Factory A1', 'DC 2') :   0.0 : 500.0 :  None : False : False :  Reals
    (2, 'Factory A1', 'DC 3') :   0.0 :   0.0 :  None : False : False :  Reals
    (2, 'Factory A1', 'DC 4') :   0.0 :   0.0 :  None : False : False :  Reals
    (2, 'Factory A2', 'DC 1') :   0.0 :   0.0 :  None : False : False :  Reals
    (2, 'Factory A2', 'DC 2') :   0.0 :   1.0 :  None : False : False :  Reals
    (2, 'Factory A2', 'DC 3') :   0.0 : 502.0 :  None : False : False :  Reals
    (2, 'Factory A2', 'DC 4') :   0.0 : 497.0 :  None : False : False :  Reals
    (2, 'Factory A3', 'DC 1') :   0.0 :   0.0 :  None : False : False :  Reals
    (2, 'Factory A3', 'DC 2') :   0.0 :   0.0 :  None : False : False :  Reals
    (2, 'Factory A3', 'DC 3') :   0.0 :   0.0 :  None : False : False :  Reals
    (2, 'Factory A3', 'DC 4') :   0.0 :   6.0 :  None : False : False :  Reals
     (2, 'Factory B', 'DC 1') :   0.0 :   0.0 :  None : False : False :  Reals
     (2, 'Factory B', 'DC 2') :   0.0 :   0.0 :  None : False : False :  Reals
     (2, 'Factory B', 'DC 3') :   0.0 :   0.0 :  None : False : False :  Reals
     (2, 'Factory B', 'DC 4') :   0.0 :   0.0 :  None : False : False :  Reals
1
​

I am not able to understand the results (I find them a bit weird):

  • if the capacities, the costs and the demand are the same for the two weeks, why are the results different, with different final costs?
  • why is the prompted cost (10042024) different from the cost I calculate from the given results (10001999)?

Solution

  • The answer returned by GLPK (and then Pyomo) is correct (insomuch as it matches the data in your example).

    The reason that the "years" are different is due to degeneracy in your model formulation: any allocation of production to factories that assigns 1000 units (of any product) to Factory A1, 1000 units to Factory A2, and 6 units to Factory A3 in each year will yield an (optimal) objective of 5021012 for that year. The solver can return any valid combination, and which combination you get back depends on the solver implementation (thing like the row/column ordering and the solver's internal simplex pivot implementation).