A mathematical formulation and heuristic approach for the heterogeneous fixed fleet vehicle routing problem with simultaneous pickup and delivery

This study considers a variant of the vehicle routing problem (VRP) called the heterogeneous VRP with simultaneous pickup and delivery (HVRPSPD). The HVRPSPD may broadly be defined as identifying the minimum cost routes and vehicle types. To solve the HVRPSPD, first, we propose a polynomial-size mixed integer programming formulation. Because the HVRPSPD is an NP-hard problem, it is difficult to determine the optimal solution in a reasonable time for moderate and large-size problem instances. Hence, we develop a hybrid metaheuristic approach based on the simulated annealing and local search algorithms called SA-LS. We conduct a computational study in three stages. First, the performance of the mathematical model and SA-LS are investigated on small and medium-size HVRPSPD instances. Second, we compare SA-LS with the constructive heuristics, nearest neighborhood and Clarke-Wright savings algorithms, adapted for the HVRPSPD. Finally, the performance of SA-LS is evaluated on the instances of the heterogeneous VRP (HVRP), which is a special case of the HVRPSPD. Computational results demonstrate that the mathematical model can solve small-size instances optimally up to 35 nodes; SA-LS provides good quality solutions for medium and large-size problems. Moreover, SA-LS is superior to simple constructive heuristics and can be a preferable solution method to solve HVRP and VRPSPD instances successfully.


1.
Introduction. At the present time, an increasing population and scarce resources emphasize the importance of logistics and distribution systems for the economy of countries, sectors, and companies, and this can be seen in dramatic reports. Toth and Vigo reported that the use of computerized procedures for distribution process planning produces substantial savings (generally from 5% to 20%) in global transportation costs [63]. A more recent report from the World Bank states that logistics activities in France, which has a higher LPI score in the world ranking, represents 10% of its national gross domestic product, e200 billion turnover, and 1.8 million jobs [3]. Coşar and Demir emphasized that efficient logistics enables countries to participate in global supply chains and exploit their comparative advantages [15]. Their findings have important developmental implications, for example, the more an industry is transportation-sensitive, the greater the reduction in transportation costs.
The vehicle routing problem (VRP) introduced by Dantzig and Ramser is a core application field of operational-level logistics [16]. The VRP, which is an NP-hard problem, can be defined as identifying the minimum cost routes to satisfy customer needs in addition to some side constraints.
The classical VRP consists of identical (homogeneous) vehicles. However, in real-life applications, the vehicles in a fleet may have different features in terms of, for example, purchasing cost, unit transportation cost, and capacity. Moreover, customer (freight) needs may require different vehicle properties. Thus, to reduce the logistics expenditure, both the decisions of fleet selection and route identification have to be examined by firms. In particular, fleet selection is a strategic decision in fleet investment. The classical VRP also assumes that customers have either a delivery or pickup demand. In particular applications, each customer has both delivery and pickup demands simultaneously. In this situation, major economic benefits can be obtained when both activities are performed by the same vehicle.
The heterogeneous VRP (HVRP) may be encountered in various real-life applications. For example, it is common practice to use vehicles with different capacities in residential and commercial waste collection activities. Other examples are shipping services for FedEx Ground and single copy newspaper delivery operations. The VRPSPD is mainly related to the reverse logistics. The management of the reverse flow of products, and raw and work-in-process materials are more challenging topics for companies. Environmental issues are another example where reverse logistics forces firms to use their distribution network more effectively, for example, in waste management, recycling, reprocessing, remanufacturing, and the evaluation of used products. In particular, The VRPSPD can be encountered in many real-life practices such as; the grocery store chains where delivery of fresh food to markets and pickup of outdated items from markets; beverage industry where delivery of soft drink to markets and pickup of empty bottles from markets; foundry industry where delivery of purified reusable sand and pickup of used sand.
In this study, we consider the heterogeneous VRPSPD (HVRPSPD), which includes more realistic features by considering the heterogeneity of the fleet, and simultaneous pickup and delivery. We propose a polynomial-size decision model for the HVRPSPD. Because of the NP-hardness of the HVRPSPD, the proposed mixed integer programming (MIP) formulation can obtain optimal solutions for small-size problems. Thus, we propose a hybrid heuristic algorithm based on simulated annealing (SA) and local search (LS) called SA-LS to solve medium and large-size instances of the problem. It should be noted that SA is a popular neighborhood search metaheuristic algorithm that can escape from the local optimum. In out implementation, while SA is used to reach different points of the solution space, LS is used to intensively search a particular subspace. To evaluate the effectiveness of the developed hybrid heuristic algorithm, we also adapt well-known and basic constructive heuristics, nearest neighborhood (NN) and Clarke-Wright savings (CWS) algorithms, for the HVRPSPD.
We conduct a three-stage experimental study to investigate the performance of the MIP formulation and SA-LS on HVRSPD instances derived from HVRP instances in the literature. The first stage evaluates the performance of the MIP formulation and SA-LS on small and medium-size HVRSPD instances, and the second stage compares SA-LS with NN and CWS, on large-size HVRSPD instances. Finally, the last stage evaluates the performance of SA-LS on HVRP instances, which is a special case of the HVRPSPD.
The main contributions of this study are as follows: 1) A rarely studied special type of VRP is considered. This problem is directly related to reverse logistics, which is one of the key concepts in green logistics, recycling, and waste management.
2) A new hybrid algorithm based on SA and LS is proposed as a solution approach for the HVRPSPD. 3) Two basic constructive heuristics known in VRP literature are adapted to quickly generate initial solutions for heuristics. 4) A new demand separation method is developed to derive HVRPSPD instances from HVRPs. Additionally, new test instances for the HVRPSPD are introduced. Especially, the relationship of HVRPSPD with several routing problems in logistics is described (see Section 3.3). Furthermore the detailed description about how to adapt the proposed formulation for the variants is provided in detail.
The remainder of the paper is organized as follows: In Section 2, we present a literature review of related works. In Section 3, we define the HVRPSPD, and then present an MIP formulation for the problem and describe the adaptation of the formulation to special cases of the HVRPSPD. In Section 4, we explain the details of the proposed hybrid metaheuristic for the HVRPSPD. We report computational results in Section 5, and discuss the conclusion and future research in Section 6.
2. Literature review. Because there are only a few studies of the HVRPSPD in the literature, first we review the literature for two sub-problems: HVRP and VRPSPD.
Two main types of the HVRP have been studied in the literature for cases in which the fleet has an unlimited or limited number of vehicles. Golden et al. was the first to propose an HVRP in which the fleet has an unlimited number of vehicles of each type [23]. This problem is referred to as the fleet size and mix VRP [23], fleet size and composition VRP [21], and vehicle fleet composition problem [52]. Taillard was the first to study an HVRP that had a limited number of vehicles of each type in the fleet [58]. This type is more realistic and referred to by several terms in the literature, such as the VRP with a heterogeneous fleet of vehicles [58] and heterogeneous fixed fleet VRP [61]. Although these two problems are very similar, their applications are different: the HVRP for an unlimited number of vehicles is more appropriate for strategic decisions when a company wants to buy a vehicle fleet and needs to define its size and composition, whereas the HVRP for a limited number of vehicles better represents the operational decisions of defining the vehicles that should be used, among those available, to serve customers [11]. Related with this paper, a more recent work is done by Juan et al. to analyze different fleet configurations and proposed several approximation methods [27]. Additionally, it is also possible to classify HVRPs depending on whether fixed and transportation costs are considered. Baldacci et al. presented a classification scheme and provided exact solution algorithms [7] and Baldacci et al. provided latest advances and challenges for the HVRP [6]. We refer interested readers to the papers of [24], [44], [26] and [34] for an extensive review and recent past of this problem and its variants.
The VRPSPD was first introduced by Min [40], who studied a book circulation problem for a public library. The author provided an MIP formulation for the problem and proposed a heuristic algorithm based on the cluster first-route second approach to obtain a solution. In the following years, studies were conducted on the VRPSPD. The mathematical models developed for the VRPSPD can be found in [18], [41] and [42]. Dell'Amico et al. proposed the first exact algorithm for the VRPSPD [17]. Exact dynamic programming and state space relaxation approaches were used for the sub-pricing problem. The algorithm can solve instances up to 40 customers optimally. Ropke and Pisinger developed a large neighborhood search heuristic for the VRPSPD and some VRP variants [50]. Bianchessi and Righini proposed a heuristic approach based on SA for the VRPSPD [9]. Ai and Kachitvichyanukul developed a heuristic algorithm based on particle swarm optimization (PSO) [1]. Gajpal and Abad presented saving-based heuristics for the VRPSPD, where the heuristics were based on merging two existing routes [19]. The feasibility of the new route was checked using a cumulative net-pickup approach. Subramanian et al. presented a parallel heuristic approach based on iterative LS and variable neighborhood search algorithms [55]. Furthermore, Subramanian et al. proposed a branch-cut-and-price algorithm for the VRPSPD [57]. Goksal et al. proposed a hybrid metaheuristic algorithm based on PSO in which the variable neighborhood descent (VND) algorithm was implemented for LS [22]. Tasan and Gen presented a genetic algorithm (GA), which implemented a permutation-based representation [62]. The initial population was constructed randomly, and genetic operators, crossover, and mutation were implemented on the members of the population. Polat et al. considered the VRPSPD with time limitations [46]. In this study, a mathematical formulation based on the model by Montane and Galvao [41] was provided. To solve the problem, the authors presented a perturbation-based neighborhood search algorithm. Avci and Topaloglu developed a hybrid adaptive LS solution approach based on SA and VND for the VRPSPD and one of its variants in which customers were visited in a mixed order [4].
To the best of the authors' knowledge, there are few studies on the HVRPSPD. Kececi et al. give the mathematical formulations and propose a simple constructive heuristic method [31]. Kececi et al. also developed a matheuristic (mathematical model based heuristic) approach [32]. Avci and Topaloglu proposed a hybrid metaheuristic algorithm for the HVRPSPD in which an LS procedure with a nonmonotone threshold adjusting strategy was integrated with the tabu search [5]. In their algorithm, the authors applied an encoding structure, which was proposed in our previous study [30]. They also provided an arc-based MIP formulation, including variables with three indices, which was adapted from Montané and Galvao's formulation [41]. The performance of the proposed hybrid metaheuristic algorithm was tested by comparing its results with those of a GA on two randomly generated test instance sets, where each set had only 14 problems. Ç etin and Gencer considered the heterogeneous fleet VRPSPD with time windows [13]. The proposed MIP formulation for the problem was based on Dethloff's model [18]. In a recent paper a variant of HVRPSPD which includes the location decision is studied in Wang [64].
Unlike the aforementioned studies, we propose a hybrid algorithm based on the SA and LS and, as observed from the literature review, this type of hybridization of SA and LS is the first application to the HVRPSPD. The rationale behind this hybridization originates from the fact that combining features of different heuristics in a complementary manner can result in more robust and effective optimization tools. Thus, in this paper, we propose an effective solution approach for the HVRP-SPD by hybridizing SA with LS, called SA-LS. SA searches for good solutions in the solution space, whereas LS improves the best solution obtained by SA for each iteration. Recently, the hybridization of SA and LS has been used to solve various combinatorial optimization problems in the literature (e.g., scheduling by Liao et al. [35] and Calleja et al. [12], and VRP by Allahyari et al. [2], Zhang et al. [67] and Masmoudi et al. [38]).
3. Mixed integer programming formulation for the HVRPSPD. In this section, after defining the HVRPSPD, we present a MIP formulation for the problem and then explain how the formulation can be revised for special cases. in the fleet. For each k ∈ B, there are T k available vehicles with the capacity of Q k and fixed cost of f k . Each arc (i, j) ∈ A is associated with a nonnegative cost, c ij = θ k l ij where l ij is the distance from node i to j and l ij = l ji . The triangular inequality holds (i.e., l ij + l jk ≥ l ik ) ∀i, j, k ∈ N . θ k is the variable transportation cost per unit distance of vehicle type k ∈ B. Each customer i ∈ N has both a pickup (p i ) and delivery (d i ) demand, and 0 ≤ d i , p i ≤ Q k , ∀k ∈ B, d 0 = p 0 = 0. The problem is to determine both the vehicle tours with minimum total cost and vehicle type on each tour under the following restrictions: • exactly one vehicle type is used on each tour; • each customer is visited by one vehicle type; • each vehicle starts and ends its tour at the depot; and • the total vehicle load is not allowed to exceed the capacity.

Proposed MIP formulation.
To solve the HVRPSPD, we propose an MIP formulation with a polynomial number of constraints and decision variables. The MIP formulation is obtained from the adaptation of Kara and Derya's formulation [28] for the VRP. Additionally, the MIP formulation of the HVRPSPD generalizes the formulations of the VRPSPD proposed by Karaoglan [29] and Montané and Galvao [41] in terms of a heterogeneous fleet. Because the proposed formulation includes auxiliary decision variables defined on each arc of the distribution network, it is also called an arc-based formulation (ABF). Waters was the first researcher to propose the ABF for the VRP, in which the subtour elimination constraints inspired by mass balance equations were used to restrict both capacity and tour length [66]. Kara and Derya strengthened the model for the same problem using tighter bounding constraints that resulted in a better linear relaxation value [28]. We adapt the subtour elimination and capacity constraints proposed by Kara and Derya [28] in our formulation considering both pickup and delivery activities. The decision variables of the formulation are given as follows: 0-1 Decision variables: x ijk equal to 1 if arc (i, j) is on the tour of vehicle type k and 0 otherwise, ∀i, j ∈ N, ∀k ∈ B.
Continuous decision variables: y k : number of type k vehicles in the fleet (k ∈ B), m: number of vehicle tours.
Auxiliary decision variables: z ij : remaining delivery load if the vehicle passes from i to j immediately after visiting node i, 0 otherwise; t ij : load picked up if the vehicle passes from i to j immediately after visiting node i, 0 otherwise. The proposed mathematical model is given as follows: The objective function is the minimization of the total transportation and vehicle utilization costs and is denoted by subject to the degree constraints for the depot given by The following constraints force the solution to have exactly one vehicle type moving in and out at any node; hence, it provides vehicle type continuity: On any arc, the constraint does not allow the violation of the vehicle capacity. The following constraints relate to the delivery and pickup load, respectively, and both prevent the subtours: The following constraints bound the delivery and pickup loads, respectively: The following constraints initially provide value to the pickup and delivery loads, respectively, because a vehicle starts and ends its tour with an empty load: The following constraint provides at most m vehicles used in the fleet: The following constraint is the vehicle type availability restriction: The following constraint ensures a feasible solution, which is that the leaving arcs from the depot associated with vehicle type k must be exactly the same as the number of vehicles of type k: The integrality and non-negativity constraints are In the ABF, the number of binary decision variables is O(n 2 b), number of positive variables is O(2n 2 ), number of positive and integer variables is O(b), and number of restrictions is O(3n 2 ).

Special cases of HVRPSPD.
It should be noted that the proposed MIP formulation of the HVRPSPD is easily adaptable to special variants, which arise in logistics problems. Figure 1 shows the connections between the HVRPSPD and its variants. An HVRP is obtained by setting either p i or d i for all i ∈ N to zero in According to the fleet size, variants of the HVRPSPD can be grouped into two classes. The proposed formulation of the HVRPSPD can easily be adapted for the limited fleet HVRPSPD with fixed costs (LHVRPSPDF) by setting θ k1 = θ k2 , ∀k 1 , k 2 ∈ B, k 1 = k 2 or θ k = 1, ∀k ∈ B. For this case, the routing costs do not depend on the vehicle type and there will only be vehicle dependent fixed costs. For the opposite case, where fixed costs do not depend on the vehicle type (or do not exist) but routing costs do, the problem becomes the limited fleet HVRPSPD with dependent routing costs (LHVRPSPDD). The proposed formulation of the HVRPSPD can easily be adapted for the LHVRPSPDD by setting In the SHVRPSPD, there is a limited number of heterogeneous vehicles available for the service with no fixed costs, and routing costs do not depend on the vehicle type. However, each customer may include restrictions on the vehicle types that may visit it. To adapt the HVRPSPD formulations to this special case, let B j , ∀j ∈ N \ {0} be the set of vehicle types that can visit the customer. Set Q k = M, ∀k ∈ B \ B j , j ∈ N \ {0} for the vehicle types that are not allowed to visit the customer; set θ k1 = θ k2 , ∀k 1 , k 2 ∈ B j , k 1 = k 2 or θ k = 1, ∀k ∈ B j , j ∈ N \ {0} for the vehicle types that are allowed to visit the customer, where M is a very large number; ignore the second term of the objective function given in (1). If there is an unlimited number of vehicles in the fleet, the problem becomes the HVRPSPD with unlimited fleet (UHVRPSPD). A HVRPSPD instance can be converted to UHVRPSPD by setting T k = +∞ or by omitting (14) in the formulation. If there is an unlimited fleet with routing costs independent of the vehicle types and with vehicle dependent fixed costs, the HVRPSPD formulation can be adapted to the unlimited fleet HVRPSPD with fixed costs (UHVRPSPDF) by omitting (14) and setting θ k = 1, ∀k ∈ B.
For the opposite case, where fixed costs do not depend on the vehicle type (or do not exist) but routing costs do with an unlimited fleet of vehicles, the HVRPSPD formulation can be adapted to the UHVRPSPDD by omitting (14) and setting if there are no fixed costs. With these reductions, all of these aforementioned special cases can be solved using the proposed MIP formulation for the HVRPSPD.
4. Hybrid simulated annealing algorithm for the HVRPSPD. Because of the NP-hardness of HVRPSPD, the MIP formulation cannot directly be applied to solve optimally large and even medium-size problems. Thus, heuristics are required to rapidly obtain reasonable solutions. In this study, to solve medium and large-size HVRPSPD instances, we propose a hybrid heuristic method, which combines SA and LS. After SA was proposed by Metropolis et al. [39], it was improved by Kirkpatrick et al. [33] to solve optimization problems. The name of the algorithm was derived from the analogy between solving optimization problems and simulating the annealing process of solids.
Studies in recent decades have considered hybrid algorithms, in which metaheuristics are used together with LS, or other metaheuristic or exact solution algorithms to take advantage of more than one approach [60]. LS is also a well-known simple heuristic that searches neighbors of the current solution until no improved solution is obtained.
In our hybrid algorithm, while SA is used to reach different points of the solution space, LS is used to intensively search a particular subspace. The proposed hybrid algorithm is called SA-LS. This successive use of SA and LS continues until stopping criteria are met.
In this section, the basic structure of SA-LS is described by considering the representation of the solution, generation of the initial solution, neighborhood structures, and LS.

Solution representation.
Representation is an important issue and it directly affects the quality of metaheuristics. Generally, data structures or representations used in different problems are not identical. In SA-LS, we use a |N | × |N | matrix which is known as a permutation encoding where every sequence is a string of numbers to represent the solutions of the HVRPSPD. Each row in the matrix provides the route and vehicle-type information of the solution. The first element in the row indicates the type of vehicle that is assigned to the route and the remaining elements indicate the sequence of customers that are visited by the vehicle. Vehicle type 3 visits customers 7, 1, and 8 for the first route, vehicle type 1 visits customers 2 and 9 for the second route, and vehicle type 3 visits customers 4, 5, 3, 10, and 6 for the third route. In this example, there are three vehicle types; however, only two are assigned to routes.

4.2.
Generating the initial solution. Initial solutions are used to start search mechanisms. There are two common approaches to generate an initial solution: one is to generate it randomly and the other is to use specific heuristics developed for the problem. Because heuristics obtain better quality solutions than a randomly generated approach, using a heuristic algorithm is mostly preferred to avoid the computational burden. In this study, we propose a heuristic approach to initialize a solution for the HVRPSPD. This approach, based on the giant tour and its partition, is a type of "route first-cluster second" approach, which is very common in VRP literature. After this approach was proposed by Beasley [8] for the VRP, it was successfully applied to variants of VRP, for example, by Golden et al. [23], Prins [47], Prins [48], and Goksal et al. [22].
The proposed initialization approach to generate the starting solution for the HVRPSPD consists of three main steps: Generate a giant tour, partition and check feasibility. These steps are described as follows: For the first step, a giant tour is generated by solving the traveling salesman problem (TSP) in which the pickup and delivery demands of customers and the capacity restrictions of vehicles are omitted. We use the Concorde TSP solver to optimally solve the TSP, that is, to obtain an ordered set of customers E = {E 1 , E 2 , ..., E |N |−1 } over the given network.
For the second step, the giant tour along E is partitioned into routes. An acyclic auxiliary graph H = (Ê, F ) is built to partition the giant tour. In this graph, is the set of nodes, which contains a dummy node 0 and the nodes 1 to |N | − 1 for customers E 1 to E |N |−1 . Arc (i, j, k) in the graph represents the vehicle of type k leaving the depot, visiting customers E i+1 to E j , and returning to the depot. Weightĉ ijk of arc (i, j, k) is calculated as in Equation (20) where M is a sufficiently large number.ĉ An optimal partition of E is obtained by a minimum cost path from nodes 0 to |N | − 1 in H, where at most T k arc (i, j, k) must be selected because there is an availability restriction for the vehicle of type k. This problem is called the shortest path problem with resource constraints [25]. Although it is an NP-hard problem, it can be solved quickly in practice using dynamic programming methods [48]. Thus, we use the following MIP formulation (SPPRC) to optimally partition E into routes. In this formulation, the binary decision variablex ijk = 1 if and only if arc (i, j, k) is in the solution, otherwisex ijk = 0. The SPPRC formulation is given as follows: The objective function minimizes the total arc weights: subject to, that satisfy exactly one arc that leaves the origin and enters the destination, respectively. The following constraint yields exactly one arc that enters and leaves the intermediate nodes: The following constraint guarantees that at most T k arc (i, j, k) are selected in the optimal solution: The following constraint defines the integrality conditions of the decision variables: x For the third step, after partitioning the giant tour into routes, every route is checked to determine whether the maximum load along the route exceeds the vehicle capacity. If a violation is observed in any route in terms of the vehicle capacity, that route is repaired by considering it as a VRPSPD with one vehicle. As was indicated in Section 3.3, the VRPSPD with one vehicle is a special case of the HVRPSPD. Thus, this problem is easily solved by the MIP formulation of the HVRPSPD considering the following settings: N = N , B = {k * }, m = 1, where N is a node set that consists of a depot and nodes from the tour to be repaired and k * is the type of vehicle on the route. The inter-route moving strategies are defined as follows and the illustrative examples for inter-route moving strategies are shown in Figure 3: Shift(1,0): Node i from tour r 1 is transferred to tour r 2 . Shift(2,0): Two adjacent nodes i and j from tour r 1 are transferred to tour r 2 . Swap(1,1): Node i from tour r 1 and node j from tour r 2 are exchanged. Swap(2,1): Two adjacent nodes, i and j, from tour r 1 are replaced by node k from tour r 2 . Swap(2,2): Two adjacent nodes, i and j, from tour r 1 are exchanged by another two adjacent nodes, k and l, from tour r 2 . Cross: Arcs between nodes i and j from tour r 1 and between k and l from tour r 2 are removed. Then two arcs connecting nodes i and l and nodes k and j are inserted. K-Shift: Subset of consecutive nodes from tour r 1 is transferred to the end of tour r 2 .
The intra-route strategies are explained as follows and Figure 4 shows the examples of intra-route moving strategies: Or-opt: One, two, or three adjacent nodes from a tour are removed and inserted into another position of the same tour. 2-opt: Two non-adjacent arcs are deleted and another two are added in such a manner that a new tour is generated. Exchange: Positions of two nodes i and j in a tour are exchanged. Reverse: Tour direction is reversed if there is a chance to reduce the maximum load along the tour.
SA-LS implements the best improvement strategy when searching within the solution space by a moving mechanism to select a neighbor of the current solution. Furthermore, it implements only feasible moves that do not exceed capacity restrictions. We conducted several tests on inter-route moving mechanisms to determine   [59]. In this study, the SA is hybridized with simple LS, called SA-LS, to obtain a good solution through the interactive use of two algorithms together within a reasonable solution time.
The pseudo code of the LS is shown in Figure 5 (lines 00-28). The LS procedure simply takes a feasible solution (S LS ) as an input and attempts to improve it. If there is no improvement, it returns the input solution (S LS ) as the output. In the beginning of LS one moving strategy was randomly chosen from the neighborhood structures of Type-1 to obtain a new solution at a certain iteration of the SA-LS algorithm (lines 01-03). If the objective function value of the new solution is better than the current solution (line 04), then the new solution is set as the current solution (line 06). Next, a random number between zero and one is generated uniformly (line 07). If the random number is greater than p, then the following two steps (lines [11][12][13][14][15][16] are repeated while the number of iterations without any improvement in the objective function is less than a certain number (this number is set as 5 for the preliminary trials): 1) A new solution is obtained by a moving strategy chosen randomly among the neighborhood structures of Type-2 and 2) if the objective value of the new solution is less than the current solution, then the new solution is set to the current solution. In this way, diversification is performed by jumping to different points in the solution space with the neighborhood structures of Type-1. Additionally, intensification is performed by searching a particular subspace with the neighborhood structures of Type-2. Finally, the neighborhood structures of Type-3 are all applied in the given order (lines [18][19][20][21][22][23][24][25]. If the objective value of the new solution is worse than the current solution, then the LS procedure returns the input solution as the output without any update (line 28).

4.5.
Steps of SA-LS. The overall steps of the methodology proposed for the HVRPSPD are shown in Figure 6 (lines 00-16) and summarized as follows: SA-LS starts with an initial solution S (line 01). The current temperature T is set to T 0 (line 02). The main steps of SA-LS between lines 04 and 14 are repeated until a stopping condition is met. At each temperature of SA-LS, the current feasible solution is submitted to the LS procedure as an input to improve the solution (line 05). Within the LS procedure, a new solution S new is generated in the neighborhood of the current solution S using a moving mechanism, which is randomly selected within a group of strategies (called Type-1 moving strategies). If the new solution S new is better than the current solution S, then the search process continues with simple LS using a moving mechanism, which is randomly chosen each time within a group of strategies (called Type-2 moving strategies), and updates the new solution S new . The LS continues until no improvement is obtained for five successive iterations. At the end of the LS, four moving strategies (called Type-3 moving strategies) are implemented in a predetermined order for S new , and S new is updated.
Eventually, the LS procedure returns a new solution if the objective function is improved, and if there is no improvement, the current solution is returned to SA-LS (S new := LS(S)). If the new solution is better than the current solution, that is,  Generate random number U[0,1]; 08.
If U  p then; 09.
Repeat while the # of iterations without improvement is less than 5. 11.

15.
Stop and report S; 16. END Figure 6. Pseudo code of SA-LS instances. The third phase attempts to determine the performance of SA-LS on the HVRP instances, which is a special case of the HVRPSPD. The fourth phase compares the performance of SA-LS on the VRPSPD instances, which is another special case of the HVRPSPD. The fifth phase compares the SA-LS with the results obtained in the study of Avci and Topaloglu [5]. And the last phase shows the convergence behavior of the SA-LS due to the change in different parameters. The section for computational experiments ends with some managerial insights for the companies in logistics sector. The OPL language and state-of-the-art LP/MIP solver engine CPLEX R (version 12.3) were used to model and solve the MIP models and their relaxations. The heuristic algorithms, which used the CPLEX Concert Technology, were coded in C++. Computational experiments were run on a PC with a 2.67 GHz Intel Core i5 750 CPU and 2 GB of RAM. 5.1. Test problems. Two HVRP problem sets were used to generate the HVRP-SPD test problems. The first HVRP problem set was derived by Taillard [58] from Golden et al.'s [23] problems. This set included four instances with 50 customers, two instances with 75 customers, and two instances with 100 customers. The second HVRP test set was derived by Liu and Shen [36] from Solomon's [54] problems. This set consisted of 80 instances with 100 customers. The HVRPSPD test problems could easily be generated from the HVRP instances using demand separation approaches.
We used two demand separation techniques to generate the pickup and delivery demands of customers in each HVRPSPD test problem. The first was proposed by Salhi and Nagy [51]. In this technique, r i = min{x i /y i , y i /x i } is calculated for each customer and the original demands are split using this ratio. For instance, let q i be the original demand of the customer i. Then the pickup and delivery demands are calculated by p i = (1 − r i )q i and d i = r i q i , respectively. These problems are called Type X. In the same manner, another problem type (referred to as Type Y ) is obtained by shifting each demand of the customer to the next customer's demand. The second separation technique is proposed in this study. This technique splits the original demands using the so-called "golden ratio". In mathematics, if we use the golden ratio to divide a line segment |AB|, this segment should be divided by point C such that |CB|/|AC| = |AB|/|CB| = ϕ, where ϕ is a constant ratio, which is an irrational number, and ϕ = (1 + √ 5)/2. Regarding these definitions, if i is odd and |AB| = q i , |AC| = d i and |CB| = p i , then the original demand is divided into pickup and delivery demands as p i = q i /ϕ , d i = q i − p i . If i is even and |AB| = q i , |AC| = p i and |CB| = d i , then the original demand is divided into the delivery and pickup demands as d i = q i /ϕ , p i = q i − d i . These problems are called Type W and we obtain Type Z problems as explained above for Type Y problems. As a result, 104 (4 ) large-size test instances were generated using 26 original HVRP problems and four separation procedures (Type X, Y, W, and Z ).
In addition to these problems, 520 (5 ) small and medium-size instances were obtained using the first 20, 25, 30, 35, and 40 customers in each of the 104 problems. The main characteristics of our test problems are presented in Table 1.

5.2.
Evaluating the performance of SA-LS with the MIP formulation. In the first phase, we investigated the performance of the proposed hybrid heuristic approach, SA-LS, on small and medium-size problems.
Based on our preliminary experiments, we set the following parameter values in SA-LS: the initial temperature is 380, in which an inferior solution (inferior by 40% relative to the current solution) is accepted with a probability of 0.90; the cooling rate is set to 0.95; the final temperature is 0.15, such that a solution that is inferior by 1% relative to the current solution is accepted with a probability of 0.001; and the LS procedure within SA-LS is used with the probability of 0.6. SA-LS stops whenever the temperature reaches the final temperature or the best solution is not improved for successive iterations. Each instance was run five times by the proposed heuristic with a different random number of seeds and the best of five runs for each instance was considered as the solution of the heuristic; the computation time of MIP formulation was limited to 2 CPU hours. Table 2 summarizes the computational results of SA-LS on small and mediumsize problems. The first two columns of the table show the number of customers in a HVRPSPD instance (n) and the demand separation strategy. The subsequent four columns show the average value of the percentage gap (Gap% ), average improvement of the initial solution (Imp% ), number of problems solved optimally (#OpSol ) with SA-LS, and average computation time (CPU ). The next three columns show the lower bound percentage gap obtained by the MIP formulation (Gap% ), number of problems solved optimally (#OpSol ) with MIP and average computation time of the MIP formulation (CPU ). Finally, the last three columns give the Gap%, #Op-Sol and CPU values for the matheuristic algorithm proposed by Kececi et al. [32]. It should be noted that the average and maximum percentage gaps of the heuristic solutions are calculated by considering optimal solutions, or if unknown, the best lower bounds, which are obtained by solving the MIP formulation with CPLEX for a maximum of 2 CPU hours.
From Table 2, we can observe that SA-LS obtained good results (8.19%, on average, for the optimal solution or the best lower bound) in a very short computation time. The average computation time of SA-LS was 5.87 s, whereas the MIP formulation required 2,493 s and MatH-LS 23.79 s, on average, to solve HVRPSPD instances. The average percentage gap of SA-LS was 8.19% and 12.78% for MatH-LS; however, this value reduced to 7.24% for the MIP formulation. Moreover, the performance of the MIP formulation quickly degenerated for instances with 35 and 40 customers. SA-LS optimally solved 45 out of 174 instances and MatH-LS solved 78 instances, which were optimally solved by the MIP formulation. Additionally, SA-LS improved upper bounds obtained by the MIP formulation for 216 instances. These results show that SA-LS is superior to the MIP formulation and MatH-LS as well, especially in terms of the solution time. MatH-LS can optimally solve more instances than SA-LS, however SA-LS can obtain better objective function value on the instances which could not be solved optimally, since it gives lower gap. Moreover, SA-LS can improve the initial solution by 5.20%, on average, and this improvement is between 2.80% and 11.70%.

5.3.
Investigating the performance of SA-LS with simple heuristics. In the second phase, we compare the performance of SA-LS with NN and CWS on large-size instances in terms of solution quality and solution time. In this study, we adapted the NN and CWS algorithms to solve HVRPSPD. The details of both heuristics are given in the Appendix. Table 3 summarizes the computational results of SA-LS in comparison with the NN and CWS algorithms. The first two columns of the table are the same as the previous tables. The next three columns show the average percentage deviation (Dev% ) of the best solution determined by SA-LS from the solutions determined by the NN and CWS algorithms and the best lower bounds, which are obtained by solving the MIP formulation with CPLEX for a maximum of 2 CPU hours, respectively. In comparison with the MIP model, the best lower bounds are used since MIP model may not provide an integer feasible solution for the large-size test instances. Finally, the last column contains the average computation time (CPU ) of SA-LS in seconds.
Observing Table 3, we conclude that SA-LS is superior to NN and CWS. SA-LS obtained, on average, 53% and 16% better solutions than those obtained by the NN and CWS algorithms, respectively. The solution time of the NN and CWS algorithms was very short (less than 1 s) because they are simple constructive heuristics; hence, they are not shown in the results. However, SA-LS outperformed the simple heuristics in terms of the solution quality within a reasonable computation time. It is remarkable to note that SA-LS achieved good solutions only within 103 s, on average. As compared to the lower bounds obtained by the MIP formulation, the algorithm performs solution quality less than 15 percent in average. The minimum and maximum deviations from the lower bounds of MIP formulation are 5.57 and 17.55. The gap value might be quite big because the MIP formulation is not as tight as it is necessary. With the addition of some valid inequalities, the formulation could be made tighter and it gives higher lower bounds. It is interesting to observe that both 75 and 100 nodes instances have similar average gap, which is approximately 14%, and also the test instance, in which the maximum gap value is obtained, is 75 nodes instance rather than 100 nodes instance. The test data might be decisive of this issue.

Comparison of SA-LS with several heuristics on HVRP instances.
As is known, the HVRPSPD is a generalization of the HVRP. Thus, the proposed solution methodology can be used to solve HVRP instances by considering the original demand of each customer as the delivery demand (d i = q i ) and the pickup demand as zero (p i = 0). To investigate the performance of SA-LS on the HVRP, we solved the well-known HVRP test instances of Golden et al. [23] and compared the results with several heuristics proposed in the HVRP literature. In these test instances, there were 20 problems with varying size 20 to 100. We excluded the eight instances (problems 1, 2, and 7 to 12) that used non-Euclidean distance. Golden et al. [23] considered the unlimited number of vehicles of each type for their vehicle fleet mix (VFM) problem. Because our algorithm is designed to solve the HVRP with a limited number of vehicles of each type, we calculated the number of available vehicles T k , ∀k ∈ B for each test problem as min k∈B {Q k }/ i∈N q i . Additionally, there was no variable cost θ k per distance unit of vehicle k ∈ B. Hence, we assumed that θ k = 1, ∀k ∈ B. Table 4 presents the results of SA-LS on HVRP instances. In this table, the first and second columns represent the name and size of the HVRP instances, respectively. Next each subsequent column contains the best objective function values obtained by the corresponding author(s). The last three columns contain the best objective function value obtained by SA-LS, percentage deviation (Dev% ) from the best objective calculated as 100(f SA−LS -f Best )/f Best , and average computation time per run in seconds (CPU ) over the five runs. The text in bold indicates the minimum objective function value (f Best ) in a row.
The results for the HVRP test instances are quite promising. SA-LS obtained, on average, 1.36% (average of the Dev % column) worse than the best solutions in the literature found so far. The average computation time was 16.56 s, and it varied between 0.46 s and 60.70 s. According to the results, SA-LS can be a preferable method to solve HVRP instances successfully.   [51] and compared the best results in the VRPSPD literature. In these test instances, there are 14 problems with varying size 50 to 199. To adapt the test problems to be solved by our proposed algorithm, we set and |B| = 1 and θ k = 1, f k = 0, T k = +∞, Q k = Q for each k ∈ B. Table 5 presents the results of SA-LS on VRPSPD instances. In this table, the first and second columns represent the name and size of the VRPSPD instances, respectively. The third column is the capacity of the vehicle in the problem. Next each subsequent column contains the best objective function values reported by the corresponding author(s). The last five columns contain the minimum, average, maximum and the best objective function value obtained by SA-LS over the five runs; the percentage deviation (Dev% ) from the best objective calculated as 100(f SA−LS -f Best )/f Best , and the average computation time per run in seconds (CPU ) over the five runs. The text in bold indicates the minimum objective function value (f Best ) in a row.
The results for the VRPSPD test instances are quite promising. SA-LS obtained, on average, 1.86% (average of the Dev% column) worse than the best solutions in the literature found so far. The average computation time was 31.16 s, and it varied between 2 s and 90 s. According to the results, SA-LS can be a preferable method to solve VRPSPD instances successfully. 5.6. Comparison of SA-LS with the results of Avci and Topaloglu [5]. Avci and Topaloglu [5] are proposed a hybrid metaheuristic algorithm to solve the HVRPSPD. In their study they compare their algorithm with several metaheuristic algorithms. They conduct their experimental computations on randomly generated test instances. To investigate the performance of SA-LS on the HVRPSPD, we also solved the test instances of Avci and Topaloglu [5] and compared with their best results.
In these test instances there are two sets, where each set has 14 problems. The first set has small and medium size problems varying in size 10 to 100. The second set has medium and large size problems varying in size 150 to 550.  Table 6 presents the results of SA-LS on set 1 instances of Avci and Topaloglu [5]. In this table, the first column represents the name of problem. The second and third columns are respectively the size and the number of vehicle types in the problem. Next three columns contain the minimum and average objective function values and average CPU times over several runs, reported by Avci and Topaloglu [5]. The last five columns contain the minimum, average, maximum and the best objective function value obtained by SA-LS over the five runs; the percentage deviation (Dev% ) from the best objective calculated as 100(f SA−LS -f Best )/f Best , and the average computation time per run in seconds (CPU ) over the five runs. The text in bold indicates the minimum objective function value (f Best ) in a row.
The results of SA-LS on small size test instances show that the SA-LS surpass the metaheuristic algorithm proposed by Avci and Topaloglu [5]. The SA-LS obtained, on average, 2.51% (average of the Dev% column) better objective value than Avci and Topaloglu [5]. The 12 out of 14 instances are improved by the SA-LS. Table 7 presents the results of SA-LS on set 2 instances of Avci and Topaloglu [5]. In this table, all columns are same as in Table 6.
The results of SA-LS on large size test instances also show that the SA-LS outperform the metaheuristic algorithm proposed by Avci and Topaloglu (2016). The SA-LS obtained, on average, 19.82% (average of the Dev% column) better objective value than Avci and Topaloglu [5]. All instances are improved by the SA-LS. In some instances the improvement is higher than 20%.    Figure 8 and Figure 9 show the convergence of SA-LS with different initial temperature, final temperature and cooling rate values, respectively. As seen from the figures, the SA-LS slowly convergences to better solutions when the parameters are set to 380, 0.15 and 0.95. This result supports that the parameter values used in the SA-LS help the algorithm to avoid local optima. 5.8. Managerial insights of the results. In today's competitive environment, logistics, transportation and distribution are at the center of managers' concern. Transportation constitutes one-third of the logistics costs in the supply chain. In this regard companies try to create competitive advantages by strengthening their transportation and distribution activities and as well as making their fleet more  As an extension of the well-known VRP, the HVRPSPD is considered in this paper. HVRPSPD is based on the optimization of several important problems such as, fleet management, vehicle assignment and routing, individually and simultaneously. The HVRPSPD presents a dashboard for decision makers where not only the transportation or distribution cost, but also the fleet investment cost is taken into account. We describe a procedure to solve the problems to be able to make right and effective decisions. To this point we give some theoretical explanations and results for the problem. But in real life implementations these results correspond to some managerial insights for the logistics and distribution companies.
The first term of objective function in the MIP formulation is the total transportation or distribution cost and the second term is the fleet investment cost. The second term also reflects the occupation of the vehicles. The size of problem instance means the size of decision problem which the company would face.
The performance indicators used for evaluating the algorithms can also be interpreted from the distribution firm's point of view. Such as, a company can recognize its current state with Gap% or Dev% by examining the current performance with the best one. In business and operations management this topic is considered as gap analysis which is a four step process and helps the organizations to make use of their resources, capital and technology to reach their full potential. With the use of Gap% or Dev% the logistics company can determine where it is today and the management team can create an action plan to move the firm forward and fill in the performance gaps. Another aforementioned percentage indicator is the Imp%. This indicator can be thought as the percentage of improvement in total transportation/distribution and fleet investment costs. On the other hand it shows the decrease in total cost as well as total occupation of the vehicles in fleet.
In decision theory, there is a phenomenon called speed-accuracy trade-off. This fact explains the negative correlation between the quality and speed of any decision made. The constructive solutions methods are quick but may be less effective in solution quality. Besides the SA-LS method (or any exact solution methodology, e.g solution of MIP model ) is slower but can be more effective in finding good quality solutions, relative to the constructive algorithms. But at this point the judgements of decision makers step in. Such that, in strategic level decisions where the consequences of any decision may affect the company in long-term and it is hard to return back, the quality of the solution is more important than the speed. Furthermore in operational level decisions, to be able to act quickly may be vital than to give best decision; since the daily problems need instantaneous solutions or in some cases there might be some opportunities which could not be missed. From the distribution company's perspective the quality of the decision can be measured by the Gap% or Imp% indicators and the speed of the decision can be measured by the CPU time of the solution.
Eventually, based on experimental studies we conclude that operational decisions can be significantly effective in a way that reduces transportation and fleet investment costs and improves the customer satisfaction. 6. Conclusions. In this study, we considered the HVRPSPD. We proposed a polynomial-size node-based MIP formulation for the HVRPSPD. We explained how the MIP formulation could be adapted to solve several variants of the HVRPSPD. Because the HVRPSPD is in the NP-hard problem class, we proposed a hybrid heuristic approach based on SA and LS algorithms called SA-LS to solve medium and large-size problems. We generated the initial solutions with the giant tour approach for SA-LS. We used seven inter-route and four intra-route moving strategies to generate neighbor solutions in SA-LS. The parameters of SA-LS were determined and enhanced based on our preliminary studies.
A series of experiments was conducted to evaluate the performance of the MIP formulation and the proposed hybrid heuristic approach using test instances derived from the literature. The computational results over 620 test instances show that the proposed hybrid heuristic algorithm is computationally efficient for solving the HVRPSPD, and good quality solutions (8.19% on average) can be obtained in a reasonable computation time (approximately 6 s). Moreover, the proposed algorithm can determine solutions better than simple constructive heuristics, the NN and CWS algorithms, up to 50%. SA-LS is also efficient and preferable for solving HVRP and VRPSPD instances, which is a special case of HVRPSPD.
There are several interesting options for further research. The overall procedure could be adapted to solve the problems which are indicated as the special cases of the HVRPSPD. User friendly computer programs, which include the solution approaches proposed in this paper, may be developed with in a decision support framework that can offer service directly to the logistics firms. More sophisticated approaches or exact algorithms such as branch-and-cut, branch-and-price, column generation, which use the proposed formulations and valid inequalities, can be developed to obtain optimal solutions for medium and large-size instances of the HVRPSPD. Another powerful approach called matheuristics, which hybridizes the complex mathematical models with clever heuristic and metaheuristic algorithms, can be developed to solve the HVRPSPD. In real-life situations, there may be conflicting objectives in the problems considered. Thus, for multi-objective HVRP-SPD, algorithms that can determine Pareto optimal solution sets can be developed. Finally, it is possible to extend the problem to consider uncertainty in demand and investigate the effects of very small capacity violations with non-linear penalty functions, as well as the effect of test data to the overall solution quality.
(i, j) are on the first or last order of different routes, then all nodes are transferred to the available route with the minimum cost vehicle type. The algorithm stops when it reaches the end of the sorted savings list and a solution for the HVRPSPD is obtained. As mentioned in the NN algorithm, the load of intermediate nodes may exceed the capacity because the vehicle load depends on the order of nodes visited. Thus, the route with an infeasible vehicle load can be repaired by considering it as a VRPSPD with one vehicle.