This research proposes a DMH framework for solving the FVRP. The core idea of this framework is to decompose the primary large-scale problem into several smaller SPs using a clustering algorithm and then solve each SP exactly and independently. In the initial step, the set of customer nodes is partitioned into multiple clusters using a hierarchical clustering method. Unlike conventional clustering approaches, this method imposes an explicit constraint on the number of nodes in each cluster to control the computational time required to solve the resulting SPs. Specifically, the maximum cluster size is set to 15 nodes, based on experimental results indicating that the exact solution for SPs of this size can be achieved in less than one minute. As cluster size increases, the solution time rises significantly, reducing the time efficiency of the overall algorithm. Following the clustering phase, each SP is solved independently using an exact solver to determine its optimal solution. Finally, the solutions obtained from the SPs are interconnected through an integration process to generate the final answer for the primary problem.

The proposed approach leverages the advantages of both exact and heuristic algorithms. On one hand, the use of heuristic methods effectively reduces the overall solution time, while on the other, the quality of the solution is nearly guaranteed by solving the SPs exactly. This framework is particularly efficient and effective for large-scale problems that cannot be solved using traditional exact methods within a reasonable timeframe.

To integrate the SP solutions, the heuristic employed in this research is the cheapest insertion algorithm proposed by Renaud et al. [13] for determining the sequence of cluster visits. In this method, a representative point is first selected for each cluster. In the present study, this point is defined as the customer node closest to the cluster centroid. The algorithm starts from an initial cluster and then sequentially inserts the remaining clusters into the route, thereby determining the order in which the clusters are visited. At each step, the next cluster is selected based on the smallest increase it causes in the total route length. Once the final cluster sequence is obtained, the optimal routes generated for the individual SPs are connected according to this sequence. Owing to its low computational cost, structured procedure, and short execution time, this method provides an efficient mechanism for integrating partial solutions into a complete final solution.

Solving a VRP can be computationally demanding due to the large number of SECs. Therefore, this study proposes a novel iterative algorithm for solving VRP instances based on Selective Subtour Elimination Constraints (SSECs).

The main concept is to avoid the significant computational complexity that arises from introducing all SECs at the beginning. Instead, only the necessary constraints required at each step are added to the model. Therefore, the model is initially solved without considering any SECs. In the next step, all existing subtours in this solution are identified, and only the constraints corresponding to the elimination of those specific subtours are added to the model. The problem is then resolved with the updated model. This iteration persists until the resulting solution provides a complete tour without any subtours.

This approach significantly reduces the number of constraints in the problem, leading to increased computational efficiency and using less time and memory compared to adding all SECs at once. A comparison of the computational results between the proposed iterative approach and the method of adding all SECs simultaneously will be presented in the numerical results section.

To better illustrate the functionality of this framework, an example is provided below. As shown in Fig. 1.a, this example consists of a set of 28 customers symmetrically dispersed around the central depot and grouped into four clusters. Due to the high number of customers, obtaining an exact solution for this problem in a logical time frame is not feasible.

Hence, the framework developed in this research initially decomposes the problem into four SPs, as shown in Fig. 1.b, each of which can be solved exactly in a short time.

Customer distribution after clustering.

Accordingly, the truck and motorcycle simultaneously enter a cluster, and after servicing the customers within that cluster, they meet to allow the motorcycle to be reloaded by the truck. The delivery process then continues in the next cluster until all customers have been serviced, and both vehicles finally return to the central depot. The final solution for the primary problem is thus obtained, as shown in Fig. 2. The interpretation of the final solution is as follows: The truck and motorcycle start their route from the central depot. The truck services customers 5, 8, and 7 in sequence. Simultaneously, the motorcycle services customers 3, 4, 2, and 6. They then rendezvous at customer node 9, where the motorcycle receives four new parcels for delivery to customers 10, 11, 15, and 14, and they continue their routes. This process persists until all customers are serviced, and both vehicles ultimately return to the depot.

To better illustrate the functionality of the developed algorithm, its application to a small-scale example is described here. As shown in Fig. 3, a sample problem consisting of 15 customer nodes was considered, and the SSEC algorithm was used for routing parcel delivery to these customers. Using this approach, constraints that eliminate the subtours generated in the preceding solution step are iteratively added to the problem, thereby preventing the recurrence of those specific subtours.

By adopting this method, the algorithm successfully reached the optimal solution in only three stages (as shown in Fig. 4) by adding just 118 constraints to the model, and doing so in 29.13 seconds. In contrast, achieving the same result when all SCEs (16,383) are included in the problem from the outset requires 371.19 seconds. These results confirm that the proposed method significantly reduces problem complexity and, consequently, the solution time, without generating the entire set of SECs.

This section presents the results of implementing the DMH-FVRP framework on standard sample problems derived from the OR-Library repository1. Specifically, the study utilizes the standard problems designed by Solomon in 1987, found within the vehicle routing section, using the VRP_06 dataset.

The 25 test problems examined range in size from 5 to 50 customers, all of whom must be serviced within an eight-hour work shift by a single motorcycle and a single truck acting as a feeder vehicle, after which both vehicles must return to the central depot. The distances between nodes are calculated using Euclidean distances.

The parameters for motorcycle capacity, service time, and reloading time are assumed to follow a uniform distribution in the ranges of [3-6], [1-4], and [2-8], respectively. The proposed algorithms were coded using the Python programming language, and the exact solution of the SPs was obtained by the CPLEX solver. The hardware used in these calculations includes 20 GB of DDR4 RAM and an Intel Core i5-8850U 2.3GHz processor. The results of these computations are summarized in Table 1.

1 https://people.brunel.ac.uk/~mastjjb/jeb/info.html Table 1 presents the computational results for instances with different numbers of customers. The first column reports the number of customers. Columns two to five provide the objective function values and solution times for the FVRP scenario, obtained using the DMH-FVRP algorithm and the exact algorithm. Similarly, columns six to nine present the objective function values and solution times for the VRP scenario, obtained using the SSEC algorithm and the exact algorithm. The tenth column reports the difference between the VRP and FVRP solutions, while the final column shows the percentage deviation of the DMH-FVRP solution from the optimal FVRP solution.

As observed across all sample problems, the objective function value for the traditional VRP delivery plan is higher than that of the FVRP delivery plan. The lower objective function value of the FVRP delivery plan is because the motorcycle has a higher speed compared to the truck, allowing it to traverse the same distance in a shorter time. It is also noted that in FVRP instances with fewer than 15 customers, the solutions obtained by the exact algorithm and the DMH-FVRP framework are identical, as the clustering process places all customers into a single cluster, thus yielding the same result for both problems. However, as the number of customers increases and clustering becomes necessary, the solution deviates from the global optimum. Furthermore, the computation time for the VRP instances is generally shorter compared to the FVRP, which is justifiable given the higher complexity of the FVRP problem and the constraints associated with the synchronization of the two vehicle types.

The results indicate that the SSEC method makes it possible to obtain the optimal VRP solution for medium-sized problems in a short time, whereas calculating this response using the exact solver is only feasible for small-scale problems. Furthermore, based on the exact solutions computed for the small-sized FVRP instances, an average difference of 8.16% is observed between the exact solution and the solution derived from the DMH-FVRP framework. This value confirms the acceptable performance of the developed algorithm in solving the FVRP.

Finally, the comparison between the VRP and FVRP solutions shows an average improvement of 32.51% in the objective function value, which demonstrates the superiority of the synchronized truck-and-motorcycle delivery strategy over the traditional method of parcel delivery. This substantial improvement is achieved while the DMH-FVRP framework requires only 5.62% more computational time to obtain the FVRP solution than is required for the VRP solution.

I: Set of all customer nodes.

N: Set of all customer nodes including the depot.

(i,j): Indices corresponding to nodes.

c: Index for the truck.

m: Index for the motorcycle.

Tij c: Truck travel time between nodes i and j.

Sic: Service time for the customer i by the truck.

Qc: Truck capacity.

T: Maximum time of one work shift.

M: A large number.

Tij m: Motorcycle travel time between nodes i and j.

Sim: Service time for the customer i by the motorcycle.

Qm: Motorcycle capacity.

cc: Motorcycle reloading time.

xij: Binary variable, equals 1 if edge ij is traveled by the truck, 0 otherwise.

wi: Binary variable, equals 1 if node i is a joint node, 0 otherwise.

fi: Motorcycle inventory level before visiting the customer i.

bi: Motorcycle cumulative working time before visiting the customer i.

yij: Binary variable, equals 1 if edge ij is traveled by the motorcycle, 0 otherwise.

wij = fi ∗ wi ∗ yij: Quantity of inventory transferred from the truck to the motorcycle at the joint node i en route to the destination.

gi: Truck inventory level before visiting the customer i.

ai: Truck cumulative working time before visiting customer i.

min∑ ∑ (tc ij + sjc) xij + ∑ ∑ (tm ij + sjm) yij + cc ∑ wi. (1)

i∈N j∈N:i≠j i∈N j∈N:i≠j i∈N

∑ xij ≤ 1 ∀ i ∈ N. (2)

j∈I:j≠i ∑ xil = ∑ xlj, ∀ j ∈ N. (3)

i∈N:i≠l j∈N:l≠j

∑ yij ≤ 1, ∀ i ∈ N. (4)

j∈I:j≠i

∑ yil = ∑ ylj , ∀ j ∈ I. (5)

i∈N:i≠l j∈N:l≠j

∑ xij ≥ wj, ∀ j ∈ N. (6)

i∈N:i≠j

∑ (xij + yij) ≤ 1 + wj, ∀ i ∈ N, j ∈ I ∶ i ≠ j. (7)

i∈N:i≠j

gi ≤ Qc, ∀ i ∈ N. (8) fi ≤ Qm, ∀ i ∈ N. (9) fj ≥ fi − (yij − wi)qi − Qm(1 − yij + wi + wj), ∀ j ∈ N. (10) fj ≤ fi − (yij − wi)qi + Qm(1 − yij + wi + wj), ∀ j ∈ N. (11) fj ≥ qjyij − Qmwj, ∀ j ∈ N. (2) fj ≤ Qm(1 − wj), ∀ j ∈ N. (13)

∑ ∑ qjxij + ∑ ∑ qjyij − ∑ qjwj = ∑ qj , ∀ j ∈ N. (14)

i∈N j∈N:i≠j i∈N j∈N:i≠j j∈N j∈N

gj ≥ gi − xijqi − ∑ wyfil − M(1 − xij), ∀ i ∈ N, j ∈ I ∶ i ≠ j. (15)

gj ≤ gi − xijqi − ∑ wyfil + M(1 − xij), ∀ i ∈ N, j ∈ I ∶ i ≠ j. (16)

M ∑(xij + yij) ≥ qj, ∀ i ∈ N, j ∈ I ∶ i ≠ j. (17)

ai ≤ T, ∀ i ∈ N. (18) aj ≥ ai + tc ij + sic − 2T(1 − xij), ∀ j ∈ N. (19) bi ≤ T, ∀ i ∈ N. (20) bj ≥ bi + tm ij + sim + wi(cc − sim) − 2T(1 − yij), ∀ i ∈ N, j ∈ I ∶ i ≠ j. (21) bj ≥ aj − 2T(1 − wj), ∀ j ∈ N. (22) aj ≥ bi + tc ij + sic − 2T(1 − xij) − 2T(1 − wi), ∀ i ∈ N, j ∈ I ∶ i ≠ j. (23) wyfij ≥ fj − Qc(1 − wi) − Qc(1 − yij), ∀ i ∈ N, j ∈ I ∶ i ≠ j (24) wyfij ≤ fj − Qc(1 − wi) + Qc(1 − yij), ∀ i ∈ N, j ∈ I ∶ i ≠ j. (25) wyfij ≤ Qm(wi), ∀ i ∈ N, j ∈ I ∶ i ≠ j (26) wyfij ≤ Qm(yij) , ∀ i ∈ N, j ∈ I ∶ i ≠ j. (27) wyfij ≥ 0, ∀ i ∈ N, j ∈ I ∶ i ≠ j. (28) fi , gi , ai , bi ≥ 0 (29) xij , yij , wi ∈ {0,1}. (30)