Modular Coordination of Vehicle Routing and Bin Packing Problems in Last Mile Logistics

Abstract

Background: Logistics and transport, core of many business processes, are continuously optimized to improve efficiency and market competitiveness. The paper describes a modular coordination of vehicle routing and bin packing problems that enables independent instances of the problems to be joined together, with the aim that the vehicle routing solution satisfies all the constraints from real-world applications. Methods: The vehicle routing algorithm is based on an adaptive memory procedure that also incorporates a simple, one-dimensional bin packing problem. This preliminary packing solution is refined by a complex, three dimensional bin packing for each vehicle to identify the infeasible packages. The method iteratively adjusts virtual volumes until reaching near-optimal routes that respect bin-packing constraints. Results: The coordination enables independent applications of an adaptive memory procedure to vehicle routing and a genetic algorithm approach to bin packing while joining them in a computationally tractable way. Such a coordinated approach is applied to a frequently used public benchmark and proven to provide commensurate costs while significantly lowering algorithm complexity. Conclusions: The proposed method is further validated on a real industrial case study and provided additional savings of 14.48% in average daily distance traveled compared to the current industrial standard.

1. Introduction

Transport and logistics play an important role in many industrial activities. At, e.g., the European level, the transport industry employs 10 million people and accounts for 5% of the gross domestic product [1]. The quality of the service ensures a sufficient quantity of products at the desired location and determines their price. There are two fundamental problems in the transport optimization process: the Vehicle Routing Problem (VRP) and the Bin Packing Problem (BPP). Both directly affect the overall transport costs. The goal of a Capacitated VRP (CVRP) is to minimize the total distance that all delivery vehicles need to traverse to deliver given products to desired locations. The CVRP treats capacity as one-dimensional, i.e., incorporates a one-dimensional Bin Packing Problem (1D-BPP). The goal of a three-dimensional BPP (3D-BPP) is to distribute the packages in as few vehicles as possible while meeting the physical and legislative constraints related to packages and vehicles. These constraints address product stability, product fragility, vehicle weight limit, balance, and more. The integrated solving of these two problems as one common problem achieves higher savings in transport but also makes the solving algorithm significantly more computationally complex.

To support a clearer understanding of the proposed approach, the following subsections outline the CVRP and 3D-BPP, followed by their integration and the key contributions of this work. While the algorithms for CVRP and 3D-BPP are not the main contributions of this paper, they are briefly presented for readers’ convenience to clearly articulate how their iterative coordination contributes to the overall approach.

1.1. The Vehicle Routing Problem

The VRP is an NP-hard combinatorial problem of finding an optimal combination of vehicle routes that minimizes the total distance traveled to deliver the goods to a given set of customers. Other optimization criteria are also applicable, such as delivery time [2], fuel consumption [3], or a combination of several factors [4]. Among different variants of VRP, two have been identified as the most related to this paper. The first is the already mentioned CVRP, which considers one dimension of the package (e.g., weight or volume) but not the height, width, length, and other characteristics of the package. The second is Heterogeneous fleet VRP (HVRP), which adds multiple different vehicles to the problem.

Methods for solving the VRP can be divided into four categories: exact methods, heuristics, meta-heuristics, and reinforcement learning. Exact methods [5] find optimal solutions but are computationally feasible only for small instances. Heuristics have lower complexity but are less accurate, while meta-heuristics such as Tabu Search [6,7], variable neighborhood search [8], adaptive large neighborhood search [9], genetic algorithms [10], and ant colony optimization [11] offer better-quality solutions at the cost of longer execution time. Reinforcement learning has recently gained attention [12,13], but its practical use remains limited due to a series of assumptions that are rarely satisfied in real-world problems due to limited datasets, different constraints, and others as reported in [14].

The adaptive memory concept was introduced and combined with Tabu Search in [6]. Adaptive memory is a pool of good solutions dynamically updated during the search process. The Adaptive Memory Procedure (AMP) from this paper still holds 12/14 best-known solutions on the benchmark from [15]. The authors of [16] proposed the Elite Parts Search Method of Solutions, which generates initial solutions through a non-stochastic population-based diversification technique. The authors report reduced suboptimality from 0.23% to 0.18% on the [15] benchmark. The authors of [17] combine the AMP and Tabu Search to solve the Fixed Fleet Open Vehicle Problem. Authors receive 1.00% lower cost in 14.2% less time than for pure Tabu Search. In [18], the authors combine AMP and Variable Neighborhood Search to solve a real practical VRP for the gas delivery industry. Adaptive memory is also implemented in [19] to find the new best solutions on all benchmark instances from [20].

1.2. The 3D Bin Packing Problem

The 3D-BPP is a combinatorial problem that aims to distribute L cuboid packages, given their length, width, height, and weight into a minimum number of cuboid bins (vehicles). Packages are not allowed to overlap, and their edges are required to be parallel to the vehicle edges. There are four approaches for solving 3D-BPP: exact algorithms, constructive heuristics, meta-heuristics, and reinforcement learning. Among them, meta-heuristics achieve the most efficient results in an acceptable time. A recent example of an exact algorithm is presented in [21]. The authors propose a linear mathematical formulation of 3D-BPP while respecting physical constraints encountered in the industry: fragility, stability, orientation, and weight distribution. Constructive heuristics [22,23] are the fastest methods to solve 3D-BPP, but the disadvantage of these methods is higher suboptimality. Therefore, they are used for large problems and as subcomponents in extended variants of VRP. Given the greatest efficiency of the results, most of the scientific papers in the field are focused on meta-heuristics. The best-known meta-heuristics include Tabu Search [24] and genetic algorithm [25,26,27,28]. Other examples are gray wolf optimization [29], hierarchical search [30], and ant colony optimization [31], or hyper-heuristics like [32]. Recently, methods based on reinforcement learning have also been developed [33,34,35]. These packing methods find their application in problems where not all packages are known in advance but appear one after the other.

1.3. Integration of CVRP and 3D-BPP

The VRP is used to optimize the order of transferring items from departure points to delivery points under different conditions. Therefore, VRP determines which items will be loaded into which vehicle and in which order they will be delivered. The BPP optimizes vehicle loading with the goal of fitting as many items as possible into the vehicle. However, most VRP variants take only one dimension to check whether all packages physically fit in the vehicle, which can lead to infeasible results in real-world scenarios. It is therefore beneficial to validate such solutions with bin packing to address physical and legal constraints often ignored in routing-only approaches. This paper addresses that limitation by proposing a method that coordinates routing and packing algorithms while maintaining their independence.

In existing scientific approaches, the integration of VRP and BPP is generally achieved through an extended variant of VRP called the three-dimensional loading capacitated vehicle routing problem (3L-CVRP). The goal of this variant is to minimize the distances that delivery vehicles need to traverse while meeting the constraints of 3D packing. A precise mathematical definition of the 3L-CVRP is given in [36], which is frequently used as a benchmark in other studies for comparison purposes. There are also variants of this problem that include different industrial constraints. These include pick-up delivery, time windows, split delivery, fuel consumption, and more.

Some of the established algorithms for solving the 3L-CVRP are presented in [3,36,37,38,39,40,41]. Authors in [37] combined Greedy Randomized Adaptive Search Procedure and Evolutionary Local Search (GRASPxELS) to achieve 3.91% lower cost than the algorithm from [36]. This approach solves only basic 3D constraints, while fragility, support, and LIFO are not solved. Also, although the algorithm had a lower cost, high computational complexity makes this approach more difficult to apply to real-world examples. In [3], authors presented an Evolutionary Local Search (ELS) for 3L-CVRP that minimizes fuel consumption. The algorithm includes common industry constraints: fragility, support, LIFO (“Last In First Out”—unloading schedule), and weight limit. This algorithm achieved 4.40–9.18% lower cost than the algorithm from [36]. Although this algorithm has improved the best-known solution in many instances and executes faster than previous approaches, the algorithm still has a very high complexity. In [38], authors developed an Adaptive Large Neighborhood Search (ALNS) for 3L-CVRP with pick-up delivery and time windows. The algorithm also includes stability, support, LIFO, and maximum weight constraints. Instead of solving the original benchmark instances from [36], the authors modified the set suitable to their use case. They introduced changes to make this set suitable for weakly heterogeneous problems. These problems have less different package sizes, which makes them easier to solve for the bin packing part. The same variant combined with split delivery is tackled in [39], which proposes a local search for routing and different heuristics for packing. The authors applied the “Packing first, routing second” principle. Common industry constraints are included in the cohesive algorithm. The algorithm has low computational complexity and achieved a 0.5% savings in travel distance when applying split delivery to the weekly heterogeneous dataset. However, the authors state that their algorithm had a 55.9% higher cost compared to [36] on the strongly heterogeneous dataset. Authors in [40] achieved 8% lower cost than [36] using a hybrid meta-heuristic approach but state that there is an overflow of time and distance on large-scale instances, even with some requiring 2 h on a 2.8 GHz CPU and in C++. The authors in [41] implemented a Column Generation (CG) approach, which, unlike the other methods mentioned, belongs to exact methods instead of meta-heuristics. Although CG theoretically provides more accurate solutions, it often comes with a significantly higher computational cost, particularly when dealing with large-scale problems. The authors reported a 1.79% higher average cost than [36] and state that the solution improvement is not impressive for more iterations of CG. Some recent papers that address 3L-CVRP on different test sets are [42,43,44,45].

The mentioned algorithms that usually solve 3L-CVRP combine CVRP and 3D-BPP into a single algorithm, which leads to the extremely high computational complexity of the problem. The result is a higher running time for such algorithms. In [36], a 4% improvement in traversed delivery path was achieved using a Tabu Search after 24 h of CPU running, compared to 10 h for relatively small instances (44–64 vehicles, 141–181 packages). The GRASPxELS hybrid algorithm from [37] solves 4/9 small instances (under 25 delivery points) in less than 0.1 s, while for bigger instances (50–100 delivery points), computation time was 2851.6 s on average. Authors in [38] report that their algorithm based on ALNS requires executing the bin packing component tens of thousands of times to find a viable solution. This number grows rapidly when the problem size increases.

1.4. Paper Contributions

The hard computational complexity of integrated problems, together with the scarcity of the literature on joining the problems, prove that these are among the most challenging problems today. On the other hand, joining the problems evidently offers additional opportunities for cost-efficiency and there is a clear motivation for further research. In addition, by solving VRP or BPP individually and ignoring the counterpart, the resulting solutions are often infeasible in the real world. To cope with this, motivated by industrial needs, we approach the problems modularly, treating each problem individually and then iteratively coordinating them towards the global optimum. For this, in this paper, we propose a Modular Iterative Algorithm for Routing and Packing (MIA-RP). With this approach, the number of algorithm runs for 3D-BPP can be significantly reduced.

The proposed algorithm is based on the virtual capacities of the packages, where a CVRP with simple 1D-BPP is executed and further verified by the complex 3D-BPP that respects all the practical constraints. If the latter is unfeasible, the virtual capacity of the packages is modified and a new iteration of problems is executed. This process is repeated until both problems provide feasible solutions.

The MIA-RP is computationally simpler than other meta-heuristic approaches because it decomposes the problem into two subproblems, VRP and BPP, instead of solving two NP-hard problems at once. By separating the problems, it theoretically lowers the complexity and then coordinates the two components using an additional algorithm with relatively lower complexity, significantly improving computational efficiency. Unlike column generation approaches, the proposed algorithm does not rely on Linear Program (LP) relaxation or complex pricing subproblems, which reduces both computational and implementation complexity. This modular and heuristic coordination of VRP and BPP components enables faster execution and easier adaptation to real-world constraints without the need for specialized optimization solvers.

From a commercial perspective, MIA-RP additionally enables the coordination of independent algorithms for VRP and BPP. It even allows independent companies to manage each part of delivery logistics, enabling any combination of different VRP and BPP variants for greater economic flexibility and easier system modernization.

We identify the following scientific contribution of this paper:

(1). Modular iterative algorithm for integrated VRP and BPP optimization of commensurate suboptimality but of significant improvement in execution speed when validated and compared to state-of-art approaches and ideal benchmark from [36],

(2). Real case study benchmark and validation of the proposed algorithm, proving significant route savings of 14.48% on average.

The algorithm is developed considering the complex, real-world problem and as such aims at the reduction of computational complexity of delivery logistics optimization. This paper is organized as follows. Section 2 elaborates on a CVRP variant formulation and the corresponding optimization solver. A genetic algorithm for 3D bin packing with industrial constraints is described in Section 3. Section 4 presents the MIA-RP algorithm for routing (CVRP) and packing (3D-BPP) coordination followed by computational experiments in Section 5. Conclusions are presented in Section 6.

2. Capacitated Vehicle Routing Problem Formulation

The CVRP can be presented as follows. Let $G=(V^G,A^G)$ be a directed complete graph, where $V^G=\{0,\dots ,N\}$ denotes the vertex set and $A^G=\{(i,j):i,j\in V^G,i\ne j\}$ is the arc set. Vertex $i=0$ corresponds to the depot (starting point), and others ($i\in V^G\backslash \left\{0\right\}$) are related to one of the N customers (delivery points). The depot, indexed by $i=0$, is a special case, and every route starts and ends with it. Every arc $\{i,j\}\in A^G$ has a corresponding non-negative travel cost $c_{ij}$ from vertex i to j. Since HVRP introduces different vehicles into the problem, the vehicle fleet contains K different vehicles with volume capacities $V_k,k\in \{1,\dots ,K\}$ and weight capacities $M_k,k\in \{1,\dots ,K\}$. Every delivery point ($i\in V^G\backslash \left\{0\right\}$) has volume demand $q_i$ and weight demand $w_i$. Parameters and variables for CVRP formulation are shown in Table 1.

Let binary decision variable $b_{kij}\in \{0,1\}$ be defined as 1 if vehicle k travels directly from node i to node j, and 0 otherwise, for all $i,j\in V^G,i\ne j,k\in \{1,\dots ,K\}$. The following applies:

(1)$\sum _{k=1}^K\sum _{\begin{array}{c}j=0\\ j\ne i\end{array}}^N{b}_{kij}=1,\forall i>0,$

(2)$\sum _{k=1}^K\sum _{\begin{array}{c}j=0\\ j\ne i\end{array}}^N{b}_{kji}=1,\forall i>0,$

where (1) and (2) ensure that every delivery point is visited exactly once by exactly one vehicle (one outgoing and one incoming arc).

The goal of the VRP algorithm is to find the minimum delivery cost $C^{*}$:

(3)$C^{*}=min\sum _{k=1}^K\sum _{i=0}^N\sum _{\begin{array}{c}j=0\\ j\ne i\end{array}}^N{c}_{ij}{b}_{kij}.$

In addition, the CVRP variant requires:

(4)$\sum _{i=1}^N{q}_i\left(\sum _{j=0}^N{b}_{kij}\right)\le V_k,\forall k,$

(5)$\sum _{i=1}^N{w}_i\left(\sum _{j=0}^N{b}_{kij}\right)\le M_k,\forall k,$

to ensure that each vehicle respects its volume and weight capacity constraints.

To avoid disconnected subtours (cycles that exclude the depot), auxiliary variables $r_i\in [1,N]$ are introduced for each customer $i\in \{1,\dots ,N\}$. The following constraints are used:

(6)$r_i-r_j+(N-1)b_{kij}\le N-2,\forall i\ne j,i,j\in \{1,\dots ,N\},\forall k,$

(7)$1\le r_i\le N,\forall i\in \{1,\dots ,N\}.$

They ensure consistent customer ordering and prevent infeasible cycles.

The algorithm that solves this problem is inspired by the algorithm from [6] and is independently proposed in [46] to solve the combination of VRP variants that do not include bin packing. In this paper, it is combined with an additional bin packing algorithm. For convenience, this algorithm is briefly presented in the sequel. It consists of three parts: initial solution generating, local search, and Adaptive Memory Procedure (AMP). The AMP includes one entire CVRP step, while initial solution generating and local search are run multiple times within one step. This is described in more detail in the subsection with the description of AMP.

2.1. Initial Solution

In this step, a solution is generated using constructive heuristics, and it is called the Randomized Constrained Clarke–Wright (RCCW) algorithm. The original Clarke–Wright algorithm [47] is a deterministic method where each delivery point is first assigned to one vehicle in a route, and then they are sequentially merged into common routes. The order of merging is determined according to the sorted sequence of savings values:

(8)$s_{ij}=c_{i0}+c_{0j}-a{c}_{ij},$

where delivery points i and j with the highest savings are connected first, etc.

The original algorithm is modified here using two hyperparameters to introduce the stochasticity needed for integration with the population-based Adaptive Memory Procedure. The first modification is obtained by changing the coefficient a, which is equal to 1 in the original algorithm. This is first introduced in [48] along with skipping randomly selected values $s_{ij}$. Although the a parameter introduces stochasticity, it often leads to repeated solution patterns. To improve diversity and avoid convergence to the same local optima, dropout parameter $p_d$ is used, which determines a percentage of the savings matrix ignored when generating the initial solution. The impact of these two hyperparameters on suboptimality and solution diversity is analyzed in detail in [46]. The RCCW is shown in Algorithm 1.

Algorithm 1: Randomized Constrained Clarke–Wright

load inputs: vehicles, locations     select route shape a     select dropout percentage $p_d$     generate solution of I vehicle routes (0, i, 0)     calculate savings matrix S of items $s_{ij}$ using (8)     sort values in savings matrix S in decreasing order into a new list, $S_{sorted}$     drop $p_d$ items from $S_{sorted}$ and all items where $s_{ij}<0$     for $s_{ij}$ in $S_{sorted}$ do         if routes (0, …, i, 0) and (0, j, …, 0) can be merged then            merge (0, …, i, 0) and (0, j, …, 0) to (0, …, i, j, …, 0)         end if     end for     compress solution to N vehicles:          sort routes by ascending number of allowed vehicles (1D volume and weight constraints)          assign each route to the allowed vehicle containing the least routes return generated initial solution

2.2. Local Search

Local search improves the initial solution as long as there is a better solution in the local neighborhood of the current one. Local neighborhood search strategies implemented here are swap nodes, move nodes, swap drives, and move drives. An example of swap nodes can be seen in Figure 1. This procedure swaps the places of two nodes in the existing solution for a full neighborhood and accepts the best swap for further optimization in terms of cost. The best solution in the neighborhood can be found with $O\left(n^2\right)$ complexity. Multiple groups containing multiple nodes can also be modified. An example of move nodes is shown in Figure 2 and can also be applied to larger groups of nodes. These strategies are executed in alternation using a greedy approach until a local minimum is reached.

2.3. Adaptive Memory Procedure

By running the initial solution and local search algorithms multiple times, different solutions are obtained. Most often, none of them are optimal, but they all contain fragments of the optimal routes. Adaptive memory procedure recognizes common parts of different solutions and fixes such parts. These parts become a new vertex in the graph, which reduces the total number of vertices and edges and finally the search space. New solutions are found in this reduced search space, and the initial solution and local search focus more on the remaining parts of it. The pseudo-code for the adaptive memory procedure is shown in Algorithm 2.

Algorithm 2: Adaptive Memory Procedure

data = vehicles, locations, items     for i_1 = 1 : iteration_limit_1 do        select new iteration_limit_2        select good solutions $percentage$        initialize solutions_list        for i_2 = 1 : iteration_limit_2 do           initial_solution = run_initial_solution(data)           solution = run_local_search(data, initial_solution)           append solution to solutions_list        end for        if i_1 < iteration_limit_1 then           keep best $percentage$ of solutions_list           vertices = select_vertices_to_connect(solutions_list)           data = reduce_data(data, vertices)        else           best_solution = select_best_solutions(solution_list, 1)        end if     end for

3. 3D Bin Packing Problem Formulation

This paper complements the common linear mathematical formulations of the BPP described as follows. A set of L rectangular packages of dimensions $l_l\times w_l\times h_l$ is given. The package l is preassigned to a specific vehicle k. A set of K available vehicles of dimensions $L_k\times W_k\times H_k$ is given. In the observed problem, each vehicle is associated with a predefined subset of packages, and the objective is to arrange the assigned packages into their respective vehicles while satisfying all physical and legislative constraints. All used indices, parameters, and variables are listed in Table 2. The conventional bin packing formulation from [21] is given in the sequel:

Find the minimum cost:

(9)$\begin{array}{c}\hfill F^{*}=\underset{u_k}{min}\sum _{k=1}^K{V}_k{u}_k,\end{array}$

s.t.

(10)$\sum _{k=1}^K{p}_{lk}=1,\forall l,$

(11)$\sum _{l=1}^L{m}_l{p}_{lk}\le M_k{u}_k,\forall k,$

(12)$\sum _{l=1}^L{v}_l{p}_{lk}\le V_k{u}_k,\forall k,$

(13)$x_l^{\prime }\le \sum _{k=1}^K{L}_k{P}_{lk},\forall l,$

(14)$y_l^{\prime }\le \sum _{k=1}^m{W}_k{P}_{lk},\forall l,$

(15)$z_l^{\prime }\le \sum _{k=1}^K{H}_k{P}_{lk},\forall l,$

(16)$x_l^{\prime }-x_l={\sigma }_{l11}{l}_l+{\sigma }_{l12}{w}_l+{\sigma }_{l13}{h}_l,\forall l,$

(17)$y_l^{\prime }-y_l={\sigma }_{l21}{l}_l+{\sigma }_{l22}{w}_l+{\sigma }_{l23}{h}_l,\forall l,$

(18)$z_l^{\prime }-z_l={\sigma }_{l31}{l}_l+{\sigma }_{l32}{w}_l+{\sigma }_{l33}{h}_l,\forall l,$

(19)$\sum _a{\sigma }_{lab}=1,\forall l,b,$

(20)$\sum _b{\sigma }_{lab}=1,\forall l,a,$

(21)$x_o^{\prime }\le x_l+(1-x_{lo})L_k,\forall l,o,$

(22)$x_l+1\le x_o^{\prime }+x_{lo}{L}_k,\forall l,o,$

(23)$y_o^{\prime }\le y_l+(1-y_{lo})W_k,\forall l,o,$

(24)$y_l+1\le y_o^{\prime }+y_{lo}{W}_k,\forall l,o,$

(25)$z_o^{\prime }\le z_l+(1-z_{lo})H_k,\forall l,o,$

(26)$z_l+1\le z_o^{\prime }+z_{lo}{H}_k,\forall l,o,$

(27)$x_{lo}+x_{ol}+y_{lo}+y_{ol}+z_{lo}+z_{ol}\ge (p_{lk}+p_{ok})-1,\forall l,o,k,$

where $l,o=1\dots L,k=1\dots K$.

Table 3 describes all elements of this mathematical formulation for BPP.

For industrial applications, different constraints are also included in the problem, such as package stability, fragility, and unloading schedule. The stability constraint requires that each package lie on a sufficiently large and flat surface. Part of its bottom area (typically 75%) must be supported either by the vehicle floor or by the top surfaces of other packages. An attribute of fragility $f_l$ is assigned to each package. If the package is fragile ($f_l=1$), only other fragile packages can be placed on its top surface, where both fragile and non-fragile packages can be stacked on a non-fragile package ($f_l=0$). To support efficient delivery at multiple stops, the unloading schedule constraint ensures that packages can be removed in order without needing to rearrange others. This is known as last-in-first-out (LIFO) in the literature. This is allowed by stacking packages in the reverse order of the planned deliveries, and with the property of the deepest bottom left with fill heuristic for stacking packages, the last stacked package can always be retrieved from the top side. Let i and j be two delivery points and i is visited before j on the route of the vehicle k. Let l and o be two packages that belong to i and j, respectively. Then, o cannot be placed between l and the rear of the vehicle or above b.

In the state of the art, genetic algorithms stand out as the most successful algorithms in solving 3D-BPP. Representative examples can be found in [25,26,27,28]. Inspired by these algorithms, we implement the genetic algorithm for 3D-BPP that incorporates the industry constraints described in the previous paragraph. The algorithm is denoted as GA-IC.

3.1. Genetic Algorithm

The solution of 3D-BPP is often represented as the sequence in which the packages will be packed, together with the corresponding positions in the vehicles. A genetic algorithm consists of two fundamental components. The first is a meta-heuristic that searches for the best packing sequence, and the second is a constructive heuristic that determines the coordinates of the packages in a vehicle based on a given sequence. The genetic algorithm is often identified with the former component. It is practical to note that constructive heuristic is important as a solution generator but also as an implicitly defined fitness function in the space of all packing sequences. In GA-IC, we apply a constructive heuristic similar to the improved deepest bottom left with fill (DBLF) presented in [26] with two differences. We refer to this modified heuristic as the Priority-Constrained DBLF. The DLBF always places the package so that the deep, bottom, and left sides are bordered by the edge of the vehicle or already stacked packages. The first difference from DBLF is that the Priority-Constrained DBLF forbids the placing of new packages behind existing packages. This modification is needed to preserve the unloading priorities given by the packing sequence. Another difference is that the location where the new package is placed depends on the shape of the package. Namely, packages in many companies are designed so that it is profitable to stack similar packages next to each other in order to leave less empty space. This algorithm therefore gives priority to stacking along the z-axis and x-axis if there is a package of the same dimensions below or next to it, and along the y-axis if it is not. This makes the algorithm more suitable for weakly heterogeneous problems. Decoding of chromosomes into 3D layouts is performed by the Priority-Constrained DBLF, which assigns feasible positions to packages based on sequence order, spatial constraints, unloading priorities, and similarity-based alignment. The resulting layout defines the fitness of the sequence. More information on the utilization of packing priorities can be found in the LIFO-correction section.

The genetic algorithm is performed iteratively until some of the stopping criteria are met. In each iteration, the whole family of packing sequences is observed. A packing sequence for one vehicle in the genetic algorithm is called an individual, and the whole family of the sequences is called a population. The term gene is reserved for the position of a given package in a sequence. Considering the population changes in an iterative process, its realization within a particular iteration is called a generation. There are four fundamental phases of the genetic algorithm: initialization, selection, crossover, and mutation. Initialization is performed only once at the beginning of the algorithm, while the remaining phases occur in each iteration.

3.1.1. Initialization

In the initialization phase, a population with a predetermined number of individuals is created. The goal of initialization is efficiency and diversity among individuals, and therefore, individuals are generated in different ways. A certain number of individuals is defined as a sequence of packages with descending volumes. The length, width, and height of the package and the property of non-fragility are also taken as the relevant values in the descending sequence of packages. Other individuals are defined as random packing sequences. In addition to the packing sequences, the corresponding package rotations are randomly selected. An example of an initialized population is shown in Figure 3. This phase ensures that the population contains both heuristic-based and random diversity to avoid premature convergence.

3.1.2. Selection

In each iteration of the genetic algorithm, the highest-quality individuals in the population are selected to create a new generation. The quality of an individual is realized as the fitness value obtained with the constructive heuristic applied to the mentioned individual. The set of selected individuals is called the inter-population, and its size corresponds to the size of the original population. In GA-IC, individuals are selected by the roulette wheel selection method, where higher-quality individuals have a higher probability of creating a new generation. The population is first sorted based on fitness, and the selection process ensures diversity while prioritizing higher-quality individuals.

3.1.3. Crossover

In this phase of the genetic algorithm, a new generation of individuals is constructed. The process of crossover is carried out on the individuals from the inter-population in pairs. Individuals from the inter-population are called parents, while constructed individuals are called children. In this process, it is important to ensure that the created individuals contain all the initially set packages. In [49,50,51,52], various crossover operators have been devised. The crossover operator used in GA-IC is inspired by order-based encoding from [49]. It is shown in Figure 4 and described in the sequel.

In the first step, a common place in the packing sequence for both parents is randomly selected. This place is called the crossover point and separates genetic material for each parent. Genetic materials up to the crossover point are carried from the first and second parent to the first and second child, respectively. The rest of the genetic material for the first child is taken from the second parent, starting from the crossover point. Packages are copied consecutively, except for already contained packages; those are omitted. Once the procedure reaches the end of the parent’s genetic material, it will continue from the beginning of the parent’s genetic material all the way to the point of intersection. By the described procedure, the first child will inherit from the second parent all packages that are not received from the first parent. For the second child, the procedure is performed analogously.

3.1.4. Mutation

Upon completion of the crossover, a random change is made to a certain percentage of individuals. This process is called mutation. In GA-IC, the mutation is realized as a position replacement of two random packages. It is performed on $5\%$ of individuals in the population.

3.1.5. Lifo-Correction

The GA-IC contains an additional phase in each iteration to address the unloading schedule. We denote it as LIFO-correction. At this phase, each packing sequence in the population is corrected in such a way that packages with a lower unloading priority take an earlier position in the sequence. Such sequences, decoded by the constructive heuristic, will generate packings that are consistent with the unloading schedule. The correction process is illustrated in Figure 5.

3.1.6. Stopping Criteria

In GA-IC, the iterative process stops when all packages are deployed or when the maximum number of iterations is reached.

To select parameters for the genetic algorithm, a grid search is implemented, and the results are available in the Supplementary Materials. In the implementation of GA-IC, the following parameters are used:

population size: 500 individuals,

elite members preserved between generations: 40,

mutants introduced per generation: 25,

crossover rate: 87% (applied to 435 individuals after elitism and mutation),

maximum number of iterations: 20.

4. Coordination of Routing and Packing

This section presents a coordination algorithm that iteratively runs algorithms for CVRP and 3D-BPP described in Section 2 and Section 3 and thus commonly optimizes these two problems. The algorithm can be applied independently of the algorithms presented in Section 2 and Section 3, serving as their wrapper. It applies to a wide range of VRP and BPP algorithms, with the option to add additional variants to each.

The MIA-RP algorithm itself does not require specific assumptions, but certain assumptions are necessary for the CVRP and 3D-BPP components:

each customer is visited exactly once,

each vehicle route starts and ends at the depot,

distances between locations are known,

all vehicles are available at the beginning with weight and volume constraints,

the cargo consists of cuboid 3D packages with defined characteristics (shape, volume, weight, fragility, support, rotation, and rotation permissibility)

unloading follows the LIFO principle where necessary.

Delivery vehicles, locations, and delivery cargo are taken as inputs into the algorithm. For each package, the variables actual volume $v_l$ and virtual volume $v_l^t$ are introduced. The actual volume is defined according to the package dimensions and does not change. The virtual volume is initially set to the value of the actual volume using (28) and is changed as needed through iterations of the algorithm.

Each iteration of the MIA-RP coordination algorithm consists of three phases:

Optimization of delivery routes by resolving the associated CVRP from (3) s.t. (1) and (2) and (4)–(7).

Validation of CVRP results by solving 3D-BPP for each individual vehicle from (9) s.t. (10)–(27), also considering shape, weight, fragility, support, rotation, and permissibility of rotation from [22].

Increasing the virtual volumes of all packages associated with vehicles that provide infeasible solutions.

The phases are described below.

4.1. Optimization of Delivery Routes

Algorithms for optimization of delivery routes consider the capacities of delivery vehicles, but due to the complexity of execution, only as one-dimensional problems. Because of this characteristic, the results of delivery routes sometimes suggest that packaging for a particular vehicle is possible, but in reality, it is not feasible due to different physical and legislative constraints on three-dimensional cargo and delivery vehicles. The algorithm for optimizing delivery routes is called at the beginning of each iteration of MIA-RP, and in the next phase, the feasibility of the obtained solutions is validated. The described procedure is denoted in Algorithm 3 as $solve\_VRP\left(input\right)$.

Algorithm 3: MIA-RP

input = vehicles, locations, items     Set initial virtual volumes with                                                                   (28)     for iteration = 1 : iteration_limit do         routes = solve_VRP (input), with                                                        (1)–(7)         stopping_ flag = True         for route in routes do             validation = solve_BPP (route), with                                             (9)–(27)             if validation $==$ False then                 increase_virtual_volumes (route), with                            (29) and (30)                 stopping_flag = False             end if         end for         if stopping_flag $==$ True then             break         end if      end for

4.2. Validation of VRP Results by Solving Packing Problems

After optimizing the delivery routes at the beginning of each MIA-RP iteration, the bin packing algorithm validates the feasibility of packing for each delivery vehicle. The algorithm seeks packing that will contain as many packages as possible. If not all packages are included in the best-found solution, the complete packing is considered infeasible. In this case, all packages for a given vehicle are added to the set of unallocated packages. To minimize the number of such situations, the algorithm is run an additional two times for vehicles that did not fit only one package. The described procedure is denoted in Algorithm 3 as $solve\_BPP\left(route\right)$.

4.3. Increasing the Virtual Volumes of All Packages Associated with Vehicles with Unfeasible Packing

If the packings for all delivery routes have not been successfully validated, a set of unallocated packages is defined as previously described. In case a package l belongs to the set of unallocated packages in the current iteration, its virtual volume is increased by a certain percentage amount denoted by ${\lambda }_l$ and passed for the next iteration of $t+1$. The percentage amount is chosen arbitrarily and usually brings the most benefit if it is in a 5–20% interval. During delivery route optimization, virtual volumes and weights are used as one-dimensional package sizes in (4) and (5). In the scenario, in which a virtual volume $v_l^t$ is increased for a given package l, the delivery route optimization algorithm will interpret the package as larger and will reserve more space for it in the next iteration as shown:

(28)$v_l^0:=v_l,$

(29)$v_l^{t+1}:=(1+{\lambda }_l)·v_l^t,\forall l.$

As a result of increases in all values of $v_l^t$, new virtual volume demands for (4) become:

(30)$q_i^{t+1}:=\sum _{l=1}^L{v}_l^{t+1}·d_{il},\forall i.$

For the next iteration of VRP, $q_i:=q_i^{t+1}$ in (4) is therefore used. When validating the results in phase 2, the bin packing algorithm uses the actual volumes and weights. This gives a real confirmation of the feasibility of the packing. By performing this procedure through a series of iterations, the possibility of delivery route results for which packing is not feasible is continuously decreasing. The procedure is denoted in Algorithm 3 as $increase\_virtual\_volumes\left(route\right)$.

4.4. Stopping Criteria

The coordination algorithm stops when all delivery routes from phase 1 are successfully validated or when the maximum number of iterations is reached. In the latter case, the validated routes are taken directly as part of the final solution, while the unallocated packages are sequentially packed into available vehicles until all are deployed. When packing unallocated packages, the delivery route optimization is not conducted. In practice, a sufficiently large maximum number of iterations is taken so that this criterion is very rarely met. Stopping criteria are verified through binary variable stopping_flag.

The output of the algorithm consists of selected routes, package locations in vehicles, and delivery costs. The delivery cost is expressed in the desired currency, and the components that determine the cost are recorded separately: the total distance traveled, the total delivery time, and the total cost of using all vehicles. The pseudo-code for the coordination algorithm is shown in Algorithm 3 and presented in Figure 6 as a flow chart. An example of a problem solved with the algorithm is shown in Figure 7. In the figure, coordination is performed in three iterations. At first, CVRP arranged all the packages for the red and blue store in one vehicle, but 3D-BPP reports that there are unplaced packages. The virtual volume increases, but a similar problem appears in the second iteration. Finally, in the third iteration, the packages are distributed in such a way that they all fit in the two given vehicles.

5. Computational Experiments

The complete MIA-RP algorithm, which encompasses both the CVRP and 3D-BPP components, was implemented in Python 3.9.13 without relying on any external solver. The algorithm was tested on two datasets. The first is a publicly available benchmark with strongly heterogeneous packages, which introduces higher computational complexity due to significant size variations and is used to demonstrate comparability with state-of-the-art methods. The second dataset, provided by our industrial partner, is more relevant to our target application and contains weakly heterogeneous packages. The algorithm is better suited for the latter, as it focuses on optimizing VRP constraints while BPP only checks constraints. However, it is also tested on publicly available strongly heterogeneous benchmark instances to highlight its versatility and efficiency across different dataset types.

5.1. Experiments on Instances from the Literature

The MIA-RP algorithm was tested on the benchmark from [36] as the most common benchmark to test the combination of CVRP and 3D-BPP. The benchmark contains 27 instances with strongly heterogeneous packages, each with 5 combinations of constraints consisting of 3D-loading, support, last in first out (LIFO), and fragility.

As part of our evaluation, we tested hyperparameters on nine representative benchmark instances with the fragility constraint to determine values for the parameter ${\lambda }_l$, which is used to increase virtual package volumes based on size and fragility. Figure 8 compares solutions when varying the base ${\lambda }_l$ and the adjustments for small and non-stackable packages. Results showed that increasing ${\lambda }_l$ for small and non-stackable packages improves solution quality. Based on these findings, we selected ${\lambda }_l=10\%$ for regular packages, ${\lambda }_l=15\%$ for small packages, ${\lambda }_l=14\%$ for non-stackable packages, and ${\lambda }_l=19\%$ for small non-stackable packages. A package is considered small if its volume is below the median volume in the instance. Testing additional variations for large packages and those with one dominant dimension (one side exceeding one-half of vehicle length) did not yield improvements and were not included in the final configuration.

The test results for all 27 benchmark instances using selected ${\lambda }_l$ values are presented in Table 4. The best of ten runs is taken for each instance, and the lowest cost is selected, where priority is given to the minimum number of vehicles.

The MIA-RP works better in simple cases with less complex constraints than in cases with more constraints. This is due to the nature of the algorithm. The algorithm primarily solves the VRP with all its constraints, while the BPP and its constraints are checked after the optimization process. Therefore, the algorithm primarily optimizes the main cost (distance) for constraints related to VRP, while with constraints related to BPP, the focus is on keeping the complexity of the algorithm. When adding more BPP constraints, especially when both support and LIFO must be satisfied, the sum of package volumes for the optimal solution is considerably lower than the volume capacity of the vehicle. Therefore, increasing the virtual capacity has a lower impact. This situation is more common in larger problems because only one vehicle without a solution is enough for MIA-RP to run a new iteration with larger virtual volumes. Therefore, for a more complex set of BPP constraints, larger problems may require more delivery vehicles, while in a simpler set of constraints, all solutions were found within the given set.

When considering only 3D loading, all instances are solved with a minimum number of vehicles, which enables comparison of the average cost. Table 5 shows this comparison with papers from the literature: the original paper [36] that published the benchmark (TS), papers [37] (GRASPxELS) and [3] (ELS), which, according to our research, have the best results on the same benchmark, an approach related to ours from [41] (CH), and with recent work using Constraint Programming (CP) with some new best solutions [42]. The MIA-RP achieved a 0.50% lower cost than the algorithm from [36], 2.30% than [41], and 0.66% compared to [42]. Algorithms from [3,37] achieved 3.44% and 3.93% lower cost than MIA-RP, respectively. These results show that MIA-RP is comparable to state-of-the-art algorithms in cost.

The main advantage of MIA-RP is lower complexity, which is reflected in the results for execution speed and shown in Table 5. The MIA-RP-generated results are on average 18.0–254.0% faster than for other algorithms, but considering that the algorithms were run on different processors (CPUs) in different programming languages, a more detailed analysis of the execution speed is needed. While Python is expected to be generally much slower than C or C++, some parts of the code are executed in a commensurate time. The algorithm from [36] is implemented in C and [3,37] are implemented in C++. To make the results of different approaches comparable, we analyze how fast each of the algorithms solves the problems for small ($K<25$), medium ($25\le K<50$), and large ($K\ge 50$) instances. This is motivated by [3,36], where the authors generate instance categories in the same way. We excluded the results from [41] because computing time is missing on the same instances and [42] because the lower limit of time is 300 s and is therefore unsuitable for this complexity analysis. Other state-of-the-art algorithms from Table 5 have extremely variable execution times for problems of the same size. In each one, there are some examples where a smaller problem is solved over 100 times slower than a larger problem. The MIA-RP has a much more stable runtime, and examples like this cannot be found in this table. To avoid excessive variability in execution time, for each of the three groups of problem sizes, we took the median execution time. In Figure 9, the increase in execution time in medium and large instances compared to small ones for each of the algorithms on a logarithmic scale is shown. The MIA-RP is a clear winner here, with 38 times slower execution for large than for small problems, while other approaches achieved 285–30,666 times slower execution. Following the analysis of the proposed algorithm implementation to ideal established benchmarks, it can be concluded that we sacrifice up to 3.93% suboptimality to obtain up to 807 times faster execution speed in these particular cases. This indicates a significant reduction in computational complexity compared to other approaches. This simplification is necessary and motivated by the goal of real-world implementation, which considers many additional constraints and significantly increased problem dimensions.

Despite being particularly adapted for weakly heterogeneous packages common to real-world problems, MIA-RP achieved excellent results on this benchmark, being the fastest among the tested approaches and ranking third in solution quality. Stable runtime and the potential for significant acceleration through implementation in a faster programming language confirm that the algorithm outperforms other approaches in execution time. The behavior of MIA-RP is consistent with established insights. In smaller instances, solution costs are similar to those of other methods.

With more BPP constraints, the solution cost grows faster than in approaches with a similar average cost. This is standard for constraint-rich problems. Although additional VRP constraints were not tested, the design of the algorithm suggests that MIA-RP would likely handle them better. As the instance size increases, MIA-RP consistently outperforms methods with comparable average costs, showing better scalability. Finally, the variability of the results remains within the expected range for heuristic methods at this level of suboptimality. In 26/27 instances, the standard deviation is lower than 2% of the average cost, with an overall average of 0.88%, which supports the stability and reliability of the proposed approach.

5.2. Experiments on Real-World Instances

The algorithm described was tested on the example of an industry partner. The partner is a logistic company that delivers products to a large number of locations, of which 118 are located in the analyzed region. Testing was performed on a five-day dataset with a number of packages between 335 and 704 delivered to 61–106 locations. In Table 6, the CVRP cost through iterations of MIA-RP and the number of unplaced packets obtained by solving 3D-BPP are monitored. Three different values of ${\lambda }_l$ are tested on the first of five days from the dataset. Those values are 5%, 10%, and 20% and are the same for all packages l and are referred to as $\lambda$ in further analysis. For each $\lambda$, it was tested how three values change through iterations. The first is the CVRP cost while ignoring the 3D-BPP. The second value is the number of unplaced packets in such a solution. Unplaced packets are all packets that have not been packed after 3D-BPP has attempted to place all packets into vehicles as specified by CVRP. The MIA-RP algorithm ends when this number reaches zero. The number of drives is the third tracked value. The higher the virtual package volumes, the more required drives are expected.

The lower the $\lambda$, the higher the expected number of required MIA-RP iterations. A combination of higher $\lambda$ and more iterations generally leads to a higher final CVRP cost. This is shown in the table where the CVRP cost has a growing trend through the iterations. This growth is faster with higher $\lambda$. It is expected that with twice the parameter $\lambda$, almost twice as many iterations are needed to find a valid solution. Therefore, with a high $\lambda$, which leads to a small number of iterations, it often happens that the CVRP cost in the last iteration is significantly higher than in the previous iteration. This leads to a higher average solution cost. On the other hand, low values of $\lambda$ give an average lower cost of CVRP, but the duration of the algorithm execution increases.

Distances (costs) obtained in the last CVRP iteration of MIA-RP were compared with three other scenarios in Table 7. The first is the distance obtained in the first iteration of MIA-RP, i.e., the distance obtained by CVRP without using 3D-BPP. The results obtained here are infeasible in practice, as CVRP considers only 1D bin packing and therefore regularly provides a solution where packages associated with certain vehicles are not all possible to fit.

The leftover, unplaced packages are then added to the additional one or more vehicles and delivered to the remaining destination after recalculating the CVRP for them. This (second) scenario is denoted as CVRP+3D-BPP in Table 7.

The third distance for comparison is the distance traveled by the industrial partner’s vehicles. The results of the industrial partner do not follow the original distance traveled by the vehicles, which is higher than the value in the table. Instead, the assumption was made that vehicles took the same (shortest) routes used for the CVRP. In this way, the influence of bad route selection between two delivery points was ignored, and only the order of visiting the points performed by the industrial partner was taken into account.

The results for the first scenario showed an average of 12.92% lower traveled distance than for MIA-RP (with $\lambda =5$), but this solution is not valid because the constraints from 3D-BPP are not taken into account. In the second scenario, which is feasible, unlike the first, the total distance traveled is much higher. With MIA-RP, the delivery was made with 15.97% lower distance traveled in the average case, 21.35% in the best, and 10.55% in the worst case. Compared to the results of the industrial partner, MIA-RP brought an average saving of 14.48%, 26.76% lower cost for the best case, and 0.35% higher cost for the worst case. Comparing our average savings with 10%, 10.74%, 8%, and 0.5% from [3,38,39,40], a significant difference can be seen in favor of our algorithm. Still, for a deeper and fairer comparison, it is necessary to put the listed approaches in the same context. Here, we provide a tentative analysis as a basis for future observations.

It should be noted that there is a difference between the nominal and real capacity of the vehicle. Our algorithm considers the nominal (lower) capacity, while the industrial partner used the real (higher) capacity of the vehicle. Also, it is possible that the partner increased the nominal volume of packages. The differences between these volumes are especially visible on days 4 and 5, where our algorithm uses more vehicles, even though the cost of MIA-RP compared to the industrial partner is equal or lower. On the other hand, on day 1, the algorithm managed to arrange all packages into fewer vehicles despite this aggravating circumstance.

Table 7 also analyzes how the value of $\lambda$ affects the final CVRP cost and the total number of iterations, similar to the analysis from Table 6. For $\lambda =5\%$, the average CVRP cost after MIA-RP is 446.62 km; for $\lambda =10\%$, the average cost is 453.83 km; and for $\lambda =20\%$, the average cost is 459.80 km. The average number of iterations for these three cases is 13.6, 7.6, and 4.4, respectively. The code was written in the Python programming language, and the average execution time on a single Intel i5-10210U CPU (1.6GHz) with 16 GB RAM was 91 min, 50 min, and 29 min, respectively. The average iteration duration is 6.62 min.

5.3. Practical Implications, Limitations, and Future Work

The MIA-RP provides a computationally efficient and modular framework for solving complex logistics problems. By coordinating individual VRP and BPP components, the approach allows effortless integration of additional constraints or, if needed, the substitution of the VRP or BPP solvers with alternative algorithms. This modularity enables organizations to adapt the system to their specific operational needs without being bound to a particular solution design. However, as with any modular approach, its effectiveness can be limited if one of the components consistently produces suboptimal solutions. Different algorithmic choices for VRP and BPP modules can be explored in future work to further enhance performance.

The MIA-RP demonstrates significant improvements in execution speed while maintaining solution quality comparable to existing state-of-the-art algorithms for 3L-CVRP, making it applicable to dynamic and operationally complex settings. Implementing the VRP and BPP algorithms in more efficient programming languages could additionally reduce computational effort. Furthermore, execution speed could be improved by reducing the number of iterations, which can be achieved by dynamically adjusting the parameter $\lambda$ based on the proportion of unsuccessfully assembled load in a given iteration.

6. Conclusions

In this paper, the modular iterative algorithm for routing and packing is proposed to coordinate vehicle routing problem (VRP) and bin packing problem (BPP) algorithms in a modular manner. This paper also describes algorithms for CVRP and 3D-BPP. The goal of this set of algorithms is to enable solving real problems in the field of logistics in a way that covers a large number of constraints. The MIA-RP was tested using the example of a partner company where all described constraints were considered. On the five-day dataset, the MIA-RP algorithm achieved a 15.97% lower cost than a simple algorithm that fills a new vehicle with excess packages and a 14.48% lower cost than the industrial partner. This shows that the analyzed algorithm achieves significant savings compared to the current state of the industry.

The MIA-RP algorithm was also tested on a widely used benchmark from the literature, where 0.66% lower to 3.93% higher costs and 7.5–807 times lower execution times are obtained compared to state-of-the-art algorithms of higher complexity. The best-known solutions were found in 6/27 instances.

Author Contributions

Conceptualization, N.P., A.K., and V.L.; methodology, N.P., A.K., and V.L.; software, N.P. and A.K.; validation, N.P. and A.K.; formal analysis, N.P., A.K., and V.L.; investigation, N.P. and A.K.; resources, V.L.; data curation, N.P.; writing—original draft preparation, N.P. and A.K.; writing—review and editing, N.P. and V.L.; visualization, N.P. and A.K.; supervision, V.L.; project administration, V.L.; funding acquisition, V.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data is unavailable due to privacy.

Acknowledgments

Subject matter disclosed herein forms an essential part of the patent application no. P20221281A. The authors would like to thank GDi LLC for their collaboration.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Swap nodes (2 groups of 1 node).

Figure 2 Move one node.

Figure 3 An example of an initialized population $P_0$ of 100 individuals. Each individual consists of a sequence of packages generated by a certain rule, with randomly selected rotations.

Figure 4 Order-based crossover, with one crossover point. The permutations in the figure represent packing sequences, i.e., individuals.

Figure 5 LIFO-correction. According to the corrected sequence, the packages with the highest priority (1) will be loaded last and unloaded first.

Figure 6 MIA-RP flow chart.

Figure 7 An example of MIA-RP. The figures show routing and packaging solutions through three iterations of the algorithm.

Figure 8 Analysis of different hyperparameter combinations. Each label shows base ${\lambda }_l$, ${\lambda }_l$ increment for small packages and ${\lambda }_l$ increment for non-stackable packages in percentages.

Figure 9 Comparison of execution time increase over problem size for MIA-RP with different state-of-the-art methods [3,34,35].

Table 1

Definitions of parameters and variables used in the CVRP formulation.

Parameters:
$c_{ij}$	Cost of traveling from node i to node j,
$q_i$	Volume demand of customer i,
$w_i$	Weight demand of customer i,
$V_k$	Volume capacity of vehicle k,
$M_k$	Weight capacity of vehicle k,
K	Number of available vehicles,
N	Number of customers,
Variables:
$b_{kij}=\left\{\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right.$	if vehicle k travels directly from node i to j,otherwise,
$r_i$	Subtour elimination variable for customer i,
$k=1\dots K,i,j=0\dots N,i\ne j.$

Table 2

Extended definitions of parameters and variables in [21].

Parameters:
L	total number of packages,
$m_l$	weight of package l,
$l_l\times w_l\times h_l$	length × width × height of package l,
$v_l$	volume of package l,
$V_k$	volume capacity of vehicle k,
$L_k\times W_k\times H_k$	length × width × height of vehicle k,
$M_k$	weight capacity of vehicle k,
Variables:
$u_k=\left\{\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right.$	if vehicle k is used,otherwise,
$p_{lk}=\left\{\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right.$	if package l is in vehicle k,otherwise,
$x_l,y_l,z_l\in \mathbb{N}$	location of the front/left/bottom corner of package l,
$x_l^{\prime },y_l^{\prime },z_l^{\prime }\in \mathbb{N}$	rear/right/top corner of package l,
${\sigma }_{iab}=\left\{\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right.$	if the side b of package l is along the a-axis,otherwise,
$x_{lo}=\left\{\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right.$	if package l is on the right of package o ($x_o^{\prime }\le x_l$),otherwise,
$y_{lo}=\left\{\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right.$	if package l is behind package o ($y_o^{\prime }\le y_l$),otherwise,
$z_{lo}=\left\{\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right.$	if package l is above package o ($z_o^{\prime }\le z_l$),otherwise,
$d_{il}=\left\{\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right.$	if package l is part of demand i,otherwise,
$f_l=\left\{\begin{array}{c}1\hfill \\ 0\hfill \end{array}\right.$	if package l is fragile (non-stackable),otherwise,
$l,o=1\dots L,k=1\dots K,i=0\dots N.$

Table 3

Description of objective and constraints.

Element	Description
Objective function (9)	Minimizes the total volume of used bins.
Constraint (10)	Ensures each package is assigned to exactly one bin.
Constraint (11)	Ensures the weight capacity of each bin is not exceeded, similarly to constraint (5).
Constraint (12)	Extends constraint (4) to account for the three-dimensional, real volume of packages. Unlike the VRP, the BPP considers actual 3D dimensions rather than scalar (1D representation) volumes. This necessitates additional spatial constraints.
Constraints (13)–(15)	Ensure that each package is placed within the physical bordersof its vehicle.
Constraints (16)–(20)	Enable orthogonal rotation of packages.
Variables $x_{ik},y_{ik},z_{ik}$	Binary variables indicating the relative position of package i with respect to package j along each spatial axis.
Constraints (21)–(26)	Define the semantics of the relative position variables.
Constraint (27)	Ensures that packages from the same vehicle do not overlap.

Table 4

The MIA-RP results on the ideal benchmark dataset from [36]. The numbers represent delivery costs and execution time.

Instance	All Constraints		No Fragility		No LIFO		No Support		3D Loading Only
	Cost	Time [s]	Cost	Time [s]	Cost	Time [s]	Cost	Time [s]	Cost	Time [s]
E016-03 m	287.75	77.0	287.75	98.4	287.75	66.4	287.40	96.4	287.40	75.4
E016-05 m	334.96	2.7	334.96	2.6	334.96	3.1	334.96	2.8	334.96	2.8
E021-04 m	380.46	93.5	380.46	104.8	364.28	14.4	362.28	15.7	362.27	16.8
E021-06 m	430.89	11.5	430.89	5.6	430.89	5.7	430.89	4.5	430.89	5.2
E022-04 g	465.01	57.2	457.50	56.5	424.35	80.7	418.95	78.3	395.64	47.0
E022-06 m	496.28	27.2	496.28	29.9	495.85	5.7	495.85	10.1	495.85	5.6
E023-03 g	789.77	72.2	747.30	39.3	742.24	39.4	732.24	61.2	732.52	27.9
E023-05 s	856.40	128.1	827.39	63.3	775.69	50.2	785.93	50.1	730.66	40.5
E026-08 m	635.93 *	86.2	643.61	106.7	630.13	67.9	630.13	86.7	630.13	73.4
E030-03 g	940.87	185.8	864.80	123.4	785.73	132.9	741.02	58.8	711.57	45.7
E030-04 s	890.27	281.5	819.37	121.6	761.54	101.2	751.87	96.3	695.71	38.2
E031-09 h	616.95	14.1	610.00	16.9	616.92	14.9	616.16	20.3	610.00	17.2
E033-03 n	2788.68	186.2	2687.24	87.6	2571.62	104.0	2606.60	142.1	2426.16	107.3
E033-04 g	1632.60	346.3	1613.02	159.4	1417.40	138.5	1453.64	208.0	1282.99	127.9
E033-05 s	1687.18	318.6	1601.61	159.1	1333.50	158.3	1339.93	156.5	1201.23	99.8
E036-11 h	708.87	223.6	698.61	31.8	702.70	17.8	705.55	52.1	698.61	14.7
E041-14 h	875.12	31.5	872.56	79.1	877.39	32.2	884.72	45.6	871.64	21.7
E045-04 f	1348.33	418.5	1284.82	421.2	1187.04	274.1	1166.33	272.6	1115.97	210.5
E051-05 e	851.13	390.6	808.89	341.8	727.90	221.1	743.80	255.0	692.00	220.4
E072-04 f	783.62	1683.6	717.68	1367.3	605.95	589.7	627.52	865.9	532.16	576.4
E076-07 s	1311.32	1404.0	1215.66	1006.2	1088.16	858.2	1070.98	1012.8	980.85	469.4
E076-08 s	1388.62	1786.3	1304.02	1442.2	1186.85	899.3	1189.18	963.1	1091.33	802.3
E076-10 e	1397.00	2363.7	1274.88	1577.8	1074.29	635.3	1135.57	872.7	1019.87	852.2
E076-14 s	1369.94	1883.4	1244.12	1295.5	1078.49	932.9	1129.94	1161.4	1054.40	968.6
E101-08 e	1755.02	2980.5	1746.31	4476.5	1412.89	3287.2	1415.82	2928.9	1278.25	1635.8
E101-10 c	2013.26	3713.8	1979.83	4120.5	1564.28	1912.7	1656.86	3265.9	1483.92	2807.1
E101-14 s	1797.93	3902.1	1734.98	4412.3	1529.47	4653.7	1535.23	2579.0	1396.57	2646.1
Average	1067.93	839.6	1025.4	805.46	926.23	566.6	935.16	569.0	871.98	442.8

* A solution for instance was found with a higher number of vehicles than suggested by [36].

Table 5

Comparison of delivery costs between MIA-RP and existing approaches in 3L-CVRP on the ideal benchmark dataset from [36].

	TS Gendreau et al. (2006) [36]		GRASPxELS Lacomme et al. (2013) [37]		ELS Zhang et al. (2015) [3]		CG Mahvash et al. (2017) [41]		CP Kucuk et al. (2022) [42]		MIA-RP Out of 10 Runs
Instance	Best Cost	Time [s]	Best Cost	Time [s]	Best Cost	Time [s]	Best Cost	Time [s] *	Best Cost	Time [s]	Avg. Cost **	std.	Best Cost	Time
E016-03 m	297.65	3.4	297.65	0.0	297.65	0.3	290.92	19.7	282.95	379	287.47	0.14	287.40	75.4
E016-05 m	334.96	0.6	334.96	0.0	334.96	0.1	334.89	5.7	334.96	308	334.96	0.00	334.96	2.8
E021-04 m	362.27	448.1	362.27	0.2	362.27	1.5	362.18	36.5	362.27	376	364.77	2.58	362.27	16.8
E021-06 m	430.89	11.1	430.89	0.0	430.89	0.1	430.78	12.44	430.88	311	435.70	3.87	430.89	5.2
E022-04 g	395.64	0.5	379.43	0.1	395.64	2.7	396.08	38.17	389.87	427	396.83	3.56	395.64	47.0
E022-06 m	495.85	14.7	495.85	0.0	495.85	0.2	495.71	18.65	495.85	311	501.25	4.47	495.85	5.6
E023-03 g	742.24	1.8	725.43	4.9	725.44	1.7	739.83	47.50	738.45	405	732.52	0.00	732.52	27.9
E023-05 s	735.14	104.9	735.14	1.1	730.66	10.1	735.03	42.79	708.32	437	734.49	4.74	730.66	40.5
E026-08 m	630.13	977.8	630.13	0.1	630.13	1.6	628.08	35.24	625.10	402	630.13	0.00	630.13	73.4
E030-03 g	717.90	410.7	687.57	32.1	706.30	218.3	756.02	128.18	689.97	433	731.78	13.26	711.57	45.7
E030-04 s	718.25	208.1	718.24	1.8	718.25	3.3	760.07	156.18	695.71	423	695.71	0.00	695.71	38.2
E031-09 h	614.60	1302.7	610.00	2.0	610.00	4.2	610.06	70.07	610.23	334	610.92	1.84	610.00	17.2
E033-03 n	2316.56	2317.3	2306.04	86.9	2306.04	121.9	2378.53	185.79	2433.45	436	2439.80	13.52	2426.16	107.3
E033-04 g	1276.60	2121.3	1184.44	3600.2	1184.27	427.5	1318.86	295.60	1231.32	412	1295.47	15.76	1282.99	127.9
E033-05 s	1196.55	2916.4	1161.11	689.3	1149.92	656.9	1361.15	330.06	1177.54	425	1242.33	38.08	1201.23	99.8
E036-11 h	698.61	863.0	698.61	0.0	698.61	3.8	698.42	98.12	698.61	312	705.19	5.56	698.61	14.7
E041-14 h	906.42	753.2	861.79	1.2	861.79	11.2	861.57	84.05	861.79	312	875.87	3.43	871.64	21.7
E045-04 f	1124.33	2198.9	1078.41	2030.8	1092.01	1180.9	1149.08	606.64	1093.21	434	1125.16	6.82	1115.97	210.5
E051-05 e	680.29	1390.3	658.34	3429.6	656.96	1216.0	694.38	1280.08	700.21	447	696.84	7.72	692.00	220.4
E072-04 f	529.00	7007.5	503.30	1469.7	503.90	2574.5	544.31	1464.27	566.70	580	534.24	3.81	532.16	576.4
E076-07 s	1004.40	6262.5	921.25	4697.4	923.74	2402.9	1033.31	1700.46	1040.67	698	989.75	8.17	980.85	469.4
E076-08 s	1068.96	2078.7	1009.45	3348.3	1001.63	2184.5	1108.33	1180.41	1090.96	721	1101.56	9.26	1091.33	802.3
E076-10 e	1012.51	4314.1	976.46	1889.1	959.49	1353.7	1034.40	1353.24	1035.70	758	1024.18	6.54	1019.87	852.2
E076-14 s	1063.61	1052.5	1047.75	682.8	1035.80	1228.9	1109.59	1090.0	1083.37	781	1069.21	11.58	1054.40	968.6
E101-08 e	1371.32	500.9	1219.77	4658.4	1166.99	3256.1	1337.55	2435.5	1404.60	908	1289.46	16.41	1278.25	1635.8
E101-10 c	1557.12	1075.0	1393.76	3066.6	1353.48	2573.5	1496.75	2815.44	1430.18	1247	1500.11	12.53	1483.92	2807.1
E101-14 s	1378.52	3983.2	1304.82	2422.3	1285.70	2610.1	1418.53	3553.39	1484.88	1085	1413.70	14.28	1396.57	2646.1
Average	876.31	1567.4	841.96	1522.6	837.72	816.5	892.01	706.82	877.69	522.3	879.98	7.70	871.98	442.8

* CG execution times are based on an extended set of constraints, which may affect direct comparability. ** for all instances indicative purpose only.

Table 6

CVRP cost in each of the MIA-RP iterations for day 1 with different virtual volume increase per items values.

Iteration Number	0	1	2	3	4	5	6	7	8
Virtual volume increase per items ($\lambda$) of 5%:
CVRP cost (km)	281.66	282.03	280.65	282.45	285.29	282.33	283.65	286.08	286.88
Number of unplaced packets	35	7	27	29	12	9	18	5	0
Number of drives	3	3	3	3	3	3	3	3	3
Virtual volume increase per items ($\lambda$) of 10%:
CVRP cost (km)	281.64	280.37	282.74	283.41	-	-	-	-	-
Number of unplaced packets	35	25	6	0	-	-	-	-	-
Number of drives	3	3	3	3	-	-	-	-	-
Virtual volume increase per items ($\lambda$) of 20%:
CVRP cost (km)	281.64	281.84	283.80	311.72	-	-	-	-	-
Number of unplaced packets	35	20	5	0	-	-	-	-	-
Number of drives	3	3	3	4	-	-	-	-	-

Table 7

Comparison of different optimization procedures on the real-world case study of 5 days.

Day No.	Num. of Packages	Num. of Delivery Points	Num. of Vehicles	Virtual Vol. Increase per Items ($\mathit{\lambda }$)	First Iteration CVRP Solution Cost	Num. of Unplaced Packets	CVRP+3D-BPP Cost	Num. of MIA-RP Iterations	CVRP Cost After MIA-RP	Industry Partner Cost
1	335	61	4	5%	281.64 km3 drives	35	331.38 km4 drives	8	286.88 km3 drives	391.69 km4 drives
1	335	61	4	10%	281.64 km3 drives	35	331.38 km4 drives	3	283.41 km3 drives	391.69 km4 drives
1	335	61	4	20%	281.64 km3 drives	35	331,38 km 4 drives	3	311.72 km4 drives	391.69 km4 drives
2	445	86	5	5%	509.29 km3 drives	96	628.37 km4 drives	14	562.08 km 5 drives	681.00 km5 drives
2	445	86	5	10%	509.29 km3 drives	96	628.37 km4 drives	8	557.64 km5 drives	681.00 km5 drives
2	445	86	5	20%	509.29 km3 drives	96	628.37 km4 drives	6	579.37 km 6 drives	681.00 km5 drives
3	437	82	5	5%	360.44 km 4 drives	113	523.32 km 5 drives	15	411.60 km 5 drives	509.26 km5 drives
3	437	82	5	10%	360.44 km 4 drives	113	523.32 km5 drives	8	425.31 km 5 drives	509.26 km5 drives
3	437	82	5	20%	360.44 km 4 drives	113	523.32 km5 drives	4	416.82 km 5 drives	509.26 km 5 drives
4	454	83	4	5%	393.27 km 4 drives	110	586.09 km 5 drives	18	480.00 km 6 drives	478.30 km4 drives
4	454	83	4	10%	393.27 km 4 drives	110	586.09 km 5 drives	11	481.09 km5 drives	478.30 km4 drives
4	454	83	4	20%	393.27 km 4 drives	110	586.09 km 5 drives	4	478.51 km 5 drives	478.30 km 4 drives
5	704	106	5	5%	399.90 km4 drives	161	589.42 km 6 drives	13	492.56 km6 drives	543.33 km 5 drives
5	704	106	5	10%	399.90 km 4 drives	161	589.42 km 6 drives	8	521.70 km 7 drives	543.33 km5 drives
5	704	106	5	20%	399.90 km 4 drives	161	589.42 km 6 drives	5	512.59 km 7 drives	543.33 km 5 drives

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/logistics9020070/s1. Table S1: The results of the hyperparameter grid search used for genetic algorithm.

References

1. Transport in the European Union—Current Trends and Issues. 2024; Available online: https://transport.ec.europa.eu/news-events/news/new-eu-transport-report-current-trends-and-issues-2024-06-27_en (accessed on 25 May 2025).

2. Seyyedhasani, H.; Dvorak, J.S. Using the Vehicle Routing Problem to reduce field completion times with multiple machines. Comput. Electron. Agric.; 2017; 134, pp. 142-150. [DOI: https://dx.doi.org/10.1016/j.compag.2016.11.010]

3. Zhang, Z.; Wei, L.; Lim, A. An evolutionary local search for the capacitated vehicle routing problem minimizing fuel consumption under three-dimensional loading constraints. Transp. Res. Part B Methodol.; 2015; 82, pp. 20-35. [DOI: https://dx.doi.org/10.1016/j.trb.2015.10.001]

4. Nolz, P.C.; Absi, N.; Feillet, D.; Seragiotto, C. The consistent electric-Vehicle routing problem with backhauls and charging management. Eur. J. Oper. Res.; 2022; 302, pp. 700-716. [DOI: https://dx.doi.org/10.1016/j.ejor.2022.01.024]

5. Laporte, H.; Mercure, G.L.; Nobert, Y. An exact algorithm for the asymmetrical capacitated vehicle routing problem. Networks; 1986; 16, pp. 33-46. [DOI: https://dx.doi.org/10.1002/net.3230160104]

6. Rochat, Y.; Taillard, E. Probabilistic diversification and intensification in local search for vehicle routing. J. Heuristics; 1995; 1, pp. 147-167. [DOI: https://dx.doi.org/10.1007/BF02430370]

7. Gmira, M.; Gendreau, M.; Lodi, A.; Potvin, J.Y. Tabu search for the time-dependent vehicle routing problem with time windows on a road network. Eur. J. Oper. Res.; 2021; 288, pp. 129-140. [DOI: https://dx.doi.org/10.1016/j.ejor.2020.05.041]

8. Sze, J.F.; Salhi, S.; Wassan, N. An adaptive variable neighbourhood search approach for the dynamic vehicle routing problem. Comput. Oper. Res.; 2024; 164, 106531. [DOI: https://dx.doi.org/10.1016/j.cor.2024.106531]

9. Saker, A.; Eltawil, A.; Ali, I. Adaptive Large Neighborhood Search Metaheuristic for the Capacitated Vehicle Routing Problem with Parcel Lockers. Logistics; 2023; 7, 72. [DOI: https://dx.doi.org/10.3390/logistics7040072]

10. Chen, C.M.; Lv, S.; Ning, J.; Wu, J.M.T. A Genetic Algorithm for the Waitable Time-Varying Multi-Depot Green Vehicle Routing Problem. Symmetry; 2023; 15, 124. [DOI: https://dx.doi.org/10.3390/sym15010124]

11. Teng, Y.; Chen, J.; Zhang, S.; Wang, J.; Zhang, Z. Solving dynamic vehicle routing problem with time windows by ant colony system with bipartite graph matching. Egypt. Inform. J.; 2024; 25, 100421. [DOI: https://dx.doi.org/10.1016/j.eij.2023.100421]

12. Vinyals, O.; Fortunato, M.; Jaitly, N. Pointer networks. Adv. Neural Inf. Process. Syst.; 2015; 28, pp. 2692-2700.

13. Pan, W.; Liu, S.Q. Deep reinforcement learning for the dynamic and uncertain vehicle routing problem. Appl. Intell.; 2023; 53, pp. 405-422. [DOI: https://dx.doi.org/10.1007/s10489-022-03456-w]

14. Dulac-Arnold, G.; Levine, N.; Mankowitz, D.J.; Li, J.; Paduraru, C.; Gowal, S.; Hester, T. Challenges of real-world reinforcement learning: Defnitions, benchmarks and analysis. Mach. Learn.; 2021; 110, pp. 2419-2468. [DOI: https://dx.doi.org/10.1007/s10994-021-05961-4]

15. Christofides, N. The vehicle routing problem. Combinatorial Optimization; Christofides, N.; Mingozzi, A.; Toth, P.; Sandi, C. Wiley: Chichester, UK, 1979; pp. 315-338.

16. Tarantilis, C.D. Solving the vehicle routing problem with adaptive memory programming methodology. Comput. Oper. Res.; 2005; 32, pp. 2309-2327. [DOI: https://dx.doi.org/10.1016/j.cor.2004.03.005]

17. Li, X.; Leung, S.C.H.; Tian, P. A multistart adaptive memory-based tabu search algorithm for the heterogeneous fixed fleet open vehicle routing problem. Expert Syst. Appl.; 2012; 39, pp. 365-374. [DOI: https://dx.doi.org/10.1016/j.eswa.2011.07.025]

18. Simeonova, L.; Wassan, N.; Salhi, S.; Nagy, G. The heterogeneous fleet vehicle routing problem with light loads and overtime: Formulation and population variable neighbourhood search with adaptive memory. Expert Syst. Appl.; 2018; 114, pp. 183-195. [DOI: https://dx.doi.org/10.1016/j.eswa.2018.07.034]

19. Nikolopoulou, A.I.; Repoussis, P.P.; Tarantilis, C.D.; Zachariadis, E.E. Adaptive memory programming for the many-to-many vehicle routing problem with cross-docking. Oper. Res.; 2019; 19, pp. 1-38. [DOI: https://dx.doi.org/10.1007/s12351-016-0278-1]

20. Morais, V.W.; Mateus, G.R.; Noronha, T.F. Iterated local search heuristics for the Vehicle Routing Problem with Cross-Docking. Expert Syst. Appl.; 2014; 41, pp. 7495-7506. [DOI: https://dx.doi.org/10.1016/j.eswa.2014.06.010]

21. Paquay, C.; Schyns, M.; Limbourg, S. A mixed integer programming formulation for the three-dimensional bin packing problem deriving from an air cargo application. Int. Trans. Oper. Res.; 2016; 23, pp. 187-213. [DOI: https://dx.doi.org/10.1111/itor.12111]

22. Paquay, C.; Schyns, M.; Limbourg, S. A tailored two-phase constructive heuristic for the three-dimensional Multiple Bin Size Bin Packing Problem with transportation constraints. Eur. J. Oper. Res.; 2018; 267, pp. 52-64. [DOI: https://dx.doi.org/10.1016/j.ejor.2017.11.010]

23. Lodi, A.; Martello, S.; Vigo, D. Heuristic algorithms for the three-dimensional bin packing problem. Eur. J. Oper. Res.; 2002; 141, pp. 410-420. [DOI: https://dx.doi.org/10.1016/S0377-2217(02)00134-0]

24. Crainic, T.G.; Perbolib, G.; Tadei, R. TS² PACK: A two-level tabu search for the three-dimensional bin packing problem. Eur. J. Oper. Res.; 2009; 195, pp. 744-760. [DOI: https://dx.doi.org/10.1016/j.ejor.2007.06.063]

25. Gonçalves, J.F.; Resende, M.G.C. A biased random key genetic algorithm for 2D and 3D bin packing problems. Int. J. Prod. Econ.; 2013; 145, pp. 500-510. [DOI: https://dx.doi.org/10.1016/j.ijpe.2013.04.019]

26. Kang, K.; Moon, I.; Wang, H. A hybrid genetic algorithm with a new packing strategy for the three-dimensional bin packing problem. Appl. Math. Comput.; 2012; 219, pp. 1287-1299. [DOI: https://dx.doi.org/10.1016/j.amc.2012.07.036]

27. Su, B.; Xie, N.; Yang, Y. Hybrid genetic algorithm based on bin packing strategy for the unrelated parallel workgroup scheduling problem. J. Intell. Manuf.; 2021; 32, pp. 957-969. [DOI: https://dx.doi.org/10.1007/s10845-020-01597-8]

28. Luo, Q.; Rao, Y.; Peng, D. GA and GWO algorithm for the special bin packing problem encountered in field of aircraft arrangement. Appl. Soft Comput.; 2022; 114, 108060. [DOI: https://dx.doi.org/10.1016/j.asoc.2021.108060]

29. Wang, P.; Rao, Y.; Luo, Q. An Effective Discrete Grey Wolf Optimization Algorithm for Solving the Packing Problem. IEEE Access; 2020; 8, pp. 115559-115571. [DOI: https://dx.doi.org/10.1109/ACCESS.2020.3004380]

30. Chen, M.; Li, K.; Zhang, D.; Zheng, L.; Fu, X. Hierarchical Search-Embedded Hybrid Heuristic Algorithm for Two-Dimensional Strip Packing Problem. IEEE Access; 2019; 7, pp. 179086-179103. [DOI: https://dx.doi.org/10.1109/ACCESS.2019.2953531]

31. Levine, J.; Ducatelle, F. Ant Colony Optimisation and Local Search for Bin Packing and Cutting Stock Problems. J. Oper. Res. Soc.; 2004; 55, pp. 705-716. [DOI: https://dx.doi.org/10.1057/palgrave.jors.2601771]

32. Guerriero, F.; Saccomanno, F.P. A hierarchical hyper-heuristic for the bin packing problem. Soft Comput.; 2023; 27, pp. 12997-13010. [DOI: https://dx.doi.org/10.1007/s00500-022-07118-4]

33. Verma, R.; Singhal, A.; Khadilkar, H.; Basumatary, A.; Nayak, S.; Singh, H.V.; Kumar, S.; Sinha, R. A Generalized Reinforcement Learning Algorithm for Online 3D Bin-Packing. arXiv; 2020; arXiv: 2007.00463

34. Zhang, J.; Zi, B.; Ge, X. Attend2Pack: Bin Packing through Deep Reinforcement Learning with Attention. arXiv; 2021; arXiv: 2107.04333

35. Tian, R.; Kang, C.; Bi, J.; Ma, Z.; Liu, Y.; Yang, S.; Li, F. Learning to multi-vehicle cooperative bin packing problem via sequence-to-sequence policy network with deep reinforcement learning model. Comput. Ind. Eng.; 2023; 177, 108998. [DOI: https://dx.doi.org/10.1016/j.cie.2023.108998]

36. Gendreau, M.; Iori, M.; Laporte, G.; Martello, S. A tabu search algorithm for a routing and container loading problem. Transp. Sci.; 2006; 40, pp. 342-350. [DOI: https://dx.doi.org/10.1287/trsc.1050.0145]

37. Lacomme, P.; Toussaint, H.; Duhamel, C. A GRASP × ELS for the vehicle routing problem with basic three-dimensional loading constraints. Eng. Appl. Artif. Intell.; 2013; 26, pp. 1795-1810. [DOI: https://dx.doi.org/10.1016/j.engappai.2013.03.012]

38. Koch, H.; Bortfeldt, A.; Wäscher, G. A hybrid algorithm for the vehicle routing problem with backhauls, time windows and three-dimensional loading constraints. Spectr.; 2018; 40, pp. 1029-1075. [DOI: https://dx.doi.org/10.1007/s00291-018-0506-6]

39. Bortfeldt, A.; Yi, J. The Split Delivery Vehicle Routing Problem with three-dimensional loading constraints. Eur. J. Oper. Res.; 2020; 282, pp. 545-558. [DOI: https://dx.doi.org/10.1016/j.ejor.2019.09.024]

40. Rodríguez, D.A.A.; Martínez, D.A.; Escobar, J.W. A hybrid matheuristic approach for the vehicle routing problem with three-dimensional loading constraints. Int. J. Ind. Eng. Comput.; 2022; 13, pp. 421-434. [DOI: https://dx.doi.org/10.5267/j.ijiec.2022.1.002]

41. Mahvash, B.; Awasthi, A.; Chauhan, S. A column generation based heuristic for the capacitated vehicle routing problem with threedimensional loading constraints. Int. J. Prod. Res.; 2017; 55, pp. 1730-1747. [DOI: https://dx.doi.org/10.1080/00207543.2016.1231940]

42. Küçük, M.; Yildiz, S.T. Constraint programming-based solution approaches for three-dimensional loading capacitated vehicle routing problems. Comput. Ind. Eng.; 2022; 171, 108505. [DOI: https://dx.doi.org/10.1016/j.cie.2022.108505]

43. Sitompul, C.; Horas, O.M. A Vehicle Routing Problem with Time Windows Subject to the Constraint of Vehicles and Good’s Dimensions. Int. J. Technol.; 2021; 12, pp. 865-875. [DOI: https://dx.doi.org/10.14716/ijtech.v12i4.4294]

44. Chi, J.; He, S. Pickup capacitated vehicle routing problem with three-dimensional loading constraints: Model and algorithms. Transp. Res. Part E Logist. Transp. Rev.; 2023; 176, 103208. [DOI: https://dx.doi.org/10.1016/j.tre.2023.103208]

45. Moura, A.; Pinto, T.; Alves, C.; de Carvalho, J.V. A Matheuristic Approach to the Integration of Three-Dimensional Bin Packing Problem and Vehicle Routing Problem with Simultaneous Delivery and Pickup. Mathematics; 2023; 11, 713. [DOI: https://dx.doi.org/10.3390/math11030713]

46. Perić, N.; Begović, S.; Lešić, V. Optimization of Heterogeneous Last-Mile Delivery of Fresh Products Considering Traffic Congestions and Other Real-World Parameters. IEEE Access; 2025; 13, pp. 87193-87207. [DOI: https://dx.doi.org/10.1109/ACCESS.2025.3569578]

47. Clarke, G.; Wright, J.W. Scheduling of Vehicles from a Central Depot to a Number of Delivery Points. Oper. Res.; 1964; 12, pp. 519-643. [DOI: https://dx.doi.org/10.1287/opre.12.4.568]

48. Gaskell, T.J. Bases for vehicle fleet scheduling. Oper. Res. Q.; 1967; 18, pp. 281-295. [DOI: https://dx.doi.org/10.1057/jors.1967.44]

49. Raidl, G.R.; Kodydek, G. Genetic Algorithms for the Multiple Container Packing Problem. Proceedings of the International Conference on Parallel Problem Solving from Nature; Amsterdam, The Netherlands, 27–30 September 1998; Volume 1498, pp. 875-884.

50. Chu, P.C.; Beasley, J.E. A Genetic Algorithm for the Multidimensional Knapsack Problem. J. Heuristics; 1998; 4, pp. 63-86. [DOI: https://dx.doi.org/10.1023/A:1009642405419]

51. Hinterding, R. Mapping, Order-Independent Genes and the Knapsack Problem. Proceedings of the 1st IEEE Conference on Evolutionary Computation; Orlando, FL, USA, 27–29 June 1994; pp. 13-17.

52. Bean, J.C. Genetic Algorithms and Random Keys for Sequencing and Optimization. ORSA J. Comput.; 1994; 6, pp. 154-160. [DOI: https://dx.doi.org/10.1287/ijoc.6.2.154]