Box Suite Recommendation

An online retailer needs a suite of $p\in\mathbb{N}$ boxes to ship its online customer orders. Moreover, $k\in\{0\}\cup\llbracket p-1\rrbracket=\left\{0:p-1\right\}$ boxes in the suite may be prescribed (or locked). For example, the online retailer may need some special boxes for shipping liquid-containing items such as liquid detergent. To save money, the online retailer wants such a suite that minimizes total shipping (shipping plus material) cost. One approach to designing this suite is to select $I\in\mathbb{N}$ historical customer shipments and $J\in\mathbb{N}$ candidate boxes, where $J>p$. Each historical customer shipment is assigned a unique index $i\in\llbracket I\rrbracket$, and each candidate box is assigned a unique index $j\in\llbracket J\rrbracket$. The set of historical customer shipments should be a small, but statistically significant, randomly sampled subset of the online retailer’s past (e.g., within the previous year) online customer shipments. Each historical customer shipment $i\in\llbracket I\rrbracket$ consists of $N_{i}\in\mathbb{N}_{0}$ 3D rectangular cartons with positive outer lengths, widths, and heights $\left\{\left(p_{in},q_{in},r_{in}\right)\right\}_{n=1}^{N_{i}}$, where $\left(p_{in},q_{in},r_{in}\right)\in\mathbb{R}_{>0}^{3}$ for $i\in\llbracket I\rrbracket$ and $n\in\llbracket N_{i}\rrbracket$, and $M_{i}\in\mathbb{N}_{0}$ foldable items with positive outer lengths, widths, and heights $\left\{\left(s_{im},t_{im},u_{im}\right)\right\}_{m=1}^{M_{i}}$, where $\left(s_{im},t_{im},u_{im}\right)\in\mathbb{R}_{>0}^{3}$ for $i\in\llbracket I\rrbracket$ and $m\in\llbracket M_{i}\rrbracket$. For the purposes of packing, foldable items are assumed to be liquid so that they may be deformed to fit into arbitrarily-shaped empty spaces inside a box. The set of $J$ candidate boxes should finely discretize the space of all possible boxes and must include the $k$ locked boxes that must be in the suite. The indices of those $k$ locked boxes are prescribed in the subset $T\subset\llbracket J\rrbracket$, where $\left|T\right|=k$. Each candidate box $j\in\llbracket J\rrbracket$ is characterized by a positive inner length, width, and height $\left(x_{j},y_{j},z_{j}\right)\in\mathbb{R}_{>0}^{3}$. Therefore, the inner volume of candidate box $j\in\llbracket J\rrbracket$ is $V_{j}=x_{j}y_{j}z_{j}\in\mathbb{R}_{>0}$. It is assumed that the candidate boxes are sorted by nondecreasing inner volume, which may be realized in $O\left(J\log J\right)$, so that $V_{j}\leq V_{j+1}$ for $j\in\llbracket J-1\rrbracket$. The optimal suite is obtained by selecting a subset $S\subset\llbracket J\rrbracket$ of the $J$ boxes, such that $\left|S\right|=p$ and $T\subset S$, that ships the $\hat{I}$ packable shipments, where $\hat{I}\leq I$ (since not all of the $I$ shipments necessarily fit in the $J$ candidate boxes), with minimum cost. Algorithm (1) gives a method for solving this optimization problem. The next few paragraphs describe the key parts of Algorithm (1).

The algorithm begins by sorting each box’s dimensions in nonincreasing order. Then, the algorithm determines whether each shipment fits into each candidate box, recording the result in the binary fitting matrix $B\in\mathbb{Z}_{2}^{I\times J}$. It is simple to solve the fitting problem when there are 0 or 1 cartons in the shipment. There are brute force algorithms for solving the fitting problem for 2 and 3 cartons in the shipment, but these are omitted from the algorithm for conciseness. For 2 or more cartons in the shipment, the algorithm first attempts to stack the cartons along each of the box’s three orthogonal axes, to see if they fit. If simple stacking does not work, then the algorithm uses the fitting MILP, a feasibility MILP described in Section 4, to solve the fitting problem. Since the fitting MILP must be solved via a third-party MILP solver, which may require a commercial license and which is computationally expensive, substantial run time improvements can be realized by utilizing the brute force fitting algorithms for 2 and 3 cartons. The source code accompanying [43] provides implementations of the brute force algorithms for 2 and 3 cartons; however, the algorithms presented in that code must be modified to permit rotations.

Also note that in order for a shipment to fit into a candidate box, the box’s inner volume must be greater than or equal to the shipment’s liquid volume and each individual carton in the shipment must fit inside the box. For efficiency, the algorithm checks that these necessary conditions are satisfied first before attempting to solve the stacking problem or fitting MILP when there are 2 or more cartons in the shipment. Moreover, if brute force fitting algorithms for 2 and 3 cartons are available, these may be used to check that all pairs and all triples of cartons fit inside the box before attempting to solve the stacking problem or fitting MILP.

Next, the algorithm removes shipments that could not be packed into any candidate box, leaving $\hat{I}\leq I$ packable shipments. For each packable shipment $i\in\llbracket\hat{I}\rrbracket$, $J_{i}\subset J$ denotes the set of candidate boxes into which packable shipment $i$ fits. Then, the algorithm constructs the nonnegative cost matrix $C\in\mathbb{R}_{\geq 0}^{\left(\hat{I}+k\right)\times J}$ which records the cost of shipping each packable shipment into each candidate box. If packable shipment $i\in\llbracket\hat{I}\rrbracket$ fits in candidate box $j\in\llbracket J\rrbracket$ (i.e., if $j\in J_{i}$), the cost $C_{ij}\in\mathbb{R}_{\geq 0}$ to ship packable shipment $i\in\llbracket\hat{I}\rrbracket$ in candidate box $j\in\llbracket J\rrbracket$ is computed. Note that in order to compute the cost for packable shipment $i\in\llbracket\hat{I}\rrbracket$, the data for shipment $W_{i}\in\llbracket I\rrbracket$ is needed; that is, $W_{i}$ is the shipment index of packable shipment $i$. The data for shipment $W_{i}\in\llbracket I\rrbracket$ may include the outer dimensions and weights of each item, the shipping carrier (e.g., USPS, FedEx, or UPS), the shipping service (e.g., 1 day, 2 day, or 3 day), and the shipping zone, which is determined by the locations of the shipping store or warehouse and the customer. The cost is computed via a detailed formula, which is omitted here, that depends on the shipping cost charged by the carrier and service combination, the cost of the corrugate used to construct the box, the cost to transport the box blank from the box manufacturer to the ratailer’s stores or warehouses, the cost of the dunnage used to fill the empty space between the packed shipment and the box’s interior, and the cost of the tape used to seal the top and bottom flaps of the box shut. More simply, if instead of minimizing cost, the online retailer wishes to minimize box outer (or inner) volume shipped or material weight shipped, then $C_{ij}$ is the outer (or inner) volume of box $j$ or the weight of the material (corrugate, dunnage, and tape) used to ship packable shipment $i$ in box $j$, respectively.

Let $\Gamma\in\mathbb{R}_{>0}$ be a sufficiently large positive real number such as $\left[\sum_{i=1}^{\hat{I}}\max_{j\in J_{i}}C_{ij}\right]+1\in\mathbb{R}_{>0}$. The constant $\Gamma$ serves as a penalty cost to impose constraints on the solution suite $S$. In order to ensure that the solution suite $S$ ships all the packable shipments, $C_{ij}$ is set to $\Gamma$ if packable shipment $i$ does not fit in candidate box $j$, i.e., if $j\in\llbracket J\rrbracket\setminus J_{i}$. In order to force the boxes in $T$ into the solution suite $S$, $k$ fake shipments are appended to the set of packable shipments and $k$ rows are added to the bottom of $C$, where each row, indexed by $i\in\left\{\hat{I}+1,\hat{I}+2,\ldots,\hat{I}+k\right\}$, stores the shipping costs for a fake shipment representing box $T_{i-\hat{I}}$. The cost of shipping the fake shipment indexed by $i\in\left\{\hat{I}+1,\hat{I}+2,\ldots,\hat{I}+k\right\}$ in box $T_{i-\hat{I}}$ is 0 and in all other boxes is $\Gamma$. That is, $C_{iT_{i-\hat{I}}}=0$ and $C_{ij}=\Gamma$ for $j\in\llbracket J\rrbracket\setminus\left\{T_{i-\hat{I}}\right\}$.

The optimization problem that must be solved is

Instead, the following optimization problem is solved:

The optimization problem (3.2) is an instance of the $p$-median problem [21, 22, 12], or more precisely the $p$-facility location problem [11]. In the literature, the $p$-median problem is sometimes also called the $k$-median problem [27, 28, 3]. Given $n\in\mathbb{N}$ customers, $m\in\mathbb{N}$ candidate facilities, $p\in\llbracket m\rrbracket$ candidate facilities to open, and nonnegative costs $d\in\mathbb{R}_{\geq 0}^{n\times m}$, where $d_{ij}$ is the cost of serving customer $i\in\llbracket n\rrbracket$ with candidate facility $j\in\llbracket m\rrbracket$, the $p$-median problem is to find (or open) a subset $S\subset\llbracket m\rrbracket$ comprised of $p$ candidate facilities that serves all the customers with minimum cost, that is, that minimizes the sum of the costs of serving each customer with its minimum cost open facility. Mathematically, the $p$-median problem is

To equate (3.2) to (3.3), the shipments are the customers, the boxes are the facilities, $\hat{I}+k=n$, $J=m$, and $C=d$.

The contrapositive of Lemma (1) gives the following corollary.

In summary, if (3.2) does not have a solution with cost less than $\Gamma$, then (3.1) does not have a solution, and if (3.2) does have a solution with cost less than $\Gamma$, then (3.1) has a solution and the solution set realized by (3.1) equals that realized by (3.2).

The $p$-median problem is NP-hard [36]. There are many methods to solve the $p$-median problem including MILP, Lagrangian relaxation, heuristics (e.g., myopic algorithm, neighborhood search, exchange algorithm, Lin–Kernighan neighborhood exchange algorithm), metaheuristics (e.g., simulated annealing, genetic algorithm, tabu search, heuristic concentration, variable neighborhood search, ant colony optimization, scatter search, and GRASP), and hyper-heuristics [12, 13, 38, 54, 45, 55].

The exchange algorithm [62] is a simple heuristic which starts with a feasible solution $S$ comprised of $p$ candidate facilities and finds a pair of facilities, $a\in S$ and $b\in S^{\mathsf{c}}$, such that the new feasible solution $\left(S\setminus\left\{a\right\}\right)\cup\left\{b\right\}$ decreases the cost of serving all the customers compared to the original feasible solution $S$. This exchange or swapping procedure is repeated iteratively until no further improvement is realized, depending on the definition of the local neighborhood about a feasible solution. There are several variations of the exchange algorithm which search over different local neighborhoods, including the first improvement [66], best improvement [23], and Lin–Kernighan [37] neighborhoods. Moreover, there are several efficient implementations of the exchange algorithm [66, 15, 23, 57] which provide significant computational savings compared to a naïve implementation. Due to its simplicity, the exchange algorithm is only effective at solving very small instances of the $p$-median problem. However, the exchange algorithm is often used as a local search algorithm by other more complicated solution methods, such as Lagrangian relaxation and metaheuristics.

Given $n\in\mathbb{N}$ customers, $m\in\mathbb{N}$ candidate facilities, $p\in\llbracket m\rrbracket$ candidate facilities to open, and nonnegative costs $d\in\mathbb{R}_{\geq 0}^{n\times m}$, where $d_{ij}$ is the cost of serving customer $i\in\llbracket n\rrbracket$ with candidate facility $j\in\llbracket m\rrbracket$, the MILP formulation of the $p$-median problem (3.3) is [13, 38]

In (3.13), $y\in\mathbb{Z}_{2}^{m}$ determines which candidate facilities are open and $x\in\mathbb{R}_{\left[0,1\right]}^{n\times m}$ determines which open facility is assigned to each customer. (3.13) formulates $x\in\mathbb{R}_{\left[0,1\right]}^{n\times m}$ instead of $x\in\mathbb{Z}_{2}^{n\times m}$, because the former imposes far fewer integrality constraints and therefore is much more computationally tractable for a MILP solver. However, in the solution of (3.13), $x_{ij}$ may not be 0 or 1 if customer $i\in\llbracket n\rrbracket$ may be served by more than one open facility with the same minimum cost. It is trivial to transform a solution $x\in\mathbb{R}_{\left[0,1\right]}^{n\times m}$ of (3.13) into an equivalent solution $\tilde{x}\in\mathbb{Z}_{2}^{n\times m}$ as follows. Let $S\equiv\left\{j\in\llbracket m\rrbracket\colon y_{j}=1\right\}\subset\llbracket m\rrbracket$ denote the open facilities in the solution of (3.13). Then, for each $i\in\llbracket n\rrbracket$, let

breaking ties arbitrarily if the minimum cost of serving customer $i$ over the open facilities in $S$ is not unique. Then for each $i\in\llbracket n\rrbracket$ and for each $j\in\llbracket m\rrbracket$, $\tilde{x}\in\mathbb{Z}_{2}^{n\times m}$ is constructed via:

A third-party MILP solver, such as Gurobi [20], CPLEX [26], or CBC [8], must be used to solve the MILP (3.13); the reader is referred to “Solving the Fitting MILP” in Section 4 for more information about the particular MILP solvers mentioned here. However, MILP solvers are only able to solve fairly small instances of the $p$-median problem to within a prescribed absolute or relative optimality gap of the globally optimal solution. Because MILP solvers rely on branch-and-bound algorithms, if an MILP solver is able to find a solution (not necessarily the globally optimal solution), it is also able to provide a lower bound of the globally optimal solution via the absolute or relative optimality gap.

The objective function $\mathcal{J}$ that is minimized in (3.13) is $\mathcal{J}\left(x\right)\equiv\sum_{i=1}^{n}\sum_{j=1}^{m}d_{ij}x_{ij}$. Construct the Lagrangian function $\mathcal{L}$ by adjoining the fourth constraint, $\sum_{j=1}^{m}x_{ij}=1\quad\forall i\in\llbracket n\rrbracket$, in (3.13) to the objective function $\mathcal{J}$ via the vector of Lagrange multipliers $\lambda\in\mathbb{R}^{n}$:

The primal form of the MILP formulation (3.13) of the $p$-median problem is

The primal problem (3.17) is equivalent to (3.13). Now construct the dual function $\mathcal{D}$:

The following are important properties of the dual function $\mathcal{D}$.

1. 1.

   For fixed $\lambda$, it is trivial to construct $x\in\mathbb{Z}_{2}^{n\times m}$ and $y\in\mathbb{Z}_{2}^{m}$ that minimize the Lagrangian $\mathcal{L}$ subject to satisfying the constraints in (3.18) [12, 13, 38]. Therefore, it is simple to solve the constrained Lagrangian minimization problem and compute $\mathcal{D}\left(\lambda\right)$.

2. 2.

   For fixed $\lambda$, the solution $x\in\mathbb{Z}_{2}^{n\times m}$ and $y\in\mathbb{Z}_{2}^{m}$ of $\mathcal{D}\left(\lambda\right)$ may be easily transformed into a corresponding primal solution, i.e., a feasible solution of the $p$-median problem [12, 13, 38].

3. 3.

   For fixed $\lambda$, $\mathcal{D}\left(\lambda\right)$ is a lower bound of the solution to the primal problem (3.17) [10].

4. 4.

   $\mathcal{D}$ is concave as a function of $\lambda$ [10].

5. 5.

   $\mathcal{D}$ is piecewise linear, and is therefore nondifferentiable, as a function of $\lambda$ [10].

The dual form of the MILP formulation (3.13) of the $p$-median problem is

Since the dual function $\mathcal{D}$ is concave and nondifferentiable, the dual problem (3.19) may be solved by subgradient optimization [12, 13, 38], which is a particular method of convex optimization. The solution to the dual problem (3.19) is a lower bound to the solution of the primal problem (3.17) [10]; in some cases, the solution to the dual problem (3.19) may even coincide with the solution to the primal problem (3.17). In the Lagrangian relaxation approach to solving the $p$-median problem (3.3), the dual problem (3.19) is solved via iterative subgradient optimization, starting from an initial guess $\lambda$ of the Lagrange multipliers. In each iteration of the subgradient optimization, 1) a lower bound is constructed by solving the constrained Lagrangian minimization problem and computing the dual function $\mathcal{D}$ in (3.18) for the current estimate $\lambda$ of the Lagrange multipliers, 2) a corresponding upper bound is constructed from the lower bound, and 3) the current estimate $\lambda$ of the Lagrange multipliers is updated. If in a given iteration, a new smallest upper bound is found, it may be further improved via a heuristic such as the exchange algorithm. As described in [13], if in a given iteration, a new smallest upper bound or a new largest lower bound is found, candidate facilities can be forced into and out of the optimal solution, which can reduce the computational time of Lagrangian relaxation. By keeping track of the smallest upper bound and the largest lower bound found over all iterations, Lagrangian relaxation provides both upper and lower bounds of the globally optimal solution to the $p$-median problem. The iterations continue until some stopping criterion is satisfied, such as executing a maximum number of iterations or realizing a prescribed absolute or relative optimality gap. Instead of adjoining the fourth constraint, $\sum_{j=1}^{m}x_{ij}=1\quad\forall i\in\llbracket n\rrbracket$, in (3.13) to the objective function $\mathcal{J}$, it is possible to adjoin the fifth constraint, $x_{ij}\leq y_{j}\quad\forall i\in\llbracket n\rrbracket\quad\forall j\in\llbracket m\rrbracket$, in (3.13) to the objective function $\mathcal{J}$ via a matrix of nonnegative Lagrange multipliers $\lambda\in\mathbb{R}_{\geq 0}^{n\times m}$; this alternative approach is discussed in [12].

POPSTAR [51] is a freely available $p$-median problem solver implemented in C++. POPSTAR solves the $p$-median problem via a hybrid metaheuristic that combines GRASP with path-relinking and the genetic algorithm [56]. Moreover, POPSTAR performs local searches via a fast implementation of the exchange algorithm [57]. Reference [45], published in 2007, comprehensively surveyed many methods, excluding hyper-heuristics, for solving the $p$-median problem and concluded that the overall algorithm implemented by POPSTAR is the best. dc2 [14] is a recent, freely available $p$-median problem solver implemented in C that is still under development. dc2 solves the $p$-median problem via several metaheuristics, including GRASP, path-relinking, and a new metaheuristic called disperse construction, and performs local searches via several fast implementations of the exchange algorithm [66, 23, 57]. Unlike POPSTAR, dc2 is multithreaded and therefore is able to exploit the parallelism offered by multicore CPUs. POPSTAR and dc2 can solve instances of the $p$-median problem with several thousand customers, several thousand candidate facilities, and $p\leq 20$ in a few minutes on a modern laptop. Several metaheuristics, including GRASP and the genetic algorithm, have been implemented on GPGPUs to solve the $p$-median problem [59, 1].

The box suite recommendation may be validated by packing the optimization shipment set and a much larger randomly sampled shipment set into the recommended suite $S^{*}$, where each shipment is assigned to the minimum cost box in the suite into which it fits. Several metrics, such as percentage of shipments packed into each box, percentage of total cost shipped by each box, percentage of all box outer volume shipped by each box, and percentage liquid void space, collected for both packings can be compared. If the metrics are similar, then the cost savings afforded by the recommended suite $S^{*}$ should be expected to hold on all shipments. An alternative validation method is to pack the optimization shipment set and a much larger randomly sampled shipment set into several suites $Q\cup S^{*}$, where each suite $S\in Q$ satisfies $T\subset S$, $\left|S\right|=p$, and $J_{i}\cap S\neq\O\,\forall i\in\llbracket\hat{I}\rrbracket$, and ensure that the cost reductions (comparing the cost of each suite $S\in Q$ to the cost of $S^{*}$) predicted by the optimization shipment set agree with those predicted by the large shipment set. If the various metrics or cost reductions are dissimilar, then Algorithm (1) should be run again using a larger set of randomly sampled historical customer shipments.

The recommended suite $S^{*}$ can be further refined (or fine-tuned) by running Algorithm (1) again on a new set of candidate boxes $J^{\prime}$ obtained by taking small variations above and below the inner lengths, widths, and heights of the unlocked boxes $S^{*}\setminus T$ in the recommended suite $S^{*}$. Like the original set of candidate boxes $J$, the new set of candidate boxes $J^{\prime}$ must also include the $k$ locked boxes.