A new two-phase heuristic for a problem of food distribution with compartmentalized trucks and trailers

Intl. Trans. in Op. Res. 0 (2021) 1–34 DOI: 10.1111/itor.13071 INTERNATIONAL TRANSACTIONS IN OPERATIONAL RESEARCH A new two-phase heuristic for a problem of food distribution with compartmentalized trucks and trailers Laura Davila-Pena , David R. Penas and Balbina Casas-Méndez∗ Group of Optimization Models, Decision, Statistics and Applications (MODESTYA), Department of Statistics, Mathematical Analysis and Optimization, Institute of Mathematics (IMAT), Faculty of Mathematics, Universidade de Santiago de Compostela, Santiago de Compostela, Spain E-mail: [email protected] [Davila-Pena]; [email protected] [R. Penas]; [email protected] [Casas-Méndez] Received 30 July 2020; received in revised form 13 September 2021; accepted 27 September 2021 Abstract This paper presents a new formulation for the routing problem in which the available fleet consists of trucks and trailers divided into compartments. Solving the model for large instances is computationally expensive. Therefore, we introduce and implemented a two-phase heuristic algorithm. In the first phase, an initial solution is generated through a constructive heuristic algorithm based on concepts from the classic Clarke–Wright algorithm. In the second phase, the initial solution is improved by an iterated tabu search metaheuristic. Our algorithm was tested on 21 instances that were converted from the classic truck and trailer routing problem. The results of our computational study prove the effectiveness of our proposal; the algorithm always finds a feasible solution, which in small-sized problems it is proven to be of good quality. In addition, the algorithm outperforms previous approaches for some truck and trailer routing problem instances. Furthermore, an application of the proposed model and heuristic is demonstrated in the field of agricultural logistics by comparing the obtained results. Keywords: truck and trailer routing problem; compartmentalized vehicles; construction heuristic algorithm; tabu search; logistics 1. Introduction In recent years, transport logistics has played a fundamental role in industry. Many public and private companies are interested in developing computational tools to design their routes, with objectives such as minimizing costs and/or maximizing the distribution of products. Thus, vehicle routing problems (VRPs) are a popular type of combinatorial optimization problem, through which transport routes for vehicles visiting a set of customers located at different places can be modeled. ∗Corresponding author. © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. 2L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 Fig. 1. Possible solution to the TTRP. Solving these mathematical optimization problems is a challenge for operational researchers. When new features from real-world applications are considered, such as capacitated vehicles, delivering in limited time windows, and stochastic behaviors, new variations of the original VRP arise, creating a need to develop new models and solution techniques. One promising modification of the VRP is the truck and trailer routing problem (TTRP) proposed by Chao (2002), which incorporates accessibility restrictions. In this variant, a fleet of trucks and trailers visits a set of customers, where some customers (vehicle customers; VCs) can be served by a complete vehicle (i.e., a truck pulling a trailer), while others are only reachable by a truck alone (truck customers; TCs). Examples of TCs are customers in inner-city areas, mountainous regions, or places where maneuvering or access with a trailer is not possible. To solve this problem, we distinguish three types of routes: pure truck routes (PTRs), which can be traveled only by trucks, pure vehicle routes (PVRs), which can be traveled entirely by a complete vehicle, and mixed vehicle routes (MVRs), which consist of a main tour traveled by a complete vehicle and one or more sub-tours traveled only by the truck part of the vehicle. Figure 1 illustrates a possible solution to the TTRP. Although this model can be very useful in many land-based logistical applications, the presence of three different types of routes makes solving the associated optimization problem more difficult, suiting it to the application of heuristics and metaheuristics, such as in Lin et al. (2009). Real-world applications include farm milk collection (Caramia and Guerriero, 2010b), delivery by feed mills (Lin et al., 2009), and the provisioning of infrastructure services in urban areas with accessibility restrictions (Parragh and Cordeau, 2017). Another interesting variation of the VRP is the multi-compartment case, where different products must be split into different storages during transport, making it challenging to maximize the use of vehicle capacity on the generated routes. Although the inclusion of compartments adds extra complexity, it can be a requirement in real logistical applications, as explained by Guitián de Frutos and Casas-Méndez (2019). Therefore, this paper proposes a novel mixed integer linear programming (MILP) approach to combine the TTRP with product compartmentalization, which we call the multi-compartment TTRP (MC-TTRP). The combination of these two features is motivated by the needs of a Spanish agricultural cooperative that produces feed for cattle, which we used to test, illustrate, and apply the proposed formulation. A tentative MC-TTRP was originally introduced in an unpublished preliminary work by Davila-Pena (2019). © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 3 Meanwhile, the effectiveness of the Clarke–Wright algorithm (Clarke and Wright, 1964) in building a solution for different VRPs and the requirement to solve MC-TTRPs that involve a relatively high number of customers served as motivation to modify this heuristic algorithm for the case of the MC-TTRP. To the best of our knowledge, only Derigs et al. (2013) reported adaptation of this algorithm to build an initial TTRP solution, although no details about such adaptation were provided. In addition to this constructive method, we propose a metaheuristic approach based on an iterated tabu search to improve the initial solution obtained. Both the constructive and improvement phases are integrated into a novel two-stage algorithm to solve the MC-TTRP. A corresponding computational study was conducted through a series of instances created from other existing ones in the literature, obtaining excellent results. However, these problems could not be benchmarked with the exact model due to the computational time required. In contrast, a series of small-sized real-world problems were solved, achieving solutions that are competitive and close to those provided by the exact method. The remainder of this paper is organized as follows. Section 2 reviews related work. Section 3 presents the current case study in detail and an in-depth description and formulation of the MCTTRP. Section 4 presents the two-stage heuristic to solve the MC-TTRP. Section 5 reports the computational results for the designed heuristic on the MC-TTRP, using a set of instances adapted from those in literature and data of a real application. Finally, Section 6 summarizes the main conclusions of our study. 2. Related work To solve the logistics of an agricultural cooperative that distributes feed for cattle to a large number of customers, many of whom have accessibility restrictions, the TTRP appears to be a satisfactory model. Although Chao (2002) introduced the term TTRP, previous works have incorporated trailers to solve similar case studies. The first approach could be that presented by Semet and Taillard (1993). These authors proposed a VRP that considered the use of trailers under accessibility restrictions. Semet (1995) proposed another example describing a new variant of the VRP formulated as an integer linear programming (ILP) problem called the partial accessibility constrained VRP (PACVRP). Despite being very similar to Chao’s TTRP, it has specific differences, such as the utilization of all available trucks. Other studies have considered a heterogeneous fleet of vehicles composed of trucks and trailers, such as the case of Gerdessen (1996), whose model is known as the VRP with trailers (VRPT). Moreover, Chao et al. (1998) studied the site-dependent VRP (SDVRP), where every customer has a specific type of vehicle assigned. Some seminal papers on the TTRP do not offer a mathematical formulation through an MILP model, although Scheuerer (2004) presented a formulation of the TTRP by Chao (2002). This turns out to be an adaptation of the proposal by Semet (1995) for the PACVRP and can be considered as the motivation for the current paper. Other researchers have built new models based on the proposal of Chao (2002) to meet various real-world requirements, such as a TTRP with time windows (TTRPTW) proposed by Lin et al. (2011). In the TTRPTW, besides its type and demand, each customer has three associated measurable times: the earliest and latest time of day at which it can be served and the service time required. Recently, Accorsi and Vigo (2020) considered a generalization of the TTRP, the extended single © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 4L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 TTRP (XSTTRP), which contains, all together, a variety of node types that were previously considered only separately: truck customers, vehicle customers with and without parking places, and parking-only locations. In the XSTTRP, a single vehicle, consisting of a truck and a detachable trailer, is used to serve a set of customers with known demand and accessibility constraints. Regarding TTRP solution methods, heuristics are popular approximation-based strategies for solving mediumto large-scale instances. In fact, heuristics have been used in the solutions of several VRP variants, with Lespay and Suchan (2021) being one of the most recent references. That study considered the problem of a food company’s distribution center. This was solved by constructing an initial solution, which was subsequently improved using a guided local search. Gerdessen (1996) proposed constructive and improvement heuristics for solving the VRPT. Semet (1995) described a two-stage heuristic method for obtaining PACVRP solutions: the first phase of the algorithm involves assigning trailers to trucks and determining the optimal allocation between customers and trucks/vehicles, and then, the second phase builds the routes. Other works, such as those by Chao (2002) and Scheuerer (2006), also proposed two-phase methods, where they first defined an initial solution by applying constructive procedures and then used improvement metaheuristics based on techniques such as tabu search (Glover and Laguna, 1998). Later, Caramia and Guerriero (2010a) combined a mathematical programming and local search approach to solve the TTRP, and they compared their results with those of Chao (2002) using a set of benchmarks. Furthermore, the TTRP can be addressed using a metaheuristic approach, as in the work by Lin et al. (2011), where a simulated annealing algorithm was designed to find approximate TTRP solutions according to given time windows, achieving improved results in 11 of the 21 instances of Chao (2002). In addition, in an original research, Derigs et al. (2013) analyzed different variants of the TTRP and proposed two-stage heuristics for solving these problems, starting by building an initial solution and then moving to an improvement phase combining techniques such as local search (LS) and large neighborhood search (LNS). The behavior of the heuristics created for the TTRPTW were compared with the heuristic proposed by Lin et al. (2011). Depending on the TTRP variant considered, the authors applied a specific construction heuristic and, among them, an adaptation of the Clarke–Wright savings algorithm stood out. In terms of exact solution methods, recent references include the paper by Parragh and Cordeau (2017), which proposes a branching and pricing algorithm for the TTRPTW. It adapts the LNS algorithm to obtain good initial columns. Compared with existing metaheuristic algorithms, such as those designed by Lin et al. (2011) and Derigs et al. (2013), they obtained highly competitive results. Some instances with up to 100 customers were optimally solved. Rothenbächer et al. (2018) also solved the TTRPTW exactly using a branching, bounding, and cutting algorithm. Their computational studies showed that their algorithm outperforms existing approaches on the TTRP and TTRPTW benchmark instances used in the literature. To solve the XSTTRP, Accorsi and Vigo (2020) developed a fast and efficient hybrid metaheuristic based on a four-phase solution approach, in which the main improvement phase consists of an iterated local search. Another challenging variation of the VRP arises when customers demand various types of products that cannot be mixed. This is the case in the multi-compartment VRP (MC-VRP), which was initially presented in Brown and Graves (1981) and Brown et al. (1987), whose objective was the distribution of petroleum products in the United States. Concerning the solving of multi-compartment problems, different approaches have been followed in recent years based on heuristics and metaheuristics. Simple constructive algorithms, such as the © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 5 Clarke–Wright algorithm, have also been successfully adapted in this context, as can be seen in the literature. El Fallahi et al. (2008) compared a constructive algorithm, memetic algorithm, and tabu search, concluding that the results provided by the tabu search were slightly better, although it required more computation time. Muyldermans and Pang (2010) used the Clarke–Wright savings algorithm to obtain a feasible initial solution. Subsequently, they performed a local search with movements taken from the literature and improved the quality of the solution previously obtained through a metaheuristic based on a guided local search. They performed a sensitivity analysis on certain parameters (number of customers and their demands, depot location, vehicle capacity, or number of products). Their computational study included a comparison with the work of El Fallahi et al. (2008). Derigs et al. (2011) considered a model with a homogeneous fleet, that is, all vehicles have the same number of compartments, all with equal capacities. This problem is a particular case of that addressed in the current paper. They implemented their own benchmarks and a collection of optimization methods capable of obtaining high-quality solutions, which covered a wide range of alternative approaches to construction, such as LS, LNS, and metaheuristics. Coelho and Laporte (2015) defined and compared four categories of multi-compartment problems. They proposed two formulations for each case and presented a branching and cutting algorithm to solve singleand multi-period cases containing up to 50 and 20 customers, respectively. Mendoza et al. (2010) extended the MC-VRP to the case in which the demands are stochastic, giving rise to the MCVRP with stochastic demands (MC-VRPSD). Mendoza et al. (2011) proposed a set of constructive heuristics to solve this problem, which included stochastic versions of the nearest neighbor, nearest insertion, and savings-based approaches, adapted to the multi-compartment scenario. Among the most recent investigations of solution methods for the MC-VRP, we highlight those by Henke et al. (2015, 2019). Starting from a real problem of collecting glass containers, a model formulation and branch-and-cut algorithm for solving the problem to optimality were presented. The performance of the proposed algorithm was evaluated through extensive numerical experiments. Furthermore, the economic benefits of introducing compartments to vehicles were investigated. Silvestrin and Ritt (2017) proposed a tabu search heuristic algorithm and integrated it with an iterated local search to solve the MC-VRP. In several experiments, they analyzed the performance of the algorithm and compared it with results in the literature, finding that it produces better solutions than those provided by other existing heuristic algorithms. They considered an initial solution obtained by the Clarke–Wright savings algorithm extended to handle multiple compartments. Metaheuristics based on iterated local searches have shown very good behavior in various VRP variants (cf. Alvarez et al., 2018, who proposed efficient metaheuristics based on iterated local search and simulated annealing). Alinaghian and Shokouhi (2018) presented a new mathematical model for the multi-depot MC-VRP. They designed a hybrid algorithm composed of adaptive large neighborhood search (ALNS) and variable neighborhood search (VNS). The results were compared to the exact solutions of small instances and compared with each other in large instances. Ostermeier and Hübner (2018) proposed an MC-VRP with a fleet of vehicles with flexible compartments. The aim of their work was to demonstrate the benefits of considering a mixed fleet consisting of both singlecompartment and compartmentalized vehicles. The problem was solved using LNS. Ostermeier et al. (2021) introduced a typology for MC-VRPs and extensively reviewed the existing literature. They also made suggestions for future research. Finally, the work of Caramia and Guerriero (2010b) should be highlighted as the first (and to the best of our knowledge, the only) to consider the TTRP with compartments, which we refer © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 6L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 to as the MC-TTRP hereinafter. They investigated a VRP in which at most one type of product could be assigned to each compartment. Furthermore, they established the additional constraint that some delivery locations were small and inaccessible by large vehicles. Due to the similarity between the problem addressed in that paper and the present one, it is considered convenient to point out the differences between the two studies. First, with regard to actual motivation, the problem analyzed by Caramia and Guerriero (2010b) was for milk collection on farms by an Italian company, while in our real-world case study, which will be described in more detail in the next section, the problem of the distribution of feed among members of a Spanish agricultural cooperative was analyzed. Regarding the model and methodology used, Caramia and Guerriero (2010b) proposed two mathematical programming models. One of them aimed to assign vehicles to farmers with the objective of minimizing the number of vehicles used, satisfying restrictions on capacity, demand, and types of milk. It should be noted that the group of farmers was divided into four zones, and an initial allocation of vehicles was made to each zone. The fleet considered was heterogeneous. The second model was used to minimize the lengths of the routes. In the proposed methodology, the possibility of serving VCs on sub-tours was not permitted. Accordingly, they used a two-phase heuristic. Such a process might result in no feasible solutions with respect to the times of work shifts. Therefore, a multiple-restart mechanism was implemented, and additional constraints, local search, and a tabu list were added to avoid cycling. Following a different approach, in our setup, there is a homogeneous truck fleet and a trailer fleet, and vehicle pre-assignments are not made to groups of customers. In this case, the sub-tours on an MVR, traveled by only a truck, can visit both TCs and VCs. In addition, what makes an important difference is that the formulation of a single model covering the whole problem is provided. We solved the model exactly for small-sized instances and then developed a two-phase heuristic, in which a generalization of the Clarke–Wright algorithm is used to find an initial solution that is then improved by a tabu search. The results of the heuristic were compared to the optimal solutions in problems where it was possible to do so, and a comprehensive computational study was conducted by creating MC-TTRP test problems for the heuristic. It is also worth mentioning that we studied the performance of our heuristic on TTRP instances, in addition to studying the scope of our model with instances built from real data. 3. Problem description and formulation 3.1. Case study The motivation for this study stemmed from the needs of a Spanish cooperative that produces and distributes feed for farm animals. The company is located in Galicia, a region in the northwestern Spain with an area of 29,565 km2spread over four provinces and 315 municipalities. The cooperative, which was created 16 years ago, currently has a total clientele of more than 1500 farmers distributed throughout the four provinces of Galicia (although not all of them order from the feed factory) and covering 60 municipalities across a large geographical area. The annual amount of feed produced exceeds 150,000 tons. The agricultural company produces different types of feed, and farmers usually place one or two orders per month. The number of daily orders is approximately 40, where each order ranges from 500 to more than 14,000 kg. The average number of annual orders per feed customer is © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 7 approximately 17. There are also occurrences such as the loss or incorporation of new customers. The roads leading to some of the farms or the farms themselves are inaccessible by large trailers. Moreover, customers sometimes request different types of feed because they have different species of animals. Naturally, goods that are not of the same type cannot be mixed. Thus, it is necessary to have compartmentalized vehicles. In addition to not being able to mix different kinds of feed in the same hopper, the same compartment cannot be used to supply two different customers because the cooperative does not have technology to measure out each customer’s supply from their vehicles. The purpose of this study was to provide a tool for the cooperative to automatically design routes for each vehicle such that their restrictions are met and the distance traveled is minimized. Each day, new orders may be received, trucks may experience breakdowns, and customers may change their demands at short notice. All these factors suggest that route planning is only useful within two or three days at most. A team of agricultural engineers designed a comprehensive global positioning system (GPS) that can monitor various vehicle routes. The GPS provides all the geographic information required to provide the data to solve the problem. We also know the capacity of each compartment, the demands of different customers, and whether a trailer can access each farm as well as its load restrictions. 3.2. Multi-compartment truck and trailer routing problem (MC-TTRP) As stated before, this paper proposes an MC-TTRP model—a novel MILP implementation of the TTRP with multi-compartmentalized vehicles. The MC-TTRP can be described as follows. Let G=(N,E) be an undirected, weighted graph consisting of a node set N={0,1,...,n}, representing the depot ({0}) and customers ({1,...,n}), and an arc set E={(i,j):i,j∈N,i= j}, representing the arcs that can be traveled between different nodes. N1and N2are subsets of Nthat contain the n1and n−n1VCs and TCs, respectively. A nonnegative cost cij,(i,j)∈E, is assigned to each arc, which represents the distance a vehicle must travel from ito j. Each node i∈Nrequires a service time si, which, in the case of the depot, refers to the time required to load the vehicles. For transportation, a set KT={1,...,mL,...,mT}of trucks and set KL={1,...,mL}of trailers are available. KT 1={1,...,mL}and KT 2={mL+1,...,mT} are the subsets of KTthat consist of trucks that can pull a trailer and pure trucks (without trailer attached), respectively. Note that |KT 1|=|KL|,thatis,therearemLcomplete vehicles. In addition, mLis the number of trailers, and mTis the number of trucks (mL≤mT). Complete vehicles and pure trucks are assumed to be homogeneous. Let QTbe the capacity of each truck and QLbe the capacity of each trailer. Hence, QT+QLis the capacity of a complete vehicle. As mentioned above and illustrated in Fig. 1, three different types of routes can appear in this variant of the TTRP. For MVRs, the complete vehicle leaves the depot and serves some VCs; this part of the MVR is known as the main tour. The main tour is entirely covered by a complete vehicle and starts and ends at the depot. During the tour of an MVR, it is possible to uncouple the trailer from the truck and leave it parked at one of the VC locations to start a sub-tour (or even at the depot, which is always a candidate for trailer parking places). VCs and TCs can be served in a subtour because they are performed by a pure truck. Sub-tours begin and end at the parking place (the depot or any of the VCs of the main tour), also known as the root of the sub-tour. There are no © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 8L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 restrictions on the number of sub-tours in an MVR or on the number of sub-tours that can start from the same VC on a given main tour, as long as the vehicle capacity restrictions are satisfied. That is, the demands transported on the MVR cannot exceed the capacity of a complete vehicle, QT+QL,andQTcannot be surpassed in a sub-tour. Another type of route is a PVR, which is fully traveled by a complete vehicle, implying that only VCs can be delivered to and their demands cannot exceed QT+QL. On the contrary, PTRs serve both types of customers because trucks travel without a trailer attached. The demand transported on a PTR cannot exceed QT. For the sake of simplicity, it is also assumed that travel costs are the same for all vehicles, regardless of whether a trailer is attached. Each trailer r∈KLis divided into a set of compartments, HL, where QH Lis the capacity of a trailer compartment. Similarly, each truck k∈KTis split into a set of compartments, HT,whereQH Tis the capacity of a truck compartment. Furthermore, let the set F={1,...,nF}of feed types be given. Each node (except for the depot) has a nonnegative demand dif (i∈N\{0},f∈F) for every feed type. The demands must be served at customers’ locations and transported from the depot without the feed types being mixed. In addition, products for different customers cannot be carried within the same compartment. The total demands of each customer must be met by the same vehicle, and it is possible to divide a customer’s demand for the same feed type among several compartments. Trucks and trailers have a maximum usage time allowed of D, an average speed of vm, and legal capacities, LTand LL, respectively, which may appear depending on regulations or laws in some specific areas. Regarding the decision variables involved in the model, xkr ij and yklv ij (both binary) are related to the construction of routes. The former involves routes covered by complete vehicles (MVRs or PVRs), and it takes a value of 1 if the complete vehicle consisting of truck k∈KT 1and trailer r∈KL travels from node ito j(i,j∈N1∪{0}); otherwise, it is 0. In contrast, yklv ij refers to routes covered only by trucks (PTRs or sub-tours of MVRs), and it takes a value of 1 if truck k∈KTtraverses the arc (i,j)onthevth route/sub-tour (v∈V={1,...,n})1with root l∈N1∪{0}; otherwise, it is 0. For l=0, the associated tour is a PTR (k∈KT 2) or a sub-tour in an MVR whose root is the depot (k∈KT 1). However, if l∈N1, then such a root refers to the VC of the main tour of an MVR working as a trailer parking place to start a sub-tour (and k∈KT 1). The remaining variables are related to the vehicle compartments. In particular, ZTk i,f,ht takes values in [0,1] and represents the proportion of compartment ht ∈HTof truck k∈KTcarrying feed f∈Ffor customer i∈N, whereas Uk i,f,ht is a binary variable equal to 1 if ZTk i,f,ht >0 and 0 otherwise. Analogously, ZLr i,f,hl takes values in [0,1] representing the proportion of compartment hl ∈HLof trailer r∈KLloaded with feed f∈Ffor customer i∈N1, while Vr i,f,hl is a binary variable equal to 1 if ZLr i,f,hl >0and 0 otherwise. The objective of the MC-TTRP is to determine a set of vehicle tours that minimizes the total cost of all edges to be traveled, that is, the total distance of the solution routes, such that all constraints are met, that is, the demands are satisfied, no vehicle capacities are exceeded, and the restrictions of access to customers and constraints related to the loading of compartments are considered. From now on, we will denote the set of VCs and the depot, N1∪{0},asN0 1, and we will denote the set of customers, N\{0},byN∗. Table 1 gives a summary of the sets, parameters, and decision variables involved in the model to facilitate better understanding of our proposal. Given this termi1Note that a root candidate l∈N1∪{0}can have as many sub-tours as there are customers, n. © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 9 Table 1 Notation of the proposed MC-TTRP Set or parameter Definition Parameter Definition {0}Depot nTotal number of customers N1SetofVCs n1Number of VCs N2SetofTCs mLNumber of trailers (and of complete vehicles) N={0}∪N1∪N2Set of nodes (customers and depot) mTNumber of trucks N∗=N\{0}Set of customers cij Distance a vehicle must travel from ito j N0 1=N1∪{0}Set of VCs and depot nFNumber of different types of feed KT 1Setoftrucksthatcanhitchatrailer dif Demand of customer ifor feed f KT 2Set of pure trucks QTCapacity of a truck KT=KT 1∪KT 2Setoftrucks QLCapacity of a trailer KLSetoftrailers QH TCapacity of a truck’s hopper FSet of different types of feed QH LCapacity of a trailer’s hopper HTSet of truck hoppers DMaximum time allowed to use a truck/vehicle HLSet of trailer hoppers vm Average speed of trucks/vehicles VSet of tours/sub-tours leaving a specific root siService time required for node i∈N LTLegal capacity of trucks LLLegal capacity of trailers Variable Definition xkr ij Binary variable equal to 1 if truck kwith trailer rpasses through arc (i,j); 0 otherwise yklv ij Binary variable equal to 1 if truck kpasses through arc (i,j) on route/sub-tour vwith parking place l;0 otherwise Uk i,f,ht Binary variable equal to 1 if compartment ht of truck kis loaded with feed ffor customer i; 0 otherwise Vr i,f,hl Binary variable equal to 1 if compartment hl of trailer ris loaded with feed ffor customer i; 0 otherwise ZTk i,f,ht Proportion of compartment ht of truck kloaded with feed ffor customer i ZLr i,f,hl Proportion of compartment hl of trailer rloaded with feed ffor customer i nology and notation, the objective function and constraints of the MC-TTRP can be formulated as follows: minimize:  i∈N0 1 j∈N0 1 k∈KT 1 r∈KL cijxkr ij + i∈N j∈N k∈KT l∈N0 1 v∈V cijyklv ij (1) subject to  i∈N0 1 k∈KT 1 r∈KL xkr ij + i∈N k∈KT l∈N0 1 l=j v∈V yklv ij =1,j∈N1;(2)  i∈N0 1 k∈KT 1 r∈KL xkr ij + i∈N k∈KT ykjv ij ≤2,j∈N1;v∈V;(3) © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 16 L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 Fig. 2. Calculation of sij for i,j∈N. Fig. 3. Calculation of ˆ sij for i∈N1and j∈N2. those tours directly connected to the depot, while ˆ Rconsiders the sub-tours whose roots are VCs. Consequently, as indicated in line 3 of Algorithm 1, we initialize both matrices as follows: R=⎛ ⎜ ⎜ ⎜ ⎝ 010 020 . . .. . .. . . 0n0 ⎞ ⎟ ⎟ ⎟ ⎠=ˆ R, assuming that routes (0,i,0) for i∈N1are traveled by complete vehicles and routes (0,j,0) for j∈N2are covered by trucks alone. Each customer’s corresponding row represents the previous and next customer in the transport network, as appropriate. Furthermore, at the beginning of the proposed algorithm, two saving matrices are created: Sand ˆ S. The entries of Sare computed as sij =ci0+c0j−cij (for i,j∈N,i= j)2, which represents the standard savings matrix. This scenario is illustrated in Fig. 2. However, because we have two different classes of customers and given that VCs can serve as trailer parking places, we also consider savings ˆ sij. These values are calculated as follows: •ˆ sij =sij for i,j∈N1or i,j∈N2. •ˆ sij =cj0+c0j−cij −cji for i∈N1and j∈N2. These savings are the result of removing route (0,j,0) and converting it into a sub-tour whose root is i. That is, we will have the route (0,...,i,j,i,...,0) instead of (0,j,0) and (0,...,i,...,0). Figure 3 illustrates this situation. •ˆ sij =ci0+c0i−cij −cji for j∈N1and i∈N2. These savings are analogous to the previous ones, swapping iand jin their type of customer. 2Moreover, sii =0foralli∈N. © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 17 Hence, matrix ˆ Swill have the following structure: ˆ S= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ s1,1··· s1,n1ˆ s1,n1+1··· ˆ s1,n . . ..... . .. . ..... . . sn1,1··· sn1,n1ˆ sn1,n1+1··· ˆ sn1,n ˆ sn1+1,1··· ˆ sn1+1,n1sn1+1,n1+1··· sn1+1,n . . ..... . .. . ..... . . ˆ sn,1··· ˆ sn,n1sn,n1+1··· sn,n ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . We must also consider the features related to compartments. As shown in line 5 of Algorithm 1, we create matrices HT rem,HL rem,andT. The former two are initialized as HTand HL, respectively, and are updated by setting the compartments that have already been used to −1. Matrix Tplays a fundamental role in this implementation because it contains all information about the vehicles’ loading procedures; it indicates the allocation details of each customer’s demand for its corresponding vehicle. Once all initial parameters have been created, our CW implementation follows an iterative process, where different routes are merged based on the maximum saving, Sm(lines 7–28 of Algorithm 1). Thereafter, riand rjare considered to be routes containing customers iand j, respectively. By checking matrix T(line 11), the algorithm knows which vehicles cover each of these routes; thus, we can determine the PTRs, PVRs, and MVRs. Therefore, our proposed method repeats a set of steps, while Sm, the maximum saving of Sand ˆ S, is positive (stop condition). Thereafter, each iteration of the algorithm has two discernible cases: Sm∈Sor its opposite, Sm∈ˆ S. The first case occurs when Smbelongs to S(lines 9–17 of the pseudocode). Our heuristic applies the classic routing merge of the CW, which considers the type of customers (VCs or TCs) involved in such saving, which is a decisive factor in the type of route generated. Thus, four different unions between routes (by connecting customer ito customer j) can be conducted, as long as iis the first customer of route riand jthe last one of route rj: 1. Merging two PVRs: If routes riand rjare PVRs covered by different complete vehicles, where the trucks are not yet loaded, the algorithm must count the number of unavailable compartments between both trailers. Two situations arise from this: (i) If the number does not exceed the number of trailer compartments, our method can move all goods to one trailer, emptying the other. (ii) If the number exceeds the number of trailer compartments, our method starts to load one of the trucks as long as the total number of available compartments is sufficient to meet the demands of both routes. Thus, it moves the goods from the other trailer to the chosen complete vehicle, giving preference to the filling of the trailer. 2. Merging two PTRs: When both routes riand rjare PTRs, our heuristic can merge them as long as the number of unavailable compartments between both trucks does not surpass the number of truck compartments. In this case, the goods are moved to one of the trucks, leaving the other free. 3. Combining an MVR with a PVR: If one of the routes is an MVR and the other is a PVR with an empty truck, say riand rj, respectively, then their union is possible when the demands already © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 18 L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 (a) Before (left) and after (right) merging when rjhas a single customer. (b) Before (left) and after (right) merging when rjhas multiple customer. Fig. 4. Merges performed when S∈Sm. met by route rjfit in the trailer of route ri. That is, when the trailer of route rihas a sufficient number of available compartments to move goods from the trailer of route rj. 4. Inserting route riin a PTR: If route ri=(0,i,0), where iis a VC, and route rjis a PTR, then merging both routes is possible when the number of available compartments in the truck is sufficient to hold i’s demands. For all the above considered fusions, updating matrices {R,ˆ R,S,ˆ S,HT rem,HL rem,T}is necessary to contemplate the changes conducted, as indicated in line 14. If merging cannot be completed, we must simply update matrices Sand ˆ Sby setting their corresponding entries to zero. The second case occurs when ˆ Scontains Sm, which means that the current saving arises from merging a PTR and PVR (obtaining a sub-tour as a result). Lines 18–27 of Algorithm 1 correspond to this situation. In such a case, it is clear that customers involved in Smhave different types. Let us see how to proceed when i∈N1and j∈N2, because the other case is analogous. Our method can only execute merges when riis a PVR with an empty truck and rjisaPTR.Withthisinmind,two possibilities arise: •Ifrj=(0,j,0), the algorithm uncouples the truck of route riand hitches the truck of route rj, loaded with j’s demands, to the trailer of route ri. This proceeds as illustrated in Fig. 4a. • Suppose jis the first customer of the route but not the only one. In this case, we must also consider the last customer of route rj,say,l. To join routes riand rj, it is necessary that c0j+cl0−cij −cli is positive. We illustrate this merging in Fig. 4b. 4.2. Refining the construction phase Once we have exited the main loop and no saving is positive, the algorithm must check whether the resulting route configuration is feasible. The following post-processing functions (lines 30–35 of the pseudocode) are applied in the order in which they are presented to verify such feasibility, correct the tours if needed, improve the quality of the results, and limit the vehicle fleet: •Add disconnected customers: This checks if there exists a route ri=(0,i,0), which means that customer ihas not yet been served. During the initialization of R, we mentioned that the starting routes should be covered by complete vehicles or trucks alone, depending on whether iisaVCor © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 19 (a) Convert a PTR into an MVR. (b) Merge a PVR into the main tour of an MVR. (c) Merge a PTR into the sub-tour of an MVR. Fig. 5. Adjusting vehicle fleet by creating or completing MVRs. TC, respectively. Considering this, we supply customer iwith its corresponding type of vehicle, contemplating the compartment load in matrix T. After that, we are able to add the route (0,i,0) to the solution. •Adjust vehicle fleet by creating or completing MVRs: The aim of this function is to use all the available trailers in the solution obtained by eliminating excess PTRs and PVRs, either by changing or merging them into MVRs, respectively. First, because all PTRs with at least a VC can be transformed into MVRs, we convert as many PTRs into MVRs as necessary to reach the required number of trailers. Subsequently, the lowest quality PVRs and PTRs are selected to be joined either to the main tour or to one of the sub-tours of an existing MVR, provided that the result is a feasible solution. This procedure is illustrated in Fig. 5. •Adjust vehicle fleet by destroying routes: It is assumed that the routes to be deleted or preserved have already been selected to adjust the fleet. Iteratively, an attempt is made to insert the first group into the second group. The method of merging these routes depends on their nature. For instance, a PTR could only be part of a sub-tour of an MVR or join another PTR. If there are no feasible insertions of a specific route to be deleted, it is split in half. The algorithm will try to add it again in the next iteration, repeating this process until there are no more residual routes or until they cannot be divided anymore. •Adjust vehicle fleet by switching residual routes: After applying the previous function, if there are still routes to be eliminated, they will consist of a single customer. This last procedure to adjust the fleet exchanges these residual customers for others that have less load on the routes to be preserved while aiming to keep the cost function from increasing as much as possible. Once the exchange is conducted, the customer to be inserted requires less space. In this manner, the algorithm aims to relocate it to one of the existing routes, repeating the process until there are no residual routes left. •Apply tour improvement: This function performs 2-opt, 3-opt, and 4-opt* moves on the resulting routes3, including the sub-tours, individually. For each route, the move that leads to the greatest cost reduction, if any, is applied. •Apply local search moves: Finally, a local search is applied to the solution obtained from Algorithm 1. This iteratively applies a set of small modifications to intensify the search on routes close to the current solution. To implement these moves, we were inspired by the work of Derigs et al. (2013), where a hybrid approach for solving the TTRP, combining local search and large neighborhood search, was presented. We adapted some of the moves implemented in their local 3The 4-opt* procedure uses a subset of potential 4-opt moves, cf. Renaud et al. (1996). © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 20 L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 (a) (b) (c) (d) (e) (f) Fig. 6. Local search moves for the MC-TTRP. (a) Replacing the single customer of a sub-tour with a VC, (b) replacing a sub-tour root with another VC, (c) replacing a main tour customer with a TC, (d) exchanging a PTR and a sub-tour and relocating the parking place, (e) moving a main tour customer to a new sub-tour, and (f) moving a sub-tour customer to the main tour and splitting the sub-tour. search to the MC-TTRP. Therefore, this function contains both specific TTRP moves as well as the standard 2-opt, 2-opt*, exchange, and relocate operators commonly used in many VRP implementations. Figure 6 shows some of the moves performed in this step. Finally, as shown in Algorithm 1, our heuristic for the MC-TTRP completes the work, returning a feasible route and the loading configuration of all the vehicles involved in the transport network. 4.3. Iterated tabu search As stated before, in the second phase of our proposal, an iterated tabu search (ITS) is implemented. We developed this metaheuristic based on previous related works, such as that by Cordeau and Maischberger (2012) for some VRP variants and the approach implemented by Silvestrin and Ritt (2017) for the MC-VRP. The ITS starts from a feasible solution, which, in our case, is the solution provided by the constructive heuristic. This procedure is usually based on two main actions: a tabu search and a perturbation. The former involves a set of local moves applied sequentially to improve the solution, with the ability to accept slightly worse solutions or to move through the infeasible region when the search is stuck. Moreover, to avoid returning to already-visited solutions, they use the short-term memory strategy known as the tabu list. In the case of the perturbation process, the objective is to escape from local optima by exploring the vicinity of the current solution through small changes in the routes. As a general overview, the ITS aims to improve the best-known solution by combining the diversification provided by the perturbation with the intensification and diversification produced by the tabu search. Algorithm 2 describes the scheme of our proposed ITS, of which we will now present a brief outline. The input parameters are the solution returned by the CW, which is the starting point, and a maximum number of iterations used as a stopping criterion. In the main loop (lines 5–14), the current solution is perturbed at each iteration after having applied the tour improvement function (presented in Subsection 4.2). This modified solution will be used as an initial guess for the tabu © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 21 Algorithm 2. Iterated tabu search for the MC-TTRP 1: procedure ITS_MC-TTRP (cw_solution,max_iterIT S) 2: best_solution ←cw_solution 3: current_solution ←cw_solution 4: iterIT S =0 5: while iterIT S <max_iterIT S do 6: current_solution_i ←Improve tours in current_solution 7: current_solution_i,φ←Perturbs the current_solution_i 8: current_solution_i,best_TS_solution ←TabuSearch(current_solution_i,φ,iterIT S,max_iterIT S) 9: if cost(best_TS_solution)<cost(best_solution)then 10: best_solution ←best_TS_solution 11: end if 12: with probability (iterIT S/max_iterIT S)2,current_solution ←best_solution 13: iterIT S =iterIT S +1 14: end while 15: return(best_solution) 16: end procedure search. In turn, as can be seen in line 8, this latter method returns its current solution, which will be the one to be perturbed in the next iteration, and its best solution, which can update the best overall solution (lines 9–11). In this way, throughout the iterations, the algorithm perturbs and searches the current solution, which is initially the best-known solution. However, as the algorithm progresses, there is a probability of working with another solution to avoid stagnation (line 12). Finally, the algorithm terminates when the stopping conditions are fulfilled, and it outputs the best solution found during the search. Regarding the perturbation procedure, this is a key point in our proposal because it guarantees diversity in the method. First, a random number, φ, is chosen between 1% and 10% of the total number of customers. Then, φcustomers are randomly selected and removed from the solution, to be subsequently reinserted. For the deletion and insertion operations, we use the generalized insertion procedure (GENI) and unstringing and stringing (US) algorithm proposed by Gendreau et al. (1992). This is followed by 3-opt and 4-opt* local search algorithms to improve the routes. Depending on the type of customer to be deleted or inserted, as well as its position in the solution, different possibilities of moves can arise. Some examples are shown in Fig. 7. Each customer is inserted into the route and position that minimizes the increase in the solution cost. If the perturbation failed to introduce the nodes that have been removed, we allow a threshold of infeasibility in the perturbed solution. After the solution is perturbed, the algorithm performs the tabu search, whose implementation requires an additional explanation. This is an iterative procedure, where the strategy is essentially to apply the best possible MC-TTRP local search moves (some of which are shown in Fig. 6) in the neighborhood of the current solution, as long as these moves are not stored in the tabu list. The approach adopted for the tabu list is as follows. As it is known, a move involves a set of routes and the customers to be deleted or inserted into them. At each iteration, the tabu list stores for each move applied to the current solution, a set of pairs {customi,rj}, with customi being the customer removed from route rj. In this way, the tabu search only accepts those moves where © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 22 L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 (a) (b) (c) (d) (e) (f) Fig. 7. Removal and insertion operations in the perturbation for the MC-TTRP. (a) Inserting a VC into a PVR, (b) inserting a TC into a PTR, (c) inserting a TC into a sub-tour, (d) converting a PVR into an MVR by creating a sub-tour with a single TC, (e) removing a VC from a PTR, and (f) removing a VC that works as parking place from an MVR. all the customers introduced in the routes are different from all the pairs {customi,rj}storedin the tabu list. The exception to this rule occurs when the result improves the best solution found so far, in which case an aspiration criterion is applied. Consequently, the tabu list prevents the insertion of customers into the same routes for a number τof iterations. Each pair {customi,rj} has an associated survival counter that controls the number of iterations that will be in the tabu list, and that requires to be updated at each iteration of the tabu search. When one of these counters reaches τ, its corresponding pair is removed from the tabu list. The parameter τis initialized at the beginning of the tabu search and is chosen randomly from the uniform distribution on the interval [1,√n·nr], where nis the number of total customers, and nr corresponds to the number of routes in the current solution. As the algorithm applies the best possible move, it may not improve the cost of the current solution or may even be infeasible. However, this is desirable because exploring new regions using worse or/and infeasible solutions prevents us from getting stuck in a local optimum. Therefore, to measure the total cost of each route, r, we consider the following objective function: F(r)=d(r)+αC(r),(50) where d(r) is the total distance traveled, αa parameter of penalty, and C(r) denotes the excess load on both vehicles and compartments. To calibrate the αpenalty, we follow the approximation proposed by Cordeau and Maischberger (2012), which we recommend referring to for further details about this mechanism. Briefly, the αpenalty is initially set to 1 and then updated throughout the tabu search, depending on the excess load in the current solution. In the case of feasibility, that is, if constraints are not violated, the penalty is decreased by a factor 1 +γ. Otherwise, if the current solution has excess load on its routes, αis increased by 1 +γ. The parameter γ is randomly selected at the beginning of the tabu search from the uniform distribution on [0,1]. Consequently, this update strategy produces an oscillatory effect between feasible and infeasible solutions. In addition, as in other tabu search methods such as Cordeau and Maischberger (2012) or Silvestrin and Ritt (2017), a table with the most frequent moves is also used to penalize the function © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 23 Algorithm 3. Tabu search for the MC-TTRP 1: procedure TabuSearch (current_solution_i,φ,iterIT S,max_iterIT S) 2: Choose τrandomly from U[1,√n·nr] and create a new tabuList ={} 3: Create a new freqPenList ={} and choose ζrandomly from U[0,1] 4: Initialize α=1 and choose γrandomly from U[0,1] 5: best_TS_solution ←current_solution_i 6: counter_iters_without_improvement =0 7: max_iters_without_improvement =(max_iterIT S −iterIT S)·φ 8: while counter_iters_without_improvement <max_iters_without_improvement do 9: Check all local search moves in the neighborhood of the current_solution_i 10: Evaluate the cost of the moves, penalizing infeasible solutions using α(see expression (50)) 11: If the moves do not produce an improvement, they are penalized using ζand freqPenList(see expression (51)) 12: Apply the best possible move in current_solution_i according to tabuList 13: if cost(current_solution_i)<cost(best_TS_solution)then 14: best_TS_solution ←current_solution_i 15: counter_iters_without_improvement =0 16: else 17: counter_iters_without_improvement =counter_iters_without_improvement +1 18: end if 19: if current_solution_i is infeasible then 20: α=α·(1 +γ) 21: else 22: α=α/(1 +γ) 23: end if 24: Save the selected move in freqPenList 25: Update tabuList 26: end while 27: return(current_solution_i,best_TS_solution) 28: end procedure cost when the search is stuck in a local minimum. Thus, when the best possible move decreases the current solution cost, we reevaluate all possible moves according to a new objective function: F(r,M)=F(r)1+ζ ci∈cMr freqPenList(ci,r)/i,(51) where Mis a specific move, ζis a penalty uniformly randomly chosen from [0,1], cMr is the set of customers involved in Mto be inserted into a route r,freqPenListis the table of frequencies, which returns how many times a customer ci entered in the route r,andiis the current iteration. Algorithm 3 shows the main scheme for the implemented tabu search. In the first lines of the pseudocode (lines 2–7), the tabu list (tabuList) and the table with the most frequent moves (freqPenList) are created, and all the parameters explained above are initialized (τ,α,γ,and ζ). Furthermore, a maximum number of iterations without improvement is set (line 7), which, in this case, will be the stopping criterion. During the main loop (lines 8–26), the tabu search evaluates all possible moves, selecting the best possible one according to the constraints imposed by the tabu list. If the current solution is a local optimum, the algorithm selects the best move © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 24 L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 considering freqPenList and expression (51). Once a new solution is created, we check if it increases the cost of the best one found during the tabu search (lines 13–18). Next, αis calibrated based on the solution feasibility (lines 19–23). Moreover, the algorithm adds the selected move information to freqPenList. The tabu list is updated in line 25 by reducing the survival counter of the stored pairs. Notably, those pairs that reached the iteration limit inside the list are then released, and new pairs involved in the last applied move are included. Finally, after the algorithm reaches the stop condition, the tabu search returns both the best solution visited during the search and the current solution. The latter is from which Algorithm 2 will continue to work, perturbing it in the next iteration of the ITS algorithm. 5. Computational results and discussion To evaluate the efficiency of our proposals, we analyzed the impact of the MC-TTRP model and two-phase heuristic using a modified version of a set of well-known benchmarks from the literature and a real-world problem. The proposed heuristic described in Section 4 was implemented in R4.0.2. To validate its performance, a set of experiments was conducted on the Finisterrae II supercomputer, provided by the Galicia Supercomputing Centre4(CESGA), which consists of 306 nodes powered by two deca-core Intel Haswell 2680v3 CPUs with 128 GB RAM connected through an Infiniband FDR network. The mathematical model presented in Subsection 3.2 was solved using the Gurobi 8.1.0 solver. The code was run on a hexa-core Intel i7-8700 CPU with 16 GB RAM. The following subsections report the computational study. Subsection 5.1 reports a detailed study of the exact solution of the proposed MC-TTRP model using the case study. To the best of our knowledge, there are currently no existing MC-TTRP benchmark problems. Hence, before showing the computational results of our algorithm, Subsection 5.2 describes the generation of a set of new MC-TTRP test problems based on the 21 well-known TTRP cases introduced by Chao (2002). Subsection 5.3 describes the solutions obtained both in the instances created by Chao for the TTRP as well as in our generated datasets for the MC-TTRP using our proposed heuristic. Finally, the real-world application described in Subsection 3.1 was used to validate our heuristic, as reported in Subsection 5.4, comparing the quality of the results achieved to those in the case of the exact method. 5.1. Exact solving of a real example As we have real data for this case study, an optimization scenario was built to assess our MCTTRP model. In particular, we know the distances between the different nodes (customers and central depot) as well as customer demands and vehicle capacities. Moreover, the drivers work 8 hours per day, and we estimated the average vehicle speed to be 60 km/h. Furthermore, because we do not have information about the service time of each customer or the time employees need to load the trucks, we assumed it to be negligible. In addition, in our real instance of the model, 4https://www.cesga.es/ © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 25 10 customers were selected, five of each type, provided by the company’s vehicles: three trucks and two trailers. The trucks have 13 hoppers each, which can carry up to 1.5 tons of loads, while the trailers have 15 hoppers each with a maximum capacity of 2 tons. The demand dof the customers (in kilograms), according to the four types of feed distributed by the cooperative, and the matrix of distances, C, (in kilometers) between the nodes are as follows: d= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1000 0 0 2300 4000 0 0 2041 1959 0 4000 0 0 951 2000 0 0 3500 1385 0 0 3003 0 0 516 0 0 2500 978 0 3500 0 2000 0 2513 900 0 3490 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ and C= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 21201765636019222460 21 0 4 6 60585515182055 20 4 0 4 59 56 53 13 8 12 53 17 6 4 0 57545211131652 65 60 59 57 0 3 7 66 69 71 6 63 58 56 54 3 0 4 64 66 69 3 60 55 53 52 7 4 0 61 64 66 2 19 15 13 11 66 64 61 0 3 5 61 22 18 8 13 69 66 64 3 0 7 64 24 20 12 16 71 69 66 5 7 0 66 60 55 53 52 6 3 2 61 64 66 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . To solve the associated mathematical problem, Gurobi5was used as an exact MILP solver. The optimization process required 31.2 hours to obtain the optimal solution. The value of the objective function was 207 km, containing an MVR of 74 km, whose main route, 0-3-2-1-0, is traveled by truck 1 with trailer 1 attached. Customer 2 serves as a trailer parking place for the sub-tour 2-8-79-2 with the truck. The remaining customers are served by truck 3 on a PTR, 0-6-5-4-10-0. It can be seen that only a trailer is required to supply these 10 customers. We studied the effect of trailers in the solution and compared the results with those obtained if only trucks were used. In such a case, and after a runtime of 6.23 hours, we obtained the following results: the value of the objective function was 232 km; truck 1 distributes feed to customers 1–3, traveling 46 km; truck 2 travels 53 km and serves customers 7–9; and truck 3 performs the route 0-4-5-10-6-0, whose length is of 133 km. As can be seen, when the company introduces trailers, a reduction is achieved not only in the total length traveled (which decreases by 25 km) but also in the number of drivers required (which decreases from 3 to 2). Moreover, we can deduce that when using only trucks, we are faced with an MC-VRP. However, because of the large amount of feed that the cooperative must distribute daily, we suggest the incorporation of trailers to accommodate the use of the MC-TTRP model. Furthermore, given the existence of access restrictions to some farms, this model seems appropriate. The purchase or rental of these additional vehicles can indeed be a significant initial investment, but the benefits usually compensate in the long term because, among other things, it is not necessary to hire more drivers. To study the increase in the computation time when the instance is slightly modified, we tried to solve this problem when adding a new VC. By having one more node, the number of variables involved increases considerably, and this causes the execution time to exceed two weeks. 5https://www.gurobi.com/ © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 32 L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 use of operational research techniques is necessary to obtain solutions. Therefore, we introduced and implemented a new heuristic algorithm to reduce computation time. Our proposal is a twostage approach: the first phase iteratively builds an initial solution, based on the savings method of Clarke and Wright, and then the second phase aims to refine the solution. To achieve this, an iterative tabu search was designed. We conducted a thorough computational study on different instances. First, the 21 benchmark problems of Chao (2002) were analyzed and the solutions reported by both Chao (2002) and Caramia and Guerriero (2010a) were compared to ours. In addition, we suitably adapted these data sets to consider compartments, leading to 21 challenging test problems. Our heuristic always generated a feasible solution to the test problems, and the results obtained showed that our method can effectively and efficiently solve the MC-TTRP. Furthermore, this algorithm was applied to the previously mentioned real-world case of a cooperative. A comparison between our results and those provided by the exact formulation showed a significant decrease in computational cost, achieving good-quality solutions to the problems. Considering this, we believe that the proposed heuristic is a promising solution approach for the MC-TTRP. There is still much room for further research on this variant of the VRP. It would be interesting to apply the model and solution algorithm to other feed producing companies similar to that considered in this study as well as other businesses, such as milk collection and fuel distribution companies. Another open research direction is to extend the model to include modifications, such as stochastic demands, time windows, or heterogeneous vehicle fleets, to address other similar realworld problems. In addition, future work could be to develop other metaheuristic methods, such as simulated annealing or evolutionary algorithms, that take our solution as a benchmark to compare with it. Finally, in view of the different ingredients that make up our model, including route design, assignment of customers to vehicles, and loading of compartments, exploring the performance of decomposition techniques for solving it could provide good results. The code and instances required to reproduce the results reported herein are available at https://github.com/LauraDavilaPena/ITS_MC-TTRP Acknowledgments We acknowledge the computational resources provided by CESGA. Laura Davila-Pena’s research was funded by the Ministry of Education, Culture and Sports of Spain (contract FPU17/02126). David R. Penas’ research was funded by the Xunta de Galicia (post-doctoral contract ED481B2019-010). This work was also supported by the ERDF (MINECO/AEI grant MTM2017-87197C3-3-P) and by the Xunta de Galicia (Competitive Reference Group ED431C 2017/38 and ED431C 2021/24). We would also like to thank the three anonymous referees for their constructive comments and suggestions, which helped us to improve the final version of this paper. References Accorsi, L., Vigo, D., 2020. A hybrid metaheuristic for single truck and trailer routing problems. Transportation Science 54, 5, 1351–1371. Alinaghian, M., Shokouhi, N., 2018. Multi-depot multi-compartment vehicle routing problem, solved by a hybrid adaptive large neighborhood search. Omega 76, 85–99. © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 33 Alvarez, A., Munari, P., Morabito, R., 2018. Iterated local search and simulated annealing algorithms for the inventory routing problem. International Transactions in Operational Research 25, 6, 1785–1809. Brown, G.G., Ellis, C.J., Graves, G.W., Ronen, D., 1987. Real-time, wide area dispatch of mobil tank trucks. Interfaces 17, 1, 107–120. Brown, G.G., Graves, G.W., 1981. Real-time dispatch of petroleum tank trucks. Management Science 27, 1, 19–32. Caramia, M., Guerriero, F., 2010a. A heuristic approach for the truck and trailer routing problem. Journal of the Operational Research Society 61, 7, 1168–1180. Caramia, M., Guerriero, F., 2010b. A milk collection problem with incompatibility constraints. Interfaces 40, 2, 130–143. Chao, I.M., 2002. A tabu search method for the truck and trailer routing problem. Computers & Operations Research 29, 1, 33–51. Chao, I.M., Golden, B.L., Wasil, E.A., 1998. A new algorithm for the site-dependent vehicle routing problem. In Woodruff, D. (ed.), Advances in Computational and Stochastic Optimization, Logic Programming, and Heuristics Search: Interfaces in Computer Science and Operations Research. Springer, Boston, pp. 301–312. Christofides, N., Mingozzi, A., Toth, P., 1979. The vehicle routing problem. In Christofides, N., Mingozzi, A., Toth, P. and Sandi, C. (eds), Combinatorial Optimization. Wiley, Chichester, pp. 315–338. Clarke, G., Wright, J.W., 1964. Scheduling of vehicles from a central depot to a number of delivery points. Operations Research 12, 4, 568–581. Coelho, L.C., Laporte, G., 2015. Classification, models and exact algorithms for multi-compartment delivery problems. European Journal of Operational Research 242, 3, 854–864. Cordeau, J.F., Maischberger, M., 2012. A parallel iterated tabu search heuristic for vehicle routing problems. Computers & Operations Research 39, 9, 2033–2050. Davila-Pena, L., 2019. Modelos y algoritmos en una clase de problemas de rutas de vehículos. Master’s thesis, University of Santiago de Compostela. Derigs, U., Gottlieb, J., Kalkoff, J., Piesche, M., Rothlauf, F., Vogel, U., 2011. Vehicle routing with compartments: applications, modelling and heuristics. OR Spectrum 33, 4, 885–914. Derigs, U., Pullmann, M., Vogel, U., 2013. Truck and trailer routing problems, heuristics and computational experience. Computers & Operations Research 40, 2, 536–546. El Fallahi, A., Prins, C., Calvo, R.W., 2008. A memetic algorithm and a tabu search for the multi-compartment vehicle routing problem. Computers & Operations Research 35, 5, 1725–1741. Gendreau, M., Hertz, A., Laporte, G., 1992. New insertion and postoptimization procedures for the traveling salesman problem. Operations Research 40, 6, 1086–1094. Gerdessen, J.C., 1996. Vehicle routing problem with trailers. European Journal of Operational Research 93, 1, 135–147. Glover, F., Laguna, M., 1998. Tabu search. In Pardalos, P. (ed.), Handbook of Combinatorial Optimization.Springer, Boston, pp. 2093–2229. Guitián de Frutos, R.M., Casas-Méndez, B., 2019. Routing problems in agricultural cooperatives: a model for optimization of transport vehicle logistics. IMA Journal of Management Mathematics 30, 4, 387–412. Henke, T., Speranza, M.G., Wäscher, G., 2015. The multi-compartment vehicle routing problem with flexible compartment sizes. European Journal of Operational Research 246, 3, 730–743. Henke, T., Speranza, M.G., Wäscher, G., 2019. A branch-and-cut algorithm for the multicompartment vehicle routing problem with flexible compartment sizes. Annals of Operations Research 275, 2, 321–338. Lespay, H., Suchan, K., 2021. A case study of consistent vehicle routing problem with time windows. International Transactions in Operational Research 28, 3, 1135–1163. Lin, S.W., Vincent, F.Y., Lu, C.C., 2011. A simulated annealing heuristic for the truck and trailer routing problem with time windows. Expert Systems with Applications 38, 12, 15244–15252. Lin, S.W., Yu, V.F.Y., Chou, S.Y., 2009. Solving the truck and trailer routing problem based on a simulated annealing heuristic. Computers & Operations Research 36, 5, 1683–1692. Mendoza, J.E., Castanier, B., Guéret, C., Medaglia, A.L., Velasco, N., 2010. A memetic algorithm for the multicompartment vehicle routing problem with stochastic demands. Computers & Operations Research 37, 11, 1886– 1898. Mendoza, J.E., Castanier, B., Guéret, C., Medaglia, A.L., Velasco, N., 2011. Constructive heuristics for the multicompartment vehicle routing problem with stochastic demands. Transportation Science 45, 3, 346–363. © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies 34 L. Davila-Pena et al. / Intl. Trans. in Op. Res. 0 (2021) 1–34 Muyldermans, L., Pang, G., 2010. On the benefits of co-collection: Experiments with a multi-compartment vehicle routing algorithm. European Journal of Operational Research 206, 1, 93–103. Ostermeier, M., Henke, T., Hübner, A., Wäscher, G., 2021. Multi-compartment vehicle routing problems: State-of-theart, modeling framework and future directions. European Journal of Operational Research 292, 3, 799–817. Ostermeier, M., Hübner, A., 2018. Vehicle selection for a multi-compartment vehicle routing problem. European Journal of Operational Research 269, 2, 682–694. Parragh, S.N., Cordeau, J.F., 2017. Branch-and-price and adaptive large neighborhood search for the truck and trailer routing problem with time windows. Computers & Operations Research 83, 28–44. Renaud, J., Boctor, F.F., Laporte, G., 1996. A fast composite heuristic for the symmetric traveling salesman problem. INFORMS Journal on Computing 8, 2, 134–143. Rothenbächer, A.K., Drexl, M., Irnich, S., 2018. Branch-and-price-and-cut for the truck-and-trailer routing problem with time windows. Transportation Science 52, 5, 1174–1190. Scheuerer, S., 2004. Neue Tabusuche-Heuristiken für die logistische Tourenplanung bei restringierendem Anhängereinsatz, mehreren Depots und Planungsperioden. PhD thesis, University of Regensburg. Scheuerer, S., 2006. A tabu search heuristic for the truck and trailer routing problem. Computers & Operations Research 33, 4, 894–909. Semet, F., 1995. A two-phase algorithm for the partial accessibility constrained vehicle routing problem. Annals of Operations Research 61, 1, 45–65. Semet, F., Taillard, E., 1993. Solving real-life vehicle routing problems efficiently using tabu search. Annals of Operations Research 41, 4, 469–488. Silvestrin, P.V., Ritt, M., 2017. An iterated tabu search for the multi-compartment vehicle routing problem. Computers & Operations Research 81, 192–202. © 2021 The Authors. International Transactions in Operational Research published by John Wiley & Sons Ltd on behalf of International Federation of Operational Research Societies