Approximate solution of the shortest path problem with resource constraints and applications to vehicle routing problems

1.   Introduction

The vehicle routing problem (VRP) is a class of problems that entails finding an optimal set of routes for a fleet of vehicles to serve a set of customers. The objective of the VRP is to minimize vehicle routes costs while originating and terminating at a depot. The VRP was first proposed by Dantzig and Ramser ^[1] before being extended with different variants: the VRP with time window (VRPTW) ^{[2,3,4,5,6,7]}, the capacitated VRP (CVRP) ^{[8,9,10,11,12]} and other VRPs ^{[13,14,15,16,17]}. In recent years, the VRPs have attracted much interest ^[18]. Exact solutions of the VRPs generally result in branch-and-price ^[13,19], branch-and-cut ^[20,21,22] and branch-cut-and-price algorithms ^[22,23]. Ben Ticha et al. ^[24,25,26] proposed well-performing branch-and-price algorithms for VRPs on multigraphs with multiple time-related attributes and a homogeneous vehicle fleet. In such approaches, the linear relaxation in each branch-and-bound node is by solved column generation, which has proved to be a powerful approach ^[9,27,28]. Column generation has been widely used to solve a variety of large mathematical programs such as vehicle routing and crew scheduling problems ^[3,29,30].

Many authors have suggested solving a shortest path problem with resource constraints (SPPRC) introduced by Desrochers ^[31] as a multi-dimensional generalization of the shortest path problem with time windows. Resolutions methods and applications of the SPPRC have been extensively discussed in the literature ^{[3,29,31,32,33]}. These strategies range from exact methods to heuristics and meta-heuristics.

Exact SPPRC resolution techniques usually use the dynamic programming which has a pseudopolynomial complexity. Desrochers and Soumis ^[32] proposed a label correcting reaching algorithm that extends the Ford-Bellman algorithm to take resource constraints into account. The algorithm has been shown to be successful for tight resource constraints. Feillet et al. ^[34] adapted the Desrochers algorithm to solve exactly the ESPPRC (SPPRC with elementary path) pricing problems in the context of the VRPTW. Since that, this algorithm has been the backbone of a number of algorithms based on column generation applied to several important problems such as vehicle routing and crew management ^[16,26]. Table 1 summarizes some reviewed literature that applied column-generation for VRPs as classified based on solution approaches used for solving the sub-problem (SPPRC or ESPPRC).

The main drawback of exact resolution techniques is the computational time which increases with the number of resources, making them suitable for small-scale problems only. To speed up the computations and obtain solutions in a reasonable time, some authors proposed approximate, heuristics and meta-heuristics techniques ^{[14,18,50,51,52,53,54,55]}. These algorithms yield near-optimal solutions and are more suitable for real-life larger-scale applications (e.g., thousands of customers from, dozens of depots with numerous vehicles subject to a variety of constraints). Nagih and Soumis ^[14] proposed a basic heuristic to quickly produce feasible solutions. But, this approach has been limited to acyclic graphs and its speed has not been proved. Some authors have used simulated annealing ^[56], genetic algorithms ^[52,53], non-dominated sorting genetic algorithms ^[54] and simplified swarm optimization ^[55]. Meta-heuristic methods allows to avoid getting trapped in local optima but this comes at a price of slower convergence.

In this paper, we propose a Lagrangian relaxation based algorithm to approximately solve the SPPRC for arbitrary, acyclic and cyclic graphs. In our approach, called "dominance on Lagrangian cost for the SPPRC" (DLC-SPPRC), a dominance is expressed on a subset of the resources while the rest of the resources are dualized in the objective. Optimized parameter update schemes are used to ensure better performances (descent steps, Lagrange multipliers, subgradients).

The rest of the paper is organized as follows. Section 2 describes formally the SPPRC. Section 3 presents our approximate solution method that speeds up the computations while keeping a good approximate solution. In Section 4, we present an application to column generation and show the results obtained for various VRP datasets. The paper ends with the final conclusions.

2.   SPPRC

2.1.   Notations and preliminaries

Consider a graph $G = (V, A)$, where $V = N \cup \{o, d\}$ is a set of nodes $N = \{v_1, \dots, v_n\}$ that would be visited from an origin $o = v_0$ to a destination $d = v_{n+1}$ and a set of arcs $A \subset V \times V$. We denote by $\mathcal{R}$ a set of resources of cardinality $|\mathcal{R}|$. An arc $(i, j) \in A$ costs $c_{ij}$ and consumes a quantity $t_{ij}^r \geq0$ of each resource $r \in \mathcal{R}$. We suppose that the resources satisfy the triangle inequality. If $(i, j) \in A$, then the node $i$ is called a predecessor of $j$ and $j$ is called a successor of $i$. A path is a finite sequence of nodes

To each node $v_{j}$ of $P$, we associate $C_{j}^P$ which is the sum of the costs of its composite arcs up to $v_{j}$. We denote by $T_{j}^{r, P}$ the amount of resource $r \in \mathcal{R}$ used to reach the node $v_{j}$. For the sake of simplification, we drop the index $P$ from these notations and write $C_{j}$ and $T_{j}^{r}$. The path $P$ is said to be feasible if it satisfies the resource constraints

A path $P$ can also be represented by $X = (x_{ij})_{(i, j) \in A} \in \{0, 1\}$, where $x_{ij} = 1$ if the nodes $v_i$ and $v_j$ belong to $P$, and $0$ otherwise.

2.2.   Problem formulation

The SPPRC seeks a feasible path from the origin to the destination with a minimal cost. Similar to ^[3], one can formulate the SPPRC as follows:

subject to

where $X = (x_{ij})_{i, j}$ represents a path solution. Eqs (2.2) and (2.3) define the flow constraints on the graph; Eq (2.4) encodes the compatibility requirements between flow and time variables, and Eq (2.5) is the time windows constraint. If the problem is feasible, then there exists an optimal solution ^[14].

2.3.   Algorithms

To solve the SPPRC, labeling algorithms have been efficiently applied by many authors, e.g., Desrochers and Soumis ^[32], Dabia et al. ^[13] and Sun et al. ^[57,58]. Labels are defined for each node to encode partial paths, along with their cost and resource consumption. Each node receives labels and then extend them toward every possible successor node. The algorithm iteratively treats the nodes until no new labels are created.

To limit the proliferation of labels and reduce the computational time, dominance rules are introduced. Dominance relation is a partial ordering applied to eliminate non optimal solutions and retain only non dominated labels. The number of retained labels and, thus, the computation time, increase with the number of resources.

In the next section, we propose an algorithm to determine approximate solutions. The size of the search space is reduced by eliminating some labels with a heuristic dominance procedure.

3.   Lagrange SPPRC label correcting algorithm

With each path $P_i$ arriving to a node $v_i$, we associate a label (or state) vector

The path associated to a label $L$ is denoted by $X^L$. The path $P$ is extended to each successor $v_j$ with

Let $P_i$ and $P_i'$ be two feasible paths arriving to a node $v_i$. We say that $P_i$ is dominated by $P_i'$ and we write $P_i \preceq P_i'$ if the label vector of $P_i$ is component-wise less than the one of $P_i'$. On the other hand, we say that $P_i \leq P_i'$ if the first component of $P_i$ is smaller than the one of $P_i'$ or by comparing recursively the following ones if the first components are equal. At least one of the inequalities should be strict. This defines a partial order relation.

3.1.   Lagrangian relaxation

To speed up the resolution of the SPPRC (2.1)–(2.6), we propose to relax a subset of resources $\mathcal{R}_1 \subset \mathcal{R}$ and apply dominance on the remaining resources $\mathcal{R}_2 = \mathcal{R} \backslash \mathcal{R}_1$ only. This is equivalent to consider Eqs (2.4) and (2.5) on $\mathcal{R}_2$ and replace Eq (2.1) by the equation

where the Lagrangian function $\mathcal{L}$ is defined by

Since the flow is supported by at most one incoming arc at each node $j$ from its predecessors $i$, the multipliers $\lambda_{i, j}^r = \lambda_{j}^r$ are independent of $i$. The vector of all multipliers $\lambda_j =$ ($\lambda_j^r$)$_{r \in \mathcal{R}_1}$ can be obtained by solving the following maximization problem

To solve this optimization sub-problem, one can use a projected sub-gradient algorithm which proved to be efficient in practice. The updating scheme is

where $g_j^k \in \partial \mathcal{L}(\lambda, X)$ is a subgradient vector, $\tau_j^k > 0$ is a step size, and the max operator is taken for each element. We choose

where $Z_{\rm{UB}}$ is the best-known-smallest (upper)-bound of the solution of the problem described by (2.1)–(2.6) and $\theta_k\in(0, 2]$ is a scalar. One iteration for the estimation of $\lambda$ proved to be sufficient (precision and speed) when integrated into column generation.

3.2.   Dominance

The dominance at each node $j$ corresponds to the determination of the Pareto optima of the multicriteria problem applied to a set of labels. With the adopted Lagrangian approach, the label vector becomes

When the obtained final solution is not feasible (due to relaxation), we need to solve again the relaxed SPPRC by dominating on the Lagrangian cost determined in the first phase and on the resources $\mathcal{R}_2$ and by checking the resource windows in the original space $\mathcal{R}$.

Assuming that all predecessors of the node $j\in N$ have been considered, the dominance at the node $j$ can be interpreted as the determination of the Pareto optima for the multicriteria problem of $|\mathcal{R}_2|+1$ functions:

3.3.   Approximate solution procedure

First, we introduce the following notations.

● $E_i$: List of labels at a node $v_i$.

● $F_i$: The list of paths associated to the labels of $E_i$.

● $S_i$: Set of the successors of the node $v_i$ indices.

● $Q$: List of unprocessed nodes indices.

● $Pareto(L, E_j)$: Add a new label $L$ to a set of non-dominated labels $E_j$ while keeping the non-dominance property.

● $Z_L$ is the cost associated with a label $L$ (first component of Eq (3.9)).

Our improved label-correcting procedure is detailed in Algorithms 1 and 2.

4.   Experiments

To solve the $VRPTW$, we use the Dantzig-Wolfe decomposition which defines $K$ independent sub-problems and a global master problem. Applying a column generation, we alternatively solve the master problem and the $K$ sub-problems. We addressed the VRPTW instances with a column generation approach by using the SPPRC as a sub-problem.

We used the Solomon data sets ^[59] and Homberger 200 customer instances. Each instance contains customers locations, resources (two for each costumer, i.e., $|\mathcal{R}| = 2$) and constraints. These data sets are classified in three categories: the r-instances (customers are located randomly), the c-instances (customers are located in clusters) and the rc-instances (mixed random and clustered structures). Furthermore, each family of instances is divided into two types. The first type (r1xx, c1xx, rc1xx) has a short scheduling horizon, i.e., small time windows, that allows only a few customers per route (approximately 5 to 10). The second type (r2xx, c2xx, rc2xx) has a long scheduling horizon permitting many customers (more than 30) to be serviced by the same vehicle. We applied our approximate DLC-SPPRC algorithm and the classical SPPRC ^[32].

We implemented our algorithm using Java programming language. For the simulation, we used a CPU Intel Core i9-9900KF (8 cores), 3.60 GHz, RAM 32 GB, running under Windows 10 (64 bit). For the SPPRC, we implemented the Desrochers and Soumis algorithm ^[32]. Linear programs for restricted master problems are solved with ILOG CPLEX 20.1.

Tables 2–5 report the iteration number ($N_i$), the lower bound ($L^b_i$), the computational time in seconds ($T_i$) and the number of generated columns ($C_i$), where $i = 1$ for the SPPRC and $i = 2$ for the DLC-SPPRC. A gap is computed as

Figure 1 shows the total number of generated columns.

For most of the cases, we obtained the optimal solutions in reasonable times. As in ^[34], the column generation process is initiated with an adaptation of the Clarke and Wright algorithm ^[60]. The sub-problem is stopped at each iteration when $500$ labels have been extended to the depot with a negative cost.

From Tables 2–5, we can draw the following conclusions. Our algorithm was faster in all cases ($228$). In $15$ instances ($7\%$), the SPPRC algorithm could not obtain a solution in reasonable time (more than $10$ hours for one iteration). Our algorithm was $10$ times faster in $37$ cases ($16\%$ of the cases) and $1.5$ times faster in $166$ case ($73\%$ of the cases). On average, we achieved a gain of $564\%$ of the computational time. For the precision, in $188$ cases ($82\%$), we obtained similar precisions (a gap less than $1\%$) whereas the difference was more than $1\%$ only in 25 cases ($11\%$) and greater than $5\%$ only in 2 cases ($0.88\%$).

5.   Conclusions

This paper proposes an approximate solution algorithm to solve the SPPRC, which relies on a Lagrangian relaxation. Our method applies to both acyclic and cyclic graphs. Dominance is applied on a subset of the resources only. Optimized parameters updates schemes are used to ensure fast convergence. When applied to the VRP, our approach offers a good compromise between precision and computational time, and thus, can be applied to large-scale practical problems. In the future, we plan to apply our algorithm to solve various problems such as branch-and-price, branch-and-cut and branch-and-price-and-cut. In addition, we will explore more advanced optimization algorithms (heuristics, meta-heuristics, etc.) that have been successfully applied in other domains, such as, online learning, scheduling, multi-objective optimization, transportation, medicine, and data classification.

Conflict of interest

The authors declare that there is no conflicts of interest.