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:
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):
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).