A Hybrid Pricing and Cutting Approach for the Multi-Shift Full Truckload Vehicle Routing Problem

Full truckload transportation (FTL) in the form of freight containers represents one of the most important transportation modes in international trade. Due to large volume and scale, in FTL, delivery time is often less critical but cost and service quality are crucial. Therefore, efficiently solving large scale multiple shift FTL problems is becoming more and more important and requires further research. In one of our earlier studies, a set covering model and a three-stage solution method were developed for a multi-shift FTL problem. This paper extends the previous work and presents a significantly more efficient approach by hybridising pricing and cutting strategies with metaheuristics (a variable neighbourhood search and a genetic algorithm). The metaheuristics were adopted to find promising columns (vehicle routes) guided by pricing and cuts are dynamically generated to eliminate infeasible flow assignments caused by incompatible commodities. Computational experiments on real-life and artificial benchmark FTL problems showed superior performance both in terms of computational time and solution quality, when compared with previous MIP based three-stage methods and two existing metaheuristics. The proposed cutting and heuristic pricing approach can efficiently solve large scale real-life FTL problems.


Introduction
In intermodal freight transportation, a large proportion of container transportation is carried out by barges, trains or ocean-going vessels Braekers et al. (2014). Container movement activities between intermodal terminals, depots and shippers are also referred to as drayage operations and such activities are usually performed by trucks. Although drayage operations represent a small fraction of the total distance of an intermodal freight transportation, they constitute a substantial share of the shipping costs Smilowitz (2006). Consequently, major port are facing intense competition and pressure to improve the efficiency of drayage operations.
Due to labour laws and other constraints related to the working time of drivers, shift based working schedules are becoming a common practice in the transportation industry (e.g. taxis and buses). The problem studied in this paper concerns the movement of containers (commodities) between a number of terminals (docks) within a short distance located in a large international port using a homogeneous truck fleet. The transportation time window of commodities usually spans from a few hours up to several days, and covers of multiple working shifts. Thus the problem that we are addressing is essentially different from the single-shift problem studied in most of the existing full truckload routing problems (e.g. Zhang et al. (2010), Braekers et al. (2014)) because the planning horizon covers several shifts and determining the transportation shift of each commodity forms part of the optimisation decision.
In our earlier study Bai et al. (2015), a set covering model and a three-stage method were proposed for this problem. However, the computational time to optimally solve large size problems was prohibitive. Chen et al. (2013) investigate a reactive shaking variable neighbourhood search (rsVNS) and a simulated annealing hyper-heuristic method (SAHH) Chen (2016) for this problem. The rsVNS extends the original VNS which utilises the systematic changes of multiple neighbourhood functions to achieve convergence and diversification.
The SAHH applies a reinforcement learning based neighbourhood selection mechanism within a simulated annealing framework. The learning mechanism aims to adapt the algorithm to different problem instances and search scenarios by dynamically adjusting the neighbourhood selection strategies. Both rsVNS and SAHH were able to obtain feasible but inferior solutions with less computational time compared with the three-stage method.
The main contributions of this paper are twofold: 1) We fully explore the advantageous features of a previously proposed indirect solution encoding scheme, leading to some insightful findings of the differences between the multi-shift FTL problems and traditional pickup and delivery problems; 2) A pricing based column generation method is investigated, in conjunction with dynamic cuts. Our method is inspired by branch-price-and-cut algorithms but differs in three major ways. Firstly, we do not solve the pricing subproblem exactly. Instead, we quickly produce approximate solutions by introducing several optimisation strategies (See Section 5.4). Secondly, we employ metaheuristics instead of traversing the entire search tree. Thirdly, the cuts are added after the column generation process because the infeasible flow assignments (See Section 4.2) rarely occurs and they take time to be evaluated and fixed in every column generation iteration. The new algorithm can significantly speed up the solution time for large size problems. Moreover, the solution quality is also improved. Two pricing methods are proposed and tested on both real-life and artificial instances.
The remainder of this paper is organised as follows: the problem is described in detail in Section 2; a literature review is given in Section 3 followed with the mathematical model of the problem in Section 4. The proposed pricing and cutting method is illustrated in Section 5 and its the computational experiments are presented in Sections 6 and 7. Finally conclusions are drawn in Section 8.

Problem Description
The multi-shift FTL problem is concerned with transporting a set of full truckload freights (containers) between a given number of terminals within multiple working shifts. Both the operational time windows of the freights and the planning horizon can span across several shifts. Although each container is transported in a single shift, its time window covering multi-shifts and determination of the shift in which this load is serviced (transported) forms part of the decisions to optimise. The objective is to minimise the total cost while satisfying various constraints.
First, each full truckload commodity (container) has an available time for pickup and a deadline for delivery. Second, during each shift, a number of unit-capacity trucks start from the deport at the start of the shift, complete a number of transportation requests and then return to the depot before a shift ends. Finally, a service time is applied during both pickup and delivery. To clarify, we refer to the request of a full truckload movement as one unit of a commodity. A commodity is a collection of requests of full truckload freights that share identical sources, destinations and time windows.
In the context of real-world applications that this research tries to address, the total quantities of all the requests within a planning horizon can be very large (more than 1000). The number of terminals is relatively small (less than 10) and the distances between these terminals are relatively short (all reachable within a shift). These features make this problem different from problems that are studied in previous work. It has been shown in Bai et al. (2015) that a model based on a set covering formulation is more promising for this problem than other node based formulations. However, the features and reasons have not yet been sufficiently analysed. For the completeness of this paper, we include the model in Section 4, along with a detailed discussion of its advantages and disadvantages.
A literature review about the real-world applications of similar problems is given in the next section.

Literature Review
The drayage operations problem is a typical case of bidirectional multi-shift full truckload vehicle routing problems. Bai et al. (2015) highlight the core features of these truckload vehicle routing problems and discuss relationships with other variants of vehicle routing problems (VRP) from three aspects, including the directions of the flow, existence of consolidation or not, and length of the planning horizon. Here, we summaries the relevant research on the drayage operation problems which we broadly classify into drayage operations with and without relocation requirements of empty containers.

Drayage problem without relocation of empty containers
Xiubin Wang (2002) model a full truckload pickup and delivery problem with time windows (FT-PDPTW) as an asymmetric multiple travelling salesman problem with time windows (m-TSPTW) and propose a time-window discretisation scheme. Jula et al. (2005) extend the m-TSPTW model with social constraints and propose an exact algorithm based on dynamic programming. Moreover, a hybrid method combining dynamic programming and genetic algorithms (GAs) is also investigated, as well as an insertion heuristic method. Chung et al. (2007) design several types of formulations for practical container road transportation problems. The basic problem is formulated as an m-TSPTW problem, which is solved by an insertion heuristic. Gendreau et al. (2015) refer to this routing problem as the one-commodity Full-Truckload Pickup-and-Delivery Problem (1-FTPDP) and present three mathematical formulations with branch-and-cut algorithms to optimally solve the model formulations. Lai et al. (2013) propose a new routing problem that can be viewed as a vehicle routing problem with clustered backhauls (VRPCB). Solutions are obtained with the Clarke-and-Wright algorithm and improved further by a neighbourhood based metaheuristic . This work is also compared in the study of a problem with single and double container loads Ghezelsoflu et al. (2018). The distribution of more-than-one container per truck by different types of trucks has also been studied in Vidović et al. (2017) and Funke and Kopfer (2016). Soares et al. (2019) study an FTL problem with multiple types of vehicle synchronisations. A MIP model and a heuristic solution method based on the fix-and-optimise principles are proposed.

Drayage problems with relocation of empty containers
Efforts to combine the planning of loaded and empty container transports are made by several authors. Coslovich et al. (2006) analyse a fleet management problem for a container transportation company by decomposing the problem into three subproblems, which are then solved using a Lagrangian relaxation. Ileri et al. (2006) present a column generation based approach for solving a daily drayage problem. Smilowitz (2006) model a drayage operation with empty repositioning choices as a multi-resource routing problem (MRRP) with flexible tasks. The solution approach is a column generation algorithm embedded in a branch-and-bound framework. Imai et al. (2007) formulate a container transportation problem as a vehicle routing problem with full container loads (VRPFC) and solve it with a subgradient heuristic based on Lagrangian relaxation. Caris and Janssens (2009) extend this work and model the problem as a FT-PDPTW. A local search heuristic is proposed. The work is further extended by using a deterministic annealing algorithm suggested in Caris and Janssens (2010). Zhang et al. (2010) improve the time window partitioning scheme used in Xiubin Wang (2002) for container transportation in a local area with multiple depots and multiple terminals. The results indicate that good performance can be achieved compared with a reactive tabu search (RTS) method demonstrated in Ruiyou Zhang (2009). Zhang et al. (2011) also investigate the single depot and terminal problem. Again, an RTS is proposed. Vidovic et al. (2011) extend the problem analysed by Zhang et al. (2010) and Imai et al. (2007) to the multi-commodity case and formulate it as a multiple matching problem. Solutions are obtained via a heuristic approach based on calculating utilities of matching nodes. Nossack and Pesch (2013) present a new formulation for the truck scheduling problem based on a FT-PDPTW and propose a two-stage heuristic solution approach. Braekers et al. (2013) investigate a sequential and an integrated approach to plan loaded and empty container drayage operations. A single-and a two-phase deterministic annealing algorithm are presented. This solution approach is further adapted in Braekers et al. (2011) to take a bi-objective optimisation function into account. The algorithms are further improved in Braekers et al. (2014). More recently, Xie et al. (2017) investigate the empty container delivery problem in an intermodal transport system composed of a sea liner firm and a rail firm. Apart from transportation cost, the difference in marginal profits between the seaport and dry port is also considered in the objective function.
Some researchers examine drayage operations problems in dynamical situations. A survey on dynamic and stochastic vehicle routing problems can be found in Ritzinger et al. (2016).
Most of the aforementioned research work has been trying to formulate the drayage problem as some forms of classical vehicle routing problems in order to exploit the time constraint structures to prune the search space. However, as discussed in Bai et al. (2015), this type of formulations does not work well for problems where time related constraints are not very tight and node-based solution representations generally lead to unnecessarily large search space, resulting to inefficient solution methods.

Hybridising exact methods and (meta)-heuristics
This paper studies an indirect solution representation for the multi-shift FTL problem that addresses these issues and contributes to the body of research with an efficient column generation method. In many vehicle routing applications solved by column generation, the subproblem is usually viewed as an elementary shortest path problem with resource constraints or one of its variants. Nowadays, an increasing number of hybridisations between heuristics and exact approaches are developed. These methods can provide a good compromise between solution quality and computational time as they adopt the advantages of both types of methods. Puchinger and Raidl (2005) classified hybridisation of exact algorithms and (meta)-heuristics into four types, we briefly introduce them and give examples for each: 1) Collaborative Combinations -sequential execution: In this type of hybridisation, either the heuristic is executed before the exact method, or vice-versa. For example, when solving a set covering problem, a heuristic is used to generate a set of feasible columns and the exact method is used to find an optimal solution from the feasible columns. Examples of this type of hybridisation can be found in Clements et al. (1997) and Vasquez et al. (2001).
2) Collaborative Combinations -parallel or intertwined execution: Instead of executing either heuristics or exact methods sequentially, this type of method implements the algorithms in a parallel or intertwined way. Clusters or multi-processors are used to deploy the parallel implementations. There are several frameworks proposed to facilitate such implementations, such as Alba et al. (2002) Vidal et al. (2014) and Lahrichi et al. (2015).

3) Integrative Combinations -incorporating exact algorithms in heuristics:
Where exact algorithms are subordinately embedded within heuristics. For example, the solution of LP-relaxation and its dual values can be utilised in heuristically guiding neighbourhood search. Applications can be found in Marino et al. (1999) and .

4) Integrative Combinations -incorporating heuristics in exact algorithms:
This type of hybridisation is analogous with the previous one, but heuristics are embedded within exact algorithms. For example, heuristics can be used to determine bounds in branch-and-bound algorithms. Heuristics can also be used to search for columns with negative costs in the branch-and-price approach. Applications of this hybridisation method can be found in  and Strandmark et al. (2020). The column generation based method proposed in this paper falls into this category.
Please refer to Blum et al. (2011) and Muthuraman and Venkatesan (2017) for more comprehensive reviews of the hybridisation approach.

Model Formulation
The problem studied here can be defined on a graph G = (N, A) where each node i ∈ N represents a physical terminal (including the depot, i = 0). An arc (i, j) between nodes i, j ∈ N is included in the arc set A if the visit of j can be performed immediately after i. A service time t i is applied to each node i to represent the loading/unloading times of truckload commodities and the travel time of arc (i, j) is denoted as µ ij . All trucks must depart from and return to node 0 (depot). Let R be the set of all feasible routes that a truck can execute in a working shift without the complication of time window requirements from commodities. Therefore, each route r ∈ R is called distance-wise feasible.
For a given shift s, the i-th node in route r ∈ R (denoted as r i ) can only be visited within a time window (e s r i , l s r i ) where e s r i is the earliest time that a truck covering route r can start Example of a routing sharing among five commodities.
a pickup or delivery operation while l s r i is the latest time that a truck must depart from the node. Let t r i be the service time at node r i . In each route r encoding, a duplicated node is inserted if the node involves both a loading and an unloading operations (i.e. this node is both the destination and source of two different commodity flows). Therefore, if the nodes is 0-indexed in a route, loading operations are always conducted at the odd indexed nodes of a route (see Eq. (3)) and unloading operations are at the even-indexed nodes. e s r i and l s r i can be calculated using the following recursive equations: Let K denote the set of commodities for delivery. Each commodity k ∈ K is defined by a tuple (o(k), d(k), Q(k), σ(k), τ (k)), which, respectively, are the origin, destination, quantity, availability time and deadline of commodity k. Denote δ ks r i the binary variable to indicate whether the i-th node of route r can be the loading node for commodity k in shift s (δ ks r i = 1) or not (δ ks r i = 0). To speed up the computation, δ ks r i can be pre-calculated by checking the following conditions: Figure 1 presents a simple example of a feasible truck route where 0 denotes the depot.
For a 0-indexed route node list, odd numbered nodes are commodity loading nodes while even numbered nodes are unloading nodes. If a node on a route is involved with both loading and unloading, a copy of it is created so that the above rules are maintained (see more details in Bai et al. (2015)). In this example, a truck departs from the depot and picks up a unit of commodity from either commodity 1, commodity 2 or commodity 3 from node 1 and unload the commodity at node 3. Then the truck picks another unit of commodity (either commodity 4 or commodity 5) at node 4 and unload at node 2 before the truck returns to the depot.
In summary, the following notations are used for the formulation:

Sets
• N : Set of nodes in the transportation network.
• S : List of time-continuous shifts in the planning horizon.
• R : Set of all feasible truck routes within a shift.
• K : Set of full truckload commodities to be delivered.

Other parameters
• d r : The cost (distance) of route r.
• n : The maximum number of trucks available for use.

Decision variables
• x ks r i : Commodity flow of the ith node of r in s for servicing commodity k ∈ K. • y s r : The number of times a given route r ∈ R is used during shift s ∈ S and y s r ∈ N + . The model for this multi-shift FTL problem can be formulated as the follows: x ks r i ≤ δ ks r i y s r ∀i ∈ r, ∀r ∈ R, ∀k ∈ K, ∀s ∈ S (12) x ks r i = Z + ∀i ∈ r, ∀r ∈ R, ∀k ∈ K, ∀s ∈ S (13) The objective function (8) minimises the aggregated distance of all routes used in a solution. Constraint (9) restricts the total number of trucks being used in a solution.
Constraint (10) ensures all the commodities are serviced (transported) entirely. Constraint (11) requires that each arc of a route in a shift can only be used y s r times. Constraint (12) makes sure that any positive x ks r i is feasible in terms of the source, destination and operation time window of commodity k. Since binary indicator δ ks r i can be pre-calculated, this constraint can be eliminated by removing the corresponding flow variables x ks r i from the model when δ ks r i takes value of 0. This is indeed how the model was implemented in our experiments because the resulting model is a lot smaller. Constraints (13) and (14) define the domains of the decision variables.

Merits of this solution encoding
One of the most helpful benefits of this solution encoding scheme is the transformation of a previous m-TSPTW based non-linear model (e.g. the model proposed by Chen (2016)) into a linear integer model, so it can be solved using various integer programming techniques.
This was done through hiding nonlinear time related constraints into the generation of the shift-independent feasible truck route set.
For some applications (e.g. FTL with a small number of terminals), pre-computing all feasible routes is possible since the time related constraints in this problem are slightly different from those in the traditional pickup and delivery problem with time windows (PDPTW). In this multi-shift FTL problem, each commodity k has an operation time window (σ(k), τ (k)) defining its availability time and the delivery deadline. Time constraints require that both the pickup and delivery operations occur within this time window for commodity k. In PDPTW problems, two separate time windows are used, one for pickup and the other for delivery. Note that for non-time critical full truckload transportation, having one time window is reasonable since all the terminals (nodes) operate all the time, and having short time windows for both pickup and delivery does not make sense, although we acknowledge it is very different for express deliveries which are mostly for household customers.
A second benefit of this solution representation is its capability to handle nonlinear cost functions. For example, the costs of routes could be a nonlinear, complex function of the distance. It also permits to include various other constraints related to drivers (e.g. maximum driving distance, time or preferred terminals). compared with a commonly used m-TSPTW formulation, in which each container is modelled as a node in the graph. For a problem instance with hundreds or even thousands of truckload sized containers, the corresponding graph in m-TSPTW formulation could be prohibitively expensive to handle. However, in the real-world problem that we are concerned with, containers often arrive in large batches with same requirements (i.e. same O-D pairs and time windows). In an m-TSPTW formulation, any swaps of positions of these nodes (i.e. containers) in the TSP tours shall result in the same objective value (i.e. many-to-one mapping from solution encoding and objective space). This leads to a significantly larger search space with a plateau. In our proposed formulations, containers with the same property are grouped as one commodity, leading to a one-to-one mapping and a much smaller search space.

Dealing with non-compatible commodities
Although for all the practical instances that we extracted from real-world problems, the FTL model in Section 4 produces solutions that satisfy practical constraints. However, it is possible to artificially generate problem instances that the proposed FTL model returns an infeasible solution. That is, the solution is feasible for the FTL model but may still violate the time window constraints for some commodities. This happens when two time non-compatible commodities are assigned to a same route and same shift. An example of such cases is illustrated in Figure 2.
Service time push back and non-compatible commodities.
In this example, we included two non-compatible commodities (k and v) that are serviced using a feasible route r at arcs (i, i + 1) and (j, j + 1), respectively. The solution becomes infeasible when the service time of the first commodity k was delayed because the commodity available time σ(k) is later than node i's earliest arrival time e r i . Because of this, the service time of commodity k at node i is pushed back by As a service node can serve multiple commodities (e.g. commodity 1, commodity 2, and commodity 3 served by node 1 in Figure 1), if more than one commodities lead to push back, then push back(i) is calculated as the maximum value of puch backs of all commodities that served by node i.
Unless there are larger push backs by other commodities in the remaining route nodes, the push back at node i is propagated in its entirety to all the remaining nodes in the route. A violation of another following commodity v's shipment deadline (τ (v)) shall occur if the following condition is satisfied: That is, if a push back caused by a previous commodity is greater than the difference between the earliest vehicle arrival time of its destination node and a commodity's deadline, the commodity assignments along this route become infeasible.
If the resulting solution sequentially assigns two non-compatible commodities, k and v, at nodes i and j, respectively, of a same route r in the same shift s, then the following constraints should be added to ensure v is not inserted at or after k in the same route and shift.
x ks where θ is an auxiliary variable taking either 0 or 1 and M is a large positive number. K r contains the commodities that σ(k) > e r i , V r contains the commodities that τ (v) < e r j+1 .
Note that constraint (17) also prevents cases of non-compatible commodity assignments at a same node. The process of generating the cuts through constraint (17) to eliminate non-compatible commodity assignments is given in Algorithm 1.
We do not want to strongly restrict the non-compatible commodities, as shown in the above example, v is not simply forbidden to be served by the node, instead, it is still allowed to be served by the route as long as the compatibility constraint of k and v is not violated. From Algorithm 1 it can be seen that the procedure keeps tracking the push back time (push back) of each commodity in K r and maximum allowed push back time (acceptable push back) of each commodity in V r . The cut will be added to the model only if any pairs of incompatible commodities were found. That means even if the service time of the commodities served between k and v, if any, are pushed back, they are not restricted by constraint (17)   for v in W (j) do, where where W (j) is set of commodities serviced by node j.   Table 1 gives a simple illustrative example how the valid constraints can be dynamically generated into the model to prune the search space. In this example, A total of 4 feasible routes are available for selection to deliver 4 commodity k 1 , k 2 , v 1 , v 2 . They are (0,1,2,3,4,0), (0,3,4,1,2,0), (0,1,2,0), (0,3,4,0) with distances of 79, 129, 64, 75, respectively. Without considering push backs, the optimal solution is to choose route 1 twice (i.e. y s 1 = 2), with flows of x k 1 s 1 1 = 1, x k 2 s 1 1 = 1, x v 1 s 1 3 = 1, x v 2 s 1 3 = 1 because it satisfies formulas (4) to (7) and requires the least travel distance (158) to delivery all commodities. However, since σ(k 1 ) > e 1 1 (13 : 40 > 8 : 15), the service of commodity k 1 at node 1 is pushed back to 13:40 by 325 minutes, which is propagated to all the remaining nodes in route 1 (the updated e and l after pushed back by k 1 are denoted as e and l in Table 2).

Dealing with a very large set R
The proposed model also has some problems. The most critical one is the size of the feasible route set R which can increase exponentially with the number of nodes (or terminals). In Bai et al. (2015), some real-life problems have certain special features to permit some of nodes being merged, and a three-stage algorithm was able to find near optimal solutions.
However, in addition to the excessive computational time of the three-stage algorithm, the method becomes invalid for problems that do not possess these features to permit node merging.
In this paper we propose to use a column generation method to address this issue. The idea is to use the pricing information to guide the generation of promising feasible routes dynamically.

A Hybrid Column Generation Method
Column Generation is an effective approach for solving large scale integer programming problems (i.e. problems with large number of columns). It is a potentially very good method Table 2 An for the problem formulation stated in Section 4, where the feasible route set R is very large, leading to a model with a huge number of columns while the optimal solution uses a very small subset of it. We propose to use the column generation approach for this problem in which the sub-problem (pricing problem) is solved to identify the variables that should enter the basis.

The proposed solution framework
The integer programming formulation presented in section 4 is also referred to as the master problem. The Restricted Master Problem (RMP) is the master problem that considers only of a subset of truck routes R that are generated by the pricing problem (subproblem) using the dual information obtained from the Linear Programming Relaxation (LPR) of the RMP. The pricing problem and the LRP will be discussed in Section 5.4 and Section 5.3, respectively. Before the RMP is solved for the first time, no dual information is available and an initial truck routes set (see Section 5.2) is thus required to start the process. Then the LPR is solved to optimality and the dual information is obtained for calculating the reduced costs of routes during the pricing subproblem. The overall solution framework is outlined in Figure 3, followed by detailed steps of the procedure.
Our intial experiments showed that majority of computing time is consumed by the RMP solving. The algorithm was thus modified to accelerate convergence by adding multiple routes with negative reduced costs at each iteration. Details of the methods for the pricing problem is given in Section 5.4. It is hoped that by doing this the total number of RMP calls can be reduced. This process is repeated until the stopping criteria are met. Finally, in order to obtain the integer solutions, relaxed constraints associated with x ks r i and y s r are set back to their original ones during the final RMP solving.
Because the pricing problem is solved repeatedly in the column generation framework, it is crucial that the solution algorithm for the pricing subproblem is as efficient as possible. Therefore, we propose two different strategies, one for problems with small-sized R and one for problems with a large R. For the former case, we propose to adapt an explicit enumerative generation of R as priori and then try to solve the pricing subproblem when no column with negative reduced cost can be found. We apply a recursive algorithm to generate all feasible routes as described in Bai et al. (2015) before the iterative procedure starts. In the case of a large R, we propose to use heuristic approaches (see 5.5) to solve the pricing problem and the stopping criteria of the heuristic is an limitation of the number of RMP cycles (denoted by Finish in Figure 3). The overall solution framework as described above is outlined in Figure 3.

Initial set of routes
Before the RMP is solved for the first time, no dual information is available and an initial set of columns is required to start the process. We apply two methods described in detail in the next two subsections to generate an initial set of columns (routes).

Simple route initialisation
A prerequisite of constructing a basic route set is to ensure that each commodity has at least one route to service it. Thus the simplest solution is to generate a dedicated route for each commodity, in which an empty truck leaves the depot and travels to the source of a commodity, loads the commodity and delivers it to its destination. After that, the truck returns to the depot. This method works fine in some cases but may of course lead to an infeasible solution in terms of the maximum number of vehicles constraint (9). Framework of the proposed column generation process 5.2.2. Insertion heuristic method The advantage of the basic route initialisation in the previous section is its simplicity. However, very rarely will these routes be used in the optimal solutions, neither do they resemble any of the routes that are present in the optimal solutions. In this study, we proposes to use constructive heuristic methods to generate these initial routes for our column generation method. In particular, we used the same insertion heuristics described in (Chen et al. 2013). To construct routes, the task that cause minimum empty load distance is inserted by following two initialisation criteria: First, the most urgent tasks that have deadlines closer to the shift start time are inserted.
The second criterion considers tasks that have earlier availability time.
There are two benefits here. First, because the constructive heuristic produces a feasible solution for the original problem, the vehicle routes extracted from the solution shall also produce a feasible solution in our column generation method, satisfying the maximum number of vehicles constraint (constraint (9)). Second, because the pricing subproblem is solved heuristically, starting from a good set of initial vehicle routes will enable the column generation method to generate high quality solutions more quickly compared to the simple route initialisation method. The proposed method will converge to a high quality solution much faster.

Linear programming relaxation (LPR)
Linear programming relaxation (LPR) relaxes the discrete variables constraints and differs from the model presented in Section 4 in the route set, which should be attained as a R (R ⊂ R) resultant from the pricing subproblem instead of all feasible routes R. Let LPR be the relaxed model, and Let LPR-(9-12) be the constraints corresponding to the constraints (9-12) of the master problem.

Pricing methods
We present three different route price estimation methods: The first method obtains solution by enumerating and examining all possible commodity assignments for each route and the best route (i.e. the route with the most negative reduced cost) is selected. This method is referred to as Pricing problem by enumeration in this paper. For efficiency, two other pricing estimation methods (Average pricing and Demand weighted average pricing) are also investigated.
Pricing problem by enumeration. Let α s , π k , β s r i , γ ks r i be the pricing variables for constraints from the LPR-(9-12), respectively. The reduced cost for route r in shift s is: However, since routes are generated independent of shifts, the following average reduced cost is computed over all shifts for each route.
where |S| is the total number of shifts in the planning. Let W = {w 1 , w 2 , ..., } be the set of all possible commodity assignments, each of which can be delivered by one instance of route r. In the example of the route in Figure 1,  assignments for a given route r, we evaluate them all and if the reduced cost of any given w ∈ W is found to be negative for route r, it is added to the RMP. The same process is repeated for the next route r + 1 until all feasible routes are evaluated. This searching process guarantee that the reduce cost of commodity assignment in each route is examined but its efficiency is low as some routes may contain thousands of possible commodity assignments.
In order to obtain good results in a reasonable time, we investigated the following steps to improve efficiency: Firstly, the constraintsf LPR-(12) are pre-processed offline. Given a route and its start time, the feasibility of commodity in a route can be determined by formula (3)- (7) offline. This allows us to reduce a large number of decision variables that have to be handled by the model. Consequently, we lost the price values (γ) associated with the feasibility of commodity and time window of the service node.
Secondly, we do not want to explicitly restrict which shift that a route belongs to, as the feasible route set is meant to be same across all shifts. This has benefits from management standpoints too because drivers proficiency can be improved if they are asked to follow a fixed set of routes repeatedly. Also after long run, the set of frequently used routes can become part of the knowledge system of the transportation planning and time consuming column generation procedure may not be required anymore. Therefore, the price values of arcs (β) in each shift is not used because the efficiency is substantially degraded by generating routes dependent of shifts. Fortunately, the price of an arc can be estimated by the price of all possible commodities (π) for all shifts that can be serviced by the arc.
Thirdly, the constraint related to truck numbers is also not involved for the reduced cost calculation, due to in real-life problems, vehicle number is not critical but the efficiency is, leading to α taking zeros for all of our instances. In the end, we came to two approximated pricing methods, illustrated below.
P1: average pricing Instead of enumerating all the commodity combinations of a route and then checking the c r for each of w i , a more efficient approach is to use the average prices to estimate c r approximately. More specifically, let J be the set of all service nodes in r (e.g. nodes {1, 4} in Figure 1). Denote V j be the set of all commodities that can be serviced by a node j in r. The reduced cost c r for route r is calculated by the following equation: P2: demand weighted average pricing Though the commodities processed by a service node in a route share the same source and destination node, the quantity of the commodities varies from one to another. The simple average pricing method P1 fails to take into account the quantity of the commodities, so that large quantity commodities may be left unpaired to improve the efficiency. Therefore, the demand weighted average method tries to give priority to large commodities at the early stage. A weight ω k that is proportional to the commodity quantity Q(k) is used. The weighted average pricing method uses the following equation to estimate the reduced cost.
5.5. Heuristic column generator for large R As can be seen from Figure 3, a heuristic column generator is used within the column generation framework. As optimally solving the pricing problem involves an expensive recursive tree search, we propose to use a variable neighbourhood search (VNS) and a genetic algorithm (GA) to tackle the pricing subproblem. The goal of the metaheuristics is to identify new columns with negative reduced costs. The idea is that, instead of generating a new column (i.e. route) from scratch, it is probably more efficient to search from the existing routes through either neighbourhood moves or route combinations (i.e. crossovers).
VNS and GA are widely adopted excellent frameworks to implement these ideas. The main difference here is that the metaheuristics are guided by an objective function that heavily relies on the pricing information obtained from the linear program relaxations.

VNS
The pseudo-code of our VNS algorithm is given in Algorithm 2 and the parameters of the algorithms are listed in Table 3. In our VNS method, the neighbourhood functions include swap, 2-opt, and relocate. These operators are very similar to those used in solving the classical VRP problems. For example, the swap operator swaps two arcs of two different routes. The 2-opt exchanges two nodes on the same route. The relocate operator relocates an arc from its current route to a different one. By exploring different neighbourhood structures, the method has an increased probability to detect more diversified routes than a single neighbourhood. The neighbourhood functions are called one by one in the order of swap, 2-opt and relocate. Once a neighbourhood function can no longer find a better set of routes, the next neighbourhood is called. If, however, a better solution (e.g. a more negative reduced cost) is found, the algorithm will restart from the first neighbourhood (i.e. swap).
Before VNS starts, the initial set of columns in z is generated by the insertion heuristic (Chen et al. (2013)). As shown in Algorithm 2, for each successive iteration, z is Table 3 Abbreviations of VNS z current solution R z a set of routes present in z i index of neighbourhood i max index of the last neighbourhood function c min minimum reduced cost of route set c min modified minimum reduced cost of route set maxIteration max number of column generation iterations maxColumns max number of routes columnP ool stores the set of best routes Algorithm 2 Pseudo-Code of VNS column generator Require: z, maxIteration 1: j ← 0 2: while j < maxIteration do 3: if c min < c min then 6: i ← 1, c min ← c min 7: columnP ool ← sortByReducedCost(R , maxColumns) 8: columnP ool ← columnP ool ∪ z 9: else 10: i ← i + 1 11: return columnP ool updated subsequently. Since our VNS not aims to solve the overall problem but find out a set of feasible routes with the most negative reduced costs to be solved by the RMP, the VNS based column generator is not conventionally implemented with a shaking process. The search is guided by the pricing methods described in Section 5.4. The neighbourhood(R,i,maxColumns) function applies the i-th neighbourhood function on all routes in R z to search for new feasible routes. It returns the maximum of maxColumns distinct routes with negative reduced cost. The constraints related with feasible route pattern (Eq. (3) to (7)) are imposed. Function minReducedCost(R ) returns the minimum reduced cost of route set R . Function sortByReducedCost(R z ∪ R , maxColumns) sorts routes in R z ∪ R by their corresponding reduced costs in an ascending order and returns the top maxColumns distinct routes. The RM P (columnP ool) is the restricted master problem (see Section 5.3) based on the route set stored in columnP ool and the solution is stored in set z.
Note that a distinctive feature of our VNS based column generation method is the joint exploitation of pricing information and the current best solution z. While most heuristic column generation methods aim to construct new routes from scratch in light of new pricing information, our VNS column generation procedure (and the GA column generator) is a perturbative based neighbourhood search starting from existing columns in the basis.
Consequently, we believe our column generation methods converge much faster than the constructive methods used in the literature.

Genetic algorithm
We also investigate a Genetic Algorithm (GA) approach to tackle the pricing subproblem. The motivations are two-fold: first, at each column generation iteration, we need to obtain a set of routes with the most negative reduced costs, which the VNS may struggle to achieve as a single point search method. The GA is potentially more powerful as it can find a population of routes through evolution. Secondly, we believe that high quality routes (i.e. most reduced costs) may share some common structures which could be evolved more efficiently through crossover operations in the genetic algorithm.
Therefore, each chromosome in our generic algorithm stands for a vehicle route, leading to a variable length chromosome. More specifically, a route(chromosome) is represented as a list of nodes(genes). For example, the parent 1 illustrated in Figure 4 simply represents route 0 → 1 → 2 → 2 → 1 → 0.
The pseudo-code for the GA search is given in algorithm 4. Similar to our VNS implementation, the initial population is generated by using the insertion heuristic by Chen et al. (2013). The size of the initial population for each RMP iteration is equal to the number of distinct vehicle routes used in the solution z but increased to a pre-defined value populationSize in the subsequent generations. Other implementation details of our GA are as follows. Two-point crossover operators were adopted. The length between two crossover points is randomly generated from 0 to 2 arcs, as larger crossover length would increase the possibility of generating infeasible routes due to the violation of routes' travel time constraint. Figures 4a, 4b and 4c illustrate examples of the two-point crossover. A standard mutation operator is used in which each chromosome is subject to an uniform 2-opt mutation with probability mutationRate. The 2-opt mutation operator is the same as the 2-opt neighbourhood moves in our VNS method. A local search stage is incorporated into our GA to ensure that local optima are reached in each generation. The local search is performed every time when new individuals have been generated. More specifically, the local search phase swaps two nodes between two different routes and returns two new routes that are local optimal with regard to the neighbourhood.
We use the tournament selection method. As the first population is obtained by the insertion heuristic, it usually has smaller population size than the predetermined constant value (populationSize). The tournament size is set to populationSize × tournamentRate so that it is population dependent. The fitnesses of individuals are calculated according to the functions in Section 5.4. Note that only feasible routes that satisfy the time constraints are considered and evaluated. If their fitnesses are better than any of the routes in the columnP ool which stores the set of best routes so far, they replace the inferior routes in the columnP ool, to allow a maximum of maxColumns columns to be stored.
Finally, the algorithm terminates when the number of RMP iterations reaches a predefined parameter,maxIterations. The pseudo-code of the proposed GA is given in Algorithm 4.
Note that although the main framework of our GA is the same as many other GA implementations, the goal is very different. Our GA here does not solve the overall problem, but rather evolves a set of vehicle routes (columns) with the most negative reduced costs.
These set of routes will then be used in solving the updated RMP problems. i ← i + 1 8: return z 6. Experiments with small R For the first round of experiments, we consider instances with relatively small R. As such, all instances in the first round of experiments have seven nodes, resulting in a total of 61365 feasible routes which is close to the limit to which our model can be solved directly.
Therefore, we can compare how our methods perform in comparison with exact methods. These instances were grouped into three sets. All the instances are generated by the same parameters except the size of the planning horizon. Five instances are generated for each problem set, referred to as I4, I6 and I8, standing for shift length of 4, 6 and 8, respectively.
The information and configuration of these problem sets is illustrated in Tables 4 and 5.
In order to test the efficiency of the column generation process in the first round of experiments, the initial route set is constructed by the simple method detailed in section updateColumnPool(r) 15: return R Table 4 Configuration of the artificial instances no. of nodes: 7 (including the depot) Commodity Time Window: 1-2 hours up to the length of planning horizon Commodity Availability Time: nearly 30% commodities are available at the start of the planning horizon Emergency tasks: 10% to 30% of total commodities (i.e. time window<10h) we permit our method to use more trucks than the limit (n), but this constraint will later be restored at the end of the column generation procedure. Gurobi 8 linear programming libraries were used in conjunction with Java 7.0. These experiments were run on a PC with an Intel i7 3.40GHZ processor and 16GB RAM.
The experimental results are given in Table 6. Since the pricing by enumeration method (see Section 5.4) takes an unrealistically long time even for the smallest instances (e.g.

3-4 hours for a 4-shift instance)
, it is not used for further experiments. Column T is the total running time of the entire process, from data parsing, solving, to the solution output.
Col. shows the total number of columns being generated during the process. Obj. gives the objective value which is the total travel distance. Hereafter, P1 and P2 are short abbreviations for column generation solution methods adopting P1 average pricing and P2 demand weighted average respectively (see Section 5.4).
Overall, the results in Table 6 show that most instances are solved in 1000s or less.
In most cases, P2 generated a larger number of columns than P1 during the column generation process. On average, P2 generates 1165 more columns than P1, resulting in longer running times, but P2 uses 3089km less distance than P1. Seemingly, this fact is due to P2 generating more columns that enlarge the search space used by the model. However, we notice that for the result of instances I4-1, I4-2 , I8-1 and I8-3, P1 obtained a larger number of columns which did not result in a smaller objective value.
The performance of both algorithms is also compared with the results from the Gurobi IP solver with the default algorithm setting in two experiments. The first experiment allows the solver to solve the problem to optimality and its objective value is denoted as Obj..
In the second experiment, Gurobi was given a limited computational time (the same time taken by the slowest of P1 and P2) and the corresponding objective value is marked as obj.*. All the results are given in Table 6.
It can be seen that although Gurobi can solve all instances to optimality, it takes more than 8000s on average and sometimes more than 10h. Two tailed paired t-tests (α = 0.01) were conducted to compare the performance between P1, P2 and Gurobi. In contrast, the proposed column generation methods (P1 and P2) use significantly less time (P1 vs.
for the remaining experiments only P2 is used. Similar to the previous section, maximum maxColumns = 1000 columns are allowed to be generated by both VNS and GA at each iteration. As maxColumns is the only parameter used in the proposed VNS based column generator, parameter tuning for VNS is omitted. We now illustrate parameter tuning for the GA.

Parameter Tuning for GA
The parameters used in the proposed GA are the population size (populationSize), the generation size (generations), the probability of mutation (mutationRate), and the tournament size i.e. the tournament rate (tournamentRate). In this experiment, the mutationRate is set to 0.02 and the tournamentRate is set to 0.1 after some initial tuning. Table 7 shows the results with the algorithm with different population sizes and different number of generations for the two most challenging problem instances LB8-1 and LB8-2. Each instance was run five times and the average result of both instances is given in column Avg.. The maxIteration is set to 5 as increasing it further gives very little further improvement. With the consideration of algorithm efficiency, we choose the combination of populationSize=500 and generations=500.  20676 20128 19999 21205 20228 20081 19879 19775 19724

Comparing VNS and GA based column generation approaches
Due to the very large amount of computational time required, we select a total of six instances, two from the real-life instance set and four from the artificial instance set used by Bai et al. (2015). The instance names starting with NP are real-life instances while those starting with LB, TB, LU and TU represent (Loose, Balanced), (Tight, Balanced), (Loose, Unbalanced) and (Tight, Unbalanced) configurations of artificial instances.
The first digit of each instance name indicates the length of the planning horizon (e.g. NP4-1 is a real-life instance with a 4-shift planning horizon).
The stopping criterion maxIteration is set to 10 for both VNS and GA. Figure  suggest that the GA is a faster converging algorithm thanks to its population based search framework and capacity to evolve a set of routes instead of a single one.
As experiments show that the VNS converges in 30 RMP iterations, for a fairer comparison, we set maxIteration to 30 for both GA and VNS for all instances. In addition, a comparison was also made with previous results obtained by the three-stage method, meta-heuristic methods and lower bound reported in Bai et al. (2015).  The novel solution coding and pricing methods limit the search space for the algorithm, so its efficiency is increased compared with the results obtained by meta-heuristics (rsVNS and SAHH ). The 3-stage method performs well for tight instances, but it does less well for large and loose instances. The reason is that it employs an integer programming solver so its solution time increases exponentially with large problem sizes. The proposed column generation methods are able to find effective columns in order to reduce the problem size, therefore, compared with the 3-stage method, the solving time of column generation method is significantly decreased for large instances. However, the advantage of column generation method may not be obvious for small problem instances (i.e. tight and small instances) as the iterative RMP solving comprises a significant proportions of the run time by the algorithm.

Conclusions
We have presented an innovative FTL routing formulation assisted by dynamic cuts and investigated column generation based approaches which are particularly effective on very large instances. Unlike traditional branch-price-and-cut, it performs an incomplete search, with the aim of finding good solutions more quickly. It efficiently solves the problem using the following strategies: 1) Infeasible flow assignments are allowed in the column generation process but will be fixed by adding cuts in the end; 2) To reduce the number of decision