We start by reviewing the green vehicle and mixed fleet routing problems. Conrad and Figliozzi (2011) conducted a bounding study for an EVRP with time windows and fixed recharging time. Fukasawa et al. (2016) proposed a branch-and-cut-and-price algorithm for an EVRP with time windows. In particular, Fukasawa et al. (2016) studied the strengths of different reformulations and developed the 2-cycle-free q-routes column generation framework. Schneider et al. (2014) considered an EVRP with time windows and truck capacities. Like Schneider et al. (2014), we also allow EVs to visit any available rechargers before running out of electricity. We assume that after each visit the battery is fully charged. Schneider et al. (2014) provided a non-linear mixed integer formulation and a hybrid meta-heuristic algorithm. We seek to show that a linear formulation can capture all the elements permitting greater computational efficiency.

It should be noted some authors consider richer options for recharging, while only focusing on a single type of vehicle. Desaulniers et al. (2016) studied the EVRP under four different recharging policies. The authors tailored their BCP algorithms under different scenarios and ran their models on several benchmark problems. Andelmin and Bartolini (2017) considered an alternative fuel charge vehicle with delayed recharging policy. They introduced tighter cycle elimination techniques which allow exact solutions for more than 100 customers in the previously unsolvable instances.

Although many EVRP works focus on constant or linearized electricity consumption, research that involves accurate energy estimations is of importance in providing guidance for the optimization models. For example, Arasu et al. (2019); Demir et al. (2012) both provided parameter estimation models, e.g. for powertrain recharging rates, nonlinear speed and electric consumption, and heuristic algorithms for solving energy-oriented VRPs. As mentioned before, the only two existing publications on similar mixed vehicle routing variants (Goeke et al. (2015); Macrina et al. (2019)) are focused on heuristic approaches. Apparently, none of the previous research provides a clear optimization formulation that can be solved by off-the-shelf solvers.

In terms of the VRP variants with common carriers, most of the research still relies on the field of meta-heuristics. The first work to address outsourcing options for VRPs is Chu C-W (2005). The author considered a single-depot routing problem with a less-than-truckload carrier option. The objective is to minimize the total network cost. Chu C-W (2005) solved the problem by a simple “modified savings” math-heuristic for up to 30 customers. The first work for outsourcing options with a time window is Moon et al. (2012). The authors assign deliveries to vehicles of different carriers in a greedy manner. A heuristic based on a genetic algorithm is proposed in Vidal et al. (2016). Gahm et al. (2017) introduced a new variant in which a heterogeneous dedicated fleet and multiple cost options are considered. This work involved only variable travel costs and three options on the common carrier costs. The problem is solved by a neighborhood search heuristic. More recently, Alcaraz et al. (2019) combined heterogeneous fleet types with last-mile outsourcing options and hours of service regulation (layovers). The authors employed three heuristics: simulated annealing, tabu search, and variable neighborhood search - in solving their problem. They examined the results on up to 100 shipments random instances and studied the stability of the three heuristics. Dang et al. (2020) considered the layover regulations and the outsourcing options in their work and solved the new problem by an ant colony based meta-heuristic approach. The resulting problem is called vehicle routing problem with common carriers and time regulations (VRPCCTR). A recent publication by Dabia et al. (2019,a) is the only one that focuses on exact methods. This work considered complex tariff structures for the outsourcing options. The authors provided two different set-partitioning reformulations, strengthened cuts and modified dominance criteria for the BCP algorithm. Moreover, they conducted relatively structural studies on the impacts of the outsourcing decisions to the routing decisions. Since VRP with common carriers is not thoroughly studied yet, their work yields many opportunities for further study and analysis.

The publications regarding hours-of-service regulations contain more diversified solution methodologies. Alcaraz et al. (2019) reviewed works that considered both layovers and outsourcing options. From their review, this hybrid setting is still mainly solved by meta-heuristic approaches. Goel and Irnich (2014) presented the first exact method for the VRP variants considering layover regulations. The authors proposed a branch-and-price algorithm with a forward labeling search procedure. New dominance rules involving the layover feature were introduced. Their numerical results on Solomon instances suggested the effectiveness of the algorithm. Later, Tilk (2016) improved the work of Goel and Irnich (2014) by adding cutting planes to the master problem and a bidirectional labeling search. The acceleration strategies discussed in this work are also useful in speeding up the pricing procedures.

We solve the MVRP to optimality using BCP algorithms. With BCP, an extensive formulation is linearly relaxed. The relaxation works as a master problem which is solved using column generation. Additional valid inequalities are added to strengthen the lower bounds. Integer solutions are finally reached by branch and bound, see Lübbecke and Desrosiers (2005). The column generation process iterates between the linear reoptimization of the master problem and the solution of the pricing subproblem. In general, column generation is one of the most successful large-scale exact methods in solving combinatorial optimization problems. Desrochers et al. (1992) is the first to apply column generation in the context of the VRPTW, resulting in a branch-and-pricing algorithm. Later, Kohl et al. (1999) improved the branch-and-price algorithm by introducing subtour elimination constraints and two-path inequalities which resulted in the branch-and-cut-and-price. An additional line of research focuses on valid inequalities for column generation algorithms based on the variables of the master problem. For example, Jepsen et al. (2008) proposed a set of valid inequalities called subset-row inequalities with the variables in the Dantzig-Wolfe master problem. SR inequalities have been observed to effectively improve the lower bounds. Petersen et al. (2008) discovered how to apply generalized Chavátal-Gomory cuts to speed up the column generation solutions in the context of VRPTW. More recently, Dabia et al. (2019,b) studied the strengths of robust and non-robust cover inequalities applied on BCP algorithms. Some column generation methods often show very slow convergence partly due to heavy degeneracy problems (Rousseau et al. (2007)). Such problems arise when multiple dual solutions are associated with each primal solution. Rousseau et al. (2007) showed that the interior point method helps in stabilizing the column generation algorithms for VRPTWs.

While the master problem for most VRP variants is standard, the pricing algorithms for the subproblems are differentiated. Back in the 1990’s and earlier, researchers treated subproblems as a shortest path problem with resource constraints. Exact dynamic programming and 2-cycle eliminations were developed to price the columns with negative reduced cost, see Desrochers et al. (1992). Irnich and Villeneuve (2006) improved the elimination capability from 2-cycles to k-cycles, which increases the lower bounds of the relaxation solutions. However, the procedure is NP-Hard. Nowadays researchers mostly regard the pricing subproblem as an elementary shortest path problem with resource constraints. This consensus is attributed to the contributions of Feillet et al. (2004); Righini and Salani (2006 2008). Their works are cornerstones for the so-called labeling algorithm. The authors also introduced various accelerating techniques such as bounding, state-relaxation, bidirectional labeling search, and decremental state space relaxation, etc. Tilk et al. (2017) further improved the bidirectional labeling algorithm by dynamically adjusting the midpoint so that the number of labels generated by each direction can be controlled to a low level. Another group of researchers are focused on developing heuristic pricing algorithms to speed up the solutions for subproblems. For example, Desrochers et al. (1992) proposed a tabu search heuristic to quickly price the elementary paths with negative reduced cost. By far the most widely used route relaxation method is ng-routes, which was introduced by Baldacci et al. (2011). The method was proved by many researchers to be effective in solving difficult VRPTW instances.

We state the assumptions and policies that support the modeling of MVRPC. All dedicated fleet routes (ICVs and EVs) must start and end at the same depot. The average speed of every truck type is assumed to be constant and the same. While we regard the ICV fleet as containing multiple truck types differentiated by heterogeneous sizes, we assume there is only one type of EVs because they are single-sized. The demands for each shipment are inseparable, which means they can either be consolidated on a route or be entirely outsourced to a common carrier. Each shipment that is not outsourced must be visited exactly once by exactly one route traveled by one type of truck. A vehicle must visit the nodes within its corresponding time window such that no waiting time or late delivery is permitted. Furthermore, we assume the time windows for the rechargers and the depot are relaxed, which means they can be visited anytime that they are needed. This setting gives decision makers the data needed to optimally schedule the limited rechargers. Whenever a recharger is visited, the battery is fully recharged and the time will be counted. Yet, such time will not be regarded as working time subjected to the hours of service regulation. To make our problem more general, we do not consider the subtle break period as regulated in Federal Motor Carrier Safety Administration (2011). Instead, we only request that drivers take a layover period between 10 and 14 hours once their uninterrupted working time reaches the 14-hour limit. We realize that some previous literature considered curb weight, road infrastructures, traffic conditions in computing power consumption on a given trip, Arasu et al. (2019); Goeke et al. (2015). In our work, however, we assume that the EV fleet has a constant unit electric consumption per mile. Then the consumption on the arc is linear to the distance.

The MVRPC is defined on a directed complete graph $G({\overline{I}}_0,A)$, where ${\overline{I}}_0=\{0,1,...,n,e_1,...,e_m\}$ denotes the set of nodes and $A=\{(i,j)\in {\overline{I}}_0\times {\overline{I}}_0:i\ne j\}$ is the set of arcs. Node 0 represents the depot at which all the vehicle routes start and end. The nodes in $I=\{1,...,n\}$ represent the customers, while $E=\{e_1,...,e_m\}$ is the set of recharging stations. The recharging stations may share the same locations because it is extremely expensive to establish and operate high-voltage rechargers in separate locations. The combination of customers and recharging stations are denoted as $\overline{I}$. Inclusion of the depot is indicated by subscripting the respective set, e.g., $I_0$, ${\overline{I}}_0$, $E_0$. Each node $i\in \overline{I}$ is associated with a demand ${\omega }_i$, a time window $[a_i,b_i]$, and a service time ${\lambda }_i$. Note that the demand and the service time are positive for $\forall i\in I$ and 0 for $\forall i\in E$. Each arc $(i,j)\in A$ is associated with a non-negative travel distance $d_{ij}$, and an estimated energy consumption ${\delta }_{ij}$ determined by models discussed in Arasu et al. (2019). We assume constant energy consumption in analyzing applications in the logistic industry in which vehicles usually carry full truckloads. When payload distribution is considered, the energy assumption term ${\delta }_{ij}$ becomes a function of total vehicle weight en arc. Due to operational requirements, a maximum-allowed intra-node distance $d^{max}$ is enforced to rule out infeasible arcs. A mixed dedicated fleet consisting of electric vehicles (EVs) and internal combustion vehicles (ICVs) is available at the depot, where ICVs are usually diesel-powered trucks. Let $V=V_D\cup \left\{e\right\}$ be the set of truck types, where $V_D$ denotes the set of heterogeneous diesel trucks and $\left\{e\right\}$ is the available homogeneous electric trucks. Each truck of type $k\in V$ has a limited capacity $Q_k$, a fixed cost $f_k$, a stop cost p for each delivery, and a finite number of trucks $n_k$ available. For all fleet types $k\in V$, we use the same averaged speed $g=55$ $m\mathit{ph}$. Furthermore, for all arcs $(i,j)\in A$ and fleet types $k\in V$, let $c_{ijk}$ be the travel cost from node i to node j, where $c_{ijk}$ is an average of drivers’ wages and fuel cost for ICVs and an average of driver’s wages and energy cost for electric trucks. Since there is no reason to operate an ICV route from customer sites to recharging stations, or vice versa, we define an infeasible set that represent these arcs. Let $A^{-}$ represent the infeasible arc set, where $A^{-}=\{(i,j,k),(j,i,k):\forall i\in {\overline{I}}_0,\forall j\in E,k\in V_D\}$.

On each visit to a recharger, electric trucks are recharged to their maximum battery capacity $\phi$. The recharging time depends on the recharging rate $\alpha$ and the remaining electrics on arrival at the recharger. Common carriers are available for which demand can be outsourced. If customer $i\in I$ is outsourced to a common carrier, the corresponding cost is retrieved from the tariff as $m_i$. According to U.S. governmental hours of service regulations, the truck drivers need to take a layover when their working time reaches the upper threshold $\xi =14$ hours, where the layover period is specified between ${\beta }^{lb}=10$ hours and ${\beta }^{ub}=14$ hours (according to industry conventions).

Thereafter, MVRPC can be formulated as a mixed integer linear program (MILP). For every arc $(i,j)\in A$ and $k\in V$, let $x_{ijk}$ be a binary variable that takes value 1 only if fleet type k travels arc $(i,j)$ and 0 otherwise with $x_{ijk}=0$ is fixed if $(i,j,k)\in A^{-}$. Let $y_i$ be a binary variable which takes the value 1 only if customer i is subcontracted to a common carrier and 0 otherwise. Continuous variables $s_{ik}$ define the arrival time at node i by truck type k and $u_i$ is the accumulated demands from the depot up to node i. Continuous variables $z_i$ trace the remaining battery level on arrival at node i. $t_i$ is defined as the remaining working time before a layover on arrival at node i. $r_{ik}$ is used to determine the hours taken by the driver for layovers immediately before node i and $l_{ik}$ is the number of layovers needed precedent i. To keep the notations consistent in the formulation, we define $l_{0k}=\{l_{10k},l_{20k},...,l_{e_{m}0k}\}$ as a $(n+m)\times 1$ vector for each $k\in V$, where $l_{i0k}$ represents the number of layovers needed on the way back to depot. Note that $r_{ik}=l_{ik}=0$ for all $k\in V_D$ and $i\in E$.

For clarity, we summarize the parameters and variables of the MILP model in Table 1.

This problem $\left(\mathcal{P}\right)$ can be represented by the following arc flow formulation:

The objective function (1) minimizes the total fixed costs, travel costs, and outsourcing costs. Constraint (2) ensures that a customer node that is not outsourced is visited exactly once by a type-k truck. Constraints (3) and (4) ensures that the number of trucks in use is less than the available trucks for each type k. Constraint (5) guarantee that if a recharger i is visited, it can be used at most once. Constraints (6) and (7) are the flow conservation constraint of the trucks. Constraints (8)-(10) guarantee the customer time windows are respected, with variable layovers and recharging time added if necessary. It is important to set M and $M^{\prime }$ as small as possible so that constraints (9) and (10) are tight enough but still valid. Here, we let $M=\frac{d_{ij}}{g}+{\lambda }_i+\alpha ·\phi +T·{\beta }^{ub}+b_i$, where T is a constant parameter that equals 1 for all $i,j\ne 0$ and equals an integer that refers to the largest possible layover numbers between depot and any node. Furthermore, we let $M^{\prime }=T·{\beta }^{ub}+b_i$. Constraints (11) and (12) specify the electric consumption on each arc, which follow the non-decreasing rule. Constraints (13) and (14) trace the moves of truck type k along each arc to add necessary layovers. Again, the layover period and work time follow the non-decreasing rule node by node. Next, constraint (15) ensures the dynamic layover time added to truck type k obeys the hours of service regulation. The subtour elimination constraints (16), (17) are extended from Dang et al. (2020) to the heterogeneous case. Constraint (18) ensures the intra-node distances are feasible. Finally, we define the domains for each decision variable in constraints (19) and (20).

To further tighten the LP bounds obtained from the SP formulation, we can add additional constraints satisfied by all feasible solutions into $\left(\mathcal{R}\right)$. One way to obtain additional constraints in the ($x_q,y_i$)-space is through the transformation of valid inequalities in the space of other variables. For examples, some well-known valid inequalities in the ($x_{ijk}$)-space derived for HVRP are also valid for MVRPC. The reason is that their three-index flow formulations share the same feasible set. However, those valid inequalities are no longer strong because of the existence of $y_i$. When the ($x_{ijk}$)-variables and the ($x_q$)-variables both refer to the same feasible route, the coupling constraints between them can be defined as in constraints (25). Then, any valid inequality in the ($x_{ijk}$)-space can be lifted and transformed to a valid inequality in the ($x_q,y_i$)-space.

Constraints (23) are the obvious projection of constraints (3) and (4) from the ($x_{ijk}$)-space to the ($x_q$)-space. Considering that constraints (3) and (4) hold and ${\sum }_{i\in I}{x}_{i0k}={\sum }_{q\in {\mathsf{\Omega }}_k}{x}_q{\sum }_{i\in I}{\eta }_{i0q}={\sum }_{q\in {\mathsf{\Omega }}_k}{x}_q$ holds, we thus have ${\sum }_{q\in {\mathsf{\Omega }}_k}{x}_q\le n_k$, which is exactly constraints 23.

In this section, we describe the pricing problem in detail. Let ${\mu }_i\in \mathbb{R}$ (${\mu }_i\ge 0$) for $i\in I$ be the dual variables associated with constraints (22) and ${\gamma }_k\in \mathbb{R}$ (${\gamma }_k\le 0$), $k\in V$ be the dual variables associated with constraints (23). Let ${\overline{c}}_q$, $q\in \mathsf{\Omega }$ be the reduced cost of a complete path q that is operated by a type k truck, i.e., a closed loop route that starts and ends at the depot. The pricing problem without strengthened cuts for the problem $\left(\mathcal{R}\right)$ can be formulated as follows.

We first write down the calculation of $c_q$ in (41) explicitly, $c_q=f_k+f_k·l_q+{\sum }_{(i,j)\in A}(c_{ijk}+\frac{1}{2}p)·{\eta }_{ijq}$. Here, we use $l_q$ to denote the number of layovers spent on path q. A feasible route in $\mathsf{\Omega }$ is an elementary path starting and ending at the depot in which every node is visited exactly once, i.e., no loops or subtours exists. However, not all elementary paths are feasible as they may violate the time windows, truck capacities, battery capacity, or the intra-node distance constraints. The recharging time and layover time should be applied on the elementary paths if necessary to check if time windows of the shipments are respected. These resource constraints are required to be defined and traced while extending the paths. The functions that record the constraint-related resources are called resource extension functions (REFs , see Desaulniers et al. (1998)).

Problem (41) aims at finding a feasible dedicated fleet route with vehicle type k with the minimum reduced cost. The arc cost $c_{ijk}$ for $(i,j,k)\in A\setminus A^{-}$ is replaced by a defined edge cost ${\overline{c}}_{ijk}=c_{ijk}+\frac{1}{2}p-{\mu }_i$, where ${\mu }_0={\mu }_{0^{\prime }}=0$. In order to include the three types of cuts, the reduced cost of the pricing problem needs to be adjusted accordingly. Firstly, let $\pi$ be the dual value associated with the 2-path cuts (39) and $\tau$ be the dual value associated with the generalized multistar cuts (40). Then, for a given shipment candidate set S, $\pi$ and $\tau$ are subtracted from the reduced cost ${\overline{c}}_{ijk}$ of all arcs $(i,j)\in (I\setminus S,S)$, such that $(i,j,k)\notin A^{-}$. Therefore, we have the modified reduced cost for columns of type-k vehicle in the RMP as follows.

Similarly, consider a violated SR inequality of the form (35), defined by a subset of three shipments $S_r\subset I$, $\varphi \le |S_r|$, and ${\chi }_r<0$ as its associated dual value. Since $\chi$ is negative, it will act as a penalty. Because the terms of SR inequalities are summed over the “rows", the penalties cannot be added directly to each arc. The value $\chi$ must be subtracted from the reduced cost ${\widehat{c}}_q$ each time a new column q is passed to the master problem but the column revisits shipments in $S_r$ that define the violated SR cuts. Therefore, an additional resource is required for each violated SR cut to indicate the regeneration times of the shipments. Note that we know such a column is regenerated only if completed, because the total demands fulfilled and other shipments on route could possibly change the left-hand side of the cuts. We have the updated reduced cost for the type-k column as below.

In this setting, the subproblem corresponds to an elementary shortest path problem with resource constraints (ESPPRC) that can be solved by dynamic programming using a labeling algorithm, see Feillet et al. (2004). This well-known method has been successfully applied to many problems in routing and scheduling. Compared with labeling algorithms for other VRP variants in the literature, solving the pricing problem for MVRPC is differentiated by the following factors: 1) multiple REFs are extended separately by truck type k; 2) additional resource constraints are extended to satisfy the working time regulations, recharging needs, and recharger information; 3) necessary layovers and recharging add difficulties to the dominance tests. Moreover, we must be careful about the complexity arisen by these additional constraints. In fact, Dror (1994) has proved that the ESPPRC is strongly NP-hard for VRPTW.

Theorem. The ESPPRC for MVRPC is strongly NP Hard, whereas the SPPRC without elementarity can be solved by route-relaxed dynamic programming in pseudo-polynomial time. Proof. The idea of NP-hardness proof is attributable to Dror (1994). The dynamic programming backward recurrence for a type-k truck of ESPPRC for MVRPC can be stated as: where $F_j^q(S,t)$ denotes the minimal reduced cost of the partial route going from the depot up to node i, S is the node set selected for route q, and t is the time ready to leave node i while $t_{ij}$ is the time spent on the arcs to travel and service. From the constructions of the recurrent function and the graph G, if there exists a feasible shortest path solution from 0 to $0^{\prime }$ for VRPTW on graph G, then there is also a feasible ESPPRC solution satisfying the constraints in (44) for MVRPC. Since the shortest path problem is NP-hard, and we found a polynomial transformation from Dror (1994) to our elementary path problem, our problem is NP-hard as well. Note that, the exact dynamic programming approach in solving ESPPRC yields an exponential number of possible states $(S,t)$.

However, if we shrink the state space by retaining only the m closest shipments to the current end point of the partial path q, which can be visited before i, and relax the elementary path requirement, then the dynamic programming algorithm only explores a pseudo-polynomial number of states between 0 and $\mathcal{O}\left(\right|N|·|t|$). This method is called ng-relaxation. Another known relaxed method is state relaxation where only the solutions with minimal cost are stored and extended at each state. Details refer to Baldacci et al. (2011); Righini and Salani (2008). Therefore, SPPRC is solvable by a route-relaxed method, which gives the lower bound to the solution. This will lead to the effective labeling algorithms we discuss in the following sections. □

We first review the basic concept of labeling algorithm (Irnich and Desaulniers (2005)). In this method, labels are generated and updated to represent partial paths that start at the depot 0. Starting from an initial label associated with node 0, paths are constructed iteratively by extending its label and its descendants forward/backward on graph G. These extensions of a label relay the updates of REFs along an arc. In each iteration, a new extended label is checked for feasibility with respect to the defined critical resource (Feillet et al. (2004)) constraints. If the label breaks any critical resource, it is discarded. To avoid enumerating all feasible partial paths, different dominance criteria are developed and applied to speed up the labeling algorithms. This concept is similar to pruning a tree while eliminating partial paths that cannot generate solutions better than the current bounded minimal reduced cost and pruning for further fathom. When labels are generated in one direction from the origin to the destination, a so-called mono-directional labeling search is performed.

Probably inspired by the method of Divide and Conquer, Righini and Salani (2008) designed a new bidirectional labeling algorithm, in which forward and backward labeling are performed simultaneously often yielding considerable computational time reduction. In this case, the critical resources have to be non-decreasing along the path. In other words, the triangle inequality must be satisfied. For MVRPC, this is indeed the fact, because time windows, truck capacities, layover time accumulations, and total required electricity are all suitable monotone resources. However, the intra-node distance constraint (18) does not satisfy the triangle inequality, i.e., even if $d_{ij}>d^{max}$ for the current partial path q with end node i, j could also be feasible for later extensions of the path. Moreover, a midpoint $M\in [M_0,M_1]$ is required to be specified. For instance, if we choose time windows as the resource, any $M\in [a_0,b_{0^{\prime }}]$ can be a candidate midpoint. Bidirectional labeling algorithm works in three steps. First, labels are extended forward from the sink to the midpoint M; second, labels are extended from destination backward to M; third, two labels are merged and checked for feasibility.

We will discuss how to develop efficient mono-directional and bidirectional labeling algorithms for MVRPC in the following sections.

In this section, we describe the mono-directional labeling algorithm for the ESPPRC, where the reduced cost of the subproblem is defined by (43) with the three strengthened cuts. The labeling algorithm seeks to find a negative reduced-cost route that is feasible with respect to truck capacities, time windows and the current U.S. hours of service regulations.

The forward label search concept discussed here involves a relaxation technique, in which the state space explored by the exact dynamic programming method is projected onto a lower dimensional space, so that certain “non-critical" nodes are allowed to be visited several times. Different from normal VRPTW or EVRPTW, the pricing problem of MVRPC requires additional resources to keep track of the accumulated work time since the last layover and the remaining time to reach the next customer j. Goel and Irnich (2014) has illustrated that it is effective to use the minimal layover time when updating the time REFs and add necessary hours to each layover later on to avoid waiting time. We adopt the same idea and incorporate the possible recharging time for EV routes. In a forward mono-directional labeling algorithm for MVRPC, each label L corresponds to a partial path q starting at the depot 0 up to a node $i\in {\overline{I}}_0$. In the following descriptions of labeling algorithms, we will keep using $0^{\prime }$ as a duplicate of the depot, which is the destination of each route. For a type-k truck, a label $L_i$ is defined by the following attributes:

$L_i\left(cost\right)$	Reduced cost of partial path q.
$L_i\left(demand\right)$	Accumulated demand fulfilled by partial path q.
$L_i\left(recharge\right)$	Number of recharges performed along path q.
$L_i\left(time\right)$	Service start time at node $i\in {\overline{I}}_0$ when reached through partial
	path q.
$L_i\left(electric\right)$	Accumulated electricity consumption (in hours) from depot or a
	recharger up to node i along the partial path q. Note this resource
	is renewable.
$L_i(visit,n)$	A dummy resource variable that records the number of times any
	critical shipment $n\in S$ is visited along path q. The path is
	elementary only if $\forall L_i(visit,n)\le 1$.
$L_i\left(elapsed\right)$	Accumulated elapsed time (hours) since the last layover (or depot
	if no layover).
$L_i\left(work\right)$	Accumulated working time (hours) since the last layover (or depot
	if no layover).
$L_i\left(layover\right)$	Number of layovers the driver has taken up to node i.

To further shrink the searching space, we introduce the idea of unreachable nodes, see Feillet et al. (2004). In the pricing problem of MVPRC, our considered critical resources obey the triangle inequality. Then, we can identify the nodes that cannot be visited by any extension of $L_i$ because the time window and the capacity constraints are violated. These nodes are called unreachable by the label. If a label is assigned by an ICV, all of the available rechargers will not be considered in the search algorithm, and the corresponding REFs are set as 0. When a node is set as unreachable by a label, the dummy resource variable $L_i(visit,n)=1$, i.e., as if the node n had been visited. Feillet et al. (2004) observed that setting $L_i(visit,n)$ to 1 will not change the structure of the label search but will reduce the computational time. Therefore, we define another attribute ${\overline{L}}_i\left(nodes\right)$ as the infeasible set that cannot be extended from the current state $L_i$. That is, ${\overline{L}}_i\left(nodes\right)=L_i\left(nodes\right)\cup \{j\in {\overline{I}}_0\setminus L_i\left(nodes\right):L_i\left(time\right)+\frac{d_{ij}}{g}+{\lambda }_i>b_j\vee L_i\left(demand\right)+{\omega }_j>Q_k\vee \forall j\in E,k\in V_D\}$. In addition, we introduce $Tab{u}_i$ as a candidate set of arcs that are feasible to be explored. Its definition will be introduced later in the section on implementation strategies.

In the initial state at node 0, all attributes of the label L are set to 0 except for $L_0\left(latest\right)$, which is equal to ∞. The extension of a label $L_i$=$(L_i\left(cost\right)$, $L_i\left(demand\right)$, $L_i\left(recharge\right)$, $L_i\left(time\right)$, $L_i\left(electric\right)$, $L_i(visit,n)$, $L_i\left(elapsed\right)$, $L_i\left(work\right)$, $L_i\left(layover\right)$, $L_i\left(latest\right)$, $L_i\left(nodes\right)$, ${\overline{L}}_i\left(nodes\right))$ along an arc $(i,j)\in A$ is performed for each vehicle type k using the following REFs:

While updating the time resource $L_i\left(time\right)$, we need to consider the possible recharging time and the viable layover time. We first define an alert indicator ${\Delta }_{ij}$ passing through the arc $(i,j)$ to check whether or not the service hours regulation is broken, where ${\Delta }_{ij}=min\{0,\xi -L_i\left(work\right)-(\frac{d_{ij}}{g}+{\lambda }_i)\}$. A direct extension from i to j can be performed only if ${\Delta }_{ij}=0$, otherwise, a layover must be added since the service hours regulation is about to be violated. Considering the two situations, we define the REF for $L_i\left(time\right)$ as:

Note that in REF (49), $L_i\left(electric\right)={\delta }_{ij}=0$ for all partial path q that travels by ICV. Next, we define a function named $unreachabl{e}_n\left(i\right)$ to measure the feasibility status of $\forall n\in S$ for the current label $L_i$. If the time window and capacity resources are not violated, then this function returns 1, otherwise, it returns 0.

REF (53) explains how prolonged layover period(s) can be added to a partial path. For every non-depot start node on a partial path, we assume it is delivered at its beginning time window $a_i$. The layover periods are specified whenever needed with the minimum hours ${\beta }^{lb}$ as shown in REFs (48) and (52). By this scheduling policy, we avoid times during which the driver is not productive. However, this may also cause early arrivals due to time window constraints at subsequent shipments and thus miss some valid routes. By lengthening the duration of the layover periods between $[{\beta }^{lb},{\beta }^{ub}]$, we allow the possibility of measuring the before-missed routes. However, the arrival time at already routed subsequent shipments may be pushed out of their respective time windows. We therefore must calculate the resource values of $L_i\left(latest\right)$ indicating the maximum amount by which the duration of a layover can be extended without violating any resource constraints. REF (54) says that the latest end time $L_j\left(latest\right)$ of the last layover can be updated by minimizing values between the current latest time $L_i\left(latest\right)$ and the latest possible end time of the last layover without violating node j’s time window. In some cases, $L_j\left(time\right)$ is allowed to be deferred, such that we use the term $[L_i\left(latest\right)+L_i\left(elapsed\right)]$ to represent the arrival time at i. Let the term $[L_j\left(time\right)-L_i\left(time\right)]$ be the work and potential layover period between i and j. As shown in (54), $L_j\left(elapsed\right)$ is the truncated working time left after an immediate layover. Then, the duration of the latest layover can be lengthened by any value less than or equal to $NewTime=L_j\left(latest\right)-[L_j\left(time\right)-L_i\left(time\right)+L_i\left(latest\right)+L_i\left(elapsed\right)-L_j\left(elapsed\right)]$. In other words, at each state, only if the following condition (56) is satisfied can we continue extending the label.

Because, if a layover is applied on the corresponding arc $(i,j),$ the layover period can only be adjusted between $[{\beta }^{lb},{\beta }^{ub}]$. For arcs $(i,j)$ that do not take any layover either because the work time limit is never hit or because layovers have been taken in previous states, based on our policy, the current time can still be adjusted to satisfy the time window constraints for subsequent nodes. Therefore, whether or not there is a layover added on the arc, we need to check the following condition to ensure adjusting the time resource $L_j\left(time\right)$ is feasible.

When checking whether or not a new label extension from i to j satisfies the time window resource constraint, we initially need to check if $L_j\left(time\right)\le b_j$ is feasible. Next, suppose (56) and (57) are not violated under the defined situations, we should update the current route time $L_j\left(time\right)$ to ensure the earliest acceptable time $a_j$ is respected.

In each label extension, the REFs defined in (45)-(58) are updated by the following sequence.