Mathematical formulations and improvements for the multi-depot open vehicle routing problem

Abstract

The multi-depot open vehicle routing problem (MDOVRP) is a hard combinatorial optimization problem belonging to the vehicle routing problem family. It considers vehicles departing from different depots to deliver goods to customers without requiring to return to the depots once all customers have been served. The aim of this work is to provide a new two-index-based mathematical formulation for the MDOVRP as well as enhancements on alternative sub-tour elimination constraints. Through numerical experiments, using traditional benchmarks for the MDOVRP, we show that the proposed formulation outperforms existing ones. In addition, we provide insights with regards to problem instance defining parameters on the models performance.

1 Introduction

The multi-depot open vehicle routing problem (MDOVRP) considers a set of vehicles departing from different depots that are not required to return to the depot once they have delivered the goods to customers. This problem was firstly proposed by Tarantilis and Kiranoudis [14] within the context of a real-world application involving the distribution of fresh meat in a real application. They solved the problem by using a list-based threshold-accepting algorithm. Later, Liu et al. [8] proposed the first mixed integer programming (MIP) formulation for the MDOVRP. Due to the extensive time and memory requirements when applying this formulation through a general-purpose solver such as CPLEX, the authors additionally proposed a hybrid genetic algorithm (HGA). Lalla-Ruiz et al. [6] proposed an improved formulation based upon adapting and lifting the Miller–Tucker–Zemlin (MTZ) sub-tour elimination constraints as well as additional constraints.

Besides the two above-mentioned formulations, no additional formulation nor exact approach has been developed for this problem. In this sense, having efficient formulations permits advancing the knowledge of the problem as well as assessing more accurately the performance of approximate approaches. With regards to approximate approaches, Soto et al. [13] proposed multiple neighborhood search, tabu search and ejection chains for the MDOVRP. Yao et al. [15] proposed an ant colony optimization approach in the context of product seafood delivery. Lahyani et al. [5] developed a hybrid adaptive large neighborhood search for the MDOVRP. Recently, Brandao [1] proposed a memory-based iterated local search algorithm (MBILSA) that uses the information on the history of the search to solve the MDOVRP.

The goal of this paper is to provide a new two-index-based mathematical formulation for the MDOVRP as well as assess the contribution of alternative constraints for handling flow conservation and sub-tours. In this sense, the comparison with the state-of-the-art formulations reported in the related literature so far indicates that our approach outperforms previous ones in terms of both solution quality and computational time. Furthermore, for some problem instances, our proposed model is capable of solving them to optimality for the first time. Finally, we analyze the relationship between the main defining parameters for this problem (i.e., number of customers, number of depots, demand) and the models’ performance.

The remainder of this paper is organized as follows. The proposed mathematical formulation for the MDOVRP as well as additional enhancements are introduced in Sect. 2. The numerical experiments carried out in this work is reported and discussed in Sect. 3. We end with conclusions in Sect. 4.

2 Two-index mathematical formulation

The MDOVRP can be defined as follows. Let $G=(V,A)$ be a directed complete graph, where $V = N \cup D$ is the vertex or node set and $A=V \times V$ is the arc set. Vertex set $D = \{1,2,\dots ,m\}$ represents the set of uncapacitated depots, whereas vertex set $N = \{m+1,m+2,\dots ,m+n\}$ represents the customers to be served. A symmetric travelling cost, $c_{ij}>0$, is defined for each arc between each pair of vertices $(i,j), i,j \in V, i \ne j$. Without loss of generality, the traveling cost can represent, according to the application environment, the distance, time, fuel consumption, etc. between each pair of vertices. In this regard, as in [6, 8], $c_{ij}=0$, $\forall i \in N, j \in D$, that is, the travelling cost from each customer to each depot is set to zero. This way the flow constraints are considered but they do not have an impact on the objective function. Also, in this work we do not consider distance limitations (i.e., maximum route length per vehicle) as already done in the work proposing models for the MDOVRP [6] and [8] in order to provide comparable objective and time values for the same scenarios. Moreover, each depot $d \in D$ stores and supplies enough goods to serve all the customers. With this goal in mind, each depot has an unlimited fleet of identical vehicles with the same positive capacity, denoted as Q. Each customer $i\in N$ has a certain demand of goods, denoted as $q_i$, where $0 < q_i \le Q$. The MDOVRP pursues to determine the routes of minimum traveling cost satisfying the following conditions:

• 1. Each customer must be visited in exactly one route.

• 2. Each customer has to be fully served when visited.

• 3. Each vehicle departs from one of the available depots and finishes at the last customer it serves.

• 4. The fleet of vehicles is homogeneous.

• 5. The total demand of the customers on any route does not exceed the capacity of the vehicles, that is, Q.

The above-mentioned assumptions are those related to the well-known open vehicle routing problem (OVRP, [12]) with the exception of assumption 3 that only applies to cases with more than one depot.

In the following, we introduce a new mathematical formulation for solving the MDOVRP where instead of using three indices, two are used. The following are the decision variables:

$x_{ij}$:

1 if a vehicle travels from node $i \in V$ to node $j \in V$, 0 otherwise

$y_i$:

1 if the minimum travelling cost to visit a customer $i \in N$ is from a depot, 0 otherwise

$u_i \ge 0$:

Upper bound on the load already assigned to a vehicle upon leaving node $i \in V$

As can be observed, $x_{ij}$ has two indices, while in the formulations of [6, 8] three indices are used, as they include the departing depot as additional index: $x_{ij}^k$, $k \in D$.

The following parameters are used:

$\overline{q_{i}}$:

For each customer $i \in N$, maximum value of $q_j$, $\forall j \in N, i \ne j$, i.e., $\max (q_j)_{j \in N, j\ne i}$

$d_{i}$:

Minimum traveling cost from any depot $k \in D$ to customer $i \in N$, i.e., $\min (c_{ji})_{j \in D}$. Note that $d_k = 0, \forall k \in D$

$r_{i}$:

Minimum traveling cost from any customer $j \in N$ to customer $i \in N$, $i \ne j$, i.e., $\min (c_{ji})_{j \in N, j \ne i}$

The mathematical model proposed in this work is defined as follows:

$$\begin{aligned} (MDOVRP_{2i}) \ \text{ minimize } \sum _{i \in V} \sum _{j \in V} c_{ij}x_{ij} \end{aligned}$$

(1)

subject to:

$$\begin{aligned}&\sum _{i \in V, i \ne j} x_{ij}=1, \forall j \in N \end{aligned}$$

(2)

$$\begin{aligned}&\sum _{i \in V, i \ne j} x_{ij} = \sum _{i \in V, i \ne j} x_{ji} , \forall j \in V \end{aligned}$$

(3)

$$\begin{aligned}&u_i - u_j + Q x_{ij} + (Q - q_i - q_j) x_{ji} \le Q - q_j, \forall i,j \in N, i \ne j \end{aligned}$$

(4)

$$\begin{aligned}&u_i \ge q_i + \sum _{j \in N, i \ne j} q_j x_{ji} , \forall i \in N \end{aligned}$$

(5)

$$\begin{aligned}&u_i \le Q - (Q-q_i-\overline{q_i}) \sum _{k \in D} x_{ki} - \sum _{j \in N, i \ne j} q_j x_{ij} , \forall i \in N \end{aligned}$$

(6)

$$\begin{aligned}&u_k = 0, \forall k \in D \end{aligned}$$

(7)

$$\begin{aligned}&d_i + M y_i \ge r_i, \forall i \in N \end{aligned}$$

(8)

$$\begin{aligned}&x_{k'i} \ge y_i, \forall i \in N, k'=argmin(c_{ki})_{k \in D} \end{aligned}$$

(9)

$$\begin{aligned}&\sum _{i \in N} \sum _{k \in D} x_{ki} \ge \left\lceil {\frac{\sum _{i \in N} q_i}{Q}}\right\rceil \end{aligned}$$

(10)

$$\begin{aligned}&x_{ij} \in \{0,1\}, \forall i, j \in V \end{aligned}$$

(11)

$$\begin{aligned}&y_i \in \{0,1\}, \forall i \in N \end{aligned}$$

(12)

$$\begin{aligned}&u_i \ge 0, \forall i \in V \end{aligned}$$

(13)

The objective function (1) aims at minimizing the total traveling cost of the routes used to deliver the goods to the customers. Constraints (2) guarantee that each customer is visited exactly once. Constraints (3) impose the degree balance of each node, including both customers and depots. Constraints (4)–(7) eliminate the sub-tours in the solutions found by considering the capacity of the vehicles and ensuring that their capacity is not exceeded. These constraints were those studied in [6] and proposed considering the improvements and extensions of the Miller–Tucker–Zemlin subtour elimination constraints (see Miller et al. [9]) to the VRP proposed by Desrochers and Laporte [3] and Kara et al. [4]. Furthermore, as investigated in [6], constraints (8) and (9) establish that if the minimum travelling cost to visit a customer is achieved directly from a depot, then a vehicle will visit the customer directly by departing from that depot. The minimum travelling cost to visit a customer $i \in N$ is determined by comparing all distances to that customer from any other node $j \in V, j\ne i$. Constraint (10) sets a lower bound over the number of vehicles required to serve the customers. Finally, constraints (11)–(13) are the integrality and non-negativity constraints for the different kinds of variables. It is worth mentioning that constraints (2)–(7) are necessary to model the MDOVRP, while constraints (8)–(10) are additional ones devoted to strengthen the formulation.

Note that one might but need not to add a redundant constraint (14) which would explicitly rule out arcs connecting pairs of depots from the solutions.

$$\begin{aligned} \sum _{i \in D} \sum _{j \in D} x_{ij} = 0 \end{aligned}$$

(14)

Proposition 1

Any solution from the MDOVRP$_{2i}$ formulation can be transformed into a solution for the model of Lalla-Ruiz et al. [6] (i.e., MDOVRP$_{3i}$) and vice-versa.

Proof

For transforming a solution from MDOVRP$_{3i}$ to MDOVRP$_{2i}$, the following expression has to be applied $x_{ij}=\sum _k x_{ij}^k$ , $\forall i,j \in V$. On the other hand, given a solution from MDOVRP$_{2i}$, variables $x_{ij}$ can be transformed into $x_{ij}^k$ variables by using the following expressions:

$$\begin{aligned} x_{kj}= & {} x_{kj}^k, \forall k \in D, \forall j \in V \end{aligned}$$

(15)

$$\begin{aligned} x_{ij}^k\le & {} \sum _{l \in V, l \ne i} x_{li}^k, \forall i \in N, \forall j \in V, \forall k \in D, i \ne j \end{aligned}$$

(16)

$$\begin{aligned} x_{ij} - \sum _{k \in D} x_{ij}^k= & {} 0, \forall i \in N, \forall j \in V, i \ne j \end{aligned}$$

(17)

Since all solutions from MDOVRP$_{2i}$ can be transformed into solutions for MDOVRP$_{3i}$ and vice-versa, thus the considered integer formulations are equivalent. $\square$

2.1 Alternative flow-constraints

In the formulations presented in [6] and [8], the flow constraints consider the return to depots, this consideration is valid for the MDOVRP bearing that both formulations consider $c_{ki}=0\ ,\forall k \in D$ and $\forall i \in V$, is set. Nevertheless, in the context of open routes, this set of constraints can be refined to avoid the return to depots. Thus, constraints (3) can be replaced by the following new set of constraints:

$$\begin{aligned}&\sum _{i \in V, i \ne j} x_{ij} - \sum _{i \in N, i \ne j} x_{ji} \ge 0 , \forall j \in N \end{aligned}$$

(18)

$$\begin{aligned}&x_{ij} + x_{ji} \le 1, \forall i,j \in V, i \ne j \end{aligned}$$

(19)

$$\begin{aligned}&\sum _{k \in D} \sum _{j \in V} x_{jk}=0 \end{aligned}$$

(20)

Constraints (18) establish that when a customer is visited, either the route can finish there or continue to another customer but no return to a depot is now contemplated. Constraints (19) forbid symmetries. Finally, constraints (20) set variables x to 0 when returning to depots. This formulation considering the alternative flow constraints is termed as MDOVRP$_{2i-f}$.

Proposition 2

Any solution of MDOVRP$_{3i}$ exists in MDOVRP$_{2i-f}$.

Proof

To proof this, consider a feasible solution for a given problem instance composed of (x, y, u) for MDOVRP$_{3i}$. Since $x_{ik}, \forall i \in V, k \in D$ are not considered in MDOVRP$_{2i-f}$, and bearing in mind that $c_{ik}=0, \forall i \in V, k \in D$ in MDOVRP$_{3i}$, then the partial solution (i.e., without considering $x_{ik}$) of MDOVRP$_{3i}$ exists in MDOVRP$_{2i-f}$ with the same objective function value. $\square$

Proposition 3

Any solution of MDOVRP$_{2i-f}$ exists in MDOVRP$_{3i}$ once that solution is extended to incorporate the corresponding $x_{ik}$ values satisfying constraints (3).

Proof

Since $x_{ik}, \forall i \in V, k \in D$ are not used in MDOVRP$_{2i-f}$, their corresponding values have to be set. Nevertheless, as $c_{ik}=0, \forall i \in V, k \in D$ in $MDOVRP_{3i}$, any value set to $x_{ik}$ from the MDOVRP$_{2i-f}$ leads to having a solution in MDOVRP$_{3i}$ with the same objective function value as in $MDOVRP_{2i-f}$ . $\square$

Proposition 4

Any solution of MDOVRP$_{2i}$ exists in MDOVRP$_{2i-f}$ and vice-versa.

Proof

Considering Propositions 3 and  4 with the procedure discussed in Proposition 1. $\square$

Proposition 5

MDOVRP$_{2i}$, MDOVRP$_{2i-f}$, and MDOVRP$_{3i}$ are equivalent.

Proof

Consider Propositions 3, 4, and 5. $\square$

2.2 Incorporating arc-based loading constraints to MDOVRP$_{2i-f}$

As can be checked in the MDOVRP$_{2i}$, the MTZ elimination constraints [i.e., constraints (4)–(7)] are adapted to hold capacities as opposed to [8]. This allows to ensure the continuity of vehicles in terms of demand. Nevertheless, similarly as done in [10, 11], remaining load on the vehicles among arcs can be explicitly be considered and, hence, used to avoid sub-tours.

$$\begin{aligned}&\sum _{i \in V, i \ne j} u_{ij} - \sum _{i \in V, i \ne j} u_{ji} \ge q_j, \forall j \in N \end{aligned}$$

(21)

$$\begin{aligned}&(Q - q_i) \cdot x_{ij} \ge u_{ij} , \forall i,j \in N \end{aligned}$$

(22)

$$\begin{aligned}&Q \cdot x_{kj} \ge u_{kj} , \forall k \in D, j \in N \end{aligned}$$

(23)

Constraints (21) set a lower-bound on the remaining quantity of the arcs. That is, if a customer j is visited, the current load of the vehicle has to be larger than the requested quantity as well as the next customer requested quantity in the route. Constraints (22–23) set an upper bound on variables u. The formulation MDOVRP$_{2i-f}$ replacing (4)–(7) by (21)–(23) is termed as MDOVRP$_{2i-flv}$.

Proposition 6

The LP relaxation of MDOVRP$_{2i-flv}$ is tighter than that of MDOVRP$_{2i}$.

Proof

Since constraints (22)–(23) and (5)–(7) establish the bounds on u variables in the corresponding formulations, it is sufficient to show that constraints (21) are stronger than those relating u for two different customers, namely constraints (4). In doing so, consider constraints (4) in the LP relaxation of MDOVRP$_{2i}$, as it includes variables x, it leads to the following cases:

• when $x_{ij} > 0$ and $x_{ji} =0$; $u_j \ge u_i+q_j- Q +Q x_{ij}$. This can be rewritten as $u_j \ge u_i+q_j-\varDelta Q$.

• Similarly, when $x_{ji} > 0$ and $x_{ij} =0$ ; $u_j \ge u_i - Q + q_j + (Q -q_i - q_j) \varLambda x_{ij}$, that is, $u_j \ge u_i -\varLambda Q + \varLambda q_j - \varLambda q_i$.

• Finally, $x_{ji},x_{ij}> 0$; $u_j \ge u_i - Q + q_j - \varDelta Q + \varLambda Q -\varLambda q_i - \varLambda q_j$.

Where $\varDelta$ and $\varLambda$ represent the fractional value due to the multiplication with $x_{ij}$ and $x_{ji}$, respectively. Notice that in the first two cases, related constants are considered.

With regards to constraints (21), it is necessary to convert $u_{ij}$ into $u_i$ as considered in MDOVRP$_{2i}$. In doing so, the following transformation has to be conducted for each customer $j \in N$: $u_j=Q- \sum _{l \in V} u_{lj} + q_j$. That can be rewritten as $\sum _{l \in V} u_{lj}=Q- u_j + q_j$. The previous permits to transform constraints (21) in: $Q - u_j \ge \sum _{l \in V} u_{jl}$ $\implies$ $Q - u_j \ge Q - u_l + q_l$ $\implies$ $u_l \ge u_j + q_l$. Where for a customer j, l represents the customer visited after j in that route.

As can be checked, the previous expression is tighter than those obtained for MDOVRP$_{2i}$ as x variables are not incorporated. This completes the proof. $\square$

Proposition 7

The LP relaxation of MDOVRP$_{2i-flv}$ is tighter than MDOVRP$_{3i}$.

Proof

Similar to that for Proposition 6. $\square$

3 Computational results

This section presents the computational experiments carried out to assess the performance of the proposed mathematical model as well as improvements proposed in this work. With this goal in mind, we compare the performance of our formulations with that of the mathematical models proposed by Liu et al. [8] and Lalla-Ruiz et al. [6]. All computational experiments were conducted on a computer equipped with an Intel Core i7 3.7 GHz and 32 GB of RAM and all models, including the ones from the literature, were implemented in a general-purpose solver, CPLEX version 12.9, with a maximum computational time of 2 hours (7,200 seconds) and limited to one-thread.

The problem instances used in this work were initially proposed by Cordeau et al. [2] for the MDVRP. Liu et al. [8] used them for assessing their MDOVRP formulation and algorithm. These and subsequent works do not consider the distance limitation. To provide a fair comparison, i.e., using comparable values and times, as well as an equal comparison with approximate approaches proposed for this problem, we do the same. Moreover, the main parameters ruling the instances are the number of customers |N|, number of depots |D|, and capacity of the vehicles Q. This as well as additional information is provided in detail in Sect. 3.2, in Table 3. This section also provides insights on how the instance parameters influence performance.

3.1 Optimization models’ performance

The computational results obtained by the linear relaxation of the models proposed by Liu et al. [8], Lalla-Ruiz et al. [6], and the ones proposed in this paper are depicted in Table 1. The first column shows the identifier of each problem instance to solve, Inst. For each formulation, the column Value provides the optimal value and t(s.) the computation time in seconds.

Table 1 Computational results obtained from the relaxation of the models

As can be seen, for all instances, the linear relaxation is improved with the exception of p12 that provides the same value also in the 3 index formulation. As later discussed in Figs. 1 and 2, it can be inferred that a larger number of constraints is beneficial at the expense of higher computational times. Nevertheless, comparing the case without additional improvements MDOVRP$_{2i}$ with MDOVRP$_{3i}$, it can be observed that the quality of the relaxation remains the same but the average computational time is reduced by 78%.

With regards to the results obtained by the non-relaxed optimization models they are reported in Table 2. The layout of this table is similar to the previous one but in this case, we also provide columns UB and LB for reporting the upper and lower bounds reported by CPLEX within the time limit set by default, respectively. Column Gap (%) reports the relative error of the bounds in percentage.

Table 2 Computational results obtained by the mathematical models proposed by Liu et al. [8] and that proposed in this paper for the instances introduced by Cordeau et al. [2]. Best values are reported in bold

From the computational results reported in Table 1, having a two-index formulation greatly accelerates the solving time when compared to the three-index formulation. Moreover, MDOVRP$_{2i-flv}$ provides better optimal values, indicating thus a better relaxation.

With regards to the non-relaxed models reported in Table 2, it can be noted that our proposed mathematical model, either with or without enhancements, provides better performance in terms of narrower linear bounds and computational times than MDOVRP$_{Liu}$ and MDOVRP$_{3i}$. The proposed models improve or at least equal the quality of the linear bounds obtained by the model MDOVRP$_{3i}$ in the majority of the cases. In terms of gap, the model that considers the enhancements reports a gap of only 1.84% on average in comparison with the gap of 3.40% on average reported by the model without enhancements or 6.09% obtained using MDOVRP$_{3i}$. This reduction of gap finds its rationale behind the improved lower bounds provided by our model. In this regard, it is worth highlighting that the lower bounds provided by MDOVRP$_{2i-flv}$ are equal or better than the best ones of the rest of formulations.

Furthermore, our mathematical model solves to optimality for the first time the problem instances p05, p07, p08, p09. p10, p11, pr08, and pr09. It is noticeable that in the cases where state-of-the-art mathematical models provide the optimal solution, our proposed models require remarkably less computational time. For example, for the problem instance p01, the models MDOVRP$_{Liu}$ and MDOVRP$_{3i}$ require 696.78 and 47.26, respectively while the new models require about 1 second. Another example from the last reported model is pr03, the state-of-the-art model requires around 5998 seconds to solve to optimality, whereas our proposed models require 26.91 and 19.84 seconds. Finally, in those cases where the optimal solution has not yet been achieved, both proposed models enhance the quality of the bounds for all cases.

Figures 1 and 2 show the structure of the mathematical models in terms of number of rows and binaries generated by the solver when dealing with the problem instances under analysis. As can be observed, MDOVRP$_{2i-flv}$ requires a larger amount of number of rows than MDOVRP$_{2i}$ and MDOVRP$_{3i}$ which is mainly caused by the additional flow and arc-load based constraints. Nevertheless, it can be seen that both, MDOVRP$_{2i}$ and MDOVRP$_{2i-flv}$, greatly reduce the number of binaries with regards to state-of-the-art formulations. Considering this information along with Table 1, we see that the more constraints provide a tighter LP relaxation. It is worth noting that this is just a comparison by means of the number of constraints, however, some constraints are different and restrict the feasible region in a different way.

Fig. 1

Number of rows processed by CPLEX for solving each instance

Fig. 2

Number of binaries processed by CPLEX for solving each instance

The comparison of the proposed model with and without enhancements leads to the following conclusions:

• Including additional constraints as in MDOVRP$_{2i-flv}$, instead of the Miller–Tucker–Zemlin subtour elimination constraints (4)–(6), that keep track of the remaining space along arcs, improves the quality of the solutions in comparison to MDOVRP$_{2i}$ and MDOVRP$_{3i}$. Furthermore, as can be observed, the average computational times are reduced while the quality of bounds is improved. Also, it overall leads to a large number of best solutions found.

• The model with and without enhancements reports better performance in terms of computational time and linear bounds. In this regard, the gap difference with respect to those already reported in the literature is substantial. For instance, the gap of the hardest instances so far, i.e., $p08-p11$, has been reduced to less than half in MDOVRP$_{2i-flv}$. On other cases, the enhancement in terms of time performance is noticeable, e.g., pr03, where the time to provide the optimal solution is reduced from 5998.46 to less than 27 seconds with MDOVRP$_{2i}$, and 17 with MDOVRP$_{2i-flv}$.

• For some instances, the inclusion of the enhancements accelerates the convergence to optimal solutions. For example, in p05, MDOVRP$_{2i-flv}$ requires 72.96 while MDOVRP$_{2i}$ requires 6891.12. Considering the information provided in Figs. 1 and 2, we see the positive influence of adding the proposed additional constraints without increasing the number of variables.

3.2 Analysis considering problem instances’ structure

In this section we provide insights into the effects of the instance parameter settings on the performance of the models. In Table 3, the following scenario parameters are reported: number of customers |N|, number of depots |D|, and capacity of the vehicles Q. Also, we report the ratio between the number of customers and depots (customer-depots=|N|/|D|), total amount of customers’ demand ($\sum q_i$), and the minimum number of vehicles needed to satisfy that demand (minimum vehicles=$\sum q_i/Q$). Finally, we also report the gap performance of all MDOVRP optimization models.

As can be concluded from the table, in those cases where the ratio customer-depots |N|/|D| is relatively low, i.e., less than or equal to 15, all formulations provide the optimal solution. When this ratio increases from 15 to 40, we see that not in all cases the optimal solutions are provided. In some cases, when they have the same ratio, the minimum number of vehicles influences, like for example, p12 and p15, where the three formulations are solved to optimality. Those scenarios have a ratio of 40. However, when the minimum number of vehicles $\sum q_i/Q$ increases like in p18, the problem becomes harder to solve for MDOVRP$_{Liu}$.

When considering only the best state-of-the-art performing model (i.e., MDOVRP$_{3i}$) and our current proposals, we observe that the hardest instances to solve in terms of gap are p08, p09, p10, and p11. Those instance have in common that they have the highest $\sum q_i/Q$. Also, they have the highest ratio customer-depots |N|/|D| (except for pr06 that has a higher value than p10 and p11). Thus, considering the previous and current observation, we see an influence of those parameters on the performance of the models.

Finally, it should be noted that these observations only take into account structural instance parameters but the model constraints, variables, and specific instance information play a relevant role on the models performance. Therefore, a more in-depth study in this research direction might be valuable.

Table 3 Problem instances’ structure and optimization models gap performance

4 Conclusions and further research

In this work, we have presented a novel formulation and enhancements for the well-known Multi-Depot Open Vehicle Routing Problem. Both contributions provide better results in terms of computational time and number of best solutions provided, and provide tighter linear bounds than those mathematical formulations previously reported in the literature. The incorporation of constraints to track vehicles remaining capacity as well as symmetry constraints to consider empty routes remarkably improves the convergence speed over those problem instances where the optimal solutions are reached. Furthermore, the numerical results indicate that our mathematical model reduces the computational burden while achieving better solutions for those instances where the optimal solutions are not yet known. In this regard, the new formulation with and without enhancements provides new best values for all those problem instances and provides optimal solutions for some instances that had not been solved to optimality so far. Also, the analysis of the gaps provided by the models and the instance parameters showed that the relation between customers, depots, and the minimum number of required vehicles have a relevant influence for the studied instances.

Even though our proposed formulation and improvements have been tested only for the MDOVRP, we believe they can contribute and improve other variants of the VRP and multi-depot VRP (e.g., [2, 7]) in future research. Finally, a more in-depth study on the influence of instance parameters on the solution approaches performance will be a topic for further research.