Assignment of colored boxes to groups

I am working on a problem that I think can be solved as an extension to the assignment problem, but am open to hear about other ways to solve.

PROBLEM

There are a lot of colored boxes. For each color there are many different boxes that contain different number of items (some may contain 0 items, in which case we ignore this box). The goal is combine the different colored boxes into groups so that each group contains 1 box of each available color and when we remove the same number of items from each box in each group (so there will be at least 1 box with 0 items in it in each group) we minimize the number of items still in the boxes.

for example, using the following example: $$\begin{bmatrix}90 & 80 & 75 & 70\\55 & 85 & 35 & 65\\125 & 95 & 90 & 95\\45 & 110 & 95 & 115\\50 & 100 & 90 & 100\end{bmatrix}$$ Where each row represents boxes of a particular color and each column represents the number of items in each of the boxes. Here we have 5 different colors and 4 different boxes in each color. While in this example there are no empty boxes, there could be empty boxes. Basically I want to permute the rows in such a way that a number of groups are formed (4 in this case) such that each group contains only 1 box of each color and for each group we subtract the number of items in each box from the box with the minimum number of items (giving us the remaining items) and sum those remaining items up. Do this for each group, then sum up the remaining items from every group. The goal is to minimize this total remaining items count.

So for the above matrix one of several optimal solutions would be: $$\begin{bmatrix}75 & 70 & 80 & 90\\55 & 65 & 35 & 85\\125 & 90 & 95 & 95\\115 & 110 & 45 & 95\\100 & 100 & 50 & 900\end{bmatrix}$$ Where the sum of each element in a column - min element of each column and the sum of those differences $$\begin{bmatrix}195 & 110 & 130 & 30\end{bmatrix} = \begin{bmatrix}465\end{bmatrix}$$

I get optimal solutions to small problems like this, but I need to be able to go up to 24x24 size matrices and when I get to larger matrices my model will run forever. I have let it run up to 48 hours with the following problem without finding an optimal solution: $$\begin{bmatrix}90 & 76 & 75 & 70 & 50 & 74 & 12 & 68 & 70\\35 & 85 & 55 & 65 & 48 & 101 & 70 & 83 & 0\\125 & 95 & 90 & 105 & 59 & 120 & 36 & 73 & 59\\45 & 110 & 95 & 115 & 104 & 83 & 37 & 71 & 0\\60 & 105 & 80 & 75 & 59 & 62 & 93 & 88 & 59\\45 & 65 & 110 & 95 & 47 & 31 & 81 & 34 & 0\\38 & 51 & 107 & 41 & 69 & 99 & 115 & 48 & 51\\47 & 85 & 57 & 71 & 92 & 77 & 109 & 36 & 0\\39 & 63 & 97 & 49 & 118 & 56 & 92 & 61 & 56\\47 & 101 & 71 & 60 & 88 & 109 & 52 & 90 & 52\end{bmatrix}$$ It will find a feasible solution within a minute, but make no progress the rest of the time. It may be that I need to run it with an early stopping routine and take the best feasible solution, but I am hoping to get some feedback on how to improve constraints to improve the solution finding.

WHAT I HAVE TRIED

I have coded this using ortools CP-SAT solver. The code is shown below, but I have let this run for over 48 hours without finding an optimal solution. It finds a good feasible solution within 1 minute and then no further updates for the remaining 48 hours. I realize that there are $$(8!)^{10}$$ combinations (my math may be wrong but you get the point), but I was hoping the solver could find an optimal solution quicker than 48 hours, so I feel I may be missing a critical constraint in my model or some symmetry breaking constraint or some hint that may speed the solver.

from ortools.sat.python import cp_model
import pandas as pd

def solve_assignment(cost, groups):
"""
Solve assignment problem of reducing cost of assigning boxwe to groups
Arguments:  cost: cost matrix of assigning item
groups: list of groups

Returns:    solution_list: list of dictionaries [Box, Color, Group, Cost] minimizing the costs,
status:  status of solver - string (OPTIMAL, FEASIBLE, INFEASIBLE)
"""
maxL = max([max(map(int, i)) for i in cost]) + 1
num_colors = len(cost)
num_boxes = len(cost[1])
num_groups = len(groups)
all_colors = range(num_colors)
all_boxes = range(num_boxes)
all_groups = range(num_groups)

# _____________________
# Create Model Object
# ---------------------
model = cp_model.CpModel()

# -----------------------
# Create Model Variables
# -----------------------
# objective is our objective variable
objective = model.NewIntVar(0, 1000000, 'objective')

# Minimum length in each group
min_length = [model.NewIntVar(0, 200000, f'min_length_{g}') for g in all_groups]

# Helper variables to find minimum non-zero value
smallest = [[[model.NewBoolVar(f'smallest_{g}{c}{b}') for b in all_boxes] for c in all_colors] for g in all_groups]
yi = [[[model.NewIntVar(0, 200000, f'yi_{g}{c}{b}') for b in all_boxes] for c in all_colors] for g in all_groups]

# min_length[g] * M (length of group[g] (used for objective)
mmin_length = [model.NewIntVar(0, 200000, f'max_length_{g}') for g in all_groups]
sum_g = [model.NewIntVar(0, 2000000, f'sum_{g}') for g in all_groups]

# Indicator Variables for color/reel in a group x[group, color, reel]
X = [[[model.NewBoolVar(f'x[{g},{c},{b}]') for b in all_boxes] for c in all_colors] for g in all_groups]

# -------------------------
# Create Model Constraints
# -------------------------
# Each group is assigned at most one box from each color (0 length reels don't need to be assigned).
for g in all_groups:
[model.Add(sum(X[g][c][b] for b in all_boxes) <= 1) for c in all_colors]

# Each box from every color is assigned to at most one group
for b in all_boxes:
[model.Add(sum(X[g][c][b] for g in all_groups) <= 1) for c in all_colors]

# Must have the correct number of boxes in each group
for g in all_groups:
model.Add(sum(X[g][c][b] for b in all_boxes for c in all_colors) == groups[g])

# Don't use boxes that have 0 items
for g in all_groups:
for b in all_boxes:
for c in all_colors:
if cost[c][b] == 0:

# Determine smallest non-zero item count in group - Because some X[g][c][r]'s will be zero
for g in all_groups:
for c in all_colors:
for b in all_boxes:
# smallest[g][c][r] == True if X[g][c][r] == 1

# smallest[g][c][r] == False if X[g][c][r] == 0

# yi[g][c][r] == cost[c][r] if smallest[g][c][r] == True

# Set yi[g][c][r] to a large enough value that it can't be min if smallest[g][c][r] == False

# min_length[g] = min(yi[g][c][r] for c in all_colors for r in all_reels)
model.AddMinEquality(min_length[g], [yi[g][c][b] for c in all_colors for b in all_boxes])

# Determine M*min_length in a group (to use in cost function)
# sum(C[g][c][b] - min_length[g]) is same as sum(sum_g[g] - group_length[g]*min_length[g])
for g in all_groups:
model.Add(sum_g[g] == sum(cost[c][b]*X[g][c][b] for c in all_colors for b in all_boxes))

# --------------
# Break symmetry
# --------------
# TODO - Find ways to break symmetry
# Can swap if y[g1][c1][r1] - y[g2][c2][r2] <= y[g2][c2][r2]-min_g[2] and y[g1][c1][r1]-min_g[1]

# Objective cost - Minimize the differences between box item counts in each group (same as min(sum_g-m*min))
model.Add(objective == sum(sum_g[g] - mmin_length[g] for g in all_groups))
model.Minimize(objective)

# ---------
# ---------

# Setup Early Stopping when solutions don't improve within time limit (s)
#solution_printer = PrintLogger([objective], early_stopping)

# -------------------------------------
# Create a solver and solve the model
# -------------------------------------
solver = cp_model.CpSolver()
solver.parameters.num_search_workers = 8

# Sets a time limit of approx 48 hours.
solver.parameters.max_time_in_seconds = 173000.0
# solver.parameters.search_branching = cp_model.AUTOMATIC_SEARCH

# Solve the Model
status = solver.Solve(model)
#status = solver.SolveWithSolutionCallback(model, solution_printer)
#solution_printer.clean()

solution_list = []
print(solver.StatusName(status))
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
for g in all_groups:
for c in all_colors:
for b in all_boxes:
if solver.Value(X[g][c][b]) == 1:
solution_row = {'Box': b, 'Color': c, 'Group': g, 'Cost': cost[c][b]}
solution_list.append(solution_row)
print('Statistics')
print("Min Length: ", [solver.Value(min_length[i]) for i in range(len(min_length))])
print('  - conflicts : %i' % solver.NumConflicts())
print('  - branches  : %i' % solver.NumBranches())
print('  - wall time : %f s' % solver.WallTime())

return solution_list, solver.StatusName()


And it can be run with the following code:

if __name__ == '__main__':

'''
cost =  [[90, 76, 75, 70, 50, 74, 12, 68, 70],
[35, 85, 55, 65, 48, 101, 70, 83, 0],
[125, 95, 90, 105, 59, 120, 36, 73, 59],
[45, 110, 95, 115, 104, 83, 37, 71, 0],
[60, 105, 80, 75, 59, 62, 93, 88, 59],
[45, 65, 110, 95, 47, 31, 81, 34, 0],
[38, 51, 107, 41, 69, 99, 115, 48, 51],
[47, 85, 57, 71, 92, 77, 109, 36, 0],
[39, 63, 97, 49, 118, 56, 92, 61, 56],
[47, 101, 71, 60, 88, 109, 52, 90, 52]]

groups = [10, 10, 10, 10, 10, 10, 10, 10, 6]
'''

cost = [
[70, 75, 80, 90],
[35, 55, 65, 85],
[90, 95, 95, 125],
[45, 95, 110, 115],
[50, 90, 100, 100],
]

groups = [5, 5, 5, 5]

solution, sol = solve_assignment(cost, groups)
d = pd.DataFrame(solution)

# Calculate the remnant associated with each reel (length - min of group)
d['Remnant'] = d['Cost'] - d.groupby('Group')['Cost'].transform('min')
print(sum(d['Remnant']))


I realize there are a lot of combinations to search through and it maybe that this is just a difficult problem, but I'm hoping someone can help me crack this.

Here's a mixed integer linear programming formulation that I used to find an optimal solution with objective value $$1506$$ for your $$10 \times 9$$ instance. Let $$(c_{ij})$$ be the given (sparse) cost matrix. Let binary decision variable $$x_{ijg}$$ indicate whether cell $$(i,j)$$ (with $$c_{ij} \not= 0$$) is assigned to group $$g$$. Let decision variable $$y_g$$ represent $$\min\limits_{i,j:x_{ijg}=1} c_{ij}$$. The problem is to minimize $$\sum_g \left(\sum_{i,j} c_{ij} x_{ijg} - n_g y_g\right)$$ subject to \begin{align} \sum_g x_{ijg} &=1 &&\text{for all i,j} \tag1\label1\\ \sum_{i,j} x_{ijg} &=n_g &&\text{for all g} \tag2\label2\\ \sum_j x_{ijg} &\le 1 &&\text{for all i,g} \tag3\label3\\ x_{ijg} = 1 &\implies y_g \le c_{ij} &&\text{for all i,j,g} \tag4\label4 \end{align} Constraint \eqref{1} assigns each cell to exactly one group. Constraint \eqref{2} assigns exactly $$n_g$$ cells to group $$g$$. Constraint \eqref{3} prevents assigning more than one cell from each row to the same group. Indicator constraint \eqref{4} enforces $$y_g \le \min\limits_{i,j:x_{ijg}=1} c_{ij}$$. For solvers that do not support indicator constraints, you can replace \eqref{4} with a big-M constraint $$y_g - c_{ij} \le M(1-x_{ijg})$$

An optimal solution is:

70  68  50  90  76  75  70  12  74
48  70  55 101  85  83  65  35
36 125 105  95 120  90  59  73  59
45  71  83 104 115  95 110  37
62  75  59  93 105  80  88  59  60
34  65  47  95 110  81  45  31
38  69  48 115 107  99  51  41  51
71  85  47  92 109  77  57  36
49  63  61  92  97 118  56  39  56
60  88  47  90 109 101  52  52  71


Note that the $$x$$ part of the objective is constant because of \eqref{1}, so you could equivalently maximize $$\sum_g n_g y_g$$.

• Thanks @RobPratt, changing the min value to implication of X[g][c][r] seemed to do the trick - I should have seen that earlier. I had to leave constraint (1) as a <= otherwise the model became infeasible (I may be forming the constraints incorrectly though). Do you know why the way I originally had the implication would cause it to run so much longer? Feb 4 at 14:04
• Regarding (1), I omit the zero-cost cells from the problem altogether, whereas you force $x=0$ for those. You can keep the constraint as $\le 1$ if you impose it only for nonzero-cost cells. I don’t use ortools, so I’ll leave it to others to comment on its behavior. Feb 4 at 15:08
• Figured out why I couldn't get constraint (1) as an equals and fixed that. Tremendous improvement in speed. Thanks! Feb 4 at 15:55