A Hybrid Pricing and Cutting Approach for the Multi-Shift Full Truckload Vehicle Routing Problem
AA Hybrid Pricing and Cutting Approach for theMulti-Shift Full Truckload Vehicle Routing Problem
Ning Xue
University of Nottingham UK, [email protected],
Ruibin Bai*
University of Nottingham Ningbo China, [email protected],
Rong Qu
University of Nottingham UK, [email protected],
Uwe Aickelin
The University of Melbourne, [email protected]
Full truckload transportation (FTL) in the form of freight containers represents one of the most importanttransportation modes in international trade. Due to large volume and scale, in FTL, delivery time is often lesscritical but cost and service quality are crucial. Therefore, efficiently solving large scale multiple shift FTLproblems 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 hybridisingpricing 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 cutsare dynamically generated to eliminate infeasible flow assignments caused by incompatible commodities.Computational experiments on real-life and artificial benchmark FTL problems showed superior performanceboth 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 canefficiently solve large scale real-life FTL problems.
Key words : full truckload transport, column generation, pricing and cutting, metaheuristics
1. Introduction
In intermodal freight transportation, a large proportion of container transportation is car-ried out by barges, trains or ocean-going vessels Braekers et al. (2014). Container movementactivities between intermodal terminals, depots and shippers are also referred to as drayageoperations and such activities are usually performed by trucks. Although drayage opera-tions represent a small fraction of the total distance of an intermodal freight transportation, * Corresponding author a r X i v : . [ c s . A I] D ec hey constitute a substantial share of the shipping costs Smilowitz (2006). Consequently,major port are facing intense competition and pressure to improve the efficiency of drayageoperations.Due to labour laws and other constraints related to the working time of drivers, shiftbased working schedules are becoming a common practice in the transportation indus-try (e.g. taxis and buses). The problem studied in this paper concerns the movement ofcontainers (commodities) between a number of terminals (docks) within a short distancelocated in a large international port using a homogeneous truck fleet. The transportationtime window of commodities usually spans from a few hours up to several days, and coversof multiple working shifts. Thus the problem that we are addressing is essentially differentfrom the single-shift problem studied in most of the existing full truckload routing prob-lems (e.g. Zhang et al. (2010), Braekers et al. (2014)) because the planning horizon coversseveral shifts and determining the transportation shift of each commodity forms part ofthe optimisation decision.In our earlier study Bai et al. (2015), a set covering model and a three-stage method wereproposed for this problem. However, the computational time to optimally solve large sizeproblems was prohibitive. Chen et al. (2013) investigate a reactive shaking variable neigh-bourhood search (rsVNS) and a simulated annealing hyper-heuristic method (SAHH) Chen(2016) for this problem. The rsVNS extends the original VNS which utilises the systematicchanges of multiple neighbourhood functions to achieve convergence and diversification.The SAHH applies a reinforcement learning based neighbourhood selection mechanismwithin a simulated annealing framework. The learning mechanism aims to adapt the algo-rithm to different problem instances and search scenarios by dynamically adjusting theneighbourhood selection strategies. Both rsVNS and SAHH were able to obtain feasiblebut 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 advanta-geous features of a previously proposed indirect solution encoding scheme, leading to someinsightful findings of the differences between the multi-shift FTL problems and traditionalpickup and delivery problems; 2) A pricing based column generation method is investi-gated, in conjunction with dynamic cuts. Our method is inspired by branch-price-and-cutalgorithms but differs in three major ways. Firstly, we do not solve the pricing subproblem xactly. Instead, we quickly produce approximate solutions by introducing several optimi-sation strategies (See Section 5.4). Secondly, we employ metaheuristics instead of traversingthe entire search tree. Thirdly, the cuts are added after the column generation processbecause the infeasible flow assignments (See Section 4.2) rarely occurs and they take timeto be evaluated and fixed in every column generation iteration. The new algorithm cansignificantly speed up the solution time for large size problems. Moreover, the solutionquality is also improved. Two pricing methods are proposed and tested on both real-lifeand artificial instances.The remainder of this paper is organised as follows: the problem is described in detail inSection 2; a literature review is given in Section 3 followed with the mathematical modelof the problem in Section 4. The proposed pricing and cutting method is illustrated inSection 5 and its the computational experiments are presented in Sections 6 and 7. Finallyconclusions are drawn in Section 8.
2. 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 theoperational time windows of the freights and the planning horizon can span across severalshifts. Although each container is transported in a single shift, its time window coveringmulti-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 whilesatisfying various constraints.First, each full truckload commodity (container) has an available time for pickup anda deadline for delivery. Second, during each shift, a number of unit-capacity trucks startfrom the deport at the start of the shift, complete a number of transportation requests andthen return to the depot before a shift ends. Finally, a service time is applied during bothpickup and delivery. To clarify, we refer to the request of a full truckload movement asone unit of a commodity. A commodity is a collection of requests of full truckload freightsthat share identical sources, destinations and time windows.In the context of real-world applications that this research tries to address, the totalquantities of all the requests within a planning horizon can be very large (more than1000). The number of terminals is relatively small (less than 10) and the distances between hese terminals are relatively short (all reachable within a shift). These features make thisproblem different from problems that are studied in previous work. It has been shownin Bai et al. (2015) that a model based on a set covering formulation is more promisingfor this problem than other node based formulations. However, the features and reasonshave not yet been sufficiently analysed. For the completeness of this paper, we include themodel 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 thenext section.
3. Literature Review
The drayage operations problem is a typical case of bidirectional multi-shift full truckloadvehicle routing problems. Bai et al. (2015) highlight the core features of these truckloadvehicle routing problems and discuss relationships with other variants of vehicle routingproblems (VRP) from three aspects, including the directions of the flow, existence of con-solidation or not, and length of the planning horizon. Here, we summaries the relevantresearch on the drayage operation problems which we broadly classify into drayage opera-tions with and without relocation requirements 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) extendthe m-TSPTW model with social constraints and propose an exact algorithm based ondynamic programming. Moreover, a hybrid method combining dynamic programming andgenetic 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 trans-portation problems. The basic problem is formulated as an m-TSPTW problem, which issolved by an insertion heuristic. Gendreau et al. (2015) refer to this routing problem asthe one-commodity Full-Truckload Pickup-and-Delivery Problem (1-FTPDP) and presentthree mathematical formulations with branch-and-cut algorithms to optimally solve themodel formulations. Lai et al. (2013) propose a new routing problem that can be viewedas a vehicle routing problem with clustered backhauls (VRPCB). Solutions are obtainedwith the Clarke-and-Wright algorithm and improved further by a neighbourhood based etaheuristic . This work is also compared in the study of a problem with single and dou-ble container loads Ghezelsoflu et al. (2018). The distribution of more-than-one containerper truck by different types of trucks has also been studied in Vidovi´c et al. (2017) andFunke and Kopfer (2016). Soares et al. (2019) study an FTL problem with multiple typesof vehicle synchronisations. A MIP model and a heuristic solution method based on thefix-and-optimise principles are proposed.
Efforts to combine the planning of loaded and empty container transports are made byseveral authors. Coslovich et al. (2006) analyse a fleet management problem for a containertransportation company by decomposing the problem into three subproblems, which arethen solved using a Lagrangian relaxation. Ileri et al. (2006) present a column generationbased approach for solving a daily drayage problem. Smilowitz (2006) model a drayageoperation with empty repositioning choices as a multi-resource routing problem (MRRP)with flexible tasks. The solution approach is a column generation algorithm embeddedin a branch-and-bound framework. Imai et al. (2007) formulate a container transporta-tion problem as a vehicle routing problem with full container loads (VRPFC) and solve itwith 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 isproposed. The work is further extended by using a deterministic annealing algorithm sug-gested in Caris and Janssens (2010). Zhang et al. (2010) improve the time window par-titioning scheme used in Xiubin Wang (2002) for container transportation in a local areawith multiple depots and multiple terminals. The results indicate that good performancecan be achieved compared with a reactive tabu search (RTS) method demonstrated inRuiyou Zhang (2009). Zhang et al. (2011) also investigate the single depot and terminalproblem. Again, an RTS is proposed. Vidovic et al. (2011) extend the problem analysedby Zhang et al. (2010) and Imai et al. (2007) to the multi-commodity case and formulateit as a multiple matching problem. Solutions are obtained via a heuristic approach basedon calculating utilities of matching nodes. Nossack and Pesch (2013) present a new formu-lation for the truck scheduling problem based on a FT-PDPTW and propose a two-stageheuristic solution approach. Braekers et al. (2013) investigate a sequential and an inte-grated approach to plan loaded and empty container drayage operations. A single- anda two-phase deterministic annealing algorithm are presented. This solution approach is urther adapted in Braekers et al. (2011) to take a bi-objective optimisation function intoaccount. The algorithms are further improved in Braekers et al. (2014). More recently, Xieet al. (2017) investigate the empty container delivery problem in an intermodal transportsystem composed of a sea liner firm and a rail firm. Apart from transportation cost, thedifference in marginal profits between the seaport and dry port is also considered in theobjective function.Some researchers examine drayage operations problems in dynamical situations. A surveyon 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 drayageproblem as some forms of classical vehicle routing problems in order to exploit the timeconstraint 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 constraintsare not very tight and node-based solution representations generally lead to unnecessarilylarge search space, resulting to inefficient solution methods.
This paper studies an indirect solution representation for the multi-shift FTL problem thataddresses these issues and contributes to the body of research with an efficient columngeneration method. In many vehicle routing applications solved by column generation, thesubproblem is usually viewed as an elementary shortest path problem with resource con-straints or one of its variants. Nowadays, an increasing number of hybridisations betweenheuristics and exact approaches are developed. These methods can provide a good com-promise between solution quality and computational time as they adopt the advantages ofboth types of methods. Puchinger and Raidl (2005) classified hybridisation of exact algo-rithms and (meta)-heuristics into four types, we briefly introduce them and give examplesfor each:
1) Collaborative Combinations - sequential execution : In this type of hybridisa-tion, 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 columnsand the exact method is used to find an optimal solution from the feasible columns. Exam-ples of this type of hybridisation can be found in Clements et al. (1997) and Vasquez et al.(2001). ) Collaborative Combinations - parallel or intertwined execution : Instead ofexecuting either heuristics or exact methods sequentially, this type of method implementsthe algorithms in a parallel or intertwined way. Clusters or multi-processors are used todeploy the parallel implementations. There are several frameworks proposed to facilitatesuch 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, thesolution of LP-relaxation and its dual values can be utilised in heuristically guiding neigh-bourhood search. Applications can be found in Marino et al. (1999) and Puchinger et al.(2004).
4) Integrative Combinations - incorporating heuristics in exact algorithms :This type of hybridisation is analogous with the previous one, but heuristics are embed-ded within exact algorithms. For example, heuristics can be used to determine boundsin branch-and-bound algorithms. Heuristics can also be used to search for columns withnegative costs in the branch-and-price approach. Applications of this hybridisation methodcan be found in Puchinger and Raidl (2004) and Strandmark et al. (2020). The columngeneration based method proposed in this paper falls into this category.Please refer to Blum et al. (2011) and Muthuraman and Venkatesan (2017) for morecomprehensive reviews of the hybridisation approach.
4. 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 oftruckload commodities and the travel time of arc ( i, j ) is denoted as µ ij . All trucks mustdepart from and return to node 0 (depot). Let R be the set of all feasible routes that atruck can execute in a working shift without the complication of time window requirementsfrom 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 withina time window ( e sr i , l sr i ) where e sr i is the earliest time that a truck covering route r can start igure 1 Example of a routing sharing among five commodities. a pickup or delivery operation while l sr i is the latest time that a truck must depart from thenode. Let t r i be the service time at node r i . In each route r encoding, a duplicated node isinserted if the node involves both a loading and an unloading operations (i.e. this node isboth the destination and source of two different commodity flows). Therefore, if the nodesis 0-indexed in a route, loading operations are always conducted at the odd indexed nodesof a route (see Eq. (3)) and unloading operations are at the even-indexed nodes. e sr i and l sr i can be calculated using the following recursive equations: e sr i = e sr i − + t r i − + µ r i − r i ∀ i ∈ r, r ∈ R (1) l sr i = l sr i +1 − t r i +1 − µ r i r i +1 ∀ i ∈ r, r ∈ R (2)Let K denote the set of commodities for delivery. Each commodity k ∈ K is defined by atuple ( o ( k ) , d ( k ) , Q ( k ) , σ ( k ) , τ ( k )), which, respectively, are the origin, destination, quantity,availability time and deadline of commodity k . Denote δ ksr i the binary variable to indicatewhether the i -th node of route r can be the loading node for commodity k in shift s ( δ ksr i = 1)or not ( δ ksr i = 0). To speed up the computation, δ ksr i can be pre-calculated by checking thefollowing conditions: i mod r i = o ( k ) (4) r i +1 = d ( k ) (5) l sr i ≥ σ ( k ) + t r i (6) e sr i +1 ≤ τ ( k ) (7) igure 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 whileeven numbered nodes are unloading nodes. If a node on a route is involved with bothloading and unloading, a copy of it is created so that the above rules are maintained (seemore details in Bai et al. (2015)). In this example, a truck departs from the depot andpicks up a unit of commodity from either commodity 1, commodity 2 or commodity 3from node 1 and unload the commodity at node 3. Then the truck picks another unit ofcommodity (either commodity 4 or commodity 5) at node 4 and unload at node 2 beforethe 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 ksr i : Commodity flow of the i th node of r in s for servicing commodity k ∈ K . • y sr : The number of times a given route r ∈ R is used during shift s ∈ S and y sr ∈ N + .The model for this multi-shift FTL problem can be formulated as the follows:min (cid:88) s (cid:88) r d r y sr (8)subject to (cid:88) r y sr ≤ n ∀ s ∈ S (9) (cid:88) s (cid:88) r (cid:88) i x ksr i = Q ( k ) ∀ k ∈ K (10) (cid:88) k x ksr i ≤ y sr ∀ i ∈ r, ∀ r ∈ R, ∀ s ∈ S (11) x ksr i ≤ δ ksr i y sr ∀ i ∈ r, ∀ r ∈ R, ∀ k ∈ K, ∀ s ∈ S (12) x ksr i = Z + ∀ i ∈ r, ∀ r ∈ R, ∀ k ∈ K, ∀ s ∈ S (13) y sr ∈ Z + ∀ r ∈ R, ∀ s ∈ S (14) he objective function (8) minimises the aggregated distance of all routes used in asolution. 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 sr times. Constraint(12) makes sure that any positive x ksr i is feasible in terms of the source, destination andoperation time window of commodity k . Since binary indicator δ ksr i can be pre-calculated,this constraint can be eliminated by removing the corresponding flow variables x ksr i fromthe model when δ ksr i takes value of 0. This is indeed how the model was implemented in ourexperiments because the resulting model is a lot smaller. Constraints (13) and (14) definethe domains of the decision variables. One of the most helpful benefits of this solution encoding scheme is the transformation of aprevious m-TSPTW based non-linear model (e.g. the model proposed by Chen (2016)) intoa 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 theshift-independent feasible truck route set.For some applications (e.g. FTL with a small number of terminals), pre-computing allfeasible routes is possible since the time related constraints in this problem are slightlydifferent 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 win-dow ( σ ( k ) , τ ( k )) defining its availability time and the delivery deadline. Time constraintsrequire that both the pickup and delivery operations occur within this time window forcommodity k . In PDPTW problems, two separate time windows are used, one for pickupand 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, althoughwe acknowledge it is very different for express deliveries which are mostly for householdcustomers.A second benefit of this solution representation is its capability to handle nonlinearcost functions. For example, the costs of routes could be a nonlinear, complex function ofthe distance. It also permits to include various other constraints related to drivers (e.g.maximum driving distance, time or preferred terminals).
10 of 34 third benefit of this solution representation is the reduced size of the search spacecompared with a commonly used m-TSPTW formulation, in which each container is mod-elled as a node in the graph. For a problem instance with hundreds or even thousands oftruckload sized containers, the corresponding graph in m-TSPTW formulation could beprohibitively expensive to handle. However, in the real-world problem that we are con-cerned with, containers often arrive in large batches with same requirements (i.e. sameO-D pairs and time windows). In an m-TSPTW formulation, any swaps of positions ofthese 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 signif-icantly larger search space with a plateau. In our proposed formulations, containers withthe same property are grouped as one commodity, leading to a one-to-one mapping anda much smaller search space.
Although for all the practical instances that we extracted from real-world problems, theFTL model in Section 4 produces solutions that satisfy practical constraints. However, itis possible to artificially generate problem instances that the proposed FTL model returnsan infeasible solution. That is, the solution is feasible for the FTL model but may stillviolate 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 exampleof such cases is illustrated in Figure 2. i i +1 j j +1 00 !" direct arc indirect arcs - . / . / !" - . . Figure 2 Service time push back and non-compatible commodities.
In this example, we included two non-compatible commodities ( k and v ) that are ser-viced using a feasible route r at arcs ( i, i + 1) and ( j, j + 1), respectively. The solution
11 of 34 ecomes infeasible when the service time of the first commodity k was delayed because thecommodity available time σ ( k ) is later than node i ’s earliest arrival time e r i . Because ofthis, the service time of commodity k at node i is pushed back by push back ( i ) = σ ( k ) − e r i (15)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 allcommodities 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 theroute. A violation of another following commodity v ’s shipment deadline ( τ ( v )) shall occurif the following condition is satisfied: push back ( j + 1) > e r j +1 − τ ( v ) (16)That is, if a push back caused by a previous commodity is greater than the differencebetween 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 followingconstraints should be added to ensure v is not inserted at or after k in the same route andshift. x ksr i ≤ M θ and x vsr j ≤ M (1 − θ ) i ≤ j ∈ r, ∀ k ∈ K r , ∀ v ∈ V r , ∀ s ∈ S (17)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 assignmentsat a same node. The process of generating the cuts through constraint (17) to eliminatenon-compatible commodity assignments is given in Algorithm 1.We do not want to strongly restrict the non-compatible commodities, as shown in theabove example, v is not simply forbidden to be served by the node, instead, it is still
12 of 34 llowed to be served by the route as long as the compatibility constraint of k and v isnot violated. From Algorithm 1 it can be seen that the procedure keeps tracking the pushback 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 onlyif any pairs of incompatible commodities were found. That means even if the service timeof the commodities served between k and v , if any, are pushed back, they are not restrictedby constraint (17) unless there are larger push back by other commodities result in a delayin shipment. Algorithm 1
Valid cuts generation for eliminating infeasible flow assignments
Require: r ∈ R x , where R x is the set of routes used in the current solution and x is the vector offlow variables x ksr i K r = ∅ , V r = ∅ , accu = 0 (cid:46) accu is the accumulated push back time along r for i in r do push back( i ) ← if ( i mod 2=1) then for k in W ( i ) do , where W ( i ) is set of commodities serviced by node i . if σ ( k ) > e r i then k .push back ← σ ( k ) − e r i − accu K r .add( k ) push back( i ) ← max( σ ( k ) − e r i , push back( i )) if push back( i ) > then Propagate puch back( i ) to all the remaining nodes in r accu + =puch back( i ) for j in r do if ( j mod 2=1) then for v in W ( j ) do , where where W ( j ) is set of commodities serviced by node j . if τ ( v ) < e r j +1 then v .acceptable push back ← e r j +1 − τ ( v ) V r .add( v ) for k in K r do for v in V r do if k .push back time > v .acceptable push back then k and v not compatible in r (cid:46) output constraint: x ksr i ≤ M θ and x vsr j ≤ M (1 − θ ) Table 1 An example: A problem with 4 commodities
Commodity Available Deadline Start node End node k k v v
13 of 34 able 1 gives a simple illustrative example how the valid constraints can be dynamicallygenerated into the model to prune the search space. In this example, A total of 4 feasibleroutes are available for selection to deliver 4 commodity k , k , v , v . 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. Withoutconsidering push backs, the optimal solution is to choose route 1 twice (i.e. y s = 2), withflows of x k s = 1 , x k s = 1 , x v s = 1 , x v s = 1 because it satisfies formulas (4) to (7) andrequires the least travel distance (158) to delivery all commodities. However, since σ ( k ) >e (13 : 40 > k at node 1 is pushed back to 13:40 by 325minutes, which is propagated to all the remaining nodes in route 1 (the updated e and l after pushed back by k are denoted as e (cid:48) and l (cid:48) in Table 2).It can be seen that the push back at node 1 by commodity k resulted in commodi-ties v or v not being serviced according to the optimal solution due to τ ( k ) < e (13 :10 <
15 : 00) , τ ( v ) < e (16 : 00 <
16 : 50) , τ ( v ) < e (15 : 50 <
16 : 50). Thus, K = { k } , V = { k , v , v } and 3 constraints ( x k s ≤ M θ and x k s ≤ M (1 − θ ), x k s ≤ M θ and x v s ≤ M (1 − θ ), x k s ≤ M θ and x v s ≤ M (1 − θ )) are subsequently added to the model. The trueoptimal objective, after adding the valid constraints, increased to 208 by choosing route 1to deliver k , v and route 2 to deliver v , k (i.e. x k s = 1 , x v s = 1 , y s = 1, x v s = 1 , x k s =1 , y s = 1). R The proposed model also has some problems. The most critical one is the size of the feasibleroute set R which can increase exponentially with the number of nodes (or terminals). InBai et al. (2015), some real-life problems have certain special features to permit some ofnodes 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, themethod becomes invalid for problems that do not possess these features to permit nodemerging.In this paper we propose to use a column generation method to address this issue. Theidea is to use the pricing information to guide the generation of promising feasible routesdynamically.
5. A Hybrid Column Generation Method
Column Generation is an effective approach for solving large scale integer programmingproblems (i.e. problems with large number of columns). It is a potentially very good method
14 of 34 able 2 An example: Time window e and l of nodes Node of route 1 0 1 2 3 4 0Index of route 1 0 1 2 3 4 5 e l e (cid:48) l (cid:48) e l e (cid:48) l (cid:48) e l e (cid:48) l (cid:48) e l e (cid:48) l (cid:48) e , l : before push back; e (cid:48) , l (cid:48) : after push back 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 verysmall subset of it. We propose to use the column generation approach for this problemin which the sub-problem (pricing problem) is solved to identify the variables that shouldenter the basis. The integer programming formulation presented in section 4 is also referred to as the masterproblem. The Restricted Master Problem (RMP) is the master problem that considers onlyof a subset of truck routes R that are generated by the pricing problem (subproblem) usingthe dual information obtained from the Linear Programming Relaxation (LPR) of theRMP. 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 availableand an initial truck routes set (see Section 5.2) is thus required to start the process. Thenthe LPR is solved to optimality and the dual information is obtained for calculating the
15 of 34 educed costs of routes during the pricing subproblem. The overall solution framework isoutlined in Figure 3, followed by detailed steps of the procedure.Our intial experiments showed that majority of computing time is consumed by the RMPsolving. The algorithm was thus modified to accelerate convergence by adding multipleroutes with negative reduced costs at each iteration. Details of the methods for the pricingproblem is given in Section 5.4. It is hoped that by doing this the total number of RMPcalls 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 ksr i and y sr areset back to their original ones during the final RMP solving.Because the pricing problem is solved repeatedly in the column generation framework, itis 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 andone for problems with a large R . For the former case, we propose to adapt an explicitenumerative generation of R as priori and then try to solve the pricing subproblem whenno column with negative reduced cost can be found. We apply a recursive algorithm togenerate all feasible routes as described in Bai et al. (2015) before the iterative procedurestarts. In the case of a large R , we propose to use heuristic approaches (see 5.5) to solve thepricing problem and the stopping criteria of the heuristic is an limitation of the number ofRMP cycles (denoted by Finish in Figure 3). The overall solution framework as describedabove is outlined in Figure 3.
Before the RMP is solved for the first time, no dual information is available and an initialset of columns is required to start the process. We apply two methods described in detailin the next two subsections to generate an initial set of columns (routes).
A prerequisite of constructing a basic route set is toensure that each commodity has at least one route to service it. Thus the simplest solutionis to generate a dedicated route for each commodity, in which an empty truck leaves thedepot and travels to the source of a commodity, loads the commodity and delivers it to itsdestination. After that, the truck returns to the depot. This method works fine in somecases but may of course lead to an infeasible solution in terms of the maximum number ofvehicles constraint (9).
16 of 34 igure 3 Framework of the proposed column generation process
The advantage of the basic route initialisation inthe previous section is its simplicity. However, very rarely will these routes be used inthe optimal solutions, neither do they resemble any of the routes that are present in theoptimal solutions. In this study, we proposes to use constructive heuristic methods togenerate these initial routes for our column generation method. In particular, we used thesame insertion heuristics described in (Chen et al. 2013). To construct routes, the taskthat 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 feasiblesolution for the original problem, the vehicle routes extracted from the solution shall alsoproduce a feasible solution in our column generation method, satisfying the maximumnumber of vehicles constraint (constraint (9)). Second, because the pricing subproblem issolved heuristically, starting from a good set of initial vehicle routes will enable the columngeneration method to generate high quality solutions more quickly compared to the simpleroute initialisation method. The proposed method will converge to a high quality solutionmuch faster.
17 of 34 .3. Linear programming relaxation (LPR)
Linear programming relaxation (LPR) relaxes the discrete variables constraints and differsfrom the model presented in Section 4 in the route set, which should be attained as a R (cid:48) ( R (cid:48) ⊂ R ) resultant from the pricing subproblem instead of all feasible routes R . LetLPR be the relaxed model, and Let LPR-(9-12) be the constraints corresponding to theconstraints (9-12) of the master problem. We present three different route price estimation methods: The first method obtainssolution by enumerating and examining all possible commodity assignments for each routeand the best route (i.e. the route with the most negative reduced cost) is selected. Thismethod is referred to as
Pricing problem by enumeration in this paper. For efficiency,two other pricing estimation methods (
Average pricing and
Demand weightedaverage pricing ) are also investigated.
Pricing problem by enumeration . Let α s , π k , β sr i , γ ksr i be the pricing variables for con-straints from the LPR-(9-12), respectively. The reduced cost for route r in shift s is: d r + α s − (cid:88) r i β sr i − (cid:88) r i (cid:88) k δ ksr i γ ksr i (18)However, since routes are generated independent of shifts, the following average reducedcost is computed over all shifts for each route. d r + 1 | S | ( (cid:88) s α s − (cid:88) s (cid:88) r i β sr i − (cid:88) s (cid:88) r i (cid:88) k δ ksr i γ ksr i ) (19)where | S | is the total number of shifts in the planning. Let W = { w , w , ..., } be the setof all possible commodity assignments, each of which can be delivered by one instanceof route r . In the example of the route in Figure 1, all possible commodity assignmentsare W = { [1 , , [1 , , [2 , , [2 , , [3 , , [3 , , ... } . Since there are many possible commodityassignments 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 processis repeated for the next route r + 1 until all feasible routes are evaluated. This searchingprocess guarantee that the reduce cost of commodity assignment in each route is examined
18 of 34 ut its efficiency is low as some routes may contain thousands of possible commodityassignments.In order to obtain good results in a reasonable time, we investigated the following stepsto improve efficiency: Firstly, the constraintsf LPR-(12) are pre-processed offline. Givena route and its start time, the feasibility of commodity in a route can be determined byformula (3)-(7) offline. This allows us to reduce a large number of decision variables thathave to be handled by the model. Consequently, we lost the price values ( γ ) associatedwith 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 thefeasible route set is meant to be same across all shifts. This has benefits from managementstandpoints too because drivers proficiency can be improved if they are asked to follow afixed set of routes repeatedly. Also after long run, the set of frequently used routes canbecome part of the knowledge system of the transportation planning and time consumingcolumn generation procedure may not be required anymore. Therefore, the price valuesof arcs ( β ) in each shift is not used because the efficiency is substantially degraded bygenerating routes dependent of shifts. Fortunately, the price of an arc can be estimated bythe 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 costcalculation, 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 approximatedpricing methods, illustrated below. P1: average pricing
Instead of enumerating all the commodity combinations of a routeand then checking the c r for each of w i , a more efficient approach is to use the averageprices to estimate c r approximately. More specifically, let J be the set of all service nodesin r (e.g. nodes { , } in Figure 1). Denote V (cid:48) j be the set of all commodities that can beserviced by a node j in r . The reduced cost c r for route r is calculated by the followingequation: c r = d r − (cid:88) j ∈ J ( 1 | V (cid:48) j | (cid:88) k ∈ V (cid:48) j π k ) (20) P2: demand weighted average pricing
Though the commodities processed by a servicenode in a route share the same source and destination node, the quantity of the commodi-ties varies from one to another. The simple average pricing method P1 fails to take into
19 of 34 ccount the quantity of the commodities, so that large quantity commodities may be leftunpaired to improve the efficiency. Therefore, the demand weighted average method triesto give priority to large commodities at the early stage. A weight ω k that is proportionalto the commodity quantity Q ( k ) is used. The weighted average pricing method uses thefollowing equation to estimate the reduced cost. c r = d r − (cid:88) j ∈ J (cid:88) k ∈ V (cid:48) j ω k π k (21) R As can be seen from Figure 3, a heuristic column generator is used within the columngeneration framework. As optimally solving the pricing problem involves an expensiverecursive tree search, we propose to use a variable neighbourhood search (VNS) and agenetic algorithm (GA) to tackle the pricing subproblem. The goal of the metaheuristics isto identify new columns with negative reduced costs. The idea is that, instead of generatinga new column (i.e. route) from scratch, it is probably more efficient to search from theexisting routes through either neighbourhood moves or route combinations (i.e. crossovers).VNS and GA are widely adopted excellent frameworks to implement these ideas. The maindifference here is that the metaheuristics are guided by an objective function that heavilyrelies on the pricing information obtained from the linear program relaxations.
The pseudo-code of our VNS algorithm is given in Algorithm 2 and theparameters of the algorithms are listed in Table 3. In our VNS method, the neighbourhoodfunctions include swap , , and relocate . These operators are very similar to those usedin solving the classical VRP problems. For example, the swap operator swaps two arcsof two different routes. The exchanges two nodes on the same route. The relocate operator relocates an arc from its current route to a different one. By exploring differentneighbourhood structures, the method has an increased probability to detect more diver-sified routes than a single neighbourhood. The neighbourhood functions are called oneby one in the order of swap , and relocate . Once a neighbourhood function can nolonger find a better set of routes, the next neighbourhood is called. If, however, a bettersolution (e.g. a more negative reduced cost) is found, the algorithm will restart from thefirst neighbourhood (i.e. swap ).Before VNS starts, the initial set of columns in z is generated by the insertion heuris-tic (Chen et al. (2013)). As shown in Algorithm 2, for each successive iteration, z is
20 of 34 able 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 (cid:48) 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 j ← while j < maxIteration do columnP ool ← V N S ( z ) (cid:46) Algorithm 3 z ← RM P ( columnP ool ), j ← j + 1 return z Algorithm 3
Pseudo-Code of VNS()
Require: z , i max , maxColumns i ←
1, update R z by z , c min ← while i ≤ i max do R (cid:48) ← neighbourhood ( R z , i, maxColumns ) c (cid:48) min ← minReducedCost ( R (cid:48) ) if c (cid:48) min < c min then i ← c min ← c (cid:48) min columnP ool ← sortByReducedCost ( R (cid:48) , maxColumns ) columnP ool ← columnP ool ∪ z else i ← i + 1 return columnP ool updated subsequently. Since our VNS not aims to solve the overall problem but findout a set of feasible routes with the most negative reduced costs to be solved by theRMP, the VNS based column generator is not conventionally implemented with a shak-ing 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 allroutes 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 pat-tern (Eq. (3) to (7)) are imposed. Function minReducedCost ( R (cid:48) ) returns the minimumreduced cost of route set R (cid:48) . Function sortByReducedCost ( R z ∪ R (cid:48) , maxColumns ) sortsroutes in R z ∪ R (cid:48) by their corresponding reduced costs in an ascending order and returnsthe top maxColumns distinct routes. The RM P ( columnP ool ) is the restricted master
21 of 34 roblem (see Section 5.3) based on the route set stored in columnP ool and the solution isstored in set z .Note that a distinctive feature of our VNS based column generation method is the jointexploitation of pricing information and the current best solution z . While most heuristiccolumn generation methods aim to construct new routes from scratch in light of newpricing 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 theconstructive methods used in the literature. We also investigate a Genetic Algorithm (GA) approach totackle the pricing subproblem. The motivations are two-fold: first, at each column genera-tion iteration, we need to obtain a set of routes with the most negative reduced costs, whichthe VNS may struggle to achieve as a single point search method. The GA is potentiallymore powerful as it can find a population of routes through evolution. Secondly, we believethat high quality routes (i.e. most reduced costs) may share some common structures whichcould be evolved more efficiently through crossover operations in the genetic algorithm.Therefore, each chromosome in our generic algorithm stands for a vehicle route, leadingto a variable length chromosome. More specifically, a route(chromosome) is represented asa list of nodes(genes). For example, the parent 1 illustrated in Figure 4 simply representsroute 0 → → → → → z but increased to a pre-definedvalue populationSize in the subsequent generations. Other implementation details of ourGA are as follows. Two-point crossover operators were adopted. The length between twocrossover points is randomly generated from 0 to 2 arcs, as larger crossover length wouldincrease the possibility of generating infeasible routes due to the violation of routes’ traveltime constraint. Figures 4a, 4b and 4c illustrate examples of the two-point crossover. Astandard mutation operator is used in which each chromosome is subject to an uniform mutation with probability mutationRate . The 2-opt mutation operator is the sameas the 2-opt neighbourhood moves in our VNS method. A local search stage is incorporated
22 of 34 nto our GA to ensure that local optima are reached in each generation. The local searchis performed every time when new individuals have been generated. More specifically, thelocal search phase swaps two nodes between two different routes and returns two newroutes that are local optimal with regard to the neighbourhood.We use the tournament selection method. As the first population is obtained by theinsertion heuristic, it usually has smaller population size than the predetermined constantvalue ( populationSize ). The tournament size is set to populationSize × tournamentRate so that it is population dependent. The fitnesses of individuals are calculated accordingto the functions in Section 5.4. Note that only feasible routes that satisfy the time con-straints are considered and evaluated. If their fitnesses are better than any of the routesin the columnP ool which stores the set of best routes so far, they replace the inferiorroutes 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 predefinedparameter, 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 GAimplementations, 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. Algorithm 4
Pseudo-code of the GA column generator
Require: maxIterations , z , generations , populationSize , columnP ool , maxColumns while i < maxIterations do R ← z , Clear columnP ool while j < generations do R ← generateN ewP olulation ( z , R, populationSize ) (cid:46) Algorithm5 j ← j + 1 z ← RM P ( columnP ool ) i ← i + 1 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 of61365 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.
23 of 34 a) Example of crossover 1 (b) Example of crossover 2(c) Example of crossover 3
Figure 4 Examples of crossover operators in the GA
A set of randomly generated instances are used in the experiments. These instances aregenerated based on characteristics of real-life instances which are obtained from historicallyscheduled container operation data of a truck company. All artificially generated instanceshave three planning horizons of 4, 6, 8, reflecting the different problem scenarios in practice.These instances were grouped into three sets. All the instances are generated by the sameparameters except the size of the planning horizon. Five instances are generated for eachproblem 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 ofexperiments, the initial route set is constructed by the simple method detailed in section
24 of 34 lgorithm 5
Pseudo-Code of generateN ewP olulation () Require: current solution z , current population R , empty new population R (cid:48) , populationSize , tournamentRate , mutationRate for r in z do if r not in R then R.add ( r ) (cid:46) Ensure R include z while R (cid:48) .size < populationSize do r ← tournamentSelection( populationSize, tournamentRate, R ) r ← tournamentSelection( populationSize, tournamentRate, R ) R (cid:48)(cid:48) ← crossover( r , r ) R (cid:48)(cid:48) ← mutation( R (cid:48)(cid:48) , mutationRate ) R (cid:48)(cid:48) ← localSearch( R (cid:48)(cid:48) ) for r in R (cid:48)(cid:48) do if fitness( r ) < then R (cid:48) .add( r ) updateColumnPool( r ) return R (cid:48) 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 horizonCommodity Availability Time: nearly 30% commodities are available at the start ofthe planning horizonEmergency tasks: 10% to 30% of total commodities (i.e. time window < Table 5 Details of the artificial instances
Instance no. of shifts no. of commodities no. of FTL unitsI4-1 4 51 360I4-2 4 56 340I4-3 4 50 266I4-4 4 87 624I4-5 4 71 305I6-1 6 77 489I6-2 6 79 564I6-3 6 94 581I6-4 6 105 783I6-5 6 99 818I8-1 8 106 888I8-2 8 120 831I8-3 8 106 939I8-4 8 124 1067I8-5 8 127 971 maxColumns ). If maxColumns is set too small, more RMP solving calls are required which are computationally veryexpensive. However, if the maxColumns is set too large, time to solve each RMP wouldalso increase (the extreme case is that all feasible columns are included in RMP and it is
25 of 34 quivalent to the original problem). Some initial experiments suggest that maxColumns =1000 provides a good trade-off. We use this value on the understanding that it may not bethe best parameter for every instance. Note that in our method, in the early search stage,we permit our method to use more trucks than the limit ( n ), but this constraint will laterbe restored at the end of the column generation procedure. Gurobi 8 linear programminglibraries were used in conjunction with Java 7.0. These experiments were run on a PC withan 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 thetotal 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. givesthe objective value which is the total travel distance. Hereafter, P1 and P2 are shortabbreviations 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 columngeneration process. On average, P2 generates 1165 more columns than P1, resulting inlonger running times, but P2 uses 3089km less distance than P1. Seemingly, this fact is dueto 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 largernumber of columns which did not result in a smaller objective value.The performance of both algorithms is also compared with the results from the GurobiIP solver with the default algorithm setting in two experiments. The first experiment allowsthe 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 timetaken 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 morethan 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.
26 of 34 able 6 Comparison of the two pricing methods (artificial data)
P1 P2 Gurobi IP Solver RandomInstance T Obj. Col. T Obj. Col. T Obj. Obj.* Obj.I4-1 23 15516 4691 21 14154 3005 1215 13746 n.a. 22530I4-2 93 16480 9448 151 15988 5976 1208 15823 n.a. 22743I4-3 3 12793 2281 6 11067 4319 155 11037 n.a. 15885I4-4 271 27557 4674 711 25642 8811 1736 25307 28819 33583I4-5 187 13407 4021 343 11435 6430 1193 11429 14624 25798I6-1 131 27566 7742 193 25540 13985 1153 24713 29542 34589I6-2 175 26719 3046 507 23374 4861 2772 21665 n.a. 31294I6-3 87 32009 3142 513 30124 5321 2604 30029 n.a. 31889I6-4 218 41301 2170 290 35935 3321 9462 33898 n.a. 54497I6-5 172 33799 2040 420 30207 3216 5406 29223 n.a. n.a.I8-1 276 53871 2863 694 50178 2724 14890 49797 n.a. 70269I8-2 323 38589 2199 958 33532 3361 17202 32668 n.a. 54667I8-3 213 44856 3539 970 39643 2701 10006 38108 n.a. n.a.I8-4 479 35850 2022 919 32307 2286 36132 31979 n.a. 46778I8-5 213 45066 2414 704 39911 3455 19170 37979 n.a. 61476Avg. 191 31025 3753 493 27936 4918 8287 27160 n.a. 38923P:pricing method; T:Total running time(s);Col.:Total columns generated; Obj.:Objective value(km).Obj*.:Objective value with a limited computational time.n.a.:Failed to find feasible solution in the given time.
Gurobi: t=-3.389, p < < < maxColumns randomly selected routes (from all possible feasible routes) of the RMP at each iteration.All the other settings were kept the same as before. Column Random in Table 6 presentsthe objective values based on the average of five runs. As can be seen, the results are signif-icantly inferior to those by P1 or P2 (P1 vs. Random: t=-6.496, p < <
7. Experiments on instances with a very large R For larger instances, the feasible route set R can become very big and therefore it becomesimpossible to enumerate them all as we did in the previous section. In this section, we inves-tigate the effectiveness and performance of the two metaheuristic approaches presented insection 5.5. As the evidence from previous experiments suggest P2 performs better than P1,
27 of 34 or 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 eachiteration. As maxColumns is the only parameter used in the proposed VNS based columngenerator, parameter tuning for VNS is omitted. We now illustrate parameter tuning forthe GA.
The parameters used in the proposed GA are the population size ( populationSize ), thegeneration size ( generations ), the probability of mutation ( mutationRate ), and the tourna-ment 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 showsthe results with the algorithm with different population sizes and different number of gen-erations for the two most challenging problem instances LB8-1 and LB8-2. Each instancewas 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. Withthe consideration of algorithm efficiency, we choose the combination of populationSize =500and generations =500.
Table 7 Experiment results for evaluating populationSize and generations populationSize 10 100 200 500 10 100 200 500 10generations 10 10 10 10 100 100 100 100 200Avg. 25043 22154 21797 21280 21930 20991 20588 20079 21647populationSize 100 200 500 10 100 200 500 1000 2000generations 200 200 200 500 500 500 500 1000 2000Avg. 20676 20128 19999 21205 20228 20081 19879 19775 19724
Due to the very large amount of computational time required, we select a total of sixinstances, two from the real-life instance set and four from the artificial instance set used byBai et al. (2015). The instance names starting with NP are real-life instances while thosestarting 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 5 liststhe average objective values of five runs of experiments. The horizontal axis defines the
28 of 34 umber of RMP iterations while the vertical axis given the objective values. The resultssuggest that the GA is a faster converging algorithm thanks to its population based searchframework and capacity to evolve a set of routes instead of a single one. (a) Instance NP4-1 (b) Instance NP8-1(c) Instance LB8-1 (d) Instance TU8-7(e) Instance TB4-4 (f) Instance LU4-6
Figure 5 Comparison of performance maxIteration=10.
As experiments show that the VNS converges in 30 RMP iterations, for a fairer com-parison, 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).
29 of 34 able 8 presents the running time and objective values obtained by the pricing methodP2 using VNS and GA column generators (
P2-VNS , P2-GA ), method (Bai et al.2015), meta-heuristic methods using simulated annealing hyper-heuristic ( SAHH ), andreactive shaking variable neighbourhood search ( rsVNS ). In terms of objective values, theaverage results in increasing order are as follows:
P2-GA < < P2-VNS < SAHH < rsVNS . The best results are highlighted in bold. The average running times show anincreasing order as follows: P2-VNS < SAHH < P2-GA < rsVNS < . Two tailedPaired t-tests ( α = 0.01) were conducted to compare P2-VNS and P2-GA with 3-stage,SAHH and rsVNS approaches. There are significant difference in the solution quality forP2-VNS vs. P2-GA(t=3.698, p < < < < < P2-GA and
P2-VNS are able to find bettersolution in less time than most of the existing algorithms.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 method performs well for tight instances, but it does less wellfor large and loose instances. The reason is that it employs an integer programming solverso its solution time increases exponentially with large problem sizes. The proposed columngeneration methods are able to find effective columns in order to reduce the problemsize, therefore, compared with the method, the solving time of column generationmethod is significantly decreased for large instances. However, the advantage of columngeneration method may not be obvious for small problem instances (i.e. tight and smallinstances) as the iterative RMP solving comprises a significant proportions of the run timeby the algorithm.
8. Conclusions
We have presented an innovative FTL routing formulation assisted by dynamic cuts andinvestigated column generation based approaches which are particularly effective on verylarge 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 usingthe following strategies: 1) Infeasible flow assignments are allowed in the column generationprocess but will be fixed by adding cuts in the end; 2) To reduce the number of decision
30 of 34 able 8 Comparisons with previous results
P2-VNS P2-GA 3-stage rsVNS SAHH
Instance T Obj. T Obj. T Obj. T Obj. T Obj.NP4-1 126 13978 405 13860 33301
186 12508 4800 13577 383 13495LU4-5 66 18242 183
79 13033 4800 13760 232 14377TU4-8 12 18920 125 17956 138
148 21338 9600 21657 2602 21689TB8-4 112 28001 193
561 28167 9600 28398 2391 28305LU8-5 73 23288 248 22453 4380
66 27406 9600 28197 434 28162Large 8355 105793 7801 n.a. n.a. 9600 142258 15848 141252Average 349 22777 847
Large instance. variables, some constraints are pre-processed offline; 3) The shift that a route belongs to isnot restricted; 4) To avoid traversing the entire search tree, metaheuristics are implementedto repeatedly identify new columns with negative reduced costs in light of both new pricinginformation and latest columns on the basis of the RMP problem.Two pricing methods and two approaches for generating initial routes were proposed andevaluated. The result indicates that the proposed solution methods improve the existingalgorithms both in terms of the computational time and the solution quality. We believethe advantageous features of the indirect solution encoding of the FTL problem haven beenfully explored in this paper and the proposed solution methods can efficiently solve real-lifedrayage container operation problem with long planning horizon covering multi-shifts.
31 of 34 cknowledgments
This work is supported by the National Natural Science Foundation of China (Grant No.71471092), Natural Science Foundation of Zhejiang Province (Grant No. LR17G010001)and Ningbo Municipal Bureau of Science and Technology (Grant No. 2017D10034).
References
Alba, E., Almeida, F., Blesa, M., Cabeza, J., Cotta, C., D´ıaz, M., Dorta, I., Gabarr´o, J., Le´on, C., Luna, J.,et al., 2002. Mallba: A library of skeletons for combinatorial optimisation, in: European Conference onParallel Processing, Springer. pp. 927–932.Bai, R., Xue, N., Chen, J., Roberts, G.W., 2015. A set-covering model for a bidirectional multi-shift fulltruckload vehicle routing problem. Transportation Research Part B: Methodological 79, 134 – 148.Blum, C., Puchinger, J., Raidl, G.R., Roli, A., 2011. Hybrid metaheuristics in combinatorial optimization:A survey. Applied soft computing 11, 4135–4151.Braekers, K., Caris, A., Janssens, G., 2013. Integrated planning of loaded and empty container movements.OR Spectrum 35, 457–478.Braekers, K., Caris, A., Janssens, G.K., 2011. A deterministic annealing algorithm for a bi-objective fulltruckload vehicle routing problem in drayage operations. Procedia - Social and Behavioral Sciences 20,344 – 353.Braekers, K., Caris, A., Janssens, G.K., 2014. Bi-objective optimization of drayage operations in the servicearea of intermodal terminals. Transportation Research Part E: Logistics and Transportation Review65, 50 – 69. Special Issue on: Modeling, Optimization and Simulation of the Logistics Systems.Caris, A., Janssens, G., 2009. A local search heuristic for the pre- and end-haulage of intermodal containerterminals. Computers & Operations Research 36, 2763 – 2772.Caris, A., Janssens, G.K., 2010. A deterministic annealing algorithm for the pre- and end-haulage of inter-modal container terminals. IJCAET 2, 340–355.Chen, J., 2016. An Intelligent Container Transportation System Using Novel Modelling and Metaheuristics.Ph.D. thesis. University of Nottingham.Chen, J., Bai, R., Qu, R., Kendall, G., 2013. A task based approach for a real-world commodity routingproblem, in: Computational Intelligence In Production And Logistics Systems (CIPLS), 2013 IEEEWorkshop on, IEEE. pp. 1–8.Chung, K.H., Ko, C.S., Shin, J.Y., Hwang, H., Kim, K.H., 2007. Development of mathematical models forthe container road transportation in korean trucking industries. Computers and Industrial Engineering53, 252 – 262.Clements, D., Crawford, J., Joslin, D., Nemhauser, G., Puttlitz, M., Savelsbergh, M., 1997. Heuristic opti-mization: A hybrid ai/or approach, in: Workshop on Industrial Constraint-Directed Scheduling.
32 of 34 oslovich, L., Pesenti, R., Ukovich, W., 2006. Minimizing fleet operating costs for a container transportationcompany. European Journal of Operational Research 171, 776 – 786. Feature Cluster: Heuristic andStochastic Methods in Optimization Feature Cluster: New Opportunities for Operations Research.Funke, J., Kopfer, H., 2016. A model for a multi-size inland container transportation problem. TransportationResearch Part E: Logistics and Transportation Review 89, 70–85.Gendreau, M., Nossack, J., Pesch, E., 2015. Mathematical formulations for a 1-full-truckload pickup-and-delivery problem. European Journal of Operational Research 242, 1008–1016.Ghezelsoflu, A., Di Francesco, M., Frangioni, A., Zuddas, P., 2018. A set-covering formulation for a drayageproblem with single and double container loads. Journal of Industrial Engineering International 14,665–676.Ileri, Y., Bazaraa, M., Gifford, T., Nemhauser, G., Sokol, J., Wikum, E., 2006. An optimization approachfor planning daily drayage operations. Central European Journal of Operations Research 14, 141–156.Imai, A., Nishimura, E., Current, J., 2007. A lagrangian relaxation-based heuristic for the vehicle routingwith full container load. European Journal of Operational Research 176, 87 – 105.Jula, H., Dessouky, M., Ioannou, P., Chassiakos, A., 2005. Container movement by trucks in metropolitannetworks: modeling and optimization. Transportation Research Part E: Logistics and TransportationReview 41, 235 – 259.Lahrichi, N., Crainic, T.G., Gendreau, M., Rei, W., Cri¸san, G.C., Vidal, T., 2015. An integrative cooperativesearch framework for multi-decision-attribute combinatorial optimization: Application to the mdpvrp.European Journal of Operational Research 246, 400–412.Lai, M., Crainic, T.G., Francesco, M.D., Zuddas, P., 2013. An heuristic search for the routing of heteroge-neous trucks with single and double container loads. Transportation Research Part E: Logistics andTransportation Review 56, 108 – 118.Marino, A., Pr¨ugel-Bennett, A., Glass, C.A., 1999. Improving graph colouring with linear programming andgenetic algorithms, in: In, Citeseer.Muthuraman, S., Venkatesan, V.P., 2017. A comprehensive study on hybrid meta-heuristic approachesused for solving combinatorial optimization problems, in: 2017 World Congress on Computing andCommunication Technologies (WCCCT), IEEE. pp. 185–190.Nossack, J., Pesch, E., 2013. A truck scheduling problem arising in intermodal container transportation.European Journal of Operational Research 230, 666 – 680.Puchinger, J., Raidl, G.R., 2004. An evolutionary algorithm for column generation in integer programming:an effective approach for 2d bin packing, in: International Conference on Parallel Problem Solving fromNature, Springer. pp. 642–651.Puchinger, J., Raidl, G.R., 2005. Combining metaheuristics and exact algorithms in combinatorial optimiza-tion: A survey and classification, in: International work-conference on the interplay between naturaland artificial computation, Springer. pp. 41–53.
33 of 34 uchinger, J., Raidl, G.R., Koller, G., 2004. Solving a real-world glass cutting problem, in: EuropeanConference on Evolutionary Computation in Combinatorial Optimization, Springer. pp. 165–176.Ritzinger, U., Puchinger, J., Hartl, R.F., 2016. A survey on dynamic and stochastic vehicle routing problems.International Journal of Production Research 54, 215–231.Ruiyou Zhang, Won Young Yun, I.M., 2009. A reactive tabu search algorithm for the multi-depot containertruck transportation problem. Transportation Research Part E: Logistics and Transportation Review45, 904 – 914.Smilowitz, K., 2006. Multi-resource routing with flexible tasks: an application in drayage operations. IieTransactions 38, 577–590.Soares, R., Marques, A., Amorim, P., Rasinm¨aki, J., 2019. Multiple vehicle synchronisation in a full truck-load pickup and delivery problem: A case-study in the biomass supply chain. European Journal ofOperational Research 277, 174–194.Strandmark, P., Qu, Y., Curtois, T., 2020. First-order linear programming in a column generation basedheuristic approach to the nurse rostering problem. Computers & Operations Research , 104945.Vasquez, M., Hao, J.K., et al., 2001. A hybrid approach for the 0-1 multidimensional knapsack problem, in:IJCAI, pp. 328–333.Vidal, T., Crainic, T.G., Gendreau, M., Prins, C., 2014. A unified solution framework for multi-attributevehicle routing problems. European Journal of Operational Research 234, 658–673.Vidovi´c, M., Popovi´c, D., Ratkovi´c, B., Radivojevi´c, G., 2017. Generalized mixed integer and vns heuristicapproach to solving the multisize containers drayage problem. International Transactions in OperationalResearch 24, 583–614.Vidovic, M., Radivojevic, G., Rakovic, B., 2011. Vehicle routing in containers pickup up and deliveryprocesses. Procedia - Social and Behavioral Sciences 20, 335 – 343.Xie, Y., Liang, X., Ma, L., Yan, H., 2017. Empty container management and coordination in intermodaltransport. European Journal of Operational Research 257, 223–232.Xiubin Wang, A.C., 2002. Local truckload pickup and delivery with hard time window constraints. Trans-portation Research Part B: Methodological 36, 97 – 112.Zhang, R., Yun, W., Kopfer, H., 2010. Heuristic-based truck scheduling for inland container transportation.OR Spectrum 32, 787–808.Zhang, R., Yun, W.Y., Moon, I.K., 2011. Modeling and optimization of a container drayage problem withresource constraints. International Journal of Production Economics 133, 351 – 359. Leading Edge ofInventory Research.uchinger, J., Raidl, G.R., Koller, G., 2004. Solving a real-world glass cutting problem, in: EuropeanConference on Evolutionary Computation in Combinatorial Optimization, Springer. pp. 165–176.Ritzinger, U., Puchinger, J., Hartl, R.F., 2016. A survey on dynamic and stochastic vehicle routing problems.International Journal of Production Research 54, 215–231.Ruiyou Zhang, Won Young Yun, I.M., 2009. A reactive tabu search algorithm for the multi-depot containertruck transportation problem. Transportation Research Part E: Logistics and Transportation Review45, 904 – 914.Smilowitz, K., 2006. Multi-resource routing with flexible tasks: an application in drayage operations. IieTransactions 38, 577–590.Soares, R., Marques, A., Amorim, P., Rasinm¨aki, J., 2019. Multiple vehicle synchronisation in a full truck-load pickup and delivery problem: A case-study in the biomass supply chain. European Journal ofOperational Research 277, 174–194.Strandmark, P., Qu, Y., Curtois, T., 2020. First-order linear programming in a column generation basedheuristic approach to the nurse rostering problem. Computers & Operations Research , 104945.Vasquez, M., Hao, J.K., et al., 2001. A hybrid approach for the 0-1 multidimensional knapsack problem, in:IJCAI, pp. 328–333.Vidal, T., Crainic, T.G., Gendreau, M., Prins, C., 2014. A unified solution framework for multi-attributevehicle routing problems. European Journal of Operational Research 234, 658–673.Vidovi´c, M., Popovi´c, D., Ratkovi´c, B., Radivojevi´c, G., 2017. Generalized mixed integer and vns heuristicapproach to solving the multisize containers drayage problem. International Transactions in OperationalResearch 24, 583–614.Vidovic, M., Radivojevic, G., Rakovic, B., 2011. Vehicle routing in containers pickup up and deliveryprocesses. Procedia - Social and Behavioral Sciences 20, 335 – 343.Xie, Y., Liang, X., Ma, L., Yan, H., 2017. Empty container management and coordination in intermodaltransport. European Journal of Operational Research 257, 223–232.Xiubin Wang, A.C., 2002. Local truckload pickup and delivery with hard time window constraints. Trans-portation Research Part B: Methodological 36, 97 – 112.Zhang, R., Yun, W., Kopfer, H., 2010. Heuristic-based truck scheduling for inland container transportation.OR Spectrum 32, 787–808.Zhang, R., Yun, W.Y., Moon, I.K., 2011. Modeling and optimization of a container drayage problem withresource constraints. International Journal of Production Economics 133, 351 – 359. Leading Edge ofInventory Research.