Advanced loading constraints for 3D vehicle routing problems

Given automated order systems, detailed characteristics of items and vehicles enable the detailed planning of deliveries including more efficient and safer loading of distribution vehicles. Many vehicle routing approaches ignore complex loading constraints. This paper focuses on the comprehensive evaluation of loading constraints in the context of combined Capacitated Vehicle Routing Problem and 3D Loading (3L-CVRP) and its extension with time windows (3L-VRPTW). To the best of our knowledge, this paper considers the currently largest number of loading constraints meeting real-world requirements and reducing unnecessary loading efforts for both problem variants. We introduce an approach for the load bearing strength of items ensuring a realistic load distribution between items. Moreover, we provide a new variant for the robust stability constraint enabling better performance and higher stability. In addition, we consider axle weights of vehicles to prevent overloaded axles for the first time for the 3L-VRPTW. Additionally, the reachability of items, balanced loading and manual unloading of items are taken into account. All loading constraints are implemented in a deepest-bottom-left-fill algorithm, which is embedded in an outer adaptive large neighbourhood search tackling the Vehicle Routing Problem. A new set of 600 instances is created, published and used to evaluate all loading constraints in terms of solution quality and performance. The efficiency of the hybrid algorithm is evaluated by three well-known instance sets. We outperform the benchmarks for most instance sets from the literature. Detailed results and the implementation of loading constraints are published online.


Introduction
In recent years, sales in online trading have risen steadily. Forecasts for the coming years predict significant growth. Therefore, efficient logistics operations are more important than ever. Through many years of research in the field of Vehicle Routing Problems (VRP), (near-) optimal tour plans can be found for many use cases. Hereby, the demand of a customer is often simplified by using a total mass or volume for the items to be delivered. In practice, solutions might be infeasible since a vehicle cannot be feasibly packed because of unbalanced loading and/ or unsafe placement of items. As more and more information on items becomes available for detailed planning, the realistic planning of transportation and of packing processes could become the key factor for cost reduction and safety, leading to an increasing interest in combined routing and loading problems.
The combined problem at hand is the Three-Dimensional Loading Capacitated Vehicle Routing Problem (3L-CVRP). It was first introduced by Gendreau et al. (2006) and assumes delivery of cuboid items laying at the depot. A homogeneous fleet of vehicles is available for transporting the items to a number of customers. Each vehicle must be equipped with a feasible packing plan considering several loading constraints. The depot and the customers have specific time windows, in which the delivery must take place. This problem variant is known as the Three-Dimensional Loading Vehicle Routing Problem with Time Windows (3L-VRPTW).
The focus of this paper is on the comprehensive examination of loading constraints. Although the consideration of different complex loading constraints leads to more realistic models, it is mainly neglected in work dealing with the 3L-CVRP problem or its variants so far. The reason is that modelling and evaluation of loading constraints are complex and require new solution approaches. We tackle this problem and integrate the current largest constraint set so far. We introduce a new variant for the robust stability constraint, which increases the stability and the performance. Moreover, we develop an approach based on the science of statics so that for the first time, the acting load on an item is distributed through the entire stack, which ensures realistic and stable packing plans. In addition, this paper considers aspects for manual unloading, reachability, the axle weights of vehicles and a balanced loading. For the latter, we introduce formulas to illustrate our implementation approach. In case of manual unloading, the items are unloaded without lifting. The reachability constraint avoids unnecessary rearrangements of items. Detailed modelling of axle weights and balanced loading prevents overloaded axles and tipping over of vehicles. The implementation of all loading constraints is published online within a solution validator written in C++ as well as in Java. The validator can be used to check the feasibility of solutions for different loading constraint sets.
All constraints are integrated in a hybrid algorithm. The hybrid algorithm consists of an inner deepest-bottom-left-fill algorithm which solves the Loading Problem and is embedded in an outer adaptive large neighbourhood search tackling the Vehicle Routing Problem. The efficiency of the hybrid algorithm is 1 3 Advanced loading constraints for 3D vehicle routing problems shown by using the instance sets by Ceschia et al. (2013), Moura and Oliveira (2009) and Zhang et al. (2017). Experiments show that the hybrid algorithm performs better than the benchmark for most instance sets.
Moreover, we have created and published an instance set consisting of 600 new instances varying systematically in the number of customers, of item types and of items. For the first time, every complex loading constraint is evaluated concerning its impact on the objective values (number of used vehicles and total travel distance) grouped by number of item types, items and customers. Our evaluations consist of over 30,000 results, and we provide all results online and in detail (e.g. routing and packing plans with the position of all items) to ensure extraordinary transparency. On this basis, we give recommendations about which constraints are reasonable based on their impact on algorithmic performance and solution quality.
The paper is organized as follows. In Sect. 2, the relevant literature is reviewed. The 3L-CVRP and the 3L-VRPTW are formulated in Sect. 3. In Sect. 4, the new definitions and the implementation variants of the loading constraints are presented. In Sect. 5, the hybrid algorithm is described, and Sect. 6 deals with the testing of the constraints. Finally, conclusions are drawn in Sect. 7.

Literature review
This paper considers the Three-Dimensional Loading Capacitated Vehicle Routing Problem (3L-CVRP) and its extension with Time Windows (3L-VRPTW), which represent a combination of the Vehicle Routing Problem (VRP) and 3D Loading constraints. As shown in Table 2, the consideration of multiple loading constraints is currently sparely researched. Thus, this paper examines the impact of different loading constraints on the results for the 3L-CVRP and the 3L-VRPTW. The current state of modelling loading constraints for both problems is analysed in the following. Gendreau et al. (2006) introduced the combined Vehicle Routing and 3D Loading Problem, namely 3L-CVRP. They solve the VRP using an "outer" tabu search, which determines customer sequences. An iteratively invoked "inner" tabu search defines the item sequence for the routes. The loading algorithms are based on the touching parameter algorithm by Lodi et al. (1999) and the bottom-left-algorithm by Baker et al. (1980). The items are packed orthogonally into the vehicle loading space (orthogonality constraint) without overlapping and respecting their dimensions (geometry constraint). The rotation of the items is only allowed along the width-length plane (rotation constraint). Each item has a mass, and the vehicle has a maximum capacity (load capacity). Moreover, a fragility flag is assigned to each item to prevent stacking fragile items on top of each other (fragility constraint). When stacking items, they must be supported by other items with a certain percentage (minimal supporting area constraint). When unloading items, it should be done by direct movements parallel to the length of the vehicle (LIFO constraint). Since 1 3 the constraints orthogonality, geometry, rotation, load capacity, fragility, minimal supporting area and LIFO are commonly considered in researches on the 3L-CVRP and its variants, this set is here defined as basic constraint set. For testing, Gendreau et al. (2006) developed 27 instances.

3L-CVRP
The 3L-CVRP has been studied intensively in recent years so that the results for this benchmark have been improved repeatedly (e.g. Tarantilis et al. 2009;Fuellerer et al. 2010;Bortfeldt 2012 andWei et al. (2014)). Tarantilis et al. (2009) used a combination of tabu search and guided local search to build the routes. For the Loading Problem, successively six packing heuristics are called until a feasible solution is found. They also present a new variant-the Capacitated Vehicle Routing Problem with Manual 3D Loading Constraints (M3L-CVRP). This variant deals with the manual handling of items, e.g. the items are small and of low mass. Therefore, the LIFO constraint is modified so that it is allowed that one item hangs over another one. This adaption of the LIFO policy, which is in this paper referred to as MLIFO, is also examined in a paper by Ceschia et al. (2013). Ceschia et al. propose a local search approach combining simulated annealing and large neighbourhood search to solve the VRP. To handle the Loading Problem, one out of nine loading heuristics based on the bottom-left-algorithm and the touching perimeter algorithm is selected. Besides the MLIFO constraint, they consider the reachability of an item for the first time within the 3L-CVRP. In the context of the Three-Dimensional Bin Packing Problem, this constraint was developed by Junqueira et al. (2013) to avoid the driver standing on items to reach other items for unloading or arranging operations. Ceschia et al. (2013) also include the item's load bearing strength (lbs), which was first mentioned by Bischoff and Ratcliff (1995) and examined in Bischoff (2003) for the Three-Dimensional Bin Packing Problem. Thus, to the best of our knowledge, Ceschia et al. currently combine the most loading constraints. Krebs and Ehmke (2021) consider detailed modelling of axle weights of vehicles for the first time for the 3L-CVRP.

3L-VRPTW
In Moura (2008) and Moura and Oliveira (2009), the VRTWLP is introduced, which corresponds to the 3L-VRPTW without the consideration of masses and stacking constraints (e.g. fragility and load capacity) and with higher stability requirements (full support) and with more rotation possibilities. Moura (2008) proposes a multiobjective genetic algorithm to generate routes (VRP). If a customer is inserted in a route, a wall-building heuristic is called to tackle the Loading Problem. This packing heuristic is also used in Moura and Oliveira (2009), where a hierarchical and a sequential approach are combined. The hierarchical one solves primarily the VRP, while the sequential one handles the VRPTW and the bin packing. 46 instances are created. The current best-known results for these instances are received by Reil et al. (2018), who solve the packing problem through a tabu search algorithm. Then, a multi-start evolutionary search minimizes the number of used vehicles while another tabu search algorithm minimizes the total travel distance. Pace et al. (2015) propose a heuristic based on simulated annealing and an iterated local search for the routing phase. Since they examine the distribution of fibre boards, a specialized loading heuristic based on a depth-first tree search and a balanced loading constraint are necessary. The latter is also adopted by Mak-Hau et al. (2018), who develop a mixedinteger linear programme model of the 3L-VRPTW with a heterogeneous fleet. Zhang et al. (2017) solve the 3L-VRPTW with a hybrid approach, consisting of a new loading heuristic and a routing heuristic based on a tabu search and an artificial bee colony algorithm. They include the basic constraint set and combine the two well-known instance sets provided by Gendreau et al. (2006) andSolomon (1987).
In this paper, we use the approach by Koch et al. (2018) proposed for the 3L-VRPTW with Backhauls, which is also used for the 3L-CVRP in Krebs and Ehmke (2021). The following Table 1 summarizes the approaches. Table 2 summarizes the related literature and highlights our contribution. As demonstrated in Table 2, this paper deals with the largest constraints set and combines the robust stability (C6b), load bearing strength (C7b) and reachability (C8) with axle weights (C9) and the balanced loading (C10) constraints. Moreover, we distribute the loads for the first time through the entire stack in the load bearing strength constraint.

Problem formulation
Following the convention by Koch et al. (2018), the 3L-VRPTW is described as follows: Let G = (N, E) be a complete, directed graph, where N is the set of n+1 nodes including the depot (node 0) and n customers to be served (node 1 to n), and E is the edge set connecting each pair of nodes. Each edge e i,j ∈ E (i ≠ j, i, j = 0, ..., n) has an associated routing distance d i,j (d i,j > 0) . The demand of customer i ∈ N ⧵ {0} consists of c i cuboid items. Let m be the total number of all demanded items. Moreover, time windows are considered by assigning three times to each node i: the ready time RT i , which is the earliest possible start time of service, the due date DD i , the latest possible start time, and the service time ST i , which specifies the needed time to (un-) load all c i items of a customer i.
Each item I i,k (k = 1, ..., c i ) is defined by mass m i,k , length l i,k , width w i,k and height h i,k . The items are delivered by at most v max available, homogenous vehicles. Each vehicle has a maximum load capacity D and a cuboid loading space defined by length L, width W and height H. It is assumed that each vehicle has a constant speed of 1 distance unit per time unit. If a vehicle arrives at an edge before its ready time, it has to wait until the ready time is reached.
Let v used be the number of used vehicles in a solution. A solution is a set of v used pairs of routes R v and packing plans PP v , whereby the route R v (v = 1, ..., v used ) is an ordered sequence of at least one customer and PP v is a packing plan containing the position within the loading space for each item included in the route.
A solution is feasible if (S1) All routes R v and packing plans PP v are feasible (see below); (S2) Each customer is visited exactly once;  A route R v must meet the following routing constraints: (R1) Each route starts and terminates at the depot and visits at least one customer; (R2) The vehicle does not arrive after the due date DD i of any location i.
Each packing plan must obey a loading set P defining a subset of the following loading constraints, which are described in detail in next Sect. 4.
(C1) Geometry: The items must be packed within the vehicle without overlapping; (C2) Orthogonality: The items can only be placed orthogonally inside a vehicle; (C3) Rotation: The items can be rotated 90 • only on the width-length plane; (C4) Load capacity: The sum of masses of all included items of a vehicle does not exceed the maximum load capacity D. (C5a) LIFO: No item is placed above or in front of item I i,k , which belongs to a customer served after customer i; (C5b) MLIFO: No item is placed on or in front of item I i,k , which belongs to a customer served after customer i; (C6a) Minimal supporting area: Each item has a supporting area of at least a percentage of its base area; (C6b) Robust stability: Each item has a supporting area of at least a percentage of its base area at any height; (C7a) Fragility: No non-fragile items are placed on top of fragile items; (C7b) Load bearing strength: The load bearing strength lbs i,k is the maximal load per area unit an item can bear. It must not be exceeded anywhere on the top face of an item; (C8) Reachability: The distance between an item and the driver must be less or equal than a certain length ; (C9) Axle weights: The loads for the front and the rear axle do not exceed the permissible axle weights FA perm and RA perm ; (C10) Balanced loading: The load of one vehicle half does not exceed a certain percentage p of D.
The 3L-CVRP and 3L-VRPTW aim at determining a feasible solution minimizing the objective values, e.g. number of used vehicles v used and the total travel distance ttd, and meeting all corresponding constraints.

3
Advanced loading constraints for 3D vehicle routing problems

Definitions and implementations of loading constraints
This section discusses the implementation details and challenges of the considered loading constraints. We introduce new realizations and implementation variants. Detailed algorithms are provided and explained in a solution validator, written in Java and C++, available via http:// github. com/ Corin naKre bs/ Solut ionVa lidat or. Table 3 gives an overview of the considered loading constraints. The loading constraints C1-C4, C5a, C6a and C7a are used as described in Gendreau et al. (2006). In Table 3, we have highlighted new developed loading constraints in bold and constraints examined for the first time for the 3L-VRPTW in italics.

Unloading sequence (C5)
The unloading sequence constraints define the order in which the items of the customers of one route should be unloaded. The purpose is to prevent costly reloading processes of items during the unloading process. In the following, two definitions, namely LIFO (C5a) and MLIFO (C5b), are shown.

LIFO (C5a)
As shown in Gendreau et al. (2006), the last-in first-out (LIFO) constraint treats the unloading sequence in the way that all items c i of a customer i are loaded and unloaded by movements parallel to the front-rear axis (x-axis) of the vehicle without moving other items. Forklifts are mostly used for this purpose, which may need to lift an item during the unloading process (cf. Ceschia et al. 2013). Therefore, no item demanded by a customer that is delivered later can be placed over I i,k or between I i,k and the rear of the vehicle.

MLIFO (C5b)
In the Manual LIFO constraint (MLIFO) introduced by Tarantilis et al. (2009), the items are (un-)loaded by manual operations without the usage of, e.g. forklifts. Consequently, the items can be (un-)loaded without lifting them. Therefore, an item demanded by a customer that is served later than customer i can hang over the item I i,k without touching its surface and without being placed between I i,k and the rear of the vehicle. The differences between the LIFO and the MLIFO constraint are visualized in Fig. 1. For both variants, it is not allowed to place an item directly on top of another item that is delivered earlier (see Fig. 1a). In contrast to the LIFO constraint, it is allowed that one item hangs over another item that is delivered earlier (see Fig. 1b).

Vertical stability (C6)
The vertical stability constraints prevent stacked items from falling on the ground. For this purpose, we show that the current definition is not sufficient and formulate the new robust stability constraint.

Minimal supporting area (C6a)
The minimal supporting area constraint ensures that a certain ratio of the base of a stacked item is supported by the upper surface of the directly underlying items (see Gendreau et al. 2006). As shown in Fig. 2, this formulation can lead to unstable, but still feasible item arrangements: When stacking several items with same This arrangement is in accordance with the minimal supporting area constraint (C6a) because for the calculation of the support for one item, only the directly underlying items are considered. According to the science of statics, this stack is not stable, because the x-value of the centre of gravity (CG) lays outside of the dimensions of the first item. Therefore, the stack would topple.

Robust stability (C6b)
As shown above, the minimal supporting area can lead to unstable stacks, since in the calculation of the item's support only the directly underlying items are considered. Therefore, we formulate the robust stability constraint as follows: For each item, the relative support of at least a percentage of needs to be guaranteed at any height from the vehicle ground to the item's bottom edge.
Multiple overhanging (C6b1): This constraint was first introduced by Ceschia et al. (2013). As the name suggests, all items of a stack are allowed to overhang. When placing an item, the minimal supporting area is checked for all underlying items: Let U be the set which includes all placed items supporting directly or indirectly the item I i,k . An item I u supports I i,k directly if the top area of item I u has direct contact with the base area of item I i,k . An item I u supports I i,k indirectly if I u directly supports any placed item which directly supports I i,k . Each coordinate for the top surface of item I a ∈ U defines a plane. Another item I b ∈ U counts to this plane if the top surface of I b is at the same level as of the plane (see items I 1,2 and I 1,3 in Fig. 3b) or if the top surface of I b is above the plane and the base area of I b is below the plane (see item I 1,3 in Fig. 3c). Each plane must obey the minimal supporting Top overhanging (C6b2): In this paper, we want to introduce another variant for the robust stability, namely the "top overhanging" constraint. In contrast to the previous approach, here, only the topmost item of a stack is allowed to hang over other items. Hence, all items of a stack must be completely supported by other items except the topmost item, which can hang over considering the minimal supporting area constraint (C6a) (see Fig. 4b). This is appropriate for high stability requirements.
Top overhanging is implemented in the following way: Let distance ceiling be the distance between the topmost item of a stack and the ceiling (see Fig. 4a). Let h min be the smallest height of any unplaced item I min of the route and I i,k be the item which should be placed on top of the stack. When stacking items, two cases can occur: 1. If distance ceiling + h i,k ≥ h min , then item I i,k as well as I min can be placed on the stack. In this case, the item I i,k must be fully supported, since I min could be placed on top of I i,k , so that I i,k is not the topmost item of the stack. 2. If distance ceiling + h i,k < h min , then no unplaced item can be stacked on top of the stack. Thus, the item I i,k is the topmost item and must therefore obey the minimal supporting area constraint (C6a).

Stacking (C7)
The stacking constraints focus on the ability of items to bear other items. In the following, different approaches are shown. The fragility constraint (C7a) as shown in Gendreau et al. (2006) is the standard approach. The load bearing strength constraint (C7b) is proposed by Bischoff (2003) for the Container Loading Problem, where each item has an additional parameter indicating the maximum load it can bear. For this load bearing strength constraint, two implementation variants are described below. The first one (C7b1) is proposed by Bischoff (2003), while another approach is developed and introduced in this paper and is based on the science of statics (C7b2).

Fragility (C7a)
As shown in Gendreau et al. (2006), a fragility flag f i,k is assigned to each item to divide them into fragile items (f i,k = 1) and non-fragile ones (f i,k = 0) . On top of a fragile item, only another fragile item can be stacked, whereas both fragile and nonfragile items can be stacked on a non-fragile item. As demonstrated in Ceschia et al. (2013), the fragility constraint (C7a) has weaknesses: It is supposed that a non-fragile item lies mostly on another non-fragile item and a very small part on a fragile one (see Fig. 5). Even if the non-fragile part on top of the fragile item would be infinitely small, the arrangement remains infeasible.

Load bearing strength (C7b)
To handle the issue described before, the actual load on the items should be considered. Therefore, the load bearing strength (LBS) constraint is introduced: Each item I i,k can support a maximum load per area described by the parameter lbs i,k . It must not be exceeded anywhere on the top face of an item. A small lbs i,k value corresponds to fragile items.
If an item I c is stacked on top of another item I u , then, a load caused by I c acts on the underlying item I u ( load c,u ). For its calculation, the percentage of support for item I c ( support c ) provided by all directly underlying items must be first determined. Then, all area units ( supportArea c,u ) between Item I c and I u must be identified.
Based on that, the support share of I u on I c is given as follows: Since the item I c could overhang, but the load must be distributed in total, the support share is increased proportionally: The load acting on item I u is given as follows: (1) support c,u = supportArea c,u l c ⋅ w c .
(2) support prop = support c,u support c . where load c is the load which has to be distributed due to item I c . Its value is explained below. When placing an item on top of another, then the load must be distributed to underlying items. There are two ways to select these items: the simplified and the complete selection.
Simplified selection (C7b1): The approach proposed by Bischoff (2003) selects all items which are underneath the base area (e.g. the footprint) of an item I c . When placing an item I c on top of a stack, then all items which are underneath the base area and which directly or indirectly support item I c are considered. In this case, not all items of the stack may contribute to the mass distribution (see item I 1,3 in Fig. 6). In this approach, load c in Eq. 3 corresponds to the mass of I c . The example in Fig. 6 shows the resulting loads for the underlying items caused only by item I 1,6 .
Complete selection (C7b2): The following approach is based on the science of statics. When placing an item I c on top of other items, all items are investigated that are located directly below item I c . Therefore, the mass of I c is distributed as load c to the directly underlying items. Then, for each of these items, the received load c is further adopted and distributed to the directly underlying items again. This is recursively repeated until the items on the ground are reached. Fig. 7 shows the same exemplary situation as Fig. 6. In this approach, all items of a stack contribute to the mass distribution. The resulting loads caused by item I 1,4 and item I 1,5 are calculated in the same way.

Reachability (C8)
When an item is (un-)loaded, then it should be guaranteed that the working equipment or the driver can reach the item when standing as close as possible to the item (cf. Junqueira et al. 2013). For this purpose, the distance r i,k of an item I i,k should be equal or less than a certain length , which represent the driver's arm length, for example.
In this paper, for the reachability of an item I i,k , all items of customers, which are served after customer i and placed above or beneath item I i,k , are considered (see Fig. 8a). The distance r i,k is defined by the front of the item which is the closest to the door (MaxFront) and the front of item I i,k .
If the distance is larger than and thus the item is not reachable, then it is tried to shift the item along the x-axis. This is achieved by searching for the maximum x-value of already placed items on the same layer (MaxShift). The new placement must obey the DBL policy. Therefore, the item I i,k is shifted until the reachability constraint is just fulfilled, which means the new distance is defined by MaxFront − (see Fig. 8b). Additionally, the new placement is tested w.r.t. the loading constraint set P. If the item is not reachable and it cannot be shifted, the placement is rejected.

Axle weights (C9)
The exceedance of the maximum axle weights of one or more axles leads to farreaching consequences with regard to vehicle safety: It increases the braking distance and, in the event of a collision, the consequences are more severe due to the increased impact energy. Therefore, the axle weights constraint respects the permissible axle weights for the front and rear axles of a vehicle. Let FA perm be the maximum load the vehicle's front axle can bear and RA perm be the maximum load for the rear axle, respectively. Both limits are given in mass units. Let L f be the length between the front axle and the loading space (see Fig. 9). The wheelbase WB is the distance between the front and the rear axle. For each placed item I i,k at the x-position x i,k , the distance s i,k between the mass centre of I i,k and the front axle must be determined.
According to the approach by Krebs and Ehmke (2021), the following formulas can be applied for a vehicle v to calculate the acting forces for the front F FA and the rear F RA axle. Hereby, g is the constant for acceleration of gravity ( g ≈ 9.81 m s 2 ). and The acting forces must be below the permissible ones, but also greater than zero to avoid uplifting: and (4)  As demonstrated in Krebs and Ehmke (2021), the constraint must be checked after each placement of an item since, an axle may become overloaded after unloading items.

Balanced loading (C10)
To prevent a reduction in vehicle stability, the load per vehicle half should be examined. Therefore, Pace et al. (2015) suggest that a percentage p of the vehicle capacity D is not exceeded. In the following, we introduce formulas for this approach.
In our implementation, the item's mass m i,k is assigned to the vehicle sides depending on its y-position ( y i,k ). If an item lays entirely on the left side of the vehicle (see Fig. 10a), its total mass is assigned to the left vehicle side. The same is true for the opposite right side (see Fig. 10b). Otherwise, the mass of the item is distributed proportionally to the vehicle sides (see Fig. 10c). The sum of all assigned masses must not exceed a certain percentage p of the load capacity D. Consequently, the following must apply: Formula 11 restricts the range to positive real values. It is used to assign the masses to the corresponding vehicle sides. The constraint for the left vehicle side is checked in Formula 12 and in 13 for the right one, respectively. Inside of the large square brackets, the position of the item with respect to the corresponding vehicle side is determined in order to assign the proportional mass.

Hybrid solution approach
We propose a hybrid solution approach consisting of a routing heuristic (adaptive large neighbourhood search) for creating routes and an embedded packing heuristic (deepest-bottom-left-fill algorithm), which optimizes the loading of the items of all customers of a route into the loading space of a vehicle. The packing heuristic generates feasible packing plans for the generated routes. This packing plan is created following a loading constraint set P, which determines the included loading constraints.

Routing heuristic
We use the routing algorithm as described in Koch et al. (2018), who modified the adaptive large neighbourhood search (ALNS) proposed by Ropke and Pisinger (2006). The algorithm by Koch et al. (2018) was developed for the 3L-VRPTW with Backhauls and is applied to the pure 3L-VRPTW in this paper, considering additional constraints. The general framework is shown in Alg. Advanced loading constraints for 3D vehicle routing problems

Initial solution
The initial solution s init is constructed [1] with the savings heuristic developed by Clarke and Wright (1964). Hereby, all routing (except S3) and loading constraints are obeyed. Based on this feasible initial set of routes, the ALNS determines other feasible improved solutions.

Iteration
In each iteration of the ALNS, one removal rem and one insertion operator inst are randomly chosen [5][6]. These are used to generate the next solution s next by removing a number of customers n rem from the solution and reinserting them again [8]. The number of customers to be removed n rem ( n min ≤ n rem ≤ n max ) is determined randomly [7]. Then, it is checked whether the generated solution meets the routing constraints [9] described in Sect. 3. The packing procedure shown in the next subsection is called here.

Evaluation function
In order to evaluate different solutions and to lead the search, the following internal evaluation function is defined. The evaluation function f for a solution s giving total routing costs is described as follows: where N miss is a set containing all customers that have not been dispatched yet, v max is the maximal number of available vehicles and v used the number of used vehicles. Each customer i, which is not yet dispatched ( i ∈ N miss ), is assigned to one vehicle (round-trip) even if this leads to an exceedence of the number of used vehicles. The penalty term pen v is used to achieve a reduction of used vehicles v used . Additionally, the total travel distance ttd(s) for a solution s is respected.

Solution acceptance
A solution is regarded better the smaller its evaluation function value is. A better and feasible solution is always accepted. A worse solution may be accepted [9] according to an acceptance probability which depends on a simulated annealing heuristic proposed by Kirkpatrick et al. (1983). In particular, the acceptance probability is adapted to the annealing process with a geometric cooling schedule. The best solution s best is updated [13] if it has a superior evaluation function value relative to the current solution s curr [12].

Inserts customers iteratively so that an increase of routing costs is minimal
Regret-2 Inserts customers iteratively so that the maximal difference of routing costs for the best and the second best insertion in different tours is achieved Regret-3 Inserts customers iteratively so that the sum of two differences of routing costs is maximal. The first difference is the routing cost for the best and the second best insertion in different tours, while the second difference results from the best and the third best insertion in different tours Table 4 shows nine removal operators and Table 5 summarises the three insertion approaches used in this paper. We use the removal and insertion operators as described and evaluated in Koch et al. (2018). After a defined number of iterations iter p , the selection probabilities for the removal and insertion operators are adjusted [16][17][18] according to their improvement of the solution. This is described in detail in the following section.

Operator selection and probability adaption
The selection of the operators is accomplished by means of the roulette wheel selection principle. Hereby, the probability to select one operator op is defined by their weighting wg op . Initially, all operators have the same selection probability ( wg op = 1).

3
Advanced loading constraints for 3D vehicle routing problems The number of iterations is counted, in which the operator op -is selected ( counter op ), -is selected and led to a new best solution ( iter best op ), -is selected and improved the current solution ( iter impr op ), -is selected and led to a worse but not yet accepted solution or a solution as good as the current solution ( iter e op ).
After a certain number of iterations iter p , the success of the operator is evaluated and described by score op , which is calculated as follows: Hereby, best , impr and e are coefficients. Then, the new weighting wg op can be calculated. A reaction factor r regulates the influence of the adaptions: Moreover, counter op , iter best op , iter impr op and iter e op are reset to zero.

Stopping criteria
If one of the following stopping criteria is met [19], the heuristic terminates, and the current best known solution is given: -Number of total iterations iter max ; -Number of iterations without improvement iter wimpr ; -Calculation time limit t max .

Packing heuristic
As packing heuristic, we use the same approach as in Krebs and Ehmke (2021), which is based on the deepest-bottom-left-fill (DBLF) algorithm proposed by Karabulut and İnceoğlu (2005). The algorithm is detailed in Alg. 2. The basic concept is to place the items as far as possible to the back (first priority), to the bottom (second priority) and to the left (third priority) of the loading space. The available free spaces in the vehicle's loading space are stored in a list.
In the following, the point of origin of a Cartesian coordinate system is assumed to be located in the deepest, bottom, leftmost point of the loading space. The driver's cab is located behind it accordingly. The length, width and height of the loading space are parallel to the x-, y-and z-axes. The placement of an item I i,k is defined by ( x i,k , y i,k , z i,k ) of the corner which is closest to the point of origin.
(15) score op = iter best op ⋅ best + iter impr op ⋅ impr + iter e op ⋅ e . Before starting the packing process, the items of each customer are sorted by means of the following priorities: 1. fragility flag f i,k (non-fragile first) 2. volume (larger volume first) 3. length l i,k (longer first) 4. width w i,k (wider first).
Then, the items are added to the packing sequence IS reversed to the customer's visiting order [1]. Let S be the set of unique cuboids representing available free spaces for placing items. Initially, the set consists of one potential space, which corresponds to the entire loading space [2]. Therefore, the first item of the packing sequence is placed in the origin. The potential spaces of the set S are always sorted based on the DBL-rule [10]. Thus, an item is placed in the deepest, bottom, leftmost point of the space. For each item I c (c = 1, ..., |IS|), a feasible placement is determined [3]. For this purpose, each space sp of the set is tested as possible item position [5] until a feasible position is found obeying all loading constraints of the loading set P [7-8].
In comparison to Karabulut and İnceoğlu (2005), the set S does not contain all available placements inside the loading space. Rather, three new spaces (front, right, top) are created based on the feasible item placement [9].
The front (right, top) space is defined by the item's front (right, top) edge and either the door (wall, ceiling) or the nearest item in front (rightmost, topmost) of the item. Then, the minimum and maximum values for the y-(x, y) and z-axis (z, x) limited by the loading space or other items are searched. Fig. 11 shows these three created spaces exemplary based on item I 3,1 . Additional three spaces are created if they are unique: Another front and right space, where the minimum z-value represents the bottom edge of item I c , and another top space, where the minimum x-value is the deepest edge of item I c . The new spaces (front, right, top) are included in the set. After each feasible placement of an item I c , the available spaces are updated [13][14], which means that all available spaces are checked w.r.t. an intersection with item I c . If one or more spaces intersect with item I c , then these spaces are decreased so that no intersection occurs. Therefore, if an item can be placed within an available space, it is guaranteed that the item does not overlap with other items or with the vehicle's walls (geometry constraint (C1)). In contrast to the approach by Karabulut and İnceoğlu (2005), an overlapping check between each item is not necessary for this approach, which improves the performance. The used space is removed from the set [11]. To increase the efficiency of the packing heuristic and to reduce the number of spaces in the set, only spaces which are large and high enough for the smallest dimensions of any unplaced item of the route are inserted in the set S. Therefore, the shortest length or width l min and height h min of any unplaced item of the route are searched [12]. Due to the permitted rotations, only the two measures l min and h min are relevant. If the length or height of any space in the set is smaller than l min or h min , the space is removed from the set [15][16][17]. Then, a placement for the next item is searched [19].
If no feasible position for the item can be found, the route is revised, and a new one must be searched by the ALNS [24][25][26].

Computational studies
In this section, we investigate the solution quality and the performance of the hybrid algorithm in the context of advanced loading constraints. We use wellknown instance sets and investigate the impact of the proposed loading constraints on the objective values by means of a new instance set. All results along with detailed packing plans are available via https:// github. com/ Corin naKre bs/ Resul ts.
The hybrid algorithm is implemented in C++ as single-core, x64-application and is compiled using the GCC version 4.8.3, compiler. The experiments were executed on a High Performance Cluster, Haswell-16-Core with 2.6 GHz.

Parameters
The parameters for the loading constraints (see Sect. 4) and for the routing heuristic (see Sect. 5) are listed in Table 6. Regarding the parameters for the routing heuristic, we performed a preliminary study to tune the parameters. As the evaluation showed, the best results were obtained by the parameters as described in Koch et al. (2018) and therefore, these parameters were set. The parameters for the loading constraints are those used in the literature so far.

Instances
For our computational study, we use the instance sets by Moura and Oliveira (2009), Ceschia et al. (2013) and Zhang et al. (2017). Moreover, a new instance set consisting of 600 instances is created. The characteristics of the instance sets are shown in Table 7. Our new instance set is available via https:// doi. org/ 10. 24352/ UB. OVGU-2020-139.
The instances vary in the number of customers, items and item types. They either have 20, 60 or 100 customers, which demand either 200 or 400 items in total. These items differ in their homogeneity: Either there are only three item types (very homogeneous), 10 item types or 100 different item types (very heterogeneous). For realistic item masses, we analysed 12,000 products from a Swedish furniture company which offers products of different categories among other housewares, decorative articles and groceries. The densities of these products vary mainly between 0.5 and 1.5 kg/dm 3 . So the densities for the items are assigned by choosing a value randomly within this interval. Thereby, it is considered that the total mass for one customer is less than the vehicle load capacity D, since a customer can only be served by a vehicle (see constraints S1 and C4).
The fragility flag is set randomly to the items, where approx. 30% are fragile. To define the parameters for the load bearing strength, the formula by Ratcliff and Bischoff (1998) is used. For each item, a value r is determined depending on the fragility flag: If an item is fragile, then the value r is randomly chosen in the interval [1.0, Advanced loading constraints for 3D vehicle routing problems  To ensure a realistic proportion between vehicle load capacity and axle weights, parameters from the two-axle truck ML180 by IVECO were chosen, which has a maximum payload of 12,595 kg. Then, a proportional factor pr was calculated on the basis of the vehicle load capacity D of an instance and the maximum payload of the IVECO truck. Thus, the following applies: pr = D 12595 . The axle weights for the front and the rear axle were then proportionally scaled on the basis of pr.

Evaluation of hybrid algorithm
This section deals with the evaluation of the hybrid algorithm concerning its solution quality and performance. Hereby, we use our instance set and the benchmark instances by Ceschia et al. (2013), Zhang et al. (2017), and Moura and Oliveira (2009). Every instance is tested five times. We present summarized results. The more detailed results are presented in the "appendix" and are available via https:// github. com/ Corin naKre bs/ Resul ts. Note again that smaller objective values ( v used and ttd) represent better results.
In Table 8, we compare the received best and average results based on the basic constraint set used in Gendreau et al. (2006). On average, the difference between best and average for the number of vehicles is 1.42%, for the total travel distance only 0.42%. The results show a tendency that the more difficult the instances are (more customers or items), the higher the deviation between best and average results for the objective values. Therefore, we see potential for improvement as the hybrid algorithm should achieve an average deviation of only 1% at most.
The following Table 9 presents the best results per instance set. Concerning the Ceschia et al. (2013) instances, some of the instances require split delivery or feature a heterogeneous vehicle fleet. Excluding these instances, seven instances remain. Ceschia et al. (2013) have a maximum time limit varying between 300 and 10000 seconds. In contrast, the calculation time limit for our computational tests is only 3600 seconds. The results are based on the basic constraint set. The hybrid algorithm achieves clearly better results than the benchmark for nearly all instances (except SD-CSS04). On average, 19.33% less vehicles are used and the total travel distance decreases by 11.87%. Moreover, a shorter calculation time is required. 1 3 Advanced loading constraints for 3D vehicle routing problems For the Zhang et al. (2017) instances tested with the basic constraint set, the hybrid algorithm achieves a reduction of 23.24% used vehicles and a reduction of the total travel distance by 19.31%, on average. Moreover, the presented hybrid algorithm achieves the results with half of the calculation times compared to the benchmark. In case of the Moura and Oliveira (2009) instances, the instances do not provide any item masses, vehicle load capacities or fragility flags. Moreover, fully support of items is required ( = 1 ) and we only rotate the items along the lengthwidth plane. The currently best-known results are received by Reil et al. (2018). In comparison to the benchmark, the number of used vehicles increases by 20.24%, while the total travel distance decreases by 2.16%. The reason for these results is that we cannot use all rotation possibilities and only a small amount of iterations are conducted due to the high number of items per vehicle. We see potential for improvements of our hybrid algorithm to be able to compete with these instances.
To summarize, the hybrid algorithm finds new best results for two of three benchmark instances. In case of a high number of items per vehicle, the hybrid algorithm is not competitive so that we need to address this in our further research.

Evaluation of loading constraints
The following subsections analyse the impact of the different loading constraints on the objective values.

Constraint sets
In order to evaluate the impact of the new loading constraints in a systematical way, one new constraint is either replaced or added based on the basic constraint set P1. The last constraint set is a combination of the most restrictive ones. Table 10 shows the loading constraints as considered in each set. Details are as follows:

3
Advanced loading constraints for 3D vehicle routing problems  1. Replacing: The constraint sets P2-P6 are created by replacing one definition or implementation variant with another. 2. Adding: The constraint sets P7-P9 are generated by adding further loading constraints to P1. 3. Combination: The constraint set P10 is a combination of replacing and adding of loading constraints.

Results
For the analysis of the loading constraints, we use our new instance set to enable comparison concerning the number of customers (n), items (m) and item types. Every instance is tested five times for each constraint set. In sum, our analysis is based on 30,000 results (600 instances, 10 constraint sets, 5 runs). In the following Table 11, we report the average results and calculate the percentage difference to the basic constraint set P1. Note again that smaller objective values ( v used and ttd) represent better results. The impact of the MLIFO constraint (C5b) is evaluated by set P2. Compared to the LIFO constraint (C5a), on average, the objective values decrease slightly by around 0.2%. The solution space of the MLIFO constraint (C5b) is larger; contrasting our expectations, a significant influence of an increasing number of customers, items or item types on the objective values is not evident, though. To be more flexible in the handling of items, we recommend to rely further on the LIFO (C5a) constraint.
The robust stability (sets P3 and P4) reduces the solution space due to its more stable definition in comparison to the minimal supporting area definition. Thus, the constraints lead to a notable increase of the objective values and the calculation time. In case of the multiple overhanging constraint, the number of used vehicles rise by 10.80%, the total travel distance by 8.27%, on average. Since the top overhanging constraint is more restrictive, the number of used vehicles increases by additional 4.16% points, the total travel distance by 2.66% points, on average. For both variants, the calculation time increases by around 60%. Moreover, the more heterogeneous the items are (more item types), the more the objective values rise, since homogeneous item stacks items do not overhang and therefore, the constraint is fulfilled. Another aspect is the calculation time. In case of instances with 200 items, the calculation time increases by around 97%. In contrast, the increase for instances with 400 items is only around 37%. The explanation is as follows: In general, the higher the number of items, the higher the calculation time and the smaller the difference to the maximum calculation time. In case of the robust stability constraints, the maximum calculation time is exploited for most instances and therefore also for those, which had a small calculation time for the basic set P1. We recommend using the top overhanging constraint since its more stable definition.
The constraint sets P5 and P6 deal with the impact of the load bearing strength constraint. In general, the simplified and the complete selection variants lead to comparable objective values. On average, the number of used vehicles increases by approx. 3.2% ( v used ), the total travel distance by approx. 2.6% for both approaches.
The objective values of some instances are even smaller since the fragility constraint (C7a) is more restrictive than the load bearing strength constraints in this case. Since the load of items is calculated for the entire stack starting from the last placed item to the vehicle floor, the calculation time increases rapidly due to the algorithmic complexity (on average by around 29%) and also the objective values increase with the number of items. Interestingly, the objective values decrease with the number of customers. An explanation could be that with a lower number of customers, a higher number of items per customer is demanded. Since all items of a customer have to be packed into a vehicle, more items are stacked on top of each other, so that the limit values for the LBS are reached. Furthermore, our results show that a higher number of item types has a positive effect on the objective values. The reason is that with homogeneous item stacks, the load is distributed over fewer items. Since both variants lead to similar results and due to the fact that the complete selection variant is more realistic, we recommend using the complete selection.
The reachability constraint (C8), evaluated by set P7, leads to an increase of the number of used vehicles by 4.06% and the total travel distance by 2.80%, on average. However, the constraint has almost no effect on the objective values for half of the instances. A higher number of items or item types leads to an increase of the objective values by some per cent points, because in case of heterogeneous items or of overall more items, it is more likely that the items block the way. As the constraint has rather small impacts and avoids unnecessary rearrangements of items during unloading, it is recommended to take it into account in practically oriented VRP computations.
The sets P8 and P9 deal with the effects of distributed masses. The axle weights constraint (C9) distributes the masses in the vehicle along the x-axis, while the load is balanced along the y-axis in the balanced loading constraint (C10). For the majority of instances, the objective values remain unchanged or increase by only a few per cent. On average, the number of used vehicles increase by around 2.9%, the total travel distance by approx. 1.8%. Moreover, there is a positive effect on the calculation time, which is reduced by around 20%. A correlation between objective values and the number of customers, items and item types is not apparent. Due to the small impact on the objective values, the positive effects on the calculation time and the great safety relevance in traffic, considering both constraints makes sense if the appropriate information is available.
Since P10 contains all new complex constraints, the objective values clearly deteriorate, whereby the number of used vehicles increases more (24.42%) than the total travel distance (17.15%), on average. A comparison with the previous results reveals that a combination of several loading constraints leads to a deterioration of the objective values but not to the extent that the sum of the deteriorations would result from individually investigated constraints. The same applies for the calculation time.
Finally, regarding the results of all loading constraints, a correlation between an increasing number of customers or items and the objective values is not evident.

3
Advanced loading constraints for 3D vehicle routing problems However, an increase of the number of item types and therefore an increase of the degree of heterogeneity tends to be correlated with an increase of the objective values except for the load bearing strength constraints, where the load can be better distributed along heterogeneous items.

Conclusions and future work
This paper continues the research on the combined Vehicle Routing Problem and 3D Loading (3L-CVRP) introduced by Gendreau et al. (2006) and the extended problem with the consideration of Time Windows (3L-VRPTW). In our implementation, the possible placement of items is represented by free spaces inside the vehicle's loading space improving the performance of the algorithm. For a more realistic modelling, new loading constraints are introduced. Since the common definition of stability leads to unstable stacks, the robust stability is investigated by means of two implementation variants. In the first one, items of a stack are allowed to overhang when respecting a minimal supporting area at any height (multiple overhanging). In the second variant, only the topmost item of a stack is allowed to overhang (top overhanging). Moreover, instead of a simple fragility flag grouping items in fragile and non-fragile ones, the load per area unit for each item is considered. For this load bearing strength constraint, also two implementation variants (simplified and complete) are investigated. Additionally, constraints regarding the reachability of items, the axle weights and the balanced loading inside the vehicle are considered as well as the Manual LIFO by Tarantilis et al. (2009). The solution quality and the performance of our hybrid algorithm are evaluated by using well-known instances by Moura and Oliveira (2009), Ceschia et al. (2013) and Zhang et al. (2017). For the latter two instance sets, the presented algorithm performs better than the benchmarks.
The impact of the loading constraints on the objective values (number of used vehicles and total travel distance) is tested by 600 new instances varying in the number of customers, items and item types. In most cases, the Manual LIFO constraint has no influence on the objective values. Both variants for the load bearing strength constraint lead to comparable results. Therefore, the usage of the realistic "complete" variant instead of the simplified one is recommended. In case of the robust stability, we recommend using the top overhanging variant since it achieves higher stability of the stacks. The axle weights and the balanced loading constraints only lead to small increases of the objective values and even decrease the calculation time. Since they increase the vehicle stability and thus the safety, we recommend to investigate these further in future research. The same applies to the reachability constraint, which prevents unnecessary rearrangements during unloading. Furthermore, our investigations showed that when combining complex constraints, the results deteriorate, but not to the extent that the sum of the deterioration would result from 1 3 Advanced loading constraints for 3D vehicle routing problems  individually investigated constraints. As future work, we suggest to improve the performance of the hybrid algorithm and of the complex loading constraints. Furthermore, we plan to determine the influence of the individual neighbourhood operators on the results as well as impact of instance features on the packing algorithm.