Search code examples
setinitializationcombinationspyomo

PYOMO LP, combination set


I would like your help with PYOMO LP. I am not sure what I am doing wrong so any feedback would be helpfull:

My data set is like bellow:

#The demand for each order
demand= {782912: 808, 782913: 3188, 782914: 2331, 782915: 847, 782916: 2163,789954:5643}

#The cost per unit produced for each order based in which factory chosen
total cost= { (782912, 'PLANT16'): 0.46, (782913, 'PLANT16'): 0.46, (782914, 'PLANT16'): 0.46, (782915, 'PLANT16'): 0.46, (782916, 'PLANT16'): 0.46}, (789954,'PLANT05'):0.90,(789954,'PLANT07'):0.91,(789954,'PLANT08'):1.13,(789954,'PLANT10'):0.12}

#The capacity of each factory
supply= {'PLANT05': 531,'PLANT07': 841,'PLANT08': 1107,'PLANT10': 981,'PLANT16': 2313}

#defining the model

model=pyo.ConcreteModel()
#Sets
model.i=pyo.Set(initialize=demand.keys())#orders
model.j=pyo.Set(initialize=supply.keys()) #factories
model.select_combos = pyo.Set(within = model.i * model.j, initialize = total_costs_per_unit.keys())
#Parameters

model.p=pyo.Param(model.i,model.j,initialize=total_costs_per_unit) # here goes the cost dictionary
p=model.p

model.d=pyo.Param(model.i,initialize=demand)
d=model.d

model.s=pyo.Param(model.j,initialize=supply)
s=model.s

#Decision variable
model.x=pyo.Var(model.i,model.j,within=pyo.NonNegativeReals)
x=model.x

#Objective function
def obj_rul(model):

  return sum(p[i,j]*x[i,j]  for i,j in model.select_combos)
#warning , not all combinations of i,j exist in my model.p, as they would not be valid solutions for the problem
model.Obj=pyo.Objective(rule=obj_rul,sense=pyo.minimize)

#Constraints

def Const1(model,i):
  return sum(x[i,j] for j in model.j)>=d[i]
model.condemand=pyo.Constraint(model.i,rule=Const1)

def Const5(model,j):
  return sum(x[i,j] for i in model.i)<=s[j]
model.consupply=pyo.Constraint(model.j,rule=Const5)

#Solving
Solver=SolverFactory('glpk')

results=Solver.solve(model)

print(results)
print('Obj funct=',model.Obj())
for i in model.i:
  for j in model.j:
    print('for order',i," from plant ",j, " is sent ",x[i,j]())

The error message I am getting is :

ERROR:pyomo.core:evaluating object as numeric value: x[782912,PLANT16]
    (object: <class 'pyomo.core.base.var._GeneralVarData'>)
No value for uninitialized NumericValue object x[782912,PLANT16]

All in all, I think if I had all combinations of 'orders' & 'plants', it would not had an issue. As it is now I dont know how to solve this problem.


Solution

  • You are on the right track here and your suspicion is correct.

    You are getting the error because some combinations of i, j are not used in your model at all, so they never become initialized, and when you try to print their value, you have problems. There are other issues in your model as well dealing with the same thing.

    So you have several options here.

    You have an "implied" constraint here that if you don't have a production cost for something, then it cannot be produced at that factory, so you could add to your data, and put in "production caps" of some kind for each item-factory combo... not real fun.

    Or you could put in a very large cost for the production of the "illegal" items so that the solver doesn't pick them.

    But you already did part of the work to what I think is the best solution, which is making a set of "legal combinations" of i, j. You just didn't bring it to the finish line. You can/should use this subset as the basis for the things that naturally comply with it like x and p. See my code. You just need to be careful for the sums that you only use i, j combos that are legal, which is pretty easy to do by screening them as I show in your updated constraints.

    Note: I multiplied all of your capacities by 10x because the initial model wasn't feasible. It is, this runs, and you should be able to modify this easy enough as needed.

    Code

    import pyomo.environ as pyo
    
    #The demand for each order
    demand= {782912: 808, 782913: 3188, 782914: 2331, 782915: 847, 782916: 2163,789954:5643}
    
    #The cost per unit produced for each order based in which factory chosen
    production_cost= {  (782912, 'PLANT16'): 0.46, 
                        (782913, 'PLANT16'): 0.46, 
                        (782914, 'PLANT16'): 0.46, 
                        (782915, 'PLANT16'): 0.46, 
                        (782916, 'PLANT16'): 0.46, 
                        (789954, 'PLANT05'): 0.90,
                        (789954, 'PLANT07'): 0.91,
                        (789954, 'PLANT08'): 1.13,
                        (789954, 'PLANT10'): 0.12}
    
    #The capacity of each factory
    supply= {'PLANT05': 5310,'PLANT07': 8410,'PLANT08': 11070,'PLANT10': 9810,'PLANT16': 23130}
    
    model=pyo.ConcreteModel()
    
    #Sets
    model.i=pyo.Set(initialize=demand.keys()) #orders
    model.j=pyo.Set(initialize=supply.keys()) #factories
    model.select_combos = pyo.Set(within = model.i * model.j, initialize = production_cost.keys())
    
    #Parameters
    model.p=pyo.Param(model.select_combos, initialize=production_cost) # here goes the cost dictionary
    # p=model.p  <-- I wouldn't rename things...confusing
    
    model.d=pyo.Param(model.i,initialize=demand)
    # d=model.d
    
    model.s=pyo.Param(model.j,initialize=supply)
    # s=model.s
    
    #Decision variable
    model.x=pyo.Var(model.select_combos, within=pyo.NonNegativeReals)
    # x=model.x
    
    #Objective function
    # def obj_rul(model):
    #   return sum(p[i,j]*x[i,j]  for i,j in model.select_combos)
    #warning , not all combinations of i,j exist in my model.p, as they would not be valid solutions for the problem
    
    # your objective is a simple sum() expression, so you don't need the complication of a rule-function combo...
    model.Obj=pyo.Objective(expr=sum(model.p[i,j]*model.x[i,j] for i,j in model.select_combos),sense=pyo.minimize)
    
    #Constraints
    
    def Const1(model,i):
      return sum(model.x[i,j] for j in model.j if (i,j) in model.select_combos)>=model.d[i]
    model.condemand=pyo.Constraint(model.i,rule=Const1)
    
    def Const5(model,j):
      return sum(model.x[i,j] for i in model.i if (i,j) in model.select_combos)<=model.s[j]
    model.consupply=pyo.Constraint(model.j,rule=Const5)
    
    #Solving
    Solver=pyo.SolverFactory('glpk')
    
    results=Solver.solve(model)
    
    print(results)
    #print('Obj funct=',model.Obj())
    for i, j in model.select_combos:
      print('for order',i," from plant ",j, " is sent ", pyo.value(model.x[i,j]))