I'm new to pyomo and not experienced with optimization tasks at all. But the problem I'm trying to solve is fairly simple, but pyomo can't find a solution to it. I have 10 entries (i index) of the same type and 3 boxes (or whatever). I'm trying to distribute my values in d (being for example kilos) between the three boxes puting them in such a way so that total weight of each box doesn't differ much from average weight, meaning that my ideal scenario is to uniformly spread elements in boxes so that weights of elements in the individual boxes are equal to average of all elements in d. Since weights of each elements different (and other factors) I would like it to be as close to the ideal case as possible (to uniform distribution by weight in the box). Unfortunately, solvers (cbc, glbk) provide only a trivial solution where it stacks everything in the first "box".
Please, help me to find my mistake or any alternative working solution. This task is crucial to me
# Create a concrete model
model = pyo.ConcreteModel()
# Set of indices I
model.I = pyo.RangeSet(1, 10)
# Set of indices J
model.J = pyo.RangeSet(1, 3)
# Parameters d and x
model.d = {1: 15, 2: 8, 3: 9, 4: 18, 5: 20, 6: 2, 7: 7, 8: 5, 9: 12, 10: 10}
model.x = pyo.Var(model.I, model.J, domain=pyo.Binary)
# average of all d's
avg = 0
for key in model.d:
avg += model.d[key]
avg /= len(model.J)
# Constraints
def constraint_rule(model, i):
return sum(model.x[i, j] for j in model.J) == 1
# Objective function
model.obj = pyo.Objective(expr=sum(sum(model.d[i] * model.x[i, j] for i in model.I)-avg for j in model.J), sense=pyo.minimize)
model.constraints = pyo.Constraint(model.I, rule=constraint_rule)
# Solve the optimization problem
solver = pyo.SolverFactory('cbc')
results = solver.solve(model)
Here's what is happening in your model: All choices are optimal.
In your objective function, you are adding up the deviations from the average. Good idea. However, all configurations lead to the same result because some boxes are heavier, some are lighter, but those are positive and negative deviations that all (algebraically) must add up to zero, because they are deviations from average. If you inspect the objective value for code as written, it is ~1.5e-15 (or zero).
So, you need a new approach. One way to do this is to just look at the deviations above the average and minimize the sum of them. Basically smushing them all down as much as possible to the average and ignoring the negative deviations. This is a formulation that does that. Realize there are other options as well, such as minimizing the maximum weight, etc.
Aside: there are 2 "shortcuts" in here. The first is the alternate way to make a constraint with a decorator, which I find more eloquent (there are many ways to make a constraint in pyomo). The second is pyo.sum_product which does what you think it does... for one variable, it is just syntactic sugar for the grand sum. See the dox for more.
import pyomo.environ as pyo
# Create a concrete model
model = pyo.ConcreteModel()
# Set of indices I
model.I = pyo.RangeSet(1, 10)
# Set of indices J
model.J = pyo.RangeSet(1, 3)
# Parameters d and x
model.d = {1: 15, 2: 8, 3: 9, 4: 18, 5: 20, 6: 2, 7: 7, 8: 5, 9: 12, 10: 10}
model.x = pyo.Var(model.I, model.J, domain=pyo.Binary)
model.amount_above_average = pyo.Var(model.J, domain=pyo.NonNegativeReals)
# average of all d's
avg = 0
for key in model.d:
avg += model.d[key]
avg /= len(model.J)
# Constraints
@model.Constraint(model.J)
def above_average_constraint(model, j):
return model.amount_above_average[j] >= sum(model.x[i, j] * model.d[i] for i in model.I) - avg
def constraint_rule(model, i):
return sum(model.x[i, j] for j in model.J) == 1
# Objective function
model.obj = pyo.Objective(expr=pyo.sum_product(model.amount_above_average))
model.constraints = pyo.Constraint(model.I, rule=constraint_rule)
# Solve the optimization problem
solver = pyo.SolverFactory('cbc')
results = solver.solve(model)
print(results) # <--- you MUST check the status == optimal
print(f'obj: {pyo.value(model.obj)}')
model.x.display()
for i, j in model.x.index_set():
if model.x[i, j].value:
print(f'put {i} in {j}')
Problem:
- Name: unknown
Lower bound: 0.66666667
Upper bound: 0.66666667
Number of objectives: 1
Number of constraints: 13
Number of variables: 33
Number of binary variables: 30
Number of integer variables: 30
Number of nonzeros: 3
Sense: minimize
Solver:
- Status: ok
User time: -1.0
System time: 0.06
Wallclock time: 0.06
Termination condition: optimal
Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Black box:
Number of iterations: 128
Error rc: 0
Time: 0.07272696495056152
Solution:
- number of solutions: 0
number of solutions displayed: 0
obj: 0.66666667
x : Size=30, Index=x_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
(1, 1) : 0 : 0.0 : 1 : False : False : Binary
(1, 2) : 0 : 1.0 : 1 : False : False : Binary
(1, 3) : 0 : 0.0 : 1 : False : False : Binary
(2, 1) : 0 : 0.0 : 1 : False : False : Binary
(2, 2) : 0 : 0.0 : 1 : False : False : Binary
(2, 3) : 0 : 1.0 : 1 : False : False : Binary
(3, 1) : 0 : 0.0 : 1 : False : False : Binary
(3, 2) : 0 : 0.0 : 1 : False : False : Binary
(3, 3) : 0 : 1.0 : 1 : False : False : Binary
(4, 1) : 0 : 1.0 : 1 : False : False : Binary
(4, 2) : 0 : 0.0 : 1 : False : False : Binary
(4, 3) : 0 : 0.0 : 1 : False : False : Binary
(5, 1) : 0 : 0.0 : 1 : False : False : Binary
(5, 2) : 0 : 1.0 : 1 : False : False : Binary
(5, 3) : 0 : 0.0 : 1 : False : False : Binary
(6, 1) : 0 : 0.0 : 1 : False : False : Binary
(6, 2) : 0 : 0.0 : 1 : False : False : Binary
(6, 3) : 0 : 1.0 : 1 : False : False : Binary
(7, 1) : 0 : 1.0 : 1 : False : False : Binary
(7, 2) : 0 : 0.0 : 1 : False : False : Binary
(7, 3) : 0 : 0.0 : 1 : False : False : Binary
(8, 1) : 0 : 0.0 : 1 : False : False : Binary
(8, 2) : 0 : 0.0 : 1 : False : False : Binary
(8, 3) : 0 : 1.0 : 1 : False : False : Binary
(9, 1) : 0 : 0.0 : 1 : False : False : Binary
(9, 2) : 0 : 0.0 : 1 : False : False : Binary
(9, 3) : 0 : 1.0 : 1 : False : False : Binary
(10, 1) : 0 : 1.0 : 1 : False : False : Binary
(10, 2) : 0 : 0.0 : 1 : False : False : Binary
(10, 3) : 0 : 0.0 : 1 : False : False : Binary
put 1 in 2
put 2 in 3
put 3 in 3
put 4 in 1
put 5 in 2
put 6 in 3
put 7 in 1
put 8 in 3
put 9 in 3
put 10 in 1