I am trying to optimize the allocation of of products inside a fictive warehouse, having a predefined number of aisles (3 in the example code below) where products can be placed. For now, the only optimization criterion that I would like to impose in the objective function is:

  • products that have been frequently part of the same order (i.e., people ordering product $A$ also ordered product $B$) should be placed in the same aisle

As input I have a list of historical orders (randomly created in the example), where each order is composed of one or more products/articles. Using the list I build a co-occurrence matrix, that shows the number of times two articles were part of the same order, for each pair of articles.

The following code (full of comments) shows how I implemented the model using Docplex (IBM Decision Optimization CPLEX Optimizer Modeling for Python):

from docplex.mp.model import Model
import numpy as np
import pandas as pd
import docplex

# docplex version: 2.22.213
# Cplex version:

################## Clean and Prepare Data ##################

# number of (historical) orders containing one or more articles
num_orders = 50
# number of different articles (products) available
num_articles = 40
# number of records, i.e., couples (o, a) of order-article, meaning that the article a is part of order o
s = (250, 1)

# aisles of the warehouse, with respective capacity (maximum number of articles that can be placed in them)
aisles = {1: 15, 2: 17, 3: 20}

# building a dataframe from the input data
order_col = np.random.randint(low=1, high=num_orders, size=s)
arts_col = np.random.randint(low=1, high=num_articles, size=s)
temp = np.array([1]*s[0]).reshape(s)
cols = np.concatenate((order_col, arts_col, temp), axis=1)
df = pd.DataFrame(cols, columns=["ORDER_NUMBER", "ARTICLE_CODE", "TEMP"])

# missing articles due to randomness in building the "ARTICLE_CODE" column of df
missing_articles_randomness = set(range(1, num_articles+1)).difference(set(df["ARTICLE_CODE"]))

to_add = pd.DataFrame()
to_add["ORDER_NUMBER"] = df["ORDER_NUMBER"].iloc[:len(missing_articles_randomness)]
to_add["ARTICLE_CODE"] = list(missing_articles_randomness)
to_add["TEMP"] = 1

# they need to be appended to the dataframe so that the corresponding row and column shows up in the co-occurrence matrix
df = df.append(to_add, ignore_index=True)

# create the co-occurrence matrix
temp_count = df.pivot_table(index="ORDER_NUMBER", columns="ARTICLE_CODE", 
                            values="TEMP", aggfunc="count").fillna(0)

coocc = temp_count.T.dot(temp_count)
np.fill_diagonal(coocc.values, 0)


################## Create Optimization Model ##################

mdl = Model("coocc_v1")

R = range(1, num_articles+1)
id_1 = [(i, j) for i in R for j in R[i:]]
# binary variable indicating whether two articles are located in the same aisle
arts_same_aisle = mdl.binary_var_dict(id_1, name="asa")

id_2 = [(c, j) for c in aisles for j in R]
# binary variable indicating whether an article a is located in an aisle c
art_aisle = mdl.binary_var_dict(id_2, name="aa")

# each article can be placed in only one aisle
constr_1 = mdl.add_constraints([mdl.sum([art_aisle[c, j] for c in aisles]) == 1 for j in R])

# each aisle has a maximum capacity
constr_2 = mdl.add_constraints([mdl.sum([art_aisle[c, j] for j in R]) <= cap for c, cap in aisles.items()])

# definition of arts_same_aisle_i_j: 1 if article i and article j are in the same aisle, 0 otherwise
for i in R:
    for j in R[i:]:
        arts_same_aisle[(i, j)] = mdl.logical_or(*[mdl.logical_and(art_aisle[c, i] == 1, art_aisle[c, j] == 1) for c in aisles])

# the 'co-occurrence score', that is, positioning articles in the same aisle if they are frequently part of the same order 
obj1 = mdl.sum([arts_same_aisle[i, j] * coocc.loc[i, j] for i in R for j in R[i:]])

# the goal is to maximize the objective function

################## Model Solving ##################

solution = mdl.solve(log_output=True)
# Total (root+branch&cut) =  759.07 sec. (477653.82 ticks)
# integer optimal solution
# objective value: 393.0

sol_df = solution.as_df()
sol_df = sol_df.loc[sol_df["name"].str.startswith("aa")]

for n, i in zip(["aisle", "article"], [1, 2]):
    sol_df[n] = sol_df["name"].apply(lambda x: x.split("_")[i])

s = sol_df.groupby("aisle")["article"].apply(list)

which gives

1                                         [15, 26, 40]
2    [1, 2, 3, 4, 5, 6, 9, 14, 17, 18, 19, 20, 23, ...
3    [7, 8, 10, 11, 12, 13, 16, 21, 22, 24, 27, 30,...
Name: article, dtype: object

Unfortunately, the solving time explodes when the number of different articles or aisles increases, and I could not find a non-trivial reformulation of the model that significantly reduces the solving time of the problem.

Is there a way to reformulate the problem to make it usable in a real-life scenario, with thousands of articles and tens / few hundreds of aisles?


2 Answers 2


This is an instance of the quadratic assignment problem, which is known to be difficult.

You have binary variable $x_{ic}$ to indicate whether article $i$ is assigned to aisle $c$ and binary variable $y_{ij}$ to indicate whether articles $i$ and $j$ are assigned to the same aisle. Your formulation maximizes the total reward $$\sum_{i<j} r_{ij} y_{ij}$$ subject to \begin{align} \sum_c x_{ic} &= 1 &&\text{for all $i$} \tag1 \\ \sum_i x_{ic} &\le b_c &&\text{for all $c$} \tag2 \\ y_{ij} &= \vee_c x_{ic} x_{jc} &&\text{for all $i<j$} \tag3 \end{align} Constraint $(1)$ assigns each article to exactly one aisle. Constraint $(2)$ imposes capacity $b_c$ on aisle $c$. Constraint $(3)$ enforces the definition of $y_{ij}$, which the solver presumably linearizes in the usual way by introducing a new variable $z_{ijc}$ to represent the product $x_{ic}x_{jc}$ and then replacing $(3)$ with linear constraints \begin{align} z_{ijc} &\le x_{ic} &&\text{for all $i<j$ and $c$} \tag4 \\ z_{ijc} &\le x_{jc} &&\text{for all $i<j$ and $c$} \tag5 \\ y_{ij} &\le \sum_c z_{ijc} &&\text{for all $i<j$} \tag6 \end{align} You can try explicitly replacing $(3)$ with $(4)$ through $(6)$. You can also try relaxing $y$ and $z$ to be nonnegative rather than binary. You can try omitting $y$ and $(6)$ and replacing the objective function with $\sum_{i<j} r_{ij} \sum_c z_{ijc}$.

An alternative approach is to use compact linearization as follows.

Introduce binary (or nonnegative) variable $w_{icjd}$ to represent $x_{ic}x_{jd}$, and maximize $$\sum_{i<j}r_{ij} \sum_c w_{icjc}$$ subject to $(1)$, $(2)$, and \begin{align} \sum_c w_{icjd} &= x_{jd} &&\text{for all $i,j,d$} \tag7 \\ w_{icjd} &\le x_{ic} &&\text{for all $i,j,c,d$} \tag8 \\ w_{icjd} &= w_{jdic} &&\text{for all $i<j,c,d$} \tag9 \\ w_{icic} &= x_{ic} &&\text{for all $i,c$} \tag{10} \end{align} Note that $(7)$ is obtained by multiplying $(1)$ by $x_{jd}$.

Yet another approach is to solve this directly as a binary quadratic programming problem: maximize $$\sum_{i<j}r_{ij} \sum_c x_{ic}x_{jc}$$ subject to $(1)$ and $(2)$.

  • 1
    $\begingroup$ Some runtime results from the reformulations, using a random co-occurrence matrix with 40 articles and three aisles as in the Python example: - original model: 4472 seconds - replace (3) with (4)->(6): 615 seconds - omit y and (6) and rewrite objective: 417 seconds - compact linearization: 3073 seconds - binary quadratic: 343 seconds $\endgroup$ Mar 2, 2022 at 13:25

I would try the reformulations suggested by @RobPratt first. If those don't scale sufficiently well, you might consider the following.

Suppose, (very) hypothetically, that we had a list of all possible aisle layouts. Let $A$ be a binary matrix with one row for every item and one column for every layout, such that $A_{ij}=1$ if and only if layout $j$ includes item $i.$ Assume that layouts were designed to respect the capacity limit. Moreover, assume you have computed a colocation score $c_j$ for each layout $j.$

Now create a binary variable $x_j$ for each layout $j$, 1 if you use it and 0 if not. Your first constraint is $$\sum_j x_j \le N$$where $N$ is the number of aisles available. . For each article $i$, you have a constraint $$\sum_j a_{ij} x_j = 1$$ forcing you to use exactly one layout containing item $i$. Your objective is to maximize $$\sum_j c_j x_j,$$ the aggregate colocation score.

Now to address the elephant in the model: the number of possible aisle layouts grows combinatorially with aisle capacity and item count, so enumerating all possible layouts will likely be prohibitive. That brings you to a couple of options.

  • For an exact solution (proven optimal), you could explore branch and price algorithms. These require specialized software (as well as a MIP solver) and there is no guarantee you get an answer in an acceptable amount of time and memory.
  • If you can live with a heuristic approach, one possibility is something along the lines of the Gilmore-Gomory approach for the one dimensional cutting stock problem. The logic is relatively straightforward.
    • Build the model I described, but with the integrality restriction relaxed ($x_j \in [0,1]$) and with a handful of aisle layouts to start, chosen so that the model is feasible (every item is in one aisle, no aisle exceeds capacity, colocation is ignored).
    • Solve the current LP relaxation and grab the dual solution.
    • Use the dual in a subproblem (which I think might be a quadratic binary knapsack in your case, but don't quote me on that) to build a new aisle layout that would have a positive reduced cost in the LP relaxation of the master problem.
    • Add that layout to the relaxed master and repeat.

You stop when the subproblem based on the latest dual values fails to find a new layout. At that point, you restore the integrality restrictions, solve the integer program and declare approximate victory.

  • 2
    $\begingroup$ This will need slight modification from the usual Gilmore-Gomory approach because here each aisle has a (potentially) different capacity. The $\sum_j x_j \le N$ constraint will instead be $\sum_j x_j \le N_k$, where $N_k$ is the number of aisles in "class" $k$. And there will be a separate subproblem for each class. $\endgroup$
    – RobPratt
    Feb 26, 2022 at 0:43
  • $\begingroup$ @RobPratt I take it that aisles having different capacities is in the Python code? In that case, you're right about separate subproblems and separate summations (where the index $j$ runs over aisles within a capacity class). It seems more realistic to me to assume that all aisles have the same capacity but items have different capacity utilization (size), in which case there would be one subproblem, still a knapsack. $\endgroup$
    – prubin
    Feb 26, 2022 at 3:44
  • $\begingroup$ Yes, the Python code has capacities $15$, $17$, and $20$ for the three aisles. $\endgroup$
    – RobPratt
    Feb 26, 2022 at 4:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.