Reformulate CPLEX optimization model of Warehouse Product Allocation

Question

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

print(docplex.__version__)
# docplex version: 2.22.213
# Cplex version: 20.1.0.0

################## 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)

coocc.head()


################## 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
mdl.maximize(obj1)

################## 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)
print(s)

which gives

aisle
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?

Answer 1

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

Answer 2

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.