A mixed integer linear programming formulation for the vehicle routing problem with backhauls

Article history: Received March 9 2018 Received in Revised Format March 16 2018 Accepted June 14 2018 Available online June 14 2018 The separate delivery and collection services of goods through different routes is an issue of current interest for some transportation companies by the need to avoid the reorganization of the loads inside the vehicles, to reduce the return of the vehicles with empty load and to give greater priority to the delivery customers. In the vehicle routing problem with backhauls (VRPB), the customers are partitioned into two subsets: linehaul (delivery) and backhaul (pickup) customers. Additionally, a precedence constraint is established: the backhaul customers in a route should be visited after all the linehaul customers. The VRPB is presented in the literature as an extension of the capacitated vehicle routing problem and is NP-hard in the strong sense. In this paper, we propose a mixed integer linear programming formulation for the VRPB, based on the generalization of the open vehicle routing problem; that eliminates the possibility of generating solutions formed by subtours using a set of new constraints focused on obtaining valid solutions formed by Hamiltonian paths and connected by tie-arcs. The proposed formulation is a generalpurpose model in the sense that it does not deserve specifically tailored algorithmic approaches for their effective solution. The computational results show that the proposed compact formulation is competitive against state-of-the-art exact methods for VRPB instances from the literature. © 2019 by the authors; licensee Growing Science, Canada


Introduction
The vehicle routing problem with backhauls (VRPB) has been extensively defined in the literature (Toth & Vigo 1999;Mingozzi et al., 1999;Osman & Wassan, 2002;Ropke & Pisinger 2006) and can be stated as a problem in determining a set of vehicle routes visiting all customer vertices, which are partitioned into two subsets.The first subset contains the vertices of the linehaul customers (LCs), where a given quantity of product is to be delivered.The second subset contains the backhaul customers (BCs), where a given quantity of product to be collected and transported to the depot.The objective is to consider the routes performed from the depot to the customers by a fleet of homogeneous vehicles, in order to satisfy the demand of the customers (products to be collected or products to be delivered).In such a case, the vehicles must first attend the customers with delivery requirements before the customers with collection requirements.For some transportation companies it is critical to avoid the reorganization of the products inside the vehicles at each delivery point.The pickups and deliveries of goods in a mixed order, or simultaneously, cause difficulties due to the rearrangements of goods on board.The VRPB adequately represents this strategic need and must satisfy the following conditions:  Each vertex must be visited exactly once by a single route.That is, each vertex is grade 2.  Each route starts and _nishes at the depot. Each customer must be fully attended when visited. All customers are serviced from a single depot. The vehicle capacity should never be exceeded in both the linehaul and backhaul route and all vehicles have the same capacity. In each circuit the linehaul vertices precede the backhaul vertices, if any.That is: o A circuit of only BCs is not allowed.o The last customer of a linehaul route is always connected with the depot or with BC who is starting a backhaul route.o The last BC of a backhaul route is always connected with the depot.o The precedence constraint is also justi_ed by the need to attend to LCs with higher priority than BCs.
Fig. 1 shows the optimal solution of a VRPB example with 20 customers; in which the first 10 customers are LCs and the other 10 are BCs.For this case, the capacity of all vehicles is the same and equal to Q = 60.The minimum number of vehicles needed to serve all the linehaul and backhaul customers is known in advance and is indicated by KL and KB, respectively.These values can be obtained by solving the binpacking problem instances associated with the corresponding customer subset, which calls for the determination of the minimum number of bins, each with capacity Q, needed to serve all customers (Toth & Vigo, 2002).To ensure feasibility, we assume that the number of vehicles needed K must be greater than or equal to the maximum value between KL and KB, i.e., K = max {KL, KB}.Thus, in this example, the minimum number of vehicles needed to serve all the linehaul and backhaul customers is KL=KB = 3.
Because the VRPB is NP-hard in the strong sense, it can be solved as a matching problem (Toth & Vigo, 2002); it has two types of routes and has a precedence constraint, so many heuristic processes are appropriate and efficient for the solution.Therefore, most existing literatures on the VRPB are related to heuristics and metaheuristics methodologies with high quality results.A comprehensive review of metaheuristic techniques for VRPB is found in Ropke and Pisinger (2006).Two literature reviews cover the main works about VRPB: the first, presented by Toth and Vigo (2002), presents the existing work up to 2002 and second, by Irnich et al. (2014a,b), focuses on complementing the review up to 2013.Concerning the exact approaches, few jobs have been proposed (Toth & Vigo, 2014).In our review, only two works were found: the first exact method is reported by Toth and Vigo (1997), in which an effective Lagrangian bound is introduced that extends the methods previously proposed for the capacitated VRP (CVRP).The resulting Branch-and-Bound algorithm is able to solve problems with up to 70 customers in total.The second exact method is proposed by Mingozzi et al. (1999), in which a set-partitioningbased approach is presented and the resulting mixed integer linear programming (MIP) is solved through a complex procedure.The results show that the approach is capable of solving undirected problems with up to 70 customers.Toth and Vigo state that no exact approaches have been proposed for VRPB during the last decade (Toth & Vigo, 2014).In our review, we have reached the same conclusion and new proposals for unified exact models of VRPB were not found, since the only two existing proposals are used to derive the relaxations on which the exact approaches are based (Toth & Vigo, 1997).
Other two problems in the literature commonly handled by exact methods, where the backhaul load is considered, are: i) the mixed vehicle routing problem with backhauls (MVRPB) and ii) simultaneous pickups and deliveries.In the first, deliveries after pickups are allowed where the linehaul and backhaul customers are mixed along the routes.In the second, customers may simultaneously receive and send goods.Although the differences between these two problems and the VRPB appear to be subtle, they are very different; direct comparisons between the problems serving pickups and deliveries in a mixed order or simultaneously with problems that deliver first and pick-up second should not be performed, since they are addressing different requirements.The VRPB is a problem with a special structure of the routes that consist of two distinct parts; a delivery and a pickup segment.A complete review of these two types of problems can be found in Ropke and Pisinger (2006), Wade and Salhi (2003) and Parragh et al. (2008).
More recently, Chávez et al. (2016) proposed a Pareto ant colony algorithm to solve a multi-objective variant of the multidepot VRPB where the aim is to minimize distance, travel time and energy consumption.A random fixed speed between 30 km/h and 90 km/h was assigned to each arc, and the function considered by Bektaş and Laporte (2011) was used to compute energy consumption.The algorithm is based on the idea of Doerner et al. (2004) which uses three matrices of pheromones for each objective function.The method was tested on new instances based on those of Salhi and Nagy (1999).In Chávez et al. (2018) a heuristic algorithm based on Tabu Search Approach for solving the VRPB is proposed.The proposed algorithm considers intensification of local information and the proper exploitation of local search procedures combined with exploratory stochastic strategies.The solution strategy considers the design of separate routes for both delivery and collection of goods.
A metaheuristic algorithm based on the ant colony optimization was presented by Chávez et al. (2015), to solve the multi-depot vehicle routing problem with delivery and collection of package.Each performed route consists of one sub-route in which only the delivery task is done, in addition to one sub-route in which only the collection process is performed.The proposed algorithm tries to find the best order to visit the customers at each performed route.
A recent survey paper with interesting conclusions and research perspectives on the VRPB, including models, exact and heuristic algorithms, variants, industrial applications and case studies, are identified in Koç and Laporte (2017).In this review, the authors highlight the importance of using matheuristic algorithms that allow the interoperation of metaheuristic and mathematical programming techniques.Additionally, they identify the need for new studies focused on developing effective and powerful exact methods to solve all available standard VRPB instances to optimality.The authors also conclude that no electric vehicle version has yet been studied for the VRPB.
In VRPB the precedence constraint, which stipulates that in each circuit the linehaul vertices precede the backhaul vertices, difficult to build an exact model from the traditional viewpoint of the CVRP.This is because traditional restrictions for the elimination of sub-tours perfectly fit into VRPs with a unique set of vertices, where the evaluation of constraints of degree or conservation of flow can be evaluated in a general way on all vertices.Adapting these restrictions to the VRPB involves special cases such as vertices at the end of a linehaul route, vertices at the start of a backhaul route and routes with linehaul customers only.
Thus, we have approached the problem from another point of view; considering a representation of each part of the problem based on a generalization of the open vehicle routing problem (OVRP) with a set of new constraints focused on maintaining the arborescence condition of both paths separately.Additionally, another set of constraints to handle tie-arcs coupling the two types of routes are proposed.
The OVRP was first proposed in the early 1980s (Schrage, 1981;Bodin et al., 1983), when there were cases where a delivery company did not own a vehicle fleet or its fleet was inadequate for fully satisfying customer demand.Therefore, contractors who were not employees of the delivery company used their own vehicles for deliveries.In these cases, the vehicles were not required to return to the central depot after their deliveries because the company was only concerned with reaching the last customer.
Compensation was not given for any driving outside of meeting this goal.Thus, the goal of the OVRP is to design a set of Hamiltonian paths to satisfying customer demand.
In practice, the OVRP formulation represents situations such as: home delivery of packages and newspapers, school bus routing, routing of coal mines material, and the shipment of hazardous materials (Braekers et al., 2015).Thus, the VRPB structure can be seen as OVRPs of linehaul and backhaul routes connected by tie-arcs.The proposed model is a general-purpose model in the sense that it does not deserve specifically tailored algorithmic approaches for their effective solution and can be solved by an integer-programming solver.The main contributions of this paper can be summarized as follows:  A unified and compact model for the VRPB is proposed, which can be a starting point for the generalization of problems shortly discussed in the literature, as are the multi-depot VRPB and the location VRPB. This paper presents a contribution to the discussion on VRPB and its feature from a new approach based on arborescence, which allows the best advantage of the structure of the problem.The proposed model can be used to derive new relaxations on which the exact approaches are based. The proposed formulation allows to solve the symmetric and asymmetric VRPB and it is able to minimize the number of used vehicles.
The rest of the paper is organized as follows: in Section 2 we first describe the problem formulation, presenting the nomenclature for the variables and parameters used in the mathematical model.We then introduce the new mixed integer linear programming (MILP) formulations based on the arborescence condition (MILP-AC) for the VRPB.We describe how the arborescence constraints operate on the different structures of the problem.In Section 3 we present a computational study performed on 142 test instances.Finally, the conclusions in Section 4 are presented.

Nomenclature
The nomenclature for the variables and parameters of the proposed model for the VRPB is summarized next.Sets:

Cij
Cost of traveling between nodes I and j.

Dj
Nonnegative quantity of product to be delivered or collected (demand) of the customer j ∈ CU.

KL, KB
Minimum number of vehicles needed to serve all the linehaul and backhaul customers, respectively.Q Capacity of the vehicles (identical vehicles).
Variables: sij Binary variable for the use of the path between nodes i, j ∈ V. ξij Binary variable for the use of the path between nodes i ∈ L and j ∈ B0 ( tie-arcs ).lij Continuous variable indicating the amount of cargo transported between nodes i and j.

The VRP with Backhauls
Note that in the optimal solution shown in Fig. 1 the linehaul routes, excluding the tie-arcs, constitute a subproblem that has a radial configuration (arborescence) starting from the depot, spanning all the linehaul vertices and ending up at a client.This subproblem we have named linehaul open vehicle routing problem (LOVRP).Similarly, the backhaul routes also have a radial configuration, entering the depot and spanning all the backhaul vertices; this sub-problem is named backhaul open vehicle routing problem (BOVRP).Thus, the arborescence characteristics allow to handle the VRPB as the solution of two open routing subproblems (ORPs) connected through a tie-arc.In the OVRP context, the necessary condition for obtaining a minimum spanning tree is that the number of arcs be equal to the number of customer nodes, which is given by the cardinality of the sets L and B for the case of the LOVRP and BOVRP, respectively.However, this condition is necessary but not sufficient because two situations can be presented: i) there may be customer nodes with a degree greater than two and ii) disconnected solutions can be obtained.To avoid the first situation, it must be considered that a spanning tree becomes a subgraph formed only by Hamiltonian paths if each customer node has a degree less than or equal to two.
Concerning the second situation, the addition of a balance condition of the demand flow by each customer node avoids getting disconnected solutions.These situations are considered in the model presented below.

Proposed Model for the VRPB
The two-index vehicle flow formulation for the VRPB is defined as follows: The objective function (1) minimizes operating costs, which correspond to the sum of the total travelling cost of the routes used to deliver and collect the goods to the customers and the total travelling cost associated with use of the tie-arcs connecting the last customer of a linehaul route with the first customer of a backhaul route or with the depot.The set of constraints ( 2)-( 7) model the LOVRP, where ( 2) and ( 3) impose the arborescent connectivity requirements.More precisely, these two constraints allow configuring one shortest spanning arborescence linehaul with fixed indegree KL at the depot vertex.In the optimal solution of the LOVRP, each route has an arborescent configuration formed by a minimum spanning tree; starting from the depot, spanning all the nodes, and ending at a customer.A comprehensive discussion on the application of concepts of graph theory for the formulation of spanning tree constraints is found in works related to optimizing the operation of distribution systems (DS) of electric power, which has been an active topic for years, with recent emphasis on smart grid initiatives.Distribution systems are most commonly operated in radial configurations in order to keep the system operation as simple as possible.In this context, a radial configuration is equivalent to a spanning tree, where there is only one path between the electrical substation and final consumers.The resulting subgraph must be connected, without cycles and with a number of arcs that is equal to the number of demand nodes.A brief review of existing approaches for imposing spanning tree constraints in the operation of the DS, can be found in Ahmadi and Martí (2015).Thus, several of the concepts of graph theory used in DS can be used in the LOVRP since both problems are quite similar.In the LOVRP context, the necessary condition for obtaining a minimum spanning tree is that the number of arcs be equal to the number of customer nodes.This necessary condition is guaranteed by the constraint (2), where the number of customer nodes is given by the cardinality of the set L. However, this constraint is necessary but not sufficient because there may be customer nodes with a degree greater than two and disconnected solutions can be obtained, as shown in Fig. 2(a).A spanning tree becomes a subgraph formed only by Hamiltonian paths if each customer node has a degree less than or equal to two.Therefore, another necessary condition is given by the set of degree constraints (4) and ( 5).The indegree constraints (4) impose that exactly one arc enters each vertex associated with a LC and, consequently, the outdegree constraints ( 5) impose that exactly one arc leaves each customer, except for those customers who are at the end of the route where the tie-arcs emerging from LC to a BC or to the depot should be considered.However, the addition of these degree constraints in directed graphs may not represent a spanning tree, because a disconnected graph can be obtained, as shown in Fig. 2(b).The addition of flow balance constraint by each customer node avoids getting disconnected solutions, since an infeasible solution is obtained when the goods leaving the depot can not reach the customers.Thus, the set of constraints reported in (3) guarantees network connectivity through the balance of the demand flow by each customer so that they are fully served when visited, and also ensures that the remaining demand of a vehicle is zero after the vehicle has serviced its last customer.The impact of these constraints is shown in Figure 2 2)-( 5) allow obtaining a linehaul arborescence structure as shown in Fig. 2(c), whose main characteristic is that each vertex is of degree 2 when tie-arcs are considered.Constraint ( 6) and ( 7) impose both the vehicle and depot capacity requirements, respectively, associated with the linehaul routes only.The first is an upper limit defined by the capacity of the vehicle to transport a given quantity of product on any linehaularc, while the second is a lower limit to the number of routes out of the depot to supply linehaul customers, which is determined by the ratio between the total demand to be collected and the vehicle capacity.Constraint (7) limits the minimum number of vehicles used on linehaul routes.When there is a choice of set partitioning or set covering as a formulation, set covering is preferred (Barnhart et al. 1998).Similarly, the set of constraints ( 8) -( 13) model the BOVRP, were (8) and ( 9) impose the antiarborescent connectivity requirements.Note that (9) guarantees the balance of demand flow in each BC, so that the product is fully collected when the customer is visited.Out-degree constraint (10) imposes that exactly one arc leave each vertex associated with a BC.In Figure 3(a), an example of the shortest spanning arborescence linehaul, with fixed indegree KL at the depot vertex, is shown in solid lines.Similarly, the shortest spanning antiarborescence backhaul, with fixed indegree KB at the depot vertex, is shown in dashed lines.In Figure 3(b), in solid lines, the impact of constraints ( 2)-( 5) is shown, allowing generating a linehaul arborescence structure.In dashed lines, the impact of constraints ( 8)-( 10) is shown, allowing generating a backhaul antiarborescence structure.
The in-degree constraint (11) imposes that exactly one arc enter each backhaul vertex, including tie-arcs coming from the LCs.Constraint ( 12) and ( 13) impose both the vehicle and depot capacity requirements, respectively, associated with the backhaul routes only.Constraint (13) limits the minimum number of vehicles used on backhaul routes.Comparing equations ( 13) and ( 7) one can see that the number of linehaul arcs leaving the depot may be different from the number of backhaul arcs arriving at the depot.This case occurs when there are tie-arcs between a linehaul route and the depot.Thus, constraint ( 14) ensures that the number of arcs leaving the depot should be equal to the number of arcs coming to depot either a tie-arc or an arc from a BC.Constraints ( 5), ( 11) and ( 14) are responsible for coupling the two problems (LOVRP and BOVRP) through the tie-arcs and to ensure the characteristic of precedence.
Because the symmetric and asymmetric versions of VRPB are considered, the uni-directional constraint (15) ensures that only one of the two variables sij or sji must be used.Finally, constraints ( 16) and ( 17) define all binary decision variables, and constraint (18) defines the real variable.Ropke & Pisinger (2006).b TV: exact algorithm proposed by Toth & Vigo (1997).Computing times in Pentium 60 MHz.Time limit of 6,000 seconds.c EHP: exact algorithm proposed by Mingozzi et al. (1999).Computing times in SGI 200 MHz.Time limit of 25,000 seconds.d MILP-AC: proposed mixed integer linear programming with radility condition formulation.Computing times in PC intel i3/3.3GHz.Time limit of 6600 seconds.e K(KLB): K is the specified number (given in advance) of vehicles to use and KLB is the number of vehicles performing routes conformed by linehaul and backhaul customers.Thus, K − KLB is the number of vehicles performing routes conformed by linehaul customers exclusively.z*: value of the best solution found by the earlier algorithms.LB(%): percentage error of the lower bound (LB) computed at the root node.It is calculated as the ratio of the LB divided by the best z* and multiplied by 100.f Gap (%): percentage gap is calculated as (CPLEX − LB)/LB.†: optimality proven for the first time.† †: new BKS.Time: overall computing time expressed in CPU seconds.Bold numbers are the new best-known solution values by exact method.
To reduce the number of binaries and to speed up computational time we apply the preconditioning valid constraints ( 19), ( 20) and ( 21) to remove unfeasible arcs.
The constraint (19) guarantees that there is no connection from one BC to a LC.The beginning of a route cannot be through a BC, because a route consisting exclusively of BCs is not allowed, which is ensured through (20).Finally, constraint (21) ensures that the tie-arcs linking a LCi with a BCj or with the depot should only be handled by the variable ξij.

Computational results for the exact algorithms
Two datasets of scenarios are used in order to show the operation and effectiveness of the proposed formulation.The first dataset, denoted as GJ dataset, was proposed by Goetschalckx & Jacobs-Blecha (1989) and contains 62 instances with a range between 20 and 150 customers.Details on how these scenarios were generated can be consulted in Toth & Vigo (2002).The second dataset, denoted as TV dataset, was proposed by Toth & Vigo (1997) and contains 33 instances between 21 and 100 customers.
The proposed model corresponds to a MIP formulation and were implemented in AMPL (Fourer et al., 2002) and solved with GUROBI 6.5 (called with the optimality gap option equal to 0%), on an Intel Core i3 computer; 3.3 GHz, 3.8 GB of RAM.
In Tables 2 and 3, the results for the VRPB using the GJ and TV datasets, respectively, are compared with those obtained by exact methods proposed by Toth and Vigo (1997) and Mingozzi et al. (1999).
Additionally, the best-known solutions (BKS) reported by the heuristic method proposed by Ropke & Pisinger (2006) is presented in the fourth column of each table, which are the state of the art of heuristic methods for VRPB instances from the literature.In Table 2, we give in columns 1-3: the problem name, the number of LCs and the number of BCs, respectively.Column 4 presents the BKS obtained by a heuristic procedure (Ropke & Pisinger, 2006).In columns 5-7 and 8-10 the results are presented for TV and EHP algorithms, respectively.
For each of these algorithms the table reports two types of information: i) the value of the best solution found by the algorithm, z*, and ii) the overall computing time expressed in CPU seconds.The Euclidean distances were rounded to one decimal and the final result was rounded to an integer.Finally in columns 9-12 the characteristics of the solution for each problem are presented, using the proposed formulation in this paper (computing times in PC Intel i3/3.3GHz and time limit of 6600 seconds).
The optimality for the GJ dataset was proved for the first time for 8 instances.One new best-known solutions was found considering both heuristic methods and exact methods and 12 new BKS were found in the tests considering a comparison between exact methods only.Optimality was achieved in 42 of the 47 scenarios.As an example, in Figure 3a the optimal solution for the instance K4 is shown.
Table 3 presents the same configuration of Table 2.For this dataset the optimality of 4 TV instances is proven for first time.Two new best-known solutions were found considering both heuristic methods and exact methods and 8 new best-known solutions values were found in the tests considering a comparison of exact methods only.Optimality was achieved in 27 of the 33 cases.It is important to note that for some instances, a better value of the objective function can be obtained with a smaller number of vehicles.This causes an increase in computational time as shown in Table 4, which reports the results for three instances where the number of vehicles can be reduced.In Table 5, the results are obtained with different rounding scheme.The Euclidean distances were rounded to integers.This rounding scheme has only been used in the literature to compare heuristic methods, Therefore the gap for all instances is reported for the first time.The exact solutions are reported and compared with the BKS obtained by a heuristic methodology.Optimality was achieved in 47 of the 62 cases, five new best-known solutions were found.As an example, in Fig. 3b the optimal solution for M4 instance is shown.The results show that for the GJ dataset the computing time for the optimal solution is less in those instances where the specified number of vehicles to use (K) is equal to the number of vehicles performing routes conformed by linehaul and backhaul customers (KLB ) and where the capacity of the vehicle is greater.This feature is obviously because the problem is less restricted under these conditions and this can be exploited in future works using heuristics techniques or to formulate new relaxations in exact approaches based in branch-and-bound techniques.For example, the solution of the BOVRP, which is simpler than the LOVRP under the above conditions, could be an interesting starting point for a more elaborate methodology.As can be seen from the computational results, the proposed model produces high quality results, obtaining equal or better upper bounds in all instances, and the final lower bounds prove stronger than those obtained by earlier methods.Ropke & Pisinger (2006).b Reference: the heuristic algorithm reporting the result is indicated in this column.TV refers to the heuristic algorithm by Toth & Vigo (1999), OW refers to the heuristic by Osman & Wassan (2002) and LNS refers to the heuristic by Ropke & Pisinger (2006).

Fig. 2 .
Fig. 2. Impact of the spanning tree constraints in the LOVRP subgraph (c), where the demands of each customer are shown with the notation (•) and the amount of goods owing through the arcs are shown with the notation [•].Thus, the constraints ( (a) Linehaul arborescence and backhaul antirborescence.(b) Hamiltonian paths when degree constraints are considered.

Table 1
Coordinates and loads for 20 customers

Table 2
Goetschalckx and Jacobs-Blecha (1989)ses fromGoetschalckx and Jacobs-Blecha (1989), the Euclidean distances were rounded to one decimal and the final result was rounded to an integer

Table 3
Toth & Vigo (1997)lts for the VRPB cases fromToth & Vigo (1997)The nomenclature of this table is the same as that presented in Table2.Time limits were 18000, 25000 and 1000 seconds for methods TV, EHP and MILP-AC, respectively.

Table 4
Examples of results minimizing the number of vehicles

Table 5
Goetschalckx & Jacobs-Blecha (1989)cases fromGoetschalckx & Jacobs-Blecha (1989)with the Euclidean distances rounded to integers BKS: best solution values obtained by heuristic algorithm and reported in a