1 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 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.

2 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.

2.1 3L-CVRP

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

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 and Wei 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.

2.2 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 multi-objective 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 mixed-integer 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) and Solomon (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 1 Summary and overview of approaches
Full size table

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.

Table 2 Summary and overview over loading constraints
Full size table

3 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} \in E\) \((i \ne j, i,j = 0, ..., n)\) has an associated routing distance \(d_{i,j}\) \((d_{i,j}> 0)\). The demand of customer \(i \in N \setminus \{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

  1. (S1)

    All routes \(R_{v}\) and packing plans \(PP_{v}\) are feasible (see below);

  2. (S2)

    Each customer is visited exactly once;

  3. (S3)

    The number of used vehicles \(v_{used}\) does not exceed the number of available vehicles \(v_{max}\);

  4. (S4)

    Each packing plan \(PP_v\) contains all \(c_{i}\) items of all customers i included in the corresponding route (\(i \in R_{v}\)).

A route \(R_{v}\) must meet the following routing constraints:

  1. (R1)

    Each route starts and terminates at the depot and visits at least one customer;

  2. (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.

  1. (C1)

    Geometry: The items must be packed within the vehicle without overlapping;

  2. (C2)

    Orthogonality: The items can only be placed orthogonally inside a vehicle;

  3. (C3)

    Rotation: The items can be rotated \(90^{\circ }\) only on the width-length plane;

  4. (C4)

    Load capacity: The sum of masses of all included items of a vehicle does not exceed the maximum load capacity D.

  5. (C5a)

    LIFO: No item is placed above or in front of item \(I_{i,k}\), which belongs to a customer served after customer i;

  6. (C5b)

    MLIFO: No item is placed on or in front of item \(I_{i,k}\), which belongs to a customer served after customer i;

  7. (C6a)

    Minimal supporting area: Each item has a supporting area of at least a percentage \(\alpha\) of its base area;

  8. (C6b)

    Robust stability: Each item has a supporting area of at least a percentage \(\alpha\) of its base area at any height;

  9. (C7a)

    Fragility: No non-fragile items are placed on top of fragile items;

  10. (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;

  11. (C8)

    Reachability: The distance between an item and the driver must be less or equal than a certain length \(\lambda\);

  12. (C9)

    Axle weights: The loads for the front and the rear axle do not exceed the permissible axle weights \(FA_{perm}\) and \(RA_{perm}\);

  13. (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.

4 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/CorinnaKrebs/SolutionValidator.

Table 3 Overview of loading constraints
Full size table

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.

4.1 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.

4.1.1 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.

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

Fig. 1
figure 1

Difference between LIFO and MLIFO

Full size image

4.2 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.

4.2.1 Minimal supporting area (C6a)

The minimal supporting area constraint ensures that a certain ratio \(\alpha\) 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 density, whereby the length or width of each item enlarge by \(\frac{1}{\alpha }\), an overhanging stack of items is created.

Fig. 2
figure 2

Unstable, feasible stack w.r.t. minimal supporting area

Full size image

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.

4.2.2 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 \(\alpha\) 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 \in U\) defines a plane. Another item \(I_b \in 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 area constraint. Otherwise, the constraint is violated and the placement of item \(I_{i,k}\) is rejected.

Fig. 3
figure 3

Determination of planes for item \(I_{1,5}\)

Full size image

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. 1.

    If \(distance_{ceiling} + h_{i,k} \ge 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. 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).

Fig. 4
figure 4

Implementation of robust stability—top overhanging

Full size image

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

4.3.1 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 non-fragile 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.

Fig. 5
figure 5

Infeasible item arrangement w.r.t. fragility constraint

Full size image

4.3.2 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:

$$\begin{aligned} support_{c, u} = \frac{supportArea_{c, u}}{l_c \cdot w_c}. \end{aligned}$$
(1)

Since the item \(I_c\) could overhang, but the load must be distributed in total, the support share is increased proportionally:

$$\begin{aligned} support_{prop} = \frac{support_{c, u}}{support_c}. \end{aligned}$$
(2)

The load acting on item \(I_u\) is given as follows:

$$\begin{aligned} load_{c, u} = support_{prop} \cdot {load_c}, \end{aligned}$$
(3)

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}\).

Fig. 6
figure 6

Mass distribution according to simplified selection based on item \(I_{1,6}\)

Full size image

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.

Fig. 7
figure 7

Mass distribution according to complete selection based on item \(I_{1,6}\)

Full size image

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.

4.4 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 \(\lambda\), 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}\).

Fig. 8
figure 8

Illustration of the distance search space for \(I_{3,1}\)

Full size image

If the distance is larger than \(\lambda\) 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 - \lambda\) (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.

4.5 Axle weights (C9)

The exceedance of the maximum axle weights of one or more axles leads to far-reaching 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.

Fig. 9
figure 9

Vehicle data

Full size image

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\approx 9.81\frac{m}{s^2}\)).

$$\begin{aligned}&s_{i,k} = L_f + x_{i,k} + l_{i,k}/2, \end{aligned}$$
(4)
$$\begin{aligned}&\frac{1}{WB} \cdot \sum _{i=1|i\in R_{v}}^{n}\sum _{k=1}^{c_i}(m_{i,k}\cdot g \cdot s_{i,k}) = F_{RA} \end{aligned}$$
(5)

and

$$\begin{aligned} \sum _{i=1|i\in R_{v}}^{n}\sum _{k=1}^{c_i}(m_{i,k}\cdot g)-F_{RA} = F_{FA}. \end{aligned}$$
(6)

The acting forces must be below the permissible ones, but also greater than zero to avoid uplifting:

$$\begin{aligned}&F_{FA} \le FA_{perm} \cdot g, \end{aligned}$$
(7)
$$\begin{aligned}&F_{RA} \le RA_{perm} \cdot g, \end{aligned}$$
(8)
$$\begin{aligned}&F_{FA} \ge 0, \end{aligned}$$
(9)

and

$$\begin{aligned} F_{RA} \ge 0. \end{aligned}$$
(10)

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.

4.6 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.

Fig. 10
figure 10

Mass distribution according to the position of \(I_{i,k}\)

Full size image

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:

$$\begin{aligned}&f(t) = {\left\{ \begin{array}{ll} t &{} t > 0 \\ 0 &{} \, {\text {else}} \end{array}\right. } \end{aligned}$$
(11)
$$\begin{aligned}&\sum _{i=1|i\in R_{v}}^{n}\sum _{k=1}^{c_i} \frac{m_{i,k}}{l_{i,k}} \cdot \bigg [f\bigg (\frac{W}{2}-y_{i,k}\bigg )-f\bigg (\frac{W}{2}-(y_{i,k}+w_{i,k}\bigg )\bigg ] \le D \cdot p \end{aligned}$$
(12)
$$\begin{aligned}&\sum _{i=1|i\in R_{v}}^{n}\sum _{k=1}^{c_i} \frac{m_{i,k}}{l_{i,k}} \cdot \bigg [f\bigg ((y_{i,k}+w_{i,k})-\frac{W}{2}\bigg )-f\bigg (y_{i,k}-\frac{W}{2}\bigg )\bigg ] \le D \cdot p \end{aligned}$$
(13)

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.

5 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.

5.1 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. 1. The corresponding line number of the algorithms are given in square brackets. In general, a solution is feasible if all loading and routing constraints are obeyed except S3 so that the used vehicles could exceed the number of available vehicles.

figure a

5.1.1 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.

5.1.2 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} \le n_{rem} \le 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.

5.1.3 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:

$$\begin{aligned} f(s) = ttd(s) + pen_v \cdot max(0, v_{used} + {|N_{miss}|} -v_{max}) + {\sum _{i \in N_{miss}}^{N_{miss}} (c_{0,i} + c_{i,0})}, \end{aligned}$$
(14)

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\in 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.

5.1.4 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].

5.1.5 Removal and insertion operators

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

Table 4 Overview removal operators
Full size table
Table 5 Overview insertion operators
Full size table

After a defined number of iterations \(iter_p\), the selection probabilities for the removal and insertion operators are adjusted [16-18] according to their improvement of the solution. This is described in detail in the following section.

5.1.6 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\)).

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:

$$\begin{aligned} score_{op} = iter_{best_{op}} \cdot \omega _{best} + iter_{impr_{op}} \cdot \omega _{impr} + iter_{e_{op}} \cdot \omega _e. \end{aligned}$$
(15)

Hereby, \(\omega _{best}\), \(\omega _{impr}\) and \(\omega _e\) are coefficients.

Then, the new weighting \(wg_{op}\) can be calculated. A reaction factor r regulates the influence of the adaptions:

$$\begin{aligned} wg_{op} = wg_{op} \cdot (1 - r) + r \cdot \frac{score_{op}}{counter_{op}} \end{aligned}$$
(16)

Moreover, \(counter_{op}\), \(iter_{best_{op}}\), \(iter_{impr_{op}}\) and \(iter_{e_{op}}\) are reset to zero.

5.1.7 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}\).

5.2 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.

Before starting the packing process, the items of each customer are sorted by means of the following priorities:

  1. 1.

    fragility flag \(f_{i,k}\) (non-fragile first)

  2. 2.

    volume (larger volume first)

  3. 3.

    length \(l_{i,k}\) (longer first)

  4. 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].

Fig. 11
figure 11

New spaces based on \(I_{3,1}\)

Full size image

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.

figure b

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-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-26].

6 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 well-known 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/CorinnaKrebs/Results.

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.

6.1 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.

Table 6 Routing and loading parameters
Full size table

6.2 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.

Table 7 Overview of instance sets
Full size table

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, 2.0]. Otherwise, r lays in the interval [1.0, 5.0]. Then, over all items, the value \(max \frac{m_{i,k}}{m_{i,k}\cdot l_{i,k}}\) is searched and for each item multiplied by r.

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 = \frac{D}{12595}\). The axle weights for the front and the rear axle were then proportionally scaled on the basis of pr.

6.3 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/CorinnaKrebs/Results. Note again that smaller objective values (\(v_{used}\) and ttd) represent better results.

Table 8 Comparison best and average results for P1, our instances
Full size table

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.

Table 9 Summarized best results for instance sets
Full size table

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.

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 (\(\alpha = 1\)) and we only rotate the items along the length-width 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.

6.4 Evaluation of loading constraints

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

6.4.1 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:

  1. 1.

    Replacing: The constraint sets P2–P6 are created by replacing one definition or implementation variant with another.

  2. 2.

    Adding: The constraint sets P7–P9 are generated by adding further loading constraints to P1.

  3. 3.

    Combination: The constraint set P10 is a combination of replacing and adding of loading constraints.

Table 10 Overview of constraints sets
Full size table

6.4.2 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.

Table 11 Deviation to P1 per constraint set, average results
Full size table

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

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