Bin packing with lexicographic objectives for loading weight- and volume-constrained trucks in a direct-shipping system

We consider a packing problem that arises in a direct-shipping system in the food and beverage industry: Trucks are the containers, and products to be distributed are the items. The packing is constrained by two independent quantities, weight (e.g., measured in kg) and volume (number of pallets). Additionally, the products are grouped into the three categories: standard, cooled, and frozen (the latter two require refrigerated trucks). Products of different categories can be transported in one truck using separated zones, but the cost of a truck depends on the transported product categories. Moreover, splitting orders of a product should be avoided so that (un-)loading is simplified. As a result, we seek for a feasible packing optimizing the following objective functions in a strictly lexicographic sense: minimize the (1) total number of trucks; (2) number of refrigerated trucks; (3) number of refrigerated trucks which contain frozen products; (4) number of refrigerated trucks which also transport standard products; (5) and minimize splitting. This is a real-world application of a bin-packing problem with cardinality constraints a.k.a. the two-dimensional vector packing problem with additional constraints. We provide a heuristic and an exact solution approach. The heuristic meta-scheme considers the multi-compartment and item fragmentation features of the problem and applies various problem-specific heuristics. The exact solution algorithm covering all five stages is based on branch-and-price using stabilization techniques exploiting dual-optimal inequalities. Computational results on real-world and difficult self-generated instances prove the applicability of our approach.


Introduction
In this paper, we present a system of bin-packing problems that arise in a directshipping system in the food and beverage industry. More than 1.2 million units of different products leave the factory of our industry partner every day. By far, the largest share of the products is transported by trucks from the factory to the distribution centers of supermarket chains (in the following denoted as warehouses). As shipping is done directly from the factory to individual warehouses with full-truck volumes, the routing of trucks plays no role. For each warehouse, the use of the utilized trucks has to be optimized by assigning shipments to trucks. Hence, the overall problem naturally decomposes by warehouse.
The transportation of products is standardized and uses euro-pallets (which are typically packed with a set of uniform boxes of a product). In our application, quantities are large so that pallets are not mixed, i.e., each pallet contains only one product. Moreover, it is legitimate to assume that products are equally distributed on pallets such that each pallet loaded with a certain product has the same weight. Let K = {1, 2, . . . , m} denote the set of products (to be shipped to a particular warehouse on a specific day). As the range of food products includes a huge variety of commodities, e.g., fruit juices, spirits, baked goods, jams, wafers, and ketchup, the weight and space requirements for the transport differ substantially by product. For each product k ∈ K , these requirements are therefore described by two attributes: • the unit weight ρ k (in kg/pallet) of a pallet loaded with product k and • the number n k of pallets that need to be shipped.
In addition, all products can be uniquely grouped into three categories: • Standard products I S require no cooling; • Cooled products I C require cooling; • Shipping of frozen products I F requires deep cooling and is the most costly transport.
As a result, the set of products is partitioned into I = I S ∪ I C ∪ I F .
The fleet consists of standard trucks and refrigerated trucks, i.e., trucks without and with a cooling system. Frozen and cooled products must be transported in refrigerated trucks, which are more costly than standard trucks. Frozen and cooled products can be mixed (using different zones), but costs are higher if frozen products are transported. Also, standard products can be transported by refrigerated trucks (again using separated zones). There are no conflicts between products, i.e., any subset of products can be packed together into a truck. All trucks are homogeneous regarding capacity. Let W be the capacity (in kg) for the physical weight and B be the capacity for the number of pallets that fit into one truck (usually B = 33 pallets).
Moreover, depending on the customer or warehouse, two policies are applied regarding the splitting of orders of the same product: forbidden: Splitting is forbidden so that all n k pallets loaded with product k ∈ K must be transported by the same truck. This is the standard case because it simplifies the handling and processing in the loading and unloading phases. allowed: Splitting is allowed so that the n k pallets of a product k ∈ K may be distributed arbitrarily over the trucks. Note that neither pallets nor packages on the pallets are split. Orders are split so that pallets of one product are transported with more than one truck. This splitting-allowed policy can help to reduce the total number of used trucks but causes inconveniences. It should not occur more often than necessary (see objective 5. below).
For a given solution, we denote a product assigned to two or more trucks as split product.
We want to present both policies, splitting-forbidden and splitting-allowed, as variants of the same packing problem. For this purpose, we introduce items which are the atomic, not-divisible objects in the bin-packing model. The set K of products and the set I of items are identical, i.e., K = I . However, we write k ∈ K to refer to products and i ∈ I to refer to items. Depending on the policy, we define two weights for each item: w i refers to the physical weight and b i to the volume expressed as number of pallets. Table 1 and Example 1 clarify the correspondence.
With these definitions of items i ∈ I with two weights (w i , b i ) i∈I and two capacities (W , B), the problem of assigning products to trucks becomes a two-dimensional bin/vector-packing problem (see Sect. 2 for references to the pertinent literature). Trivially, a feasible packing exists only if w i ≤ W and b i ≤ B holds for all items i ∈ I . Table 1 Relationship between products and items for the splitting-forbidden and splitting-allowed policies What also makes our real-world problem interesting is the non-standard objective. We are seeking for a feasible product-truck assignment optimizing the following five objectives in a strictly lexicographic sense: 1. Minimize the total number of trucks; 2. Minimize the number of refrigerated trucks; 3. Minimize the number of refrigerated trucks transporting frozen products; 4. Minimize the number of refrigerated trucks which also transport standard products; and 5. Minimize order splittings, which is either 5a. Minimize the overall number of splits; or 5b. Minimize the number of products which are split (over any number of trucks).
The fifth objective is only relevant for the splitting-allowed policy. Altogether, we show that this lexicographic objective can be optimized by solving a sequence of extended bin-packing problems. Note that it is not a-priori given how to measure the impact of order splitting. The corresponding objective will depend on the actual handling of split orders after their delivery at the destination warehouse. According to our industry partner, there are two different strategies observed with its customers: One customer handles products separately for every truck. This means that for every truck the processing effort depends on the number of different products (but also on other factors, such as the number of pallets). Therefore, it is cost-efficient to minimize the overall number of splits, i.e., packing one product separated into two trucks is preferred over packing it separated into three trucks, as stated in objective 5a. A second customer collects all pallets carrying the same product, but arriving on different trucks, in the unloading zone and then inserts them into the warehouse in one batch. Thus, the main effort of splitting is due to the necessity of reserving a collection area, which happens whenever a product order is split. Whether pallets are collected from two or more trucks hardly plays a role, as stated in objective 5b.
Note that for both splitting policies, objectives 1.-4. have to be observed in a lexicographic sense. This means that, e.g., a solution with 20 trucks and many refrigerated trucks is preferred over a solution with 21 trucks and less refrigerated trucks. An example of three different optimal solutions (splitting-forbidden, splitting-allowed and minimize either the number of splits or the number of split products) is given in Fig. 1.
The methodological contribution of our paper is the development of an exact and heuristic solution approach. The exact approach is based on branch-and-price (BaP) Lübbecke and Desrosiers 2005). For different objectives and corresponding extended bin-packing problems, we develop a unified description of the respective column-generation subproblems. We show that, in turn, each subproblem can be solved by one or several modified 2-dimensional knapsack problems (Martello and Toth 2003;Kellerer et al. 2004, ch. 9.6). These are especially involved in the case of the splitting-allowed policy. Moreover, an important and here non-trivial algorithmic component is the branching scheme. We adopt the branching scheme originally developed and successfully applied by Heßler et al. (2018) to the vector-packing problem. Stabilization with dual-optimal inequalities (Ben Amor et al. 2006;Gschwind and Irnich 2016) helps to increase the performance of the BaP approach.
The heuristic approach works in two levels: In the first level, several constructive heuristics are employed to reach a packing of items from the same category. They strive for a balanced utilization of weight and volume constraints based on weightper-pallet, i.e., density considerations, and on the surrogate weights as introduced by Caprara and Toth (2001). In the second level, a more complicated meta-scheme considers the multi-compartment and product-fragmentation features of the problem at hand and applies various problem-specific heuristics to search for good solutions. They are based on decomposing and regrouping the best solution obtained from the first level to improve the lexicographic objectives for different product categories. At first, only the splitting-forbidden policy is considered. Then, for the splitting-allowed policy we try to improve the given solution by separating certain pallets of the same product and possibly reduce the number of trucks by reassigning these items. Throughout the heuristic process, we compare the current solutions to lower bounds as a stopping criterion. It should be noted that the heuristics were also incorporated in the SAP ERP software system of our industrial partner, a European producer of foodstuffs and beverages, and are used in the daily planning process.
We conducted a series of computational experiments with the exact and heuristic algorithms based on real-world data sets with up to 430 items and on more difficult test instances derived from these. We compare the performance of our heuristic to the previously employed strategy of the company and to the proven optimal solutions. It turns out that the heuristic saves on average one truck per day compared to the previously used approach for real-world instances, with an even higher advantage for more difficult instances. At the same time, running times are moderate, often less than one second. The gap to the optimal solution is fairly moderate for real-world instances but increases for the self-generated hard test cases. Regarding the performance of the exact branch-and-price algorithm, most of the test instances can be solved to optimality within the time limit of 10 min for every subproblem. Notably exceptions arise for some of the generated, difficult instances in case of the splitting-allowed policy. We can also illustrate the highly positive effect of strengthening the ILP-models by lower bounds and of solving pricing problems on a reduced graph. Finally, combining the two approaches it turns out that adding the best heuristic solution to the set of columns of the branch-and-price approach has hardly any effect for the splitting-forbidden policy, but it gives a major improvement for the splitting-allowed policy.
The remainder of the paper is organized as follows: In Sect. 2, a brief overview of the literature on two-dimensional bin-packing problems with and without item fragmentation and multi-compartment vehicle-routing problems is given. Analytically computed lower bounds are presented in Sect. 3. Details on the exact BaP algorithm are provided in Sect. 4, including the sequence of extensive formulations used for different objectives, the unified column-generation algorithm, and stabilization techniques. The heuristic approach is presented in Sect. 5. In Sect. 6, the results of our computational experiments are presented and discussed. Final conclusions are drawn in Sect. 7.

Literature review
The classical one-dimensional bin-packing problem (BP) is one of the most fundamental problems in combinatorial optimization: A given set of items has to be packed into a minimum number of bins such that the total weight of the items in each bin does not exceed the bin capacity. Classical heuristics for BP, such as first-fit (FF), best-fit (BF) and worst-fit (WF), and their variants based on a decreasing ordering of item weights, i.e., FFD, BFD, and WFD (see, e.g., Coffman et al. 2013), will be assumed common knowledge throughout this paper.
For the exact solution of BP, two alternative and competing solution principles can be observed in the recent literature (see Delorme et al. 2016, for an overview of BP formulations). The first type of solution approach relies on pseudo-polynomial arc-flow formulations (Valério de Carvalho 1999) solved with general-purpose mixedinteger linear programming (MIP) solvers. Important refinements are network/graph compression techniques (Brandão and Pedroso 2016) and the use of reflect networks modeling packings with only half of the bin capacity (Côté and Iori 2018). The first half and the second half of a packed bin are represented and feasibly combined in the reflect network. The second type of solution approach uses a pattern-based formulation a.k.a. the Gilmore-Gomory model (Gilmore and Gomory 1961), which is solved by column generation techniques. While early works focused on the solution of the linear relaxation in order to quickly compute a strong lower bound, the necessary theory providing a complete branching scheme leading to a fully fledged BaP algorithm has been provided with the work of Vanderbeck (1999). A state-of-the-art branchprice-and-cut approach for BP is presented by Wei et al. (2020b). A hybrid of both approaches, denoted as reflect+, has been recently coined by Delorme and Iori (2020) and seems to be a very competitive approach making use of the excellent dual bounds provided by column generation, heuristically and exactly truncated reflect networks, and the power of modern MIP solvers.
Due to the high practical relevance of packing problems, many variations in the classical one-dimensional BP have been investigated. For a comprehensive overview, we refer to the surveys by Christensen et al. (2017) and Coffman et al. (2013). In the context of our real-world bin-packing variant, especially publications dealing with two-dimensional, multi-compartment, and item fragmentation aspects are of particular relevance.
Two-dimensional BPs represent a well-studied generalization and, hence, a large number of publications in this field appeared over the last decades. We do not consider packing problems in which geometric items have to be placed into bins, known as rectangle packing (Lodi et al. 2002). Relevant for our application are p-dimensional BPs with p = 2 independent dimensions (such as physical weight and volume, or value and weight). In the literature, these problems are given different names, e.g., (2-dimensional) vector (bin) packing problem (VPP, 2D-VPP) and two-constraint binpacking problem. Our problem is a variant of the 2D-VPP for which an early exact approach is given by Spieksma (1994). Recent exact approaches are based on BaP and are described in Heßler et al. (2018) and Wei et al. (2020a). General approximation results are given in Bansal et al. (2016). The case, where weights and volume are strictly ordered, is considered in Caprara et al. (2003). Simple greedy-type heuristics are presented in Aringhieri et al. (2018). A comprehensive survey of vector packing problems can be found in Christensen et al. (2017).
Multi-compartment problems arise when different item types can be packed into the same bin but have to be separated from each other (incompatibility constraints). Those settings are mainly studied in vehicle-routing problems, see the surveys (Pollaris et al. 2014;Henke 2018). In our application, there are no incompatibilities between products/items. However, the use of compartments for cooled and frozen products imposes additional costs.
The introduction of item fragmentation generalizes the classical BP by allowing items to be split among several bins at a cost. Item fragmentation for BPs is introduced and proved NP-hard by Mandal et al. (1998). More recent publications provide exact solution algorithms based on column generation and dual cuts techniques as well as employing heuristics and implicit enumeration (Casazza and Ceselli 2016). LeCun et al. (2015) presents models for practical applications of item fragmentation like minimizing the number of money transfers after a group trip and presents approximation algorithms for identical as well as for the more generic case of non-equal-sized bins. The generalization to bins that differ in both cost and capacity is tackled in the contribution Casazza (2019) by proposing new mathematical programming models that even avoid the use of fractional variables. A worst-case analysis for BP with item fragmentation is performed by Bertazzi et al. (2019). General models of valuations for fragmented items are considered for the closely related knapsack problem in the recent work Malaguti et al. (2019). Recall from Sect. 1 that in our application, products might be divisible (depending on the splitting policy) but that by definition our items are never split. Our definition of an item has implicitly solved issues related to product fragmentation by discretely splitting demand along the volume (=pallets) dimension.

Lower bounds
Various lower bounds for subsets of items or all items are used in different parts of the algorithms. For convenience, we introduce these lower bounds in this brief section.
A straightforward lower bound on the number of bins needed to pack a subset I ⊆ I can be obtained by dividing the problem along the two dimensions. For each dimension, one can compute a bound on the respective one-dimensional bin-packing problems and then use the better bound of the two: Caprara (1998) proved that this lower bound equals the rounded up value of the linear relaxation to the standard assigment-type ILP model for the two-dimensional binpacking problem, see e.g.
(1)-(6) in Caprara and Toth (2001). Another lower bound is inspired by the lower bounding procedure of Martello and Toth (1990) for bin-packing problems. We use two distinct subsets of I, namely the subset of items consuming more than half of a bin's capacity, and the disjoint subset of items that do not fit together into a bin with any item of the subset I 1 . We obtain the lower bound LB 2 (I) = |I 1 | + LB 1 (I 2 ).
In the following, the lower bound LB(I) = max{L B 1 (I), L B 2 (I)} is used for various subsets I. For the sake of brevity, we define the following shorthand notation for the lower bounds LB = LB(I ),

Branch-and-price algorithm
We propose the exact solution of the overall lexicographic minimization problem using a stage-wise modified BaP algorithm Lübbecke and Desrosiers 2005). While the first objective is minimized by solving a standard 2D-VPP, the following four stages first restrict the solution space to previously computed optimal value(s) and then apply BaP to the restricted formulation. Since the initial formulation is of the Gilmore-Gomory type (Gilmore and Gomory 1961;Heßler et al. 2018), the linear programs can be characterized as extensive, meaning that they typically exhibit an enormous number of variables. BaP is tailored to solving such extensive optimization problems. Technically, BaP is a branch-and-bound algorithm in which the linear relaxation at each node of the search tree is solved by column generation (CG). CG is an iterative method that solves a linear program by decomposing it into a restricted master problem (RMP) and a pricing subproblem (SP). In Sect. 4.1, we start by describing the set of all feasible packing patterns and define some subsets of patterns which are later used for presenting the restricted formulations that address the objectives 2.-4., 5a., and 5b. at subsequent stages. The unifying approach for the pricing subproblems presented in Sect. 4.2 is followed by the discussion of stabilization, branching schemes, and acceleration techniques in Sects. 4.3-4.5.

Pattern-based formulations
The models of all stages are variations in the covering formulation of Gilmore and Gomory (1961) that are based on patterns. A pattern describes a feasible packing of a bin with items. For each item i ∈ I , the value is an upper bound on the number of times that item i may occur in a pattern. The set of all feasible patterns is therefore Recall that for the standard case of the splitting-forbidden policy, we have q i = 1 and thus u i = 1 so that all feasible patterns are binary, i.e., (a 1 , . . . , a m ) ∈ {0, 1} m . Thus, the resulting VPP is a binary VPP (01-VPP). For the splitting-allowed policy, the coefficients a i are upper bounded by integers u i so that the resulting VPP is a bounded VPP (B-VPP). For Stages 2.-5., pattern subsets taking different product types into account have to be considered. We use the superscripts S, C, and F to refer to standard, cooled, and frozen products/items, respectively. A combination of superscripts refers to patterns that only include items of the respective categories. If a superscript is supplemented with > 0, at least one item of the category must be included. Table 2 exactly defines the relevant pattern subsets.

Stage 1.: Minimize the total number of trucks
At Stage 1., the minimization of the number of trucks is equivalent to the minimization of the total number of bins that are used. Therefore, the problem is a standard 2D-VPP. The Gilmore-Gomory formulation comprises non-negative integer variables x p for all patterns p ∈ P and is given by Standard i∈I S a i Cooled i∈I C a i Frozen i∈I F a i P S = 0 a n d = 0 ≥ 1 a n d = 0 ≥ 1 a n d ≥ 1 a n d = 0 The objective (1a) minimizes the number of bins, i.e., patterns used. Constraints (1b) ensure that the demand of all items is covered. Note that for an RMP, i.e., the linear relaxation of a model with a restricted variable/pattern set P ⊂ P, we also show the associated dual variables (π i ) i∈I of constraints (1b) that we later use in Sect. 4.2 to describe the CG process. The domain of the pattern variables is defined by (1c). The first stage can be exactly solved with the BaP algorithm presented by Heßler et al. (2018) without any adaptation.

Stage 2.: Minimize the number of refrigerated trucks
The secondary objective of minimizing the number of refrigerated trucks is equivalent to minimizing the number of patterns of the subset P S,C∪F>0 (see Table 2). The first objective imposes the additional condition that the total number of all trucks/bins is restricted to z 1. . It can be incorporated into model (1) with a single additional constraint, leading to: The objective (2a) minimizes the number of refrigerated trucks. Constraint (2b) restricts the total number of trucks to the minimum number resulting from the solution of formulation (1). A lower bound on the number of refrigerated trucks is imposed by constraint (2c). Note that constraint (2c) is not mandatory for the correctness of the model, but its addition often yields a better lower bound of the linear relaxation.

Stage 3.: Minimize the number of refrigerated trucks which contain frozen products
The tertiary objective of minimizing the number of refrigerated trucks that contain frozen products is equivalent to minimizing the number of patterns of the subset P S,C,F>0 (see Table 2). At Stage 3., we additionally have to bound the number of refrigerated trucks by z 2. , which yields to the model The objective (3a) is the minimization of the number of trucks transporting frozen products. Constraint (3b) fixes the total number of refrigerated trucks. Since this fixation makes constraint (2c) redundant, it is omitted. We add the constraint (3c) imposing a lower bound on the number of frozen trucks to strengthen the linear relaxation.

Stage 4.: Minimize the number of refrigerated trucks which also contain standard products
The fourth rate objective of minimizing the number of refrigerated trucks which also contain standard products is equivalent to the minimization of patterns of the subset P S>0,C∪F>0 . Optimal values from the Stages 1.-3. are again incorporated as additional constraints: The objective (4a) minimizes the number of refrigerated trucks which also transport at least one standard product. In the following, such a truck is denoted as mixed truck.
A lower bound on the total number of mixed trucks is LB S − (z 1. − z 2. ). The difference z 1. − z 2. gives the number of trucks containing only standard products. Subtracting this number from the lower bound for the number of trucks required to pack all standard products gives a lower bound on the number of mixed trucks. Constraint (4b) restricts the number of frozen trucks. As constraint (3c), constraint (4c) is not mandatory but helps to strengthen the linear relaxation.

Stage 5.: Minimize splitting
Recall from Sect. 1 that the last objective of minimizing splitting applies only to the splitting-allowed policy and can be put into effect in two different ways using objective 5a. or 5b.
In the case of objective 5a., minimizing the overall number of splits is equivalent to minimizing the number of different products packed into all bins. To this end, let s p = {a p i > 0 : i ∈ I } denote this number of different products/items in a pattern a p i ∈ P. Then, the actual number of splits is given by the difference of the total number of different products in all used bins and the total number m of products. Therefore, minimizing the number of splits can be formulated as Constraint (5b) restricts the number of mixed trucks as determined at Stage 4. In the case of objective 5b., it only counts whether a product is split or not, but the number of trucks a split product occupies is irrelevant. Minimizing the number of split products is equivalent to maximizing the number of integrally packed products. A product is integrally packed if all its pallets are packed on the same truck. Accordingly, for a pattern p ∈ P, we define for our minimization problem the negative number of integrally packed products as d p = −| i ∈ I : a p i = q i |. Hence, our objective is expressed by the minimization of the total value of d p over all patterns used: Constraint (5d) ensures that each item i ∈ I is packed at most q i times. Note that the latter constraint is indeed necessary because otherwise one might artificially improve the objective function by adding additional integrally packed products more often than demanded. In combination with the covering constraint (1b), both constraints (5d) and (1b) are always binding. Therefore, adding disaggregated constraints of the form p∈P a p i x p = q i for all i ∈ I does not yield a tighter formulation.

Column generation
With the models (1), (2), (3), (4), (5a)-(5b), and (5c)-(5d), we have presented six different extensive formulations that all require a CG-based solution approach. In principle, pattern generation is algorithmically simple to manage, since we will show that the subproblem solution finally boils down to solving binary or bounded 2-dimensional knapsack problems, see (Martello and Toth 2003), (Kellerer et al. 2004, ch. 9.6). In the following, we assume that a linear relaxation of the RMP is solved and has produced dual prices (π i ) i∈I and (δ con ) for the respective constraint(s) con. The task of the SP is to determine at least one negative reduced-cost pattern or to prove that no negative reduced-cost pattern exists. A first difficulty for precisely describing the associated SPs lies in the fact that each of the six SPs associates different dual prices to groups of patterns depending on whether they include standard, cooled, or frozen products/items. Handling these multiple cases is indeed confusing. A fundamental step toward a concise and unified handling is the idea to further divide each SP into 2D-KPs for subsets of patterns. Solving the SP then amounts to solving the associated 2D-KPs. The decisive point is that the reduced cost of a pattern p = (a i ) can be described as (at least for the first four stages, see below)c where the coefficients δ = δ X and σ i = σ X i for i ∈ I depend of the particular subcase X . Table 3 formalizes different (sub)cases: Depending on the stage/objective, there are between one and five different pattern subsets (see Table 2 for their formal definition) for which independent 2D-KPs need to be solved. For each pattern subset, the columns σ i and δ describe how these coefficients can be computed from the given dual solution (π i ) i∈I and (δ con ) (using only dual prices of constraints con that exist at the actual stage). Moreover, the 2D-KP can be restricted to an item subset I ⊂ I . Finally, some pattern subsets require the presence of one or two subsets of items. For example, P S>0,C,F>0 requires that at least one item of a standard product and one item of a frozen product are present in the pattern. For this purpose, let R be the set of required subsets of items, i.e., for P S>0,C,F>0 there is R = {I S , I F }. The elements of the respective sets R are displayed in the last column of Table 3. A second difficulty stems from the fact that, for the objectives 5a. and 5b., the reduced cost of a pattern is no longer completely linear in the coefficients of the pattern. However, all 2D-KPs can be formalized as follows. Nonnegative bounded integer (binary for u i = 1) variables y i for i ∈ I describe the unknown pattern coefficients a i . Moreover, for these variables y = (y i ) i∈I , the function f (y) captures the non-linear part of the reduced cost (described in detail in the next paragraphs). For a given subset of admissible items I (in the sense of the second last column of Table 3), the extended 2D-KP can now be described as The objective (6a) is the minimization of the reduced cost of the pattern. The constraints (6b) and (6c) are the capacity constraints. The defining property of some pattern subsets to contain at least one item of one or two particular categories is enforced by condition(s) (6d). The domain of the decision variables is defined by (6e). At Stages 1.-4., the function f vanishes, i.e., f ≡ 0, so that KP(I, σ i , δ, R) reduces to a standard binary or bounded 2D-KP with the additional constraints (6d). A straightforward solution approach could be to solve several standard binary 2D-KP that previously pack one item of each set J ∈ R. Instead of enumerating all combinations, the selection of previously packed items could also be embedded into a branch-and-bound algorithm. However, we formulate and solve this subproblem as a variant of the shortest-path problem with resource constraints (see next paragraph).
For the first splitting variant that minimizes the total number of splits, i.e., objective 5a., the function f must count the number of products the constructed pattern contains. This is done with For the second splitting variant that minimizes the number of split products, i.e., objective 5b., the function f must count the (negative) number of integrally packed gives the correct contribution. As an example, the two 2D-KPs arising at Stage 2. are specified in the following example.
Example 2 At Stage 2., we consider two 2D-KPs that generate patterns of the subsets P S and P S,C∪F>0 . Recall that δ 2b and δ 2c are the dual prices of constraints (2b) and (2c), respectively. Both 2D-KPs have σ i = π i for all i ∈ I and f ≡ 0. The first 2D-KP generates patterns of set P S and is defined by δ = 1 − δ 2b , I = I S , and R = {I S }, while the second generates patterns of the set P S,C∪F>0 and is defined by SPPRC-based solution Different problems KP(I, σ i , δ, R) can be formulated as shortest path problems with resource constraints (SPPRCs, Irnich and Desaulniers 2005). SPPRCs can model intricate relationships between attributes and many variants of the SPPRC can be effectively solved by dynamic-programming labeling algorithms. Such an approach has been successfully chosen for the one-dimensional BP by Wei et al. (2020b) and for the VPP by Heßler et al. (2018). The advantage of a labelingbased approach is that also subsets R of required items and the objectives 5a. and 5b. defining the function f can be considered as well as involved branching decisions (discussed later in Sect. 4.4).
For the sake of simplicity, we start with the case that any item can be contained in a pattern, i.e., I = I = {1, 2, . . . , m}. The underlying digraph G = (V , E) then consists of m + 1 vertices V = {0, . . . , m}. The arc set E consists of m + i∈I u i = i∈I (u i + 1) arcs partitioned into m subsets E(i) associated with item i ∈ I . The set E(i) contains u i + 1 parallel arcs connecting vertex i − 1 with vertex i. Thus, we have an arc set E(i) consisting of arcs e j (i) with j ∈ {0, 1, . . . , u i + 1}.
Each arc e j (i) has three attributes, two for the weights and one for the profit. For the j-th arc e j (i) of E(i), these are defined as w(e j (i)) = j · w i , b(e j (i)) = j · b i , and σ (e j (i)) = j · σ i , respectively.
An example of the digraph is sketched in Fig. 2 Partial paths that do not fulfill these conditions can be directly discarded. If I I , i.e., some categories of items are not in I, the corresponding digraph can be obtained by reducing the arc set so that for each i ∈ I \I only the null-arc (0, 0, 0) is in E(i). Subsequently, the digraph can be shrunk by merging vertices i − 1 and i for all i ∈ I \I.
At Stages 2.-5., three additional boolean attributes v = (v S , v C , v F ) for the three categories standard, cooled, and frozen are added to indicate whether a specific product category has been packed. For an item i ∈ I , we define v(e j (i)) = (1, 0, 0), (0, 1, 0), or (0, 0, 1) for all j according to the category product/item i belongs to.
The SPPRC is solved by a forward dynamic-programming labeling algorithm that works as follows: The starting point is the trivial partial path (0) with the initial label (0, 0, 0). All labels at a vertex i − 1 are extended to vertex i itera- where the maximum operator is applied component-wise. The labels at m and the associated 0-m-paths have to be filtered at the end: Only labels (w, b, v, σ ) with correct v-values provide feasible patterns, e.g., v C + v F ≥ 1 for the pattern set P S,C∪F>0 used at Stage 2. in the second subcase. Moreover, the reduced cost of the associated pattern is δ − σ .
The work of Heßler et al. (2018) has pointed out that different dominance principles between labels are possible and worth to be investigated. There clearly exists a tradeoff between the strength of a dominance rule and its computational burden resulting from the pairwise comparison of labels. A strong dominance rule can help to drastically reduce the overall number of labels that remain at termination of the labeling algorithm. The result is also fewer labels at intermediate vertices, which reduces the computational effort associated with the label extension step and later dominance comparisons. On the downside, the application of the following strong dominance rule may require the pairwise comparison of a possibly huge number of labels.
The following weaker dominance rule allows a direct O(1) comparison of labels, immediately at the time when a label is created, if labels are stored in a lookup table.
As in Heßler et al. (2018), also our pretests have confirmed that the weak dominance Rule 2 performs better than the strong Rule 1 for most benchmark instance used in the computational experiments. Therefore, we do not use Rule 1 and also omit its proof. The validity of the weak dominance Rule 2 is obvious.
Finally, we discuss the incorporation of the functions f at Stage 5. into the SPPRC modeling and solution approach. For the first objective 5a., the function f counts the number of products the pattern contains, i.e., f ( y) = |{i ∈ I : y i > 0}|. Consequently, the profit attribute on the arcs of the SPPRC digraph must consider the number of different products in the resulting pattern. Thus, a cost of 1 is subtracted for every item that is packed, i.e., every non-zero arc. An example is given in Fig. 3.
For the second objective 5b., the function f counts the (negative) number of integrally packed products the pattern contains, i.e., f ( y) = −|{i ∈ I : y i = q i }|. Hence, arc profits are modified if an item i is packed completely, i.e., the additional profit of 1 is added to the q i th arc of item i (present only if q i ≤ u i ). An example is given in Fig. 4. Note that the dual price γ of constraint (5d) is already included in the definition of σ i for all i ∈ I .

Stabilization
The CG process can be stabilized by using dual inequalities (Valério de Carvalho 2005;Ben Amor et al. 2006;Gschwind and Irnich 2016). In the following, we restrict ourselves to inequalities that are formulated only in the dual variables π i that correspond to the demand fulfillment constraints (1b) of the model solved at Stage 1. (for the subsequent stages we consider the projection onto the m-dimensional space defined by the π -variables). Let D * be the set of dual optimal solutions to the linear relaxation. For t ∈ Z m and t ∈ Z, the dual inequality (DI) t π ≤ t is a dual-optimal inequality (DOI) if D * ⊆ {π : t π ≤ t}. A set of DIs is a set of deep dual-optimal inequalities (DDOIs) if at least one optimal dual solution π * ∈ D * fulfills all inequalities of this set.
We first focus on Stage 1. and the B-VPP and 01-VPP. According to Heßler et al. (2018), pair inequalities (PIs) π h ≥ π i are DDOIs for all h, i ∈ I with w h ≥ w i and b h ≥ b i . In the primal formulation (1), these additional columns stand for the exchange of the larger item h by the smaller item i, i.e., a zero-cost column in the RMP with entries −1 and 1 for h and i, respectively. Moreover, subset inequalities   3. and 4.
5b. I S , I C , I F 1 (SIs) −π h + i∈S π i ≤ 0 are defined for any item h ∈ I and subset S ⊆ I \{h} with i∈S w i ≤ w h and i∈S b i ≤ b h . We refer to a SI as (S, h). In general, SIs are not necessarily DOIs or DDOIs for the B-VPP and the 01-VPP. Nevertheless, their associated primal columns can be added to stabilize the RMP at the risk of a possible over-stabilization. Over-stabilization can be restored with a recovery procedure of low computational effort as described in detail in (Gschwind and Irnich 2016). At Stages 2.-4., the use of DIs is still possible when categories of items are kept separate, i.e., a SI for (S, h) can only be used if all items i ∈ S and also item h all belong to the same, stage-dependent subset of exchangeable items. An overview of the applicable DIs at different stages is provided in Table 4. For example, at Stage 2., the PIs and SIs exclusively exchange either items of standard products I S or of refrigerated products I C ∪ I F . At Stages 3. and 4., exchanges are only allowed among items of the same category.
At Stage 5., we have to take into account the following speciality. The transformation of a solution with packing and DI columns into a pure packing column solution may change the number of packed products in a pattern (to be taken into account at Stage 5a.). Likewise, an integrally packed product may be transformed into a split product and vice versa (to be taken into account at Stage 5b.). In order to ensure a valid objective value after the replacement, we have to modify the right-hand side of a DI from 0 to some integer value c, i.e., add a cost c to the corresponding primal DI column. Table 4 summarizes different cases. For example, let p ∈ P be a pattern and let i∈I t i π i ≤ c be a DI. If at Stage 5a. in the primal model both corresponding variables are positive, the following replacement can be performed: Exchange item h by items {i ∈ I : t i > 0} in pattern p with original associated cost s p (i.e., the number of packed products) yields a new pattern p with associated cost s p ≤ s p + |{i ∈ I : t i > 0}| − one(u h ), where one(u) is 1 for u = 1, and 0 otherwise. The new pattern p might consist of added items {i ∈ I : t i > 0} compared to pattern p, and if u h = 1, then item h vanishes in pattern p .
Similarly, at Stage 5b., the (negative) number of integrally packed items for the new pattern p can be estimated by d p ≤ d p + 1 because item h might change from an integrally packed item to a split item.
As discussed in Gschwind and Irnich (2016) in more detail, one may experiment with different strategies which DIs and at what point in time the DI columns are added to the RMP. On the one hand, DI columns can be added at the very beginning to the initial RMP and stay there during the entire CG process (static). On the other hand, only violated DIs can be added during the CG process (dynamic). We combine these two strategies which yields a mixed approach that adds some DI columns at the start and adds violated DIs later on.
We employ this mixed strategy in the following form: For the initialization of the first RMP, for each item i, the PI π i ≥ π j with j = arg min k∈I \{i} {|w i −w k |+|b i −b k | : w i ≥ w k , b i ≥ b k } is added and all SIs with |S| = 2 are added. While solving the root node, violated PIs are added dynamically. At Stages 2.-5., we only keep DIs belonging to the subset of exchangeable items according to Table 4 (and remove the others). In the same spirit, we only dynamically add violated PIs belonging to these subsets.

Branching
We apply a single or a two-level branching scheme depending on the stage. At Stage 1., there is only one level where we apply the branching scheme suggested by Vanderbeck (1999). The idea is to formulate the 2D-VPP as a pure binary problem so that all pattern coefficients are also binary. Branching is then performed choosing two disjoint subsets I 0 , I 1 of items. In any fractional solution there exists a pair (I 0 , I 1 ) for which the number of patterns having coefficients a i equal to 0 (1) for all i ∈ I 0 (∈ I 1 ) is fractional. For a detailed explanation and examples, we refer to Heßler et al. (2018).
At Stages 2.-4., there are two branching levels. First, we branch on the number of refrigerated trucks (2c), the number of frozen trucks (3c), and the number of mixed trucks (4c), respectively. Second, we apply the above described branching rules of Vanderbeck (1999).
At Stage 5., there are also two branching levels. First, we branch on the number of integrally packed products, i.e., when the value is fractional for an item i ∈ I . Note that this is a special case of Vanderbeck branching, which we apply subsequently at the second level. At all levels, the branching variable/term is selected as one with fractional part closest to 0.5. Ties are resolved randomly. Note that all first-level and second-level branching decisions can be implemented by either adding some linear constraint(s) to the RMP or by removing certain pattern subsets. Since the branching rule of Vanderbeck (1999) ensures integrality, the overall branching scheme is complete for all stages.
Finally, we use a depth-first tree exploration strategy. The intention is to find an integer solution as soon as possible. Since the lower bounds of the linear relaxations of formulations (1a)-(5) are often tight, we quickly obtain good and sometimes optimal solutions.
Reduced SPPRC digraph for objective 1. with otherwise identical attributes as depicted in Figure 2.
The arcs e ∈ E are shown with their weights and profit (w(e), b(e), σ (e))

Acceleration techniques
When applying the splitting-allowed policy, patterns with fewer split products are preferred over patterns with more split products. Therefore, it is beneficial to already generate patterns with very few split products (if any) at the Stages 1.-4. so that the initial RMP at Stage 5. mainly consists of variables representing patterns with many integrally packed products. To foster that predominantly patterns with only integrally packed products are generated, the SP is always first solved heuristically over a reduced SPPRC digraph. In this reduced digraph, each item i ∈ I has only two associated parallel arcs (0, 0, 0) and (u i w i , u i b i , u i π i ) so that for u i = q i the item is either integrally packed or not packed at all. An example of a reduced graph is given in Fig. 5. Only if no pattern with negative reduced cost is found in the reduced digraph, SPPRC SPs are solved again using the complete digraph. This approach has two advantages: First, the 2D-KPs at the Stages 1.-4. are solved faster because in most iterations the reduced network provides negative reduced cost patterns. Secondly, Stage 5. is solved faster because the RMP already consists of beneficial columns. Computational results show that this approach significantly accelerates the CG process. For a computational analysis, we refer to Sect. 6.4.

Heuristics
For real-world applications, user acceptance is a crucial success factor, especially for constructive heuristics, where the key users want to get some insights into the decision mechanism. Therefore, we developed various heuristic ideas and presented their underlying rationals as well as their strengths and weaknesses to the decision makers of the affected company. For comparison, we also re-implemented the strategy previously employed by the dispatchers, the so-called Business Today logic.
For our own development, we started with the well-known BFD, FFD, and WFD heuristics, all of them employing the concept of surrogate weight (Sect. 5.1). Then five additional heuristic approaches were developed in accordance to the real-world situation of our industrial partner (Sect. 5.2). Comparisons to results of the Business Today logic are found in Sect. 6.

Basic single-class heuristics
At the beginning, we will introduce a number of heuristics which consider only a single class of items and do not distinguish between products of different categories. Also splitting products is not taken into account. At the end of this section, we will have obtained one combined heuristic which consists of running all heuristics described so far and taking the best solution. This combined heuristic will serve as building block for a more complicated heuristic framework described in Sect. 5.4.
As pointed out in the introduction, we are facing a highly heterogeneous product range where the unit weight ρ k of product k ∈ K varies considerably. It is obvious that truck loads containing mainly heavy goods (ρ k W /B) will leave excess pallet spaces even if the weight capacity is exhausted. On the other hand, truck loads containing mainly lightweight goods (ρ k W /B) still provide excess weight capacity even if all pallet spaces are used. Thus, a "good" loading pattern should utilize both capacity dimensions by mixing heavy and lightweight products in a suitable way.
Based on this simple observation, the transport planners previously created loading patterns filling one truck after the other in a FF manner. Each truck is manually assigned its load following a strict rule that alternately packs the heaviest and the lightest unpacked item as described in Business Today Algorithm 1.

Algorithm 1 Business Today
Sort item set I by increasing ρ i . Open truck 1 and set t ← 1. while I = ∅ do Take the last (first) item i from I in an odd (even) iteration. if i fits into truck t then Insert i into truck t. else Open a new truck t + 1 and insert i.
Practical experience showed: If each product fills only one pallet, i.e., b i = 1 for all i ∈ I , then this procedure amounts to a pairwise combination of items with large and small unit weight ρ i and shows reasonably good results. However, this is a rare case, and for all other, more general instances, inefficient solutions may occur. In the worst case, many large items (high w i and b i ) with medium density ρ i remain until the end, which leads to unnecessary additional trucks.
Naturally, Algorithm 1 suffers from neglecting the absolute weights of items and from the FF strategy of handling trucks. However, sorting items by decreasing weight w i or decreasing number of pallets b i would only be useful for instances with a clear dominance in one dimension. Our real-world test instances showed that both constraints, weight and pallet capacity, can become active, i.e., there is no clear dominance of one of the two dimensions.
To overcome this dilemma, we suggest to apply the concept of a surrogate weight proposed by Caprara and Toth (2001). The surrogate weight of an item i ∈ I constitutes a convex combination of relative weight and relative pallet number. For each i ∈ I , it is defined as Based on these surrogate weights, we can (almost) transform our packing problem into a classical one-dimensional bin-packing problem and apply well-known standard heuristics. In particular, we take FFD, BFD, and WFD, and adapt them for the solution of our truck loading problem by using the surrogate weight s i . In the following, we present further heuristic approaches based on the concept of a surrogate weight.

Block-based single-class heuristics
In this section, we introduce three heuristics called A, B, and C, which are all based on a fixed pattern of assigning blocks of items to trucks. Therefore, the lower bound of the number of trucks LB (see Sect. 3) is computed and LB trucks are opened. Then, the item set is partitioned into blocks of equal size LB. In each iteration, one block is selected and matched to the set of trucks, i.e., every item of the block is assigned to one of the trucks. The selection of blocks and the matching patterns differ between the three heuristics. If an item does not fit into the preferred truck, the item is moved to the next truck in the sequence (following a FF logic) or, if necessary, a new truck is opened, see Algorithm 2. Note that any new trucks opened during the execution of the heuristic only serve as a backup for loading items exceeding the capacity, while the matching of subsets remains restricted to the originally opened LB trucks.
All three heuristics are based on sorting the items by increasing surrogate weights s i . Their details are informally described below followed by pseudocode. An illustrative example is given in Table 5 further below.
Heuristic A assigns the LB items with largest surrogate weights to the LB trucks and then adds the LB items with smallest surrogate weights in a reversed order, such that the items with s max and s min end up in the same truck number 1. These two steps are repeated iteratively for the remaining items as described in Algorithm 3.
Heuristic B also starts by assigning the LB items with largest surrogate weights to the LB trucks, but then proceeds by assigning the next largest (w.r.t. s i ) LB items to the trucks in a mirrored order such that items with LBth and (LB + 1)-largest surrogate weight end up being assigned to the same truck number LB, see Algorithm 4.
Finally, Heuristic C repeats the first step of Heuristic A and iteratively assigns the items with largest surrogate weights to the LB trucks, always maintaining the same order of trucks, see Algorithm 5.

Algorithm 2 Subprocedure Assign item i to truck t
Try to insert i into an open truck by a FF strategy, i.e., try all trucks t, t + 1, . . . , LB and then trucks 1, 2, . . . , t − 1 until item i fits. If no such truck is found, insert i into a truck LB + 1, LB + 2, . . . by a FF strategy, opening a new truck if necessary. To illustrate the different heuristics, an example with 20 items is given in Table 5. Obviously, Heuristic A exhibits more similarities to the Business Today logic, including the feature of not keeping the items with lowest s i until the end. Heuristics B and C overcome that issue. Besides that, Heuristic B aims at uniformly filling all trucks. According to the basic intuition of how good solutions might look like, this should reduce the probability that any smaller item does not fit into its preferred truck. On the other hand, Heuristic C provides gaps of different sizes in the LB different trucks, which may offer more room to maneuver if an item does not fit into its assigned truck.

Density-based single-class heuristics
The third group of heuristics for items of the same class is based on the unit weight ρ considerations discussed at the beginning of Sect. 5.1. Recall that "good" loading patterns, in general, utilize both capacity dimensions and do not leave large excess   capacity of pallet space if weight capacity is used to the limit, and vice versa. Therefore, we can define a desired average density for an ideal load of a truck, namely the ratio optavg ρ = W /B. Clearly, a loading pattern P t of a truck t with an average density avgr t := j∈P t w j / j∈P t b j close to optavg ρ and total weight close to the capacity W also fills the available pallet space close to the limit B, and vice versa.
Based on this simple observation, we introduce two heuristics D and E. In contrast to heuristics A, B, and C, they both consider trucks one after the other and for each considered truck t they add the item yielding a new average density avgr t as close as possible to the target value optavg ρ . Of course, the item also has to fit into the truck at hand. Ties are broken by choosing the item with the largest weight, because smaller items can be expected to fit more easily in later iterations.
Heuristics D and E differ from each other only by the weight measure for the initial sorting step. Note that the sorting only determines the first item put into every truck. Heuristic D sorts by the surrogate weight s i ensuring that "large" items are assigned first and will not give rise to trucks with low utilization toward the end of the packing procedure. On the contrary, heuristic E focuses on outliers with an unusual density that are potentially hard to assign. Therefore, items are sorted in decreasing order of unit weight ρ i . By the same rationale also the reverse ordering could be a reasonable choice, but for our test data, pretests showed that this ordering is less effective.

Algorithm 6 Heuristic D (E)
Sort items I by decreasing s i (ρ i ). Let P t ← ∅ be the loading of all trucks t.
Take the first i from I and insert i into P t . Set avgr t = ρ i and remove i from I . while items exist that fit into P t do i ← arg min i∈I Insert i into P t and remove i from I . end while end while Heuristics D and E could be immediately extended to consider all trucks instead of only the current truck t for inserting items. In that case, the computation of i would have to be extended to taking the arg min over all currently open trucks, which means that an item is inserted into a truck where the minimal deviation from optavg ρ can be realized. Indeed, we also implemented this extension with the initialization that the first LB items (according to the given ordering) are packed on LB trucks. However, since improvements were only marginal and occurred only for a small subset of instances, we decided to stick to the simpler standard version described above.

Multi-class heuristic scheme of splitting-forbidden policy
All heuristics described in the previous sections focus on creating good patterns utilizing both capacity dimensions, but they take neither items of different product categories packed into different compartments of a truck nor item fragmentation into account, while the latter is considered in the following Sect. 5.5. Accordingly, we now present a heuristic scheme which is based on the single-class heuristics presented before. These are embedded into a more complicated algorithmic framework consisting of various decomposition and re-merging steps. Since the running times of the singleclass heuristics are rather low, we perform the task of solving a single-class problem for a certain subset of item I ⊆ I as a subroutine by running all nine heuristics from Sects. 5.1, 5.2, and 5.3 and taking the best solution found. More precisely, we run the Business Today algorithm, FFD, BFD, WFD, and Heuristics A to E. The result (= number of trucks) derived by this combined single-class heuristic is denoted by z heu (I). Additionally, this combined heuristic is executed with a given starting solution of prepacked trucks, i.e., some trucks have reduced residual capacities. We refrain from listing all the details of adapting the single-class heuristics to this situation.
In Sect. 3, the lower bounds LB, LB S , LB C,F , and LB F for the corresponding subsets of items are introduced. Applying a single-class heuristic only for standard items and then separately only for refrigerated items yielding heuristic solution values z S heu and z C,F heu , respectively, we can conclude the following: If z 1. sep := z S heu + z C,F heu = L B, then we have reached an optimal solution w.r.t. the objectives of Stages 1. and 2. with z 1. = z 1. sep and z 4. = 0. Otherwise, if z S heu + z C,F heu > L B, we cannot guarantee optimality of the merged solution without mixed trucks. In this case, it makes sense to gradually introduce mixed trucks while trying to reduce the total number of trucks. As already used in Sect. 4, constraint (4c), the lower bound for mixed trucks (z 4. ) is given by LB mixed = LB S + LB C,F − LB for solutions that are optimal w.r.t. Stages 1. to 3. Therefore, LB mixed = 0, if LB S + LB C,F = LB and LB mixed > 0, otherwise.
The heuristic scheme striving for a good solution of the lexicographic optimization problem defined by Stages 1.-4. with splitting forbidden works in three phases which are described in Algorithm 7. Phase 1 builds up a solution starting with frozen products. After packing I F by the combined single-class heuristic (following the Stage 4. objective), the cooled items are added which might require additional trucks (Stage 3. objective). Then a separate solution for the standard item set I S is determined. Taking the union of the two sets of trucks gives a first solution, where standard and refrigerated items are packed in separate trucks (Stage 2. objective). If the resulting number of trucks z best equals the overall lower bound LB we stop (STOP 1), since we have found a solution which is optimal w.r.t. Stages 1. and 2.
Otherwise, we try to improve the currently best solution in Phase 2: Therefore, we allow mixed trucks and pack standard products also into refrigerated trucks. To reach a better starting position for adding also large standard items to the refrigerated trucks given from Phase 1, but keeping the number of mixed trucks small, the load of the refrigerated trucks is regrouped at the beginning of Phase 2 with the goal of obtaining some refrigerated trucks with a low load and others which are almost fully packed with cooled products. This is done by subprocedure Regroup, which tries to empty the least loaded refrigerated truck by moving its items (in decreasing order of surrogate weight) to other refrigerated trucks selected by a BF strategy. If the number of frozen products is small, we try to keep them together in the regrouping in the spirit of the Stage 3. objective. Therefore, all frozen products contained in the current truck are considered as a single product and moved to a new truck together. Only if this artificial product does not fit, the frozen products are reassigned individually. This regrouping is applied to all trucks in an increasing order of surrogate load. Then the standard items I S are added to the regrouped refrigerated trucks by the combined single-class heuristic. In this way, we aim to find a packing with a lower number of trucks (Stage 1.), but without increasing the number of refrigerated trucks (Stage 2.). Again, we stop if the new solution matches the lower bound LB I (STOP 1).
A further reduction in the total number of trucks is sought in Phase 3 at the cost of increasing the number of refrigerated trucks. Again we try to improve the possibility of adding standard products to refrigerated trucks. Therefore, we take the refrigerated truck with lowest load w.r.t. the surrogate weight and split its content on two empty trucks, thereby incrementing the number of refrigerated trucks by one. The splitting is performed in a balanced way by assigning the items by a WFD strategy (w.r.t. surrogate weights). Again, all frozen items are treated as one product, if the truck contains both frozen and cooled products. Then the standard items are again added by the single-class heuristic. If the total number of trucks reaches the lower bound LB we stop the procedure. Otherwise, we unpack all standard items again and iterate, taking the refrigerated truck with lowest load that was not split before. If all refrigerated trucks given at the beginning of Phase 3 were split, another iteration is started permitting a second split operation for each refrigerated truck given at the end of the first iteration. Further iterations along this rule are possible, until at some point the number of refrigerated trucks exceeds the lower bound for all items (STOP 2). In this case, we also run the single-class heuristic for all items I to reach a possible improvement in the Stage 1. objective.

Multi-class heuristic scheme of splitting-allowed policy
If the multi-class heuristic of Sect. 5.4 reaches a solution where all lower bounds for the objectives of Stages 1.-4. are reached, obviously the result cannot be improved by splitting products. However, if this is not the case, item fragmentation offers potential for improvements.
In the following, we consider the two different splitting objectives introduced in the Introduction. z current ← z refrig for all z current refrigerated trucks t opened so far in increasing order of load w.r.t. surrogate weights do Open a new (empty) truck z refrig + 1 z refrig ← z refrig + 1 Take all (refrigerated) items currently loaded on truck t and pack them on the two trucks t and z refrig by a WFD strategy. {gives a balanced bipartition} Compute z S heu starting from the current packing of the z refrig refrigerated trucks. z best ← min{z best , z S heu } if z best = LB then STOP 1. end if if z refrig ≥ LB then Compute z heu {for empty trucks} z best ← min{z best , z heu } STOP 2. end if Remove all standard items I S from the trucks end for until false Subprocedure Regroup(P 1 , . . . , P T ) Sort the T trucks in increasing order of load w.r.t. surrogate weights for t = 1 to T do for all items i in P t in decreasing order of surrogate weights do Try to pack item i in trucks {P t+1 , . . . , P T } by a BF strategy (w.r.t. surrogate weights) end for end for Algorithm 8 Multi-class heuristic: Stage 5a., minimizing the number of splits Compute LB 1 Execute Algorithm Multi-class heuristic returning T := z best trucks with loads P 1 , . . . , P T Sort the T trucks in increasing order of load w.r.t. surrogate weights t ← 1 repeat for all items i in P t in decreasing order of surrogate weights do Try to pack item i into another truck by a BF strategy (w.r.t. surrogate weights) end for while P t = ∅ do Let i be the item in P t with largest surrogate weight {find the truck t where the largest number of pallets of i can be added}

cannot be emptied}
Return the solution with z best trucks (as given at the start of the current iteration of repeat)

STOP. end if
Move h pallets of product i from P t to P t end while Remove P t z best ← z best − 1 if z best = LB I 1 then Return current solution with z best trucks STOP. end if t ← t + 1 until false Algorithm 9 Multi-class heuristic: Stage 5b., minimizing the number of split products In this setting, it is not a-priori clear which items should be split and in which way to gain the largest improvement in the solution. We have chosen to incorporate item splitting as a post-processing step applied to a solution derived from the multi-class heuristic of Sect. 5.4. In every iteration of the main loop of our Algorithm 8, we try to empty the truck with lowest surrogate load and thus reduce the total number of trucks by one. If this turns out to be impossible, the procedure is stopped. Otherwise we proceed to the next iteration, unless we have reached a solution with LB 1 trucks, which proves optimality of objective 1. and also lets us stop the procedure.
Every iteration for the truck t with minimal surrogate weight consists of two parts: At first (for-loop) we try to further reduce the load of truck t by moving items (without splitting them) to other trucks by a BF strategy. Then (in the while-loop) we iteratively take the largest (w.r.t. surrogate weight) remaining item i in truck t and split it, i.e., we pack the largest possible fractional part of it, expressed as number of pallets, to one of the other trucks. Formally, the largest pallet number h is calculated, such that h pallets out of the b i pallets constituting i can be loaded on some truck t . Ties among trucks are broken by a BF rule. The remaining part of item i consisting of b i − h pallets is treated like an individual item in t. By moving the largest possible number of pallets, it can be expected that the total number of splits remains small. If no single pallet can be moved from truck t (i.e., h = 0), the whole procedure is stopped. Otherwise, truck t is completely empty and can be removed.

Stage 5b.: Minimizing the number of split products
In this case, it does not matter whether a large or small item i is split and whether i is split into two fragments or the maximum number of b i fragments. It is intuitively clear that the largest potential for improvement is given by splitting large products completely into separate pallets. Therefore, in our Algorithm 9 we iterate over all items in decreasing order of b i . Every considered item is split into separate pallets and then the multi-class heuristic from Sect. 5.4 is executed. The iterative process is stopped in any of the following three cases: (1) the multi-class heuristic stops with a solution matching all lower bounds for the objectives of Stages 1.-4. Since we consider split products, we can only utilize the volume based lower bound LB 1 ; (2) all products consisting of more than one pallet were split; (3) all products were considered and split.

Computational results
In this section we will report computational results for all algorithms described before. The experiments are based on real-world instances and on more difficult instances generated from these by reducing their data to more difficult item sets (see Sect. 6.2). We first report in Sect. 6.3 results for the heuristic approaches described above in Sect. 5. To analyze the performance of the algorithm we compare the heuristic with the exact solutions in detail. Then we will illustrate the performance of the exact BaP algorithm in Sect. 6.4. Furthermore, we also analyze in that section the impact of adding the heuristic solutions as columns to the initial RMP.

Details of the implementation
All heuristic algorithms presented in Sect. 5 are implemented and tested in Python 3.3.7 on a standard PC equipped with an Intel ® Core™ i5-3210M processor with 2.5 GHz and 4 GB of RAM. Note that the heuristics are also used in the daily planning task of our industrial partner. For this application they were originally implemented within the SAP ERP system of the company in the SAP-internal programming language ABAP. However, for test purposes we re-implemented the heuristics in Python.
The BaP algorithm described in Sect. 4 is implemented in C++ and compiled in release mode into a single-thread code under MS Visual Studio 2015. The experiments are conducted on a standard PC with an Intel ® Core™ i7-5930k processor clocked at 3.5 GHz and 64 GB of RAM. The RMP is solved at each column-generation iteration by means of CPLEX 12.9. Moreover, CPLEX is called after the solution of each branchand-bound node as a primal MIP-based heuristic solver using the so far generated columns. CPLEX's default values are kept for all parameters. The time limit to solve each stage is set to 10 min (600 s). In case a stage cannot be solved within the time limit, the computation is stopped and following stages are not considered.

Benchmark instances
We consider 80 real-world instances with 35 to 430 items provided by our partner from the European food and beverage industry. The cardinality bound of each truck is 33 pallets (a common industry standard), and the weight capacity is 22.8 or 24.5 tons. The weight of each item is rounded up with an accuracy of 10 kg. Preliminary tests show that both constraints, weight and cardinality, can become active, i.e., there is no clear dominance of one of the two dimensions. Not all instances contain all types of products; especially, frozen products often appear in rather small quantities or they may even be completely absent. A detailed summary of our test instances is found in Table 6.
During the computational study it turned out that for most of the instances the optimal solutions for Stages 1.-4. do not require any splitting, and thus, the objective of Stage 5. would be optimized by default. The reason for this behavior is that the instances also contain rather small items which fill up trucks without being split. Therefore, we distilled more difficult new instances from the given real-world instances avoiding this effect in the following way. The 80 real-world instances are divided into 10 groups each with 8 instances. For each group, we selected 5, 10 or 15 items with the highest weight and/or the highest number of used pallets of each instance forming a new instance, respectively. In total, we generated 30 new instances with up to 80, 160  or 240 items, respectively. All instances are available on the website http://logistik. bwl.uni-mainz.de/benchmarks.php.

Results of the heuristic approach
In this section, we report results of the heuristic approach. Table 7 shows the key results for the pure heuristic approach. For 56 out of the 80 real-world instances, all lower bounds for objectives 1.-4. are already achieved without splitting. The 24 remaining instances provide at least a theoretical improvement potential, which is exploited in 18 cases. For 13 of them, even the lower bounds for objectives 1.-4. are achieved if items are allowed to be split. While for the real-world instances exactly 70 % could be solved to the lower bounds for objectives 1.-4., this was possible only for 10 % (3 out of 30) of the generated instances, all without splitting. However, the share of solutions that could be improved by applying splitting is far higher (25 out of 27) and the same is valid for the number of instances where that improvements even led to the lower bounds for objectives 1.-4. This higher improvement is not surprising, as the average pallet size of items is more than 75 % larger for generated instances, and thus, the effect of splitting items into separate pallets is much stronger. The higher difficulty of generated instances is also reflected in the running times of the heuristics.
For none of the 110 instances, the original Business Today logic outperformed the heuristic approach, but for 89 instances, the heuristics performed better. For the other 21 (all real-world) instances, both the heuristics and the Business Today logic reached the lower bounds.
For a further performance analysis, we compare the objective values of the heuristic approach with the optimal objective values computed by the BaP approach. Detailed results are found in Table 8. The table entries have the following meaning: #opt: number of instances solved to proven optimality at the specific stage; z H = z * : number of instances in which the heuristic approach finds the optimal solution; instances with non-optimal heuristic solutions in former stages are not considered; z H > z * : number of instances in which the heuristic approach does not find the optimal solution; instances with non-optimal heuristic solutions in former stages are not considered; gap H : average gap between heuristic and exact solution, we use the gap z H − z * instead of the relative gap because the optimal solution can be z * = 0 in Stages 4.-5.; recall that solution values for Stages 1.-4. represent number of trucks, but split parameters for Stage 5.; #implicit H : number of instances for which the heuristic and optimal solution implicitly coincide, or where the current stage is not active (e.g., no Stage 3. for instances without any frozen items); -z H = z * : total number of instances in which the heuristic approach finds the optimal solution up to the specific stage: ( -z H = z * ) := (z H = z * ) + (#implicit H ); -#opt: number of instances solved to proven optimality at the current stage; in contrast to #opt, implicitly optimally solved instances are also included.
For a better understanding, we describe the entries in the second row of Table 8 in detail. At Stage 2., the branch-and-price algorithm solves 64 instances to prove optimality. For 60 (3) instances, the heuristic approach could (not) find the optimal solution. For one instance, we cannot compare the heuristic and exact approach at Stage 2. because the solution at Stage 1. differs for this instance. For 6 instances without any refrigerated items, the heuristic approach finds the optimal solution at Stage 1. so that Stage 2. is solved implicitly. Including instances without refrigerated products, the heuristic and exact solution coincides for 66 instances in total.
For the splitting-forbidden policy, the heuristic solution is optimal for 58 + 6 of 68 + 18 instances. Interestingly, in most of the cases, the heuristic approach fails to find the optimal solution at the first stage. Especially for the generated instances, only 15 of 30 instances are solved to optimality. For the remaining 12 + 12 instances, where no optimal solutions are achieved by the BaP algorithm either, the heuristic approach fails to reach the lower bounds in 7 + 9 cases at the first stage. However, the gaps are rather small. Only for 0 + 2 instances, absolute gaps larger than 1 truck occur. For the splitting-allowed policy, the heuristic approach performs well with 67 + 24 of 67 + 25 instances solved optimally up to Stage 4. The absolute deviance gap H is low up to Stage 4. with a maximum of 0.2, but the values are big for Stage 5. With 6/6 + 4/5 of 62/63 + 7/8 non-optimal instances, the heuristics fail for most of the instances at Stage 5.

Results of the branch-and-price algorithm
In this section, we report the results of our BaP algorithm. In general, reaching a stage of the lexicographic optimization task triggers the execution of a BaP algorithm solving the respective model (1)-(4), and possibly (5a)-(5b) or (5c)-(5d). However, for some instances, not all stages become active, e.g., if there are no frozen products in the instance or if a stage is solved implicitly during the preceding stage. Therefore, the total number of instances per stage to be reported in the computational experiments may be different for every stage. Moreover, the different BaP algorithms of the stages may have different success rates and thus will lead to a different number of instances remaining for the successive stages. This makes the comparison of different algorithmic approaches a delicate task. The following analyses may, therefore, seem somewhat less transparent, but our focus is on the correct interpretation of the results we obtained. The construction of the initial RMP works as follows. As the first stage is a vector packing problem without additional constraints, we start with the same initial RMP as described in (Heßler et al. 2018) by adding unit vectors, columns with only one non-negative coefficient a p i = u i , and 200 additional columns resulting from variants of the FF and BF heuristics. At Stage 2., we add other 200 columns resulting from the FF and BF heuristics that first assign items to bins that already contain the same product type. Moreover, the same FF and BF variants used in Stage 1. are performed on restricted item sets by considering standard or refrigerated products separately. Analogously, Stage 3. applies the FF and BF variants for cooled or frozen products separately. At Stage 4.-5., no additional columns are added initially to the RMPs.
At all stages, after solving the linear relaxation of a branch-and-bound node, CPLEX is invoked as a primal heuristic using the integer model and all patterns that have been generated up to this point. If the integrality gap is closed with the CPLEX solution, the branch-and-price is immediately terminated (providing a provably optimal solution). To allow for a fair comparison of computation times, CPLEX is run in single-thread mode and its computation times are counted for the overall computation time and time limit.
For each stage, the computation time is restricted to a maximum of 10 min (600 s). If a stage cannot be solved within the time limit, we cannot solve the instance exactly. Therefore, the computation is stopped at this stage and the following stages are not considered. If an instance has no frozen or no refrigerated products at all, Stage 3. or Stages 2.-4. are skipped, respectively.
At first, the BaP algorithm is tested employing the results of the heuristic approach, i.e., the corresponding columns are added to the initial RMP. We refer to this variant as base. To analyze the impact of the heuristic solution, we compare the base variant to a variant without using the heuristic solution. Moreover, the impact of the lower bounds enforced by constraints (2c), (3c), and (4c), and the graph reduction presented in Sect. 4.5 are tested. We refer to these variants as without heuristic solution, without lb, and without reduced graph, respectively. For the base(long) variant, the computation time is extended to a maximum of 20 min (1200 s). The results are summarized in Table 9. The table entries have the following meaning: Class: class of instances; real-world or self-generated; for details see Sect. 6.2; Splitting: allowed or forbidden; in case of splitting-allowed, either the number of splits (objective 5a.) or the number of split products (objective 5b.) is minimized, see Sect. 4.1; #inst: number of active instances; #inst_cons: number of active instances considered at the current stage; instances that are not solved to proven optimality in a previous stage are not considered; #opt: number of instances (out of #inst_cons) solved to proven optimality at the specific stage within 10 min (600 s); #implicit: number of instances for which the current stage is already implicitly solved by the result of a previous stage, or where the current stage is not active (e.g., at Stage 3. for instances without any frozen items); -#opt: number of instances solved to proven optimality at the current stage: -#opt:= #opt + #implicit T : average computation time in seconds over #inst_cons instances; unsolved instances are taken into account with the time limit of 10 min (600 s); Table 9 Results of the exact BaP-based approach The largest #opt and -#opt entries are bold, respectively (not considering base ( Allowed T L P : average computation time of the linear relaxation in seconds; unsolved linear relaxations are taken into account with the time limit of 10 min (600 s); Gap: average gap between the upper and lower bounds; we use the gap UB − LB instead of the relative gap because the objective value can be zero at Stages 4.-5.
We summarize and interpret the results as follows: For the splitting-forbidden policy, 68 of 80 real-world and 18 of 30 self-generated instances are solved to optimality. In particular, the algorithm performs well at Stage 3. with an average computation time of around 1 min. In comparison, the splitting-allowed policy is more difficult than splitting-forbidden. The performance at Stage 1. is similar for both policies with 101 solved instances (splitting-allowed) instead of 106 (splitting-forbidden). Up to Stage 4., 67 of the 80 real-world and 25 of the 30 generated instances are solved which is (−1) + 3 instances more compared to the splitting-forbidden policy. The performance of the algorithm at Stage 5., where 62/63 + 7/8 instances are solved, is almost the same for both objectives 5a. and 5b. with similar average solution time. As can be expected, Stage 5. is particularly difficult for the generated instances. Especially, the average deviation (Gap) at Stage 5a. is very high because in many cases the root node cannot be solved or the solution of the preceding stage is poor.
When the base algorithm is stopped after solving the root node, only four instances less are solved to optimality. This results from the strong linear relaxation of the model in all stages, the tight upper bounds computed by the heuristic approach, and the effectiveness of the MIP solver.
We also analyze the solution quality for instances not solved within the time limit: For the splitting-forbidden policy, the base algorithm cannot solve the root node of the branch-and-bound tree in 2 cases. Moreover, the absolute gap UB − LB is 1 for all non-solved instances at Stages 1. and 4., while the average gap is 12.9 for all nonsolved instances at Stage 2. Similar observations can be made for the splitting-allowed policy: At Stages 1. and 4., the absolute gap UB − LB is almost always 1, while at Stages 2. and 5., the average gap is 8 and > 50, respectively. For the splitting-allowed policy, approximately half of the non-solved instances (18 instances) are so difficult for our approach that the root node remained unsolved in one of the stages.
To analyze the impact of the heuristic solution, we compare the variants base and without heuristic solution. For the splitting-forbidden policy, the results hardly change when the initial columns from the heuristic are omitted. In total, only one instance less can be solved and the average running times differ by less than 5 s. Therefore, the use of the heuristic solution provides only a small advantage for the splitting-forbidden policy. In contrast, for the splitting-allowed policy, the effect of adding the heuristic solution is more significant, as it helps to solve 6/7+2/3 additional instances, respectively.
Adding the additional constraints (2c), (3c), and (4c) to ensure lower bounds on the number of refrigerated, frozen, and mixed trucks, respectively, has a strong impact on the performance of the BaP algorithm. Without adding these lower bounds, only 45+13 instances of the splitting-forbidden policy and 38/39+3 instances of the splittingallowed policy can be solved. In total, around one-third of the instances cannot be solved without these bounds. Moreover, the running time increases significantly for all stages and problem variants.
For the splitting-allowed policy, using a reduced graph as described in Sect. 4.5 has a positive impact for the Stages 1.-4. Already at Stage 1., 10+9 instances less are solved if the graph reduction is omitted. Overall, 8/9+1/2 instances less are solved.
When the maximal computation time is extended to 20 min (instead of 10 min), the number of solved instances slightly increases. In total, 68+20 instances of the splitting-forbidden policy and 63+8 instances of the splitting-allowed policy are solved to optimality which is four instances more compared to the time limit of 10 min.

Conclusion
Packing pallets of products from industrial producers into trucks for shipment to wholesalers is an important daily task in supply chain management. Reducing the number of trucks, even if only by one, or rearranging the package layout to diminish the number of the more costly trucks carrying cooled or frozen products, has a significant impact on costs. It also helps to reduce exhaust emissions of trucks, since the supply operation is usually performed daily or several times a week. Thus, the effort spent on optimizing the packing operation yields savings that will add up over a year to numbers of considerable economic relevance.
In this contribution, we consider a real-word packing problem where products should be packed into trucks such that a five-level lexicographic objective function is minimized. Feasibility is defined by a weight and a volume (=number of pallets) constraint for every truck. This special variant of a multi-objective optimization problem considers the total number of required trucks as main objective dominating all other goals due to the high cost of labor. However, after fixing the overall number of trucks (and thus the number of required personal) secondary objectives come into play, which concern the costly operations of cooling and freezing devices of a truck. Finally, for practical handling it is clearly convenient for the customer to receive all pallets carrying the same product in the same truck. But if the total number of trucks could be reduced, the customer will agree to a splitting of products onto several trucks. Nevertheless, such a splitting should happen as rarely as possible, which opens the question how to measure the inconvenience of splitting. We believe that such a lexicographic setting has not be considered for one-or two-dimensional bin-packing problems before.
The first part of this paper gives a branch-and-price framework for computing exact solutions of our problem. While the upper-level problem of minimizing the number of trucks can be modeled as a two-dimensional vector packing problem, with a special structure in one dimension, the lexicographic objective requires five iterative steps of successive optimization problems, in each of which columns are generated by solving appropriate instances of a shortest path problem with resource constraints (SPPRC). It is a major contribution of the paper to put the five different problems arising for different levels into a uniform framework of column generation. Furthermore, advanced concepts of stabilization and some acceleration techniques are employed.
In the second part, we describe the heuristic solution method previously employed in practice and then develop a number of new constructive heuristics for solving a single level of the packing problem. These try to balance the two constraints and generate solutions with a favorable mix of weight and volume. Then, the single-level heuristics are integrated into a more complex heuristic framework for building solutions of high quality with respect to the multi-level lexicographic objective function.
Computational experiments for real-world benchmark instances as well as structurally more difficult instances extracted from these show that the branch-and-price algorithm is highly effective in computing optimal solutions. Within 10 min of computation time it reaches proven optimality for 68 out of 80 real-world scenarios without splitting. Allowing splitting leads to more difficult problems, but still 62 of 80 are solved to optimality.
It also turns out that the heuristic framework performs much better than the approach previously applied in practice and matches the optimal solutions in 86 % of all subproblems. Note that it is also very helpful to use the heuristic solution as a starting column for the branch-and-price algorithm. Based on this evaluation, our heuristic framework has been successfully incorporated into the SAP ERP system of our industrial partner and is currently used for the daily planning of direct shipments.
In the future, it would be interesting to extend our single-source single-sink shipping problem to a more complex supply network serving multiple customers. This would open up the possibility of delivering less-than-truckload amounts to one customer and using the residual capacity of a truck for serving a second (or third) customer, if the routing costs support such a decision.