An MIP formulation for the open location‐routing problem considering the topological characteristic of the solution‐paths

In the open location‐routing problem (OLRP), one has a set of candidate depots to be installed and the vehicles start from the depot, visit all customers, and are not required to return to the depot after completing their service. Thus, the OLRP involves the problems of facility location and open vehicle routing. In this paper, a new mixed integer programming formulation for the OLRP is presented by proposing a set of constraints to obtain valid solutions formed by a graph consisting of a spanning tree in each connected component of the graph. This approach results in an alternative way of avoiding generating subtours, which significantly simplifies the set of constraints associated with the connectivity of the solution and the vehicle capacity requirements. The computational results show that the proposed formulation is competitive against state of‐the‐art methods for this type of problems.


INTRODUCTION
The open vehicle routing problem (OVRP) was first proposed in the early 1980s when there were cases where a delivery company did not own a vehicle fleet or its private fleet was inadequate for fully satisfying customer demands [8,39]. 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. Thus, the goal of the OVRP is to obtain a spanning tree to satisfying customer demands.
However, when considering that the vehicles depart from a set of available depots and finish their deliveries on their last customer, the OVRP is generalized to the multidepot OVRP (MDOVRP), where the depots are already installed and have no opening costs. Additionally, when the facility location is considered, the resulting problem is known as the open location-routing problem (OLRP), where one has a set of candidate depots to be installed and where two decisions must be made: (i) which facilities should be opened, and (ii) which vehicle routes should be built. Figure 1 shows the solution of an OLRP with six customers and two depots. Depots are denoted with rectangles and customers with circles and only the depot 1 is opened in this solution.
Thus, the OLRP consists of finding a collection of paths of minimum traveling cost connected to a set of depots with minimum setup cost to satisfy the following conditions: • Each customer has to be fully served when visited.
• The sum of the demands of the customers visited by a route cannot exceed the vehicle capacity.
Among the techniques traditionally proposed to solve the OLRP, it is a common practice to adapt models designed to obtain Hamiltonian circuits to obtain Hamiltonian paths by handling the input data of the model. This is achieved, for example, by using the strategy of adding dummy nodes or by setting to zero the cost of all arcs that have the depot as the endpoint. These data adaptation strategies have demonstrated that traditional models are able to satisfactorily solve the OLRP. However, in this paper it is shown that the solution can be obtained more efficiently when the topological characteristic of the solution-paths is taken into account from the very formulation of the problem.
The main contributions of this paper are: • Provide a new mixed integer linear formulation for the OLRP that efficiently eliminates the possibility of generating solutions formed by subtours in order to provide more competitive results. This is achieved by proposing a set of new constraints focused on obtaining valid solutions formed by a spanning tree in each connected component of the graph. • For purposes of analyzing the impact of the new restrictions on other problems, computational results are presented on the MDOVRP. Additionally, a set of constraints that allows finding more balanced routes with respect to the length traveled is presented. Thus, a distance-constrained MDOVRP is formulated, where the difference between the length of the maximum and the minimum paths should not exceed a given threshold.
The rest of the paper is organized as follows: Section 2 briefly discusses research work on specific problems that are most closely related to the issues addressed. Section 3 describes the mathematical models of the problems addressed, where the OLRP is presented first due to its generality; subsequently, this model is adapted for the MDOVRP and, finally, a distance-constrained MDOVRP (DCMDOVRP) is formulated. In Section 4, the computational results to compare the performance of the proposed models for different test instances are depicted. Finally, Section 5 presents our conclusions.

LITERATURE REVIEW
The multidepot vehicle routing problem (MDVRP) is an extension of the classical capacitated vehicle routing problem (CVRP) in which there are several depots where vehicles can start and end their routes. The literature on the exact approaches to the MDVRP is sparse. In the last decade, to the best of our knowledge, five works concentrate on the main proposals. The first of these works is that of Kek et al. [22], who proposed a linear programming (LP) model and a branch-and-bound procedure for distance-constrained CVRPs with flexible assignment of the start and end depots. The objective function is to minimize the total cost of the routes.
In [4], Baldacci and Mingozzi proposed a mathematical formulation for solving several classes of vehicle routing problems, including the MDVRP, based on the set partitioning formulation. The exact algorithm uses three types of bounding procedures based on LP-relaxation and on the Lagrangean relaxation of the mathematical formulation. A mixed integer LP model for the multidepot petrol station replenishment problem with time windows has been proposed by Cornillier et al. [16], in which a heterogeneous fleet of vehicles is available with the maximization of total net revenue as the objective function, with maximum and minimum demands.
Braekers et al. [9] propose an exact approach based on the branch-and-cut algorithm for a general heterogeneous dial-a-ride problem with multiple depots. In this work, benchmark instances are presented for the multidepot version of the heterogeneous dial-a-ride problem. The most recent exact method reporting results on the MDVRP is that of Contardo and Martinelli [14], in which the capacitated MDVRP with route length constraints is studied. The MDVRP is formulated using a vehicle-flow and a set-partitioning formulation, both of which are exploited at different stages of the algorithm. A comprehensive literature review on the MDVRP is presented in [31].
When in the MDVRP a set of candidate depots to be installed is considered; the resulting problem is known as the capacitated location-routing problem (CLRP). From the viewpoint of the exact methods, Belenguer et al. (2011) present a branch-and-cut algorithm using a two-index formulation. The authors introduce a number of valid inequalities to strengthen the LP relaxation and the largest instances solved to optimality have 5 facilities and 50 customers. Baldacci et al. [5] present an exact algorithm that exploits the following three facts: (i) the number of facility subsets to be considered can be drastically reduced, (ii) the MDVRP is much easier to solve than the CLRP, and (iii) a good lower bound for the MDVRP can be used. Contardo et al. [12] present a computational comparison of four different flow formulations for the CLRP and introduce new branch-and-cut algorithms for each of the formulations. Additionally, an exact algorithm based on cut-and-column generation is proposed by Contardo et al. in [13]. The success of both the algorithm by Baldacci et al. [5] and the one by Contardo et al. [13] depends on the quality of the computed bounds and the ability to reduce the number of facility subsets that have to be considered. Locating and routing problems are important issues and many papers have been published to address them, mainly using heuristic methods for their solutions. Several high-quality survey articles on this problem can be consulted [2,17,29,35]. The most recent survey article, which includes an analysis of all the above, is presented by Schneider et al. in [37]. One of the most recent paper in the field of the CLRP corresponds to the work of Boccia et al. [7], where an integer LP model for a new multicommodity location-routing problem is formulated, considering three decision problems: (i) location, (ii) assignment, and (iii) routing.
In the OVRP, a vehicle does not return to the depot after servicing the last customer on a route. In practice, the OVRP formulation represents situations, such as home delivery of packages and newspapers, school bus routing, routing of coal mines material, or shipment of hazardous materials [10]. The OVRP has recently received increasing attention in the literature and has focused mainly on the development of heuristic methods to find good quality solutions quickly.
Regarding the exact methods, a branch-and-cut algorithm is proposed by Letchford et al. [24] for the open version of the CVRP, addressing the capacitated problem with no distance constraints. In [32], Pessoa et al. present several branch-cut-and-price algorithms on a number of vehicle routing problem variants, among which is the capacitated OVRP, which is addressed by setting the cost of all arcs that have the depot as the endpoint to zero. Salari et al. [36] proposed a heuristic improvement procedure for the OVRP based on integer LP techniques to improve a feasible solution of a combinatorial optimization problem. Alinaghian et al. [3] proposed a mathematical model in which open paths are used into the problem of cross-docks. To model the open path, a dummy node is defined, whose distance to other nodes is considered zero, and from which the route starts. A comprehensive literature review on the OVRP is presented in [25].
Regarding the MDOVRP, Tarantilis and Kiranoudis [40] propose it for the first time within the context of food distribution. Since this problem is a generalization of the OVRP, the MDOVRP is of type NP-hard and, therefore, heuristic methods are commonly used to find good quality solutions quickly [18,21,27,30,44]. In [27], Liu et al. proposed a hybrid genetic algorithm for the MDOVRP; however, the authors additionally proposed a mixed integer programming (MIP) formulation, which is the first formal mathematical definition of the problem. The most recent MIP formulation was reported by Lalla-Ruiz et al. [23] in which a new MIP formulation was presented and compared with the unique MIP formulation given in the literature at that time [27]. The authors focused on the improvements and extensions of the Miller-Tucker-Zemlin subtour elimination constraints, always considering that each vehicle departs from one of the available depots and finishes at the last customer it serves.
To the best of our knowledge, the MDOVRP taking into account distance constraints has not been reported in the literature. However, the OVRP with distance constraints has been addressed by several authors using heuristic techniques [11,20,25,33]. The traditional approach of distance constraints is to consider that the distance traveled by any vehicle may not exceed a prespecified upper limit [19]. However, when additionally a lower limit on the length of the route is considered, then it is possible to say that the tours found may become more balanced from the point of view of the distance traveled. In Li and Fu [26] an open routing problem of school buses is proposed where the main objective is to find the set of school bus routes that satisfies the bus capacity, to minimize the total number of buses required, the total travel time, and to balance the loads and travel times between buses.
The OLRP was introduced by Yu and Lin in [45], who also proposed an exact mathematical model to solve instances with up to 50 clients and 5 depots with gaps close to zero. This model can optimally solve small OLRP instances and verify the performances of heuristic solution approaches for the problem. Additionally, they proposed a heuristic algorithm based on simulated annealing (SA) approach, which consists of two stages. First, an initial solution is constructed with a greedy approach, and second the initial solution is improved with a SA procedure with three search mechanisms: exchange, insertion and 2-opt movements. The proposal is validated with adapted cases from Barreto et al. [6], Prins et al. [34], and Tuzun and Burke [43] and the computational results indicate that the proposed heuristic method solves the OLRP efficiently.
Schneider and Löffler [38] introduce a tree-based search algorithm (TBSA), to tackle large-scale instances of the standard location routing problem. The strategy is to use granular searches in order to reduce the arc set for each facility and customer nodes instead of the shortest arcs in the entire arc set. One of the variants of the LRP addressed in this work is the OLRP, for which it provides the state-of-the-art review on this problem addressed with metaheuristic techniques.
From this review, it can be concluded that the contributions in the literature on the OLRP through exact techniques are scarce and focus mainly on the adaptation of techniques traditionally proposed to solve the OLRP, it is a common practice to adapt models designed to obtain Hamiltonian circuits to obtain Hamiltonian paths by handling the input data of the model. Additionally, the literature on exact approaches for the MDOVRP is sparse and has received very little attention from researchers. Finally, the MDOVRP with distance constraints has not been reported in the literature. These three problems are strongly related by their topological characteristics based on arborescences. To fill this gap, this work presents a new formulation focused on obtaining valid solutions formed by a spanning tree in each connected component of the graph. This formulation is used transversely in the different types of problems addressed.

PROPOSED MODEL FOR THE OLRP
The OLRP can be defined as the following graph theoretic problem. Let G = (V, A) be a complete and directed graph, where V = I ∪ J is the vertex set and A is the arc set. Vertex set I = {1, 2, …, m} represents the set of candidate capacitated depots to be installed. Vertex set J = {m + 1, m + 2, …, m + n} represents the customers to be served. Each customer j ∈ J is associated with a known nonnegative demand of goods D j to be delivered. Each candidate depot has a fictitious demand equal to zero and has an unlimited fleet of identical vehicles with the same positive capacity, denoted as Q. Note that for feasibility, the vehicle must meet the criterion that 0 < D j ≤ Q for all customers j. In addition, each candidate depot i has a limited storage capacity of goods and a set-up cost, respectively denoted as W i and O i . A nonnegative traveling cost C ij is associated with each arc (i, j) ∈ A and a nonnegative fixed cost F is associated with the use of a vehicle, which is a unique value because they are identical vehicles. C ij = U ⋅ T ij where U is the cost per unit of distance traveled and T ij is the distance between the vertices i and j.
The model is a two-index vehicle flow formulation that uses two binary decision variables: ∈ A belongs to the optimal solution and takes value 0 otherwise.
Additionally, a real variable l ij that represents the flow of goods transported by a vehicle through the arc (i, j) ∈ A is used. The two-index vehicle flow formulation for the OLRP is defined as follows: The objective function (1) minimizes the operating costs, which correspond to the sum of the setup cost of depots, the total traveling cost of the routes used to deliver the goods to the costumers, and the fixed cost associated with the use of each vehicle.
In the optimal solution of the OLRP a spanning forest is obtained, where each connected component is a subgraph that is a minimum spanning tree; starting from the depot and covering all the nodes of the subgraph. 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. DS 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 [1]. Thus, several of the concepts of graph theory used in DS can be used in the OLRP since both problems are quite similar.
In the vehicle routing problem 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 J. 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. The set of degree constraints (4) impose that exactly one arc enters to each customer node. The outdegree constraints (5) impose that exactly one arc leaves each customer, except for those customers who are at the end of the route. This set of constraints strengthen the formulation, but the addition of they, by themselves, does not guarantee obtaining a spanning tree because a disconnected graph can be obtained. 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, since the vehicle does not return to the depot. Thus, in the case of a single central depot, the constraints (2)-(5) allow to obtain a subgraph that is a spanning tree and in the case of multiple depots they allow to obtain spanning forests consisting of a spanning tree in each connected component of the subgraph.
Constraints (6) and (7) impose both the depot and vehicle capacity requirements, respectively. Constraints (8) and (9) impose capacity requirements of flow of goods transported. Finally, constraints (10) and (11) define all binary decision variables and the set of constraints (12) defines the real variables.
From constraints (10) and (11) and considering that I = {1, 2, …, m}, J = {m + 1, m + 2, …, m + n} and V = I ∪ J, the total number of binary variables Bn, used in the model, can be calculated as Bn = (m + n) 2 + m. Similarly, considering the constraints (2)-(9), the number of rows Rn can be calculated as Rn = 3 + 3n + m + 2 ⋅ (m + n) 2 . A comparison of our model with the exact model presented by Yu and Lin in [45], in terms of these two characteristics (binary variables and rows) and based on 20 OLRP instances adopted from the literature [6], is shown in Figure 2. As can be seen, our model uses a similar number of binary variables but considers fewer rows than the model proposed earlier.

Multidepot OVRP
In the MDOVRP, the objective function only aims at computing the total traveling cost of the routes used because there are no setup costs for the opening depot, and the MDOVRP cases reported in the literature do not consider a cost associated with the use of vehicles. The depots are assumed to have unlimited capacity, so constraint (6) is also removed. Finally, because there are no setup depots, the variable y i is not necessary. Thus, the resulting MDOVRP model is as follows: subject to ∶ (2)-(5), (7), (8), (11), (12).
From the previous model, the number of binary variables Bn can be calculated as Bn = (m + n) 2 + m, and the number of rows Rn can be calculated as Rn = m + n + 2 ⋅ (m + n) 2 . Lalla-Ruiz et al. [23] present a comparison of the model presented by Liu et al. [27] in terms of these two characteristics. In this work, we also carried out the same comparison in respect to the previous two works, as shown in Figure 3A,B. As can be seen, our mathematical model requires a lower number of binary variables than the models proposed earlier. From the point of view of the number of rows, our model oscillates between the maximum and minimum limits given by Liu et al. [27] and Lalla-Ruiz et al. [23], respectively. It is interesting to note that the proposed model is more stable against variations in the number of depots. For example, instances p05, p06, and p07 they all have 100 clients but each has 2, 3, and 4 depots, respectively. Both the number of binary variables and the number of rows for these instances remain almost constant, as seen in Figure 3.

Distance-constrained MDOVRP
The constraints on the total distance traveled by vehicles can represent interesting features, such as the vehicle range and the distance taken per driver. In this subsection, a mathematical model that considers the minimization of the total traveling distance is presented, with the condition that the length of each path is between a minimum and a maximum value. The MDOVRP model (13) and (14) can be adjusted to consider these distance constraints as described below.
The total distance of an open path may be calculated as the distance accumulated in the last node of the path. Thus, two new set of variables must be defined. The first is the binary variable e j , which takes the value of 1 when the radial path ends at node j, for which the inequality (5) is rewritten as the following equality: The second variable is a i ≥ 0 ∀ i ∈ V, which is defined in order to calculate the distance accumulated from the beginning of the path to the node i ∈ J. Thus, the distance balance constraint (16) ensures that for every arc (i, j), which is part of the optimal solution, the difference between the accumulated distances a j and a i is equal to the distance T ij .
Thus, a path ending at node i implies that: e i = 1, the total length of the path is a i and, therefore, the constraint (17) limits the total length of each path solution between a minimum value L and a maximum value H. The length of the longest and shortest path of the optimal solution can be defined, respectively, as: Finally, constraint (18) defines the depots as reference nodes in which the accumulated distance traveled is zero for all the routes. Because the absolute value function is nonlinear, the following equivalent linear formulation is used to replace the expression (16): The parameter M is a sufficiently large constant, whose value was calculated as the total distance of a giant tour, covering all vertices, using arcs equivalents to the greater arc of the problem, that is, The final MIP model of the DCMDOVRP is: min. (13) subject to (2)

COMPUTATIONAL RESULTS
Four datasets of instances are used in order to show the operation and effectiveness of the proposed formulation. The first three sets of data correspond to OLRP cases, where the first dataset contains 30 instances adopted from the CLRP cases designed by Prins et al. [34], the second dataset consists of 36 instances for the OLRP, adopted from Tuzun and Burke's instances designed for the CLRP [43] and the third dataset consists of 20 instances adopted from Barreto's CLRP benchmark instances [6]. From the viewpoint of exact algorithms, the results for the OLRP are compared with those published by Yu and Lin [45]. For comparative purposes, the best-known solutions (BKS) obtained by the metaheuristic TBSA are presented, as cited in Schneider and Löffler [38]. The fourth dataset contains 24 instances used by Liu et al. [27] and Lalla-Ruiz et al. [23] for the MDOVRP and introduced by Cordeau et al. [15]. This dataset is adapted and used to test the efficiency of the model for the MDOVRP considering routes with distance limits.

Results on OLRP cases
In Table 1 we present comparative results against state-of-art exact methods and metaheuristics algorithms for the OLRP. For each mathematical model, the table gives the best known upper bound value, the lower bound obtained by LP relaxation and the relative gap (in %) between the lower and upper bounds. The value of the solution obtained by all the algorithms was computed by using the Euclidian distances multiplied by 100 and rounding the final distance value to the next highest integer. The earlier results are presented as cited in Yu and Lin [45]. The results reported previously and those presented herein consider a time limit of 14 400 seconds, are run on computers with similar characteristics (3.4 GHz Core i7 and 4GB of RAM) and use the CPLEX 12.0 (processing on one thread only).
For the proposed exact model, the bold numbers under the column CPLEX show the instances that have achieved the heuristic solution BKS a . In column LB, the lower bound (LB) obtained by LP relaxation is presented for each instance and the improvements with respect to the previous method are highlighted in bold. In the column Gap (%) c , the bold numbers indicate the best gaps reported using a mathematical solution. The column Gap (%) d corresponds to the gap between the UB of the proposed model and the heuristic solution BKS a . The bold numbers under column Gap (%) e show the best known solution obtained using MIP, wherein further negative values indicate the instances that have improved the BKS b . Finally, the last two columns report the LB and the gap obtained in the root node. In these two columns, the bold values indicate the solutions that improve the final LB b and final Gap b obtained in the earlier method, respectively.
As can be seen in Table 1, the proposed model produces high-quality results, obtaining the best exact solutions in all the test problems. The performance of the proposed formulation can be compared only in cases of 100 customers or fewer because for  In all open problems, except two (100-5-1) and (100-5-2), the lower bound LB was stronger than the previous method. In six instances, the LB root and the Gap root were significantly better than the final values of LB b and Gap b , respectively. Table 2 presents the results for the OLRP instances from Tuzun and Burke's dataset [43], using the same configuration of Table 1. For seven instances with solution previously reported, the Gap c and the final LB obtained by the proposed model are of better quality. In six instances, the Gap root turns out to be stronger than the final gap reported by Yu and Lin and in four instances, the LB root is of better quality than the final LB b reported. Two new BKS ware found, but in no instance has it been possible to achieve optimality from this dataset, either in this or other studies, so it would be very interesting to study the performance of different exact techniques in this type of instances. The abbreviations used in Table 2 are the same as those defined in Table 1.
In Table 3 we present comparative results for 20 OLRP instances adopted from Barreto et al. [6]. Results from Yu and Lin achieved optimality in 10 of them, while with our model was achieved optimality in 17 of them. Compared to the BKS obtained by the metaheuristic TBSA proposed by Schneider and Löffler [38], it should be noted that our results cover all 20 instances while they cover only 13 of them. Thus, in the second column, the results of metaheuristic algorithms are presented, corresponding to the mixture of the BKS obtained with the metaheuristic TBSA proposed in [38] and with SA-based heuristic method presented in [45].
Optimality was proven for the first time for nine instances, three new BKS were found, new lower bounds were obtained for all open problems, and five instances reached optimality in the root node. Detailed routes of the new optimal/BKS, where Gap c is zero and Gap d is negative, are presented in Table A1 of the appendix, where each route r i begins at the vertex corresponding to the depot. In all the sets of data evaluated for the OLRP the proposed model found a better average LB and the performance in the root node was quite satisfactory. The abbreviations used in Table 3 are the same as those defined in Table 1. a Optimality proven for the first time.  Table 4 presents a comparison of the node of the search tree in which the solution of the instance with optimality proven was achieved for each methodology. In all cases, except one (20-5-1-b), the proposed model explored a smaller number of nodes than the one reported in the previous method.

Results on MDOVRP cases
Two types of results are presented. The first corresponds to the case without distance limits, which are compared to those presented by Liu et al. [28] and by Lalla-Ruiz et al. [23]. These results are presented here as cited in the original works. The second type of results considers routes with distance limits. Because, in our review, the results of the second type were not found, then the presented results are a reference for future comparisons. In Table 5 we present comparative results against state-of-art exact methods for the MDOVRP. The earlier results are presented as cited in Lalla-Ruiz et al. [23]. The results reported previously and those presented herein consider a time limit of 7200 seconds, are run on computers with similar characteristics (Intel Dual Core 3.5 GHz and 4 GB) and use the CPLEX 12.3 (processing on one thread only).
As can be seen from the computational results shown in Table 5, in the MDOVRP without distance limit, 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. The average gap is 1.95% in comparison to the gaps of 62.07% and 7.77% reported by the earlier methods. Five new optimal solutions are reported for instances p05, p06, p07, pr03, and pr08. Detailed routes of the new optimal/BKS, where a Gap c equal to zero was reached and the Gap d was negative, are presented in Table A2 of the appendix. In 14 of the 25 instances of the MDOVRP, an LB root of better quality than the final LB reported in the previous models was obtained. Thus, the average of the LB root of the proposed model is considerably better than the final average LB reported by Lalla-Ruiz et al.
When the distance limits are considered (Table 6), the results show that it was possible to achieve the optimal solution in the same instances as in the case without limits, but usually with larger computational time. The results for cases with distance limits were solved with CPLEX 12.5, on a computer Intel Core i7-4770 3.4 GHz and 4 GB of RAM.
The value of H for each instance was calculated from the length of the longest path, which is obtained in the model solution without distance limits. For example, for instance p05 the length of the longest path is 97.1912. From this value we define H p05 = 97. As the H value decreases, the problem becomes more restricted and the complexity of its solution increases, for example, in the case p05 a value of H < 60 does not allow us to achieve a feasible solution in 45 000 seconds. A special case is the instance p07, whose optimal solution is 614.09, but only reached gap = 0 in 139 838 seconds. Abbreviations: l max , length of the longest route of the optimal solution; l min , length of the shortest route of the optimal solution.
Note that a perfectly balanced solution may not be possible or may contain highly inefficient routes. Therefore, the aim is that the difference between the maximum and minimum route length should not exceed a given threshold. Thus, balanced walking distances are obtained because the differences are bound by constraint (17). Figure A1A-C, presented in the appendix, show the optimal solutions obtained for the p05 instance with and without distance limits.

CONCLUSIONS AND PERSPECTIVES
A new approach for modeling the OLRP by MIP is presented, where the constraints related to subtours are defined more efficiently than the proposals presented for the same problem in previous models. Three types of open problems were addressed: (i) the OLRP, (ii) the MDOVRP, and (iii) the DCMDOVRP. The computational results were obtained from tests performed on 134 instances from the literature. Therefore, the proposed models are well suited to this type of open problems and the results are of high quality and performance.
Both the OLRP and the MDOVRP obtained equal or better upper bounds and the final lower bounds prove stronger than those obtained by earlier methods. The solution of the DCMDOVRP allows us to obtain balanced walking distances, which is an important characteristic in the realistic operation of transport vehicles.
An interesting aspect for future research is the improvement of the branch-and-bound method on the open type models that we propose. The incorporation of basic relaxations used to compute lower bounds could help to improve the effectiveness of the flow based formulation. The problem can be relaxed until a K-shortest spanning arborescence problem (KSSA) (or antiarborescence). These combinatorial relaxations were never used within branch-and-bound algorithms, and the preliminary computational results show that its quality is generally poor on problems closely related to CVRP [42]. However, in problems with topological characteristics based on arborescences, antiarborescence, and trees, Toth and Vigo [41] successfully used a Lagrangian relaxation based on the solution of KSSA.