Applied Optimization with Julia
Question: What are current trends in e-commerce?
After this lecture, you will be able to:
Question: Why did the shop send two?
Question: Why might they occur?
Question: What are the consequences?
Question: What are possible mitigations?
Question: What could be our objective?
We aim to improve the SKU2-warehouse allocation to minimize the number of split parcels resulting from SKUs being stored in different warehouses.
Question: What could be the sets here?
Question: What are possible parameters?
Question: What could the transactional data look like?
| \(t_{m,i}\) | A | B | C | D |
|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 0 |
| 2 | 1 | 1 | 1 | 0 |
| 3 | 1 | 1 | 0 | 0 |
| 4 | 1 | 0 | 0 | 1 |
| 5 | 1 | 0 | 0 | 1 |
| 6 | 1 | 0 | 0 | 1 |
| 7 | 1 | 0 | 0 | 1 |
| 8 | 0 | 0 | 1 | 1 |
Question: What is your opinion on the assumption?
Question: What could be our decision variable/s?
Question: Does anyone have an idea what this could mean?
10,000 SKUs, 1,000,000 orders and 2 warehouses already lead to \(2 \times 10^{10}\) binary variables \(Y_{m,i,k}\)!
Focus on the warehouses and the coappearance of SKUs! Discard the exact information about the customer orders.
Question: What could be the objective?
Maximize the coappearance of products that are often part of the same customer orders.
T = [
1 1 1 0;
1 1 1 0;
1 1 0 0;
1 0 0 1;
1 0 0 1;
1 0 0 1;
1 0 0 1;
0 0 1 1
]
# Create the coappearance matrix
Q = T' * T
println("Coappearance matrix Q:")
display(Q)Coappearance matrix Q:4×4 Matrix{Int64}:
7 3 2 4
3 3 2 0
2 2 3 1
4 0 1 5Question: What do the principal diagonal values tell us?
Question: What could be the sets?
We can focus on the SKUs and the warehouses, making the problem much smaller!
Question: What are possible parameters?
Instead of the transactional data, we just use the coappearance matrix in our model!
Maximize the coappearance of products that are often part of the same customer orders. In more mathematical terms: Maximize the sum of all unique pair-wise values \(q_{i,j}\) of all SKUs stored in the same warehouse.
Question: What could be our decision variable/s?
As we don’t need the customer order information, we only need to make a decision for each SKU and warehouse pair!
Question: How could we formulate the variable in Julia?
using JuMP, SCIP # SCIP is a non-commercial solver for quadratic models
warehouses = ["Hamburg", "Berlin"] # Add warehouses as a vector
skus = ["A", "B", "C", "D"] # The four SKUs from our transactional data
warehouse_model = Model(SCIP.Optimizer)@variable(warehouse_model, X[i in skus, k in warehouses], Bin)2-dimensional DenseAxisArray{JuMP.VariableRef,2,...} with index sets:
Dimension 1, ["A", "B", "C", "D"]
Dimension 2, ["Hamburg", "Berlin"]
And data, a 4×2 Matrix{JuMP.VariableRef}:
X[A,Hamburg] X[A,Berlin]
X[B,Hamburg] X[B,Berlin]
X[C,Hamburg] X[C,Berlin]
X[D,Hamburg] X[D,Berlin]Maximize the sum of all unique pair-wise values \(q_{i,j}\) of all SKUs stored in the same warehouse. Note, that this is a quadratic objective function!
Question: What could the objective function look like?
\[\text{Maximize} \quad \sum_{i=2}^{|\mathcal{I}|} \sum_{j=1}^{i-1} \sum_{k \in \mathcal{K}} X_{i,k}\times X_{j,k} \times q_{i,j}\]
The quadratic terms are \(X_{i,k}\times X_{j,k}\). This objective function is based on the Quadratic Multiple Knapsack Problem (QMKP), formulated by Hiley and Julstrom (2006).
Question: How could we formulate this in Julia?
# We reuse the matrix Q = T' * T computed from our transactional data
@objective(warehouse_model,
Max,
sum(
X[skus[i], warehouses[k]] * X[skus[j], warehouses[k]] * Q[i,j]
for i in 2:length(skus)
for j in 1:i-1
for k in 1:length(warehouses)
)
)3 X[B,Hamburg]*X[A,Hamburg] + 3 X[B,Berlin]*X[A,Berlin] + 2 X[C,Hamburg]*X[A,Hamburg] + 2 X[C,Berlin]*X[A,Berlin] + 2 X[C,Hamburg]*X[B,Hamburg] + 2 X[C,Berlin]*X[B,Berlin] + 4 X[D,Hamburg]*X[A,Hamburg] + 4 X[D,Berlin]*X[A,Berlin] + X[D,Hamburg]*X[C,Hamburg] + X[D,Berlin]*X[C,Berlin]
Question: What constraints?
Ensure that each SKU is allocated at least once.
Question: What could the constraint look like?
\[\sum_{k \in \mathcal{K}} X_{i,k} \geq 1 \quad \forall i \in \mathcal{I}\]
Question: How could we change the constraint to ensure that each SKU is allocated only once?
Question: How could we add the constraint in Julia?
@constraint(warehouse_model, single_allocation[i in skus],
sum(X[i, k] for k in warehouses) >= 1
)1-dimensional DenseAxisArray{JuMP.ConstraintRef{JuMP.Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.GreaterThan{Float64}}, JuMP.ScalarShape},1,...} with index sets:
Dimension 1, ["A", "B", "C", "D"]
And data, a 4-element Vector{JuMP.ConstraintRef{JuMP.Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.GreaterThan{Float64}}, JuMP.ScalarShape}}:
single_allocation[A] : X[A,Hamburg] + X[A,Berlin] ≥ 1
single_allocation[B] : X[B,Hamburg] + X[B,Berlin] ≥ 1
single_allocation[C] : X[C,Hamburg] + X[C,Berlin] ≥ 1
single_allocation[D] : X[D,Hamburg] + X[D,Berlin] ≥ 1Ensure that the capacity of each warehouse is not exceeded.
Question: What could the second constraint be?
\[\sum_{i \in \mathcal{I}} X_{i,k} \leq c_k \quad \forall k \in \mathcal{K}\]
And that’s basically it!
Question: How could we add the second constraint in Julia?
capacities = Dict("Hamburg" => 2, "Berlin" => 2) # Add capacities
@constraint(warehouse_model, capacity[k in warehouses],
sum(X[i, k] for i in skus) <= capacities[k]
)1-dimensional DenseAxisArray{JuMP.ConstraintRef{JuMP.Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, JuMP.ScalarShape},1,...} with index sets:
Dimension 1, ["Hamburg", "Berlin"]
And data, a 2-element Vector{JuMP.ConstraintRef{JuMP.Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, JuMP.ScalarShape}}:
capacity[Hamburg] : X[A,Hamburg] + X[B,Hamburg] + X[C,Hamburg] + X[D,Hamburg] ≤ 2
capacity[Berlin] : X[A,Berlin] + X[B,Berlin] + X[C,Berlin] + X[D,Berlin] ≤ 2\[\text{Maximize} \quad \sum_{i=2}^{|\mathcal{I}|} \sum_{j=1}^{i-1} \sum_{k \in \mathcal{K}} X_{i,k}\times X_{j,k} \times q_{i,j}\]
subject to:
\[ \begin{align*} & \sum_{k \in \mathcal{K}} X_{i,k} \geq 1 && \forall i \in \mathcal{I}\\ & \sum_{i \in \mathcal{I}} X_{i,k} \leq c_{k} && \forall k \in \mathcal{K}\\ & X_{i,k} \in \{0,1\} && \forall i \in \mathcal{I}, k \in \mathcal{K} \end{align*} \]
set_attribute(warehouse_model, "display/verblevel", 0) # Hide solver output
optimize!(warehouse_model)
println("The optimal objective value is: ", objective_value(warehouse_model))
println("The optimal solution is: ", value.(X))The optimal objective value is: 6.0
The optimal solution is: 2-dimensional DenseAxisArray{Float64,2,...} with index sets:
Dimension 1, ["A", "B", "C", "D"]
Dimension 2, ["Hamburg", "Berlin"]
And data, a 4×2 Matrix{Float64}:
1.0 0.0
0.0 1.0
0.0 1.0
1.0 0.0Question: What does this value tell us?
It sums the coappearances of SKUs sharing a warehouse: \(q_{A,D} + q_{B,C} = 4 + 2 = 6\). It does not directly count the avoided split orders!
Commercial solvers are faster and more robust than open source solvers but also more expensive. During your studies, you can use most of them for free though! Nonetheless, we will only use open source solvers in this course.
Our model is a MIQP, as only the objective is quadratic. It is not a MIQCP, which would also allow quadratic constraints. MIQCP solvers, e.g. SCIP, subsume MIQPs and can thus solve our model.
| SKUs | SCIP | HiGHS | Gurobi |
|---|---|---|---|
| 100 | 118s | X | ~400s |
| 1,000 | 1,011s (18%) | X | ~200s (2%) |
| 10,000 | X | X | X |
Questions: On model assumptions
With slack capacity, the model may store one SKU in several warehouses. A shared pair then counts once per warehouse, raising the objective without avoiding any further splits.
Allocation changes:
Workload imbalance:
Detect dependencies between products and allocate them accordingly, as products within orders can have dependencies and products are bought with different frequencies!
Question: When are two SKUs independent?
If purchasing one doesn’t affect purchasing the other:
\[P(A \cap B) = P(A) \times P(B)\]
Example: If \(A\) appears in 10% of orders and \(B\) in 20%:
Use chi-square tests to detect if SKUs are dependent:
Phase 1: Build allocation using dependencies
Phase 2: Refine with local search
Result: Works on large instances, 50,000 SKUs in 10 minutes!
Question: What challenges might arise in practice?
Question: Which of 100,000+ SKUs goes to which warehouse?
In the following case study (Vlćek and Voigt 2024), we:
In the case study, SKUs were grouped by category-brand combinations to reduce the problem size while maintaining allocation quality. New SKUs inherit the allocation of their cluster!
Problem characteristics:
Heuristics compared7:
Retailer’s status quo: 6.95% of all orders were split
Theoretical reduction of these splits:
We now have covered the Quadratic Multiple Knapsack Problem and are ready to start solving some tasks in the upcoming tutorial.
References for this lecture:
For more interesting literature to learn more about Julia, take a look at the literature list of this course.
Forecast made in 2020, not an actual number↩︎
SKU: Stock Keeping Unit↩︎
Greedy Randomized Adaptive Search Procedure, Zhu et al. (2021)↩︎
Based on numerical experiments with QMK instances. Times are averages over the instances solved to optimality; the percentage shows the share of instances solved within one hour, X = none solved within one hour.↩︎
QMKP is not applicable for instance in case study↩︎