ABSTRACT
This paper is based on a practical project jointly conducted by a major trucking company and a renowned operations research consulting firm. It studies a largescale, realtime truckload pickup and delivery problem. A number of cost factors are carefully measured such as loaded/empty travel distance, travel time, crew labor, equipment rental or operational cost, and revenue for completing the movements. This paper proposes a generalized decomposition algorithm that is capable of considering sophisticated business rules. The goal is to recommend executable and efficient truck routing decisions to minimize operating costs. Numerical tests are conducted with operational data from J.B.HUNT. A fleet of 5,000 trucks is considered in this experiment. The test result not only shows significant cost savings but also demonstrates computational efficiency for realtime application.
Keywords:
Truck routing; decomposition algorithm; column generation
1 INTRODUCTION
Research on vehicle routing has seen wide applications in the transportation industry such as truck, railway, airline and pipeline, which saves on costs while covers more demands. The objective is to improve the fleet operational efficiency. This research is important within the context of increasingly automated (or computerized) systems to support decision making for the routing/scheduling routines. It can be easily extended to solving other problems that may look seemingly different, but belong to the same family of NPcompleteness (^{Li et al., 2010}15 LI Y, MIAO Q & WANG X. 2010. A Column Generation Method for the U.S. Army Logistics Air Fleet Scheduling. In: Transportation Research Record: Journal of the Transportation Research Board , 2197: 3642., ^{2012}16 LI Y, MIAO Q & WANG X. 2012. U.S. Postal Airmail Routing Optimization. In: Transportation Research Record: Journal of the Transportation Research Board , 2234: 2128. and ^{2014}14 LI Y, MIAO Q & WANG X. 2014. HighSpeed Train Network Routing with Column Generation. Transportation Research Record: Journal of the Transportation Research Board, 2466: 5867.). The social economic impact of this problem cannot be overestimated. In the past century, especially during the last 60 years, numerous research efforts have been made specially to solve the trucking industry’s routing problems.
The main contribution of this paper is to introduce a decomposition algorithm based on a practical project. The goal is to recommend executable and efficient truck routing decisions to minimize operating costs. Numerical tests are conducted with operational data from J.B.HUNT. The test result not only shows significant cost savings but also demonstrates computational efficiency for realtime application.
Within the company, a fleet of vehicles are to be scheduled and routed to serve a known set of loads, and each load is considered as a full truckload with an origin and a destination (OD) associated with time window constraints for pickup and delivery. There are a finite number of drivers (or, crew in general) available to be assigned. Each driver has to return home within a period of 14 days according to the company rule as well as union regulation to retain drivers (for example, maximum consecutive hours of driving). A route consists of a sequence of moves to be fulfilled by a single vehicle. A typical route needs to have a starting location the same as termination location, which makes a tour to the associated vehicle. Depending on the service type, certain vehicles or drivers may not be eligible. For example, Dedicated Contract Service is a service primarily utilizes semitrailer trucks to transport cargo across the country; Intermodal Service partners with railways, commercial airlines or port authorities to move containers. Furthermore, once the drivers are committed to certain loads, diversion is typically not allowed with the exception from management approval. For the convenient purpose, this study does not consider diversion in the problem formulation and presentation.
The primary objective is to minimize the operation cost. In addition to cost of labor, other major costs include fuel rate, travel distance, equipment rental and length of operation. Revenue from serving each load can be considered as a positive profit or negative cost whenever a load is covered or service is completed. Empty truck moves and empty driver moves (also called deadhead moves) do not generate revenue, thus are only considered costs to the company. It is expected that the developed algorithm is built into a decision making system to recommend routings automatically.
2 LITERATURE REVIEW
In contrast to the lessthantruckload where each vehicle carries multiple customer demands, truckload only allows a truck to serve a single customer demand each time. An example route of a vehicle starts from depot (or source), loads goods at a factory or warehouse, moves to the load destination facility, unloads the load before goes to another loading location and repeats the same truckload movement. In the end, the vehicle returns to the depot (^{Labadie & Prins, 2012}11 LABADIE N & PRINS C. 2012. Vehicle Routing Nowadays: Compact Review and Emerging Problems. In: Production System and Supply Chain Management in Emerging Countries: Best Practices, Selected Papers From the International Conference on Production Research (ICPR). Springer, Berlin, pp. 141166.). This is a typical vehicle routing/scheduling problem in literature.
There is a tremendous literature dealing with mathematical modeling of vehicle routing problems (VRPs), all of which imply that it is impossible to enumerate exhaustively all the possible routes for the vehicles. Earlier efforts to solving VRPs are summarized in extensive reviews as ^{Desrochers et al. (1990}9 DESROCHERS M, LENSTRA JK & SAVELSBERGH MWP. 1990. A Classification Scheme for Vehicle Routing and Scheduling Problems. European Journal of Operational Research , 46: 322332.) and ^{Solomon (1987}19 SOLOMON MM. 1987. Algorithms for the Vehicle Routing and Scheduling Problem with Time Window Constraints. Operations Research , 35: 254265.). ^{Bramel et al. (1992}5 BRAMEL J, COFFMAN JR EG, SHOR P & SIMCHILEVI D. 1992. Probabilistic Analysis of Algorithms for the Capacitated Vehicle Routing Problem with Unsplit Demands. Operations Research, 40: 10951106., ^{1994}6 BRAMEL J, LI CL & SIMCHILEVI D. 1994. Probabilistic Analysis of the Vehicle Routing Problem with Time Windows. American Journal of Mathematical and Management Science, 13: 267322.) presented some probabilistic analyses of earlier heuristics for the deterministic version of the problem. Later, VRPs with stochastic features drew more and more attention from the research community. These features included but were not limited to stochastic load distributions (^{Golden & Stewart, 1978}10 GOLDEN BL & STEWART WR. 1978. Vehicle Routing with Probabilistic Demands. HOGBEN D & FIFE D. (eds.). Computer Science and Statistics: Tenth Annual Symposium on the Interface. NBSS Special Publication, National Book Service, Toronto, Canada, pp. 252259.; ^{Stewart & Golden, 1983}20 STEWART WR & GOLDEN BL. 1983. Stochastic Vehicle Routing: A Comprehensive Approach. European Journal of Operational Research , 14: 371385.; and ^{Bastian & Rinnooy Kan, 1992}1 ASTIAN C & RINNOOY KAN AHG. 1992. The Stochastic Vehicle Routing Problem Revisited. European Journal of Operational Research, 56: 407412.), stochastic travel time (^{Cook & Russell, 1978}7 COOK TM & RUSSELL RA. 1978. A Simulation and Statistical Analysis of Stochastic Vehicle Routing with Timing Constraints. Decision Science, 9: 673687.; and ^{Berman & SimchiLevi, 1989}2 BERMAN O & SIMCHILEVI D. 1989. The Traveling Salesman Location Problem on Stochastic Networks. Transportation Science, 23: 5457.), or stochastic locations (^{Laporte et al., 1994}12 LAPORTE G, LOUVEAUX FV & MERCURE H. 1994. A Priori Optimization of the Probabilistic Traveling Salesman Problem. Operations Research , 42: 543549.; and ^{Bertsimas & Howell, 1993}3 BERTSIMAS DJ & HOWELL LH. 1993. Further Results on the Probabilistic Traveling Salesman Problem. European Journal of Operational Research , 65: 6895.). As the information technologies came into play, recent realtime and dynamic VRP problems became increasingly important (^{Yang et al., 2004}21 YANG J, JAILLET P & MAHMASSANI H. 2004. RealTime Multivehicle Truckload Pickup and Delivery Problems. Transportation Science , 38(2): 135148.). ^{Powell et al. (1995}17 POWELL WB, JAILLET P & ODONI A. 1995. Stochastic and Dynamic Networks and Routing. BALL MO, MAGNANTI TL, MONMA CL & NEMHAUSER GL. (eds.). Handbooks in Operations Research and Management Science, 8. Network Routing. Elsevier (NorthHolland), Amsterdam, pp. 141296.) presented a survey of dynamic fleet optimizations dealing with some general issues. Later work of ^{Powell et al. (2000}18 POWELL WB, SNOW W & CHEUNG RK. 2000. Adaptive Labeling Algorithms for the Dynamic Assignment Problem. Transportation Science , 34: 5066.) developed a practical model to consider dynamic assignment of drivers to known demands, which provided significant insights to our study problem here. Reoptimization policies are further introduced and tested in ^{Yang, et al. (1998}22 YANG J, JAILLET P & MAHMASSANI HS. 1998. Online Algorithms for Truck Fleet Assignment and Scheduling Under Realtime Information. Transportation Research Record, 1667: 107113.). The most recent reviews summarizing the stateofart techniques are available in ^{Laporte et al. (2013}13 LAPORTE G, TOTH P & VIGO D. 2013. Vehicle Routing: Historical Perspective and Recent Contributions. European Journal of Operational Research , 1: 14.), ^{Derigs et al. (2013}8 DERIGS U, PULLMANN M & VOGEL U. 2013. Truck and Trailer RoutingProblems, Heuristics and Computational Experience. Computers and Operations Research , 40: 536546.), and ^{Braekers et al. (2016}4 BRAEKERS K, RAMAEKERS K & VAN NIEUWENHUYSE I. 2016. The vehicle routing problem: State of the art classification and review. Computers & Industrial Engineering, 99: 300313.).
Published articles on implementation, however, are less popular compared with the counterpart on theoretical studies. When it comes to practical VRP projects, people often look for details about how algorithms are developed and implemented. This paper is originated from industry projects within a major trucking company and focuses on proven practical techniques. Implementation details are revealed so that the readers can have a better understanding of the business logic. It aims at utilizing a combined optimization method and advanced information technology to develop a realtime dispatching system. The computational time invested in searching for better decisions in terms of shorter routes and more revenue should be cautiously balanced with the needs of coming up with a decision in a timely manner.
3 SOLUTION APPROACH
This decision support system requires a list of components to function properly. Information is updated via continuous data feed, stored in a centralized data warehouse. Forecast module provides projections on the existing truck moves and the estimation on the future demand. Prepossessing module turns the raw data into the format that are favored by the optimizer. Optimization component is the core module that handles the mathematical formulation and sophisticated algorithm.
3.1 Information Updating
Although the information arrives all the time, it is not required to run optimizer over the entire fleet all the time. Actually, there are two types of optimization jobs. First type of job is for planning purpose and runs daily. This is normally offline and is scheduled to run overnight or at meal breaks, so that the results are immediately available to dispatchers when they start work or resume work. Second type of job is for realtime optimization. This type of job only considers a smaller set of the fleet because most vehicles are already committed to their loads/route and do not need reoptimization. This is normally triggered from managing team whenever there is a need.
3.2 Forecast Engine
Certain loads are visible before they become available for pickup. For example, the containers can be moved by trains with a projected arrival time at the pickup location. The forecast engine would gather the information from partnered carriers to make reasonable projections on the availability of their own loads. This helps the optimization engine to look beyond the current demands and make informed decision for the immediate future. Beyond the visible horizon, long term forecast is also necessary to plan ahead in terms of fleet sizing and infrastructure change at strategic level.
3.3 Preprocessing
The preprocessor needs to selectively package the raw data into the network format that the downstream optimizer requires. The basic elements in this network are link (edge) and node (vertex), which are explained in section 5. Examining the feasibility of the links can reduce the burden on the optimizer. For example, it needs to filter out an infeasible link that tries to merry a wrong type of truck to a load. Each node represents an activity such as getting a truck, load or unload. After excluding infeasible links, certain business rules will further reduce the number of links by checking the distance, time or any other resource consumed between two nodes. If certain link violates resource limit, then this link is also excluded. The ultimate goal here is to allow the optimization module find reliable routes through the network easily.
3.4 Optimization Strategy
Due to the complicated resource constraints and business rules, it is inconvenient to formulate this truck routing problem into a network flow model. The alternative approach is to apply partition or setcovering model, where the objective function minimizes the combined cost of routes being selected.
An obvious difficulty here is the attempt to explicitly enumerate all the feasible routes, which are typically required in a partition model. Given the total number of loads and available drivers/ trucks, the number of combination is extremely large and increases exponentially with the problem size. Additionally, the business rules are not easy to formulate, such as the maximum hours each driver can take in a tour. Each tour must start and end at driver’s home terminal, and there are specified time windows for pickup and delivery. Here we propose a decomposition method to bypass the necessity of evaluating all possible combinations. During the iteration process, a master problem picks the best set of routes to minimize the total cost and price each load in terms of dual variables for the subproblem. The subproblem then updates the network and generates better routes for master problem to consider with. A schema that describes the entire framework of this column generation procedure is presented below in Figure 1. Here in this problem, one can assume that a column refers to a route.
4 MASTER PROBLEM FORMULATION
where:
x_{j} decision variable, 1 if route j is selected, 0 otherwise
Ω_{k} the set of routes in the kth iteration of the column generation procedure
c_{j} the cost associated with route j (operating cost  load revenue)
N the set of loads
${\delta}_{i}^{j}$ if route j covers load i, 0 otherwise
π_{i} dual associated ith constraint (for load i)
The objective is to minimize the total cost to cover a known set of load demands. The revenue generated from moving a load is considered a negative cost when building the routes, and it usually causes the total route cost to be negative, otherwise the route would not be profitable. Equation (2) requires all the load movement to be assigned exactly once. Each load has a dual value π _{i} associated with it, which is the shadow price the load given the current master problem.
It is important to note that during the iteration, the master problem is solved as linear programming relaxation to get dual values. The side constraints are conditional and are not formulated into the master. For example:
Although the company prefers its own salaried drivers over the contracted drivers due to the lowered operating cost. DR is the total number of both. When the total number of routes being selected exceeds limit DR, it is intuitive to sort the route costs in an ascending order to pick the first DR. The uncovered loads can be either rejected or outsourced. This is to facilitate the construction of the subproblem, where all the dual values associated with Eq. (2) are utilized to reflect what master problem desires.
5 CONSTRAINED SHORTEST PATH SUBPROBLEM
Our subproblem tries to find an optimal path go through network G that does not consume more than limited resources, such as time window, duration of path, distance, etc. Because of the additional constraints applied to the path, the subproblem is therefore a Constrained Shortest Path Problem (CSPP). A path here represents a truck route in reality.
5.1 Subproblem Formulation
where
x_{ij} are the decision variables, 1 if node i is followed by node j, 0 otherwise
c_{ij} the cost to proceed from node i to node j
A the set of eligible connections, link (i, j) ∈ A
π_{j} dual associated jth constraint (for load j)
The construction of the subproblem is essential to the usefulness of the generated routes and the overall solution time. A route is a collection of links that are connected by nodes, it is also referred as a path through network. Since it has a single objective function to generate the most profitable route, different cost elements need to be normalized within the network. For example, the overall cost on link i to j is based on the load/empty factor, distance (i, j), the revenue of moving load j, and dual π_{j} for load j. Given the same distance, an emptytruck move has a baseline cost and zero revenue. A loadedtruck move may double the baseline cost but gain revenue. The changing value of π_{j} is passed from master problem to adapt to the current need.
Each route begins from the source S and ends at sink T. In most cases, the source and sink are the same physical location but have different time windows. This is due to the business rule that the driver has to return to home terminal when the route is completed. Starting from source S, each route has to connect to an artificial node where a driver is “picked up”. Assigning a reasonable setup cost for such a node would encourage the model to minimize the total number of routes/drivers used.
Other resource constraints include (1) the maximum duration of each route, which is limited to 14 days; (2) the maximum combined travel distance for each route; (3) the time window for the load, which represents the earliest and latest time to pick up or deliver. In order to solve this multiple resource constrained shortest path problem with time windows, a specialized Label Setting Algorithm is applied.
5.2 The Algorithm for Shortest Path with Resource Constraints
Let ( ${D}_{ik}^{*}$, C_{ik} ) and ( ${D}_{jk}^{*}$, C _{jk} ) be the labels representing two different paths to node k. Then the first label dominates the latter, if and only if ( ${D}_{ik}^{*}$, C _{ik} ) − ( ${D}_{jk}^{*}$, C _{jk} ) ≥ (0, 0). The first label is smaller than the latter, i.e., ( ${D}_{ik}^{*}$, C _{ik} ) $\frac{L}{<}$ ( ${D}_{jk}^{*}$, C _{jk} ), if and only if ( ${D}_{ik}^{*}$, C _{ik} ) ≠ ( ${D}_{jk}^{*}$, C _{jk} ), and the first nonnull element of (( ${D}_{ik}^{*}$, C _{ik} ) − ( ${D}_{jk}^{*}$, C _{jk} ))is positive. A label ( ${D}_{k}^{*}$, C _{k} ) is efficient if none of the labels at node k can dominate it. The path corresponding to an efficient label is defined as an efficient path. Only efficient labels and paths are kept.
Let Q _{k} be the set of labels associated with the cost lower bound of path ending at node k ∈ V , and P _{k} be the set of labels associated with feasible paths. The P _{k} defines the primal function. The primal function provides an upper bound on the cost of efficient solutions at node k. When primal and dual functions have same value for a given stage, the labels in P and Q are associated to an efficient path for the current stage.
Then the algorithm used to solve the subproblem is presented as follows:
Step 1. (Initialization).
Find $\left({d}_{j}^{*},{c}_{j}\right)=\underset{\left(i,j\right)\in A}{\mathrm{min}lex}\left\{\left({d}_{ij}^{*},{c}_{ij}\right)\right\}$, for each node $j\in VS\cdot O={\displaystyle \underset{j\in VS}{\cup}\left({Q}_{j}\right)}$.
Step 2. (Calculation of the lexicographically smallest label in O).
If $O=\varnothing $, stop; [The current label(s) at the sink node T is (are) the shortest path(s)]. Otherwise, find the lexicographically smallest label $F\left(O\right)=\underset{\left({D}_{j}^{*},{C}_{j}\right)\in O}{\mathrm{min}lex}\left\{\left({D}_{j}^{*},{C}_{j}\right)\right\}=\left({D}_{j}^{*\text{'}},{C}_{j}\right)$.
Step 3. (Look for a label $\left({D}_{k}^{*},{C}_{k}\right)\in B\left(O\right)$ to be treated).
Step 4. (Uncertainty zone).
Define the uncertainty zone for node k:
Step 5. (Replacement of the label $\left({D}_{j}^{*},{C}_{j}\right)$.
(A) Calculation of the efficient labels defining the dual function:
(B) Calculation of the feasible paths defining the primal function:
(C) Update sets: ${P}_{k},{Q}_{k}:{P}_{k}\leftarrow {P}_{k}\cup ({R}_{k}\cap R{P}_{k}),{Q}_{k}\leftarrow {Q}_{k}({D}_{j}^{*\text{'}},{C}_{j}^{\text{'}})\cup {R}_{k}$.
(D) $O\leftarrow O({D}_{j}^{*\text{'}},{C}_{j}^{\text{'}})\cup [{R}_{k}({R}_{k}\cap R{P}_{k})]$.
(E) $B(O)\leftarrow B(O)\left({D}_{j}^{*\text{'}},{C}_{j}^{\text{'}}\right)$. If $B\left(O\right)=\varnothing $, return to Step 2; otherwise return to Step 3.
A simple graph with two resource constraints (time, distance) and cost is presented in Figure 2 and Table 1, followed by an example to show how the algorithm works. Figure 2 shows the node and link connections in the network. Table 1(a) shows the drivers’ profile. Table 1(b) shows loads’ profile. Table 1(c) shows the constraint for the drivers and loads. In this example, the constraints are time window and mileage. Table 1(d) shows the resources being consumed on each link, as well as the cost on each link. Note that the actual route, which has multiple links, can be 14days long and undertake more than just two or three loads.
Where a _{0} and b _{0} are the beginning and ending pickup times, and a _{1} and b _{1} are the minimum and maximum working miles, respectively.
Deadhead is typically discouraged because moving an empty truck comes with a cost but there is no (direct) gain in revenue. Simply speaking, driver D _{0} is originally located at Philadelphia and is eligible for Load L _{0} and L _{2} at the beginning. Once the load is delivered, the driver may become available again at the load destination (New York or Boston, depending on the load). Preprocessor skips the ineligible connections and adds eligible connections in term of links to the network. This preprocessing is a necessary step to reduce the problem size of the subproblem.
Step 1. Initialization.
Step 2. $F\left(O\right)={(\mathrm{7,0,}\infty )}_{D}{}_{0}$.
Step 3. $B\left(O\right)=\{({D}_{k}^{*},{C}_{k})\in O(\mathrm{7,0,}\infty )\frac{L}{\le}({D}_{k}^{*},{C}_{k})\frac{L}{<}\left(\mathrm{9,0,}\infty \right)\}$. We treat label (7,0,−∞).
Step 4. Uncertainty zone $T\text{}{R}_{D}{}_{0}=(\mathrm{8,500})$.
Step 5. Replacement of (7, 0, −∞). $\left(A\right){R}_{D}{}_{0}=R\text{}{P}_{D}{}_{0}=(\mathrm{7,0,50});\text{}\left(C\right){P}_{D}{}_{0}={Q}_{D}{}_{0}=(\mathrm{7,0,50}).\text{}O=\{{(\mathrm{7,0,}\infty )}_{D}{}_{1},{(\mathrm{7,0,}\infty )}_{L}{}_{1},{(\mathrm{8,0,}\infty )}_{L}{}_{2},{(\mathrm{8,0,}\infty )}_{L}{}_{0},{(\mathrm{8,0,}\infty )}_{L}{}_{3},{(\mathrm{12,0,}\infty )}_{T}\}$ . We treat label (7, 0 − ∞).
Step 4. Uncertainty zone $T\text{}{R}_{D}{}_{1}=\left(\mathrm{8,500}\right)$.
Step 5. Replacement of (7, 0, −∞). $\left(A\right){R}_{D}{}_{1}=R{P}_{D}{}_{1}=\left(\mathrm{7,0,50}\right);\text{}\left(C\right){P}_{D}{}_{1}={Q}_{D}{}_{1}=\left(\mathrm{7,0,50}\right).\text{}O=\{{\left(\mathrm{7,0,}\infty \right)}_{L}{}_{1},{\left(\mathrm{7,0,}\infty \right)}_{L}{}_{2},{\left(\mathrm{8,0,}\infty \right)}_{L}{}_{0},{\left(\mathrm{8,0,}\infty \right)}_{L}{}_{3},{\left(\mathrm{12,0,}\infty \right)}_{T}\}$, We treat label (7, 0 − ∞).
Step 4. Uncertainty zone $T\text{}{R}_{L}{}_{1}=(\mathrm{12,500})$.
Step 5. Replacement of (7, 0, −∞). $\left(A\right){R}_{L}{}_{1}=R{P}_{L}{}_{1}=\left(\mathrm{11,200,}60\right);\text{}\left(C\right){P}_{L}{}_{1}={Q}_{L}{}_{1}=\left(\mathrm{11,200,}60\right).\text{}O=\{{\left(\mathrm{7,0,}\infty \right)}_{L}{}_{2},{\left(\mathrm{8,0,}\infty \right)}_{L}{}_{0},{\left(\mathrm{8,0,}\infty \right)}_{L}{}_{3},{\left(\mathrm{12,0,}\infty \right)}_{T}\}$, We treat label (7, 0 − ∞).
Step 4. Uncertainty zone $T\text{}{R}_{L}{}_{2}=\left(\mathrm{12,500}\right)$.
Step 5. Replacement of (7, 0, −∞). $\left(A\right){R}_{L}{}_{2}=R{P}_{L}{}_{2}=\left(\mathrm{11,200,}50\right);\text{}\left(C\right){P}_{L}{}_{2}={Q}_{L}{}_{2}=\left(\mathrm{11,200,}50\right).\text{}O=\{{\left(\mathrm{8,0,}\infty \right)}_{L}{}_{0},{\left(\mathrm{8,0,}\infty \right)}_{L}{}_{3},{\left(\mathrm{12,0,}\infty \right)}_{T}\}$, We treat label (8, 0 − ∞).
Step 4. Uncertainty zone $T\text{}{R}_{L}{}_{0}=\left(\mathrm{18,500}\right)$.
Step 5. Replacement of (8, 0, −∞). $\left(A\right){R}_{L}{}_{0}=R{P}_{L}{}_{0}=\left\{\left(\mathrm{9,100,}60\right),\left(\mathrm{13,300,}170\right)\right\};\text{}\left(C\right){P}_{L}{}_{0}={Q}_{L}{}_{0}=\left\{\left(\mathrm{9,100,}60\right),\left(\mathrm{13,300,}170\right)\right\}.\text{}O=\{{\left(\mathrm{8,0,}\infty \right)}_{L}{}_{3},{\left(\mathrm{12,0,}\infty \right)}_{T}\}$ , We treat label (8, 0 − ∞).
Step 4. Uncertainty zone $T\text{}{R}_{L}{}_{3}=\left(\mathrm{18,500}\right)$.
Step 5. Replacement of (8, 0, −∞). $\left(A\right){R}_{L}{}_{3}=R{P}_{L}{}_{3}=\left\{\left(\mathrm{10,150,}70\right),\left(\mathrm{14,350,}170\right)\right\};\text{}\left(C\right){P}_{L}{}_{3}={Q}_{L}{}_{3}=\left\{\left(\mathrm{10,150,}70\right),\left(\mathrm{14,350,}170\right)\right\}.\text{}O=\left\{{\left(\mathrm{12,0,}\infty \right)}_{T}\right\}$, We treat label (12, 0 − ∞).
Step 4. Uncertainty zone $T\text{}{R}_{T}=\left(\mathrm{18,500}\right)$.
Step 5. Replacement of (12, 0, −∞). $\left(A\right){R}_{T}=R{P}_{T}=\{\left(\mathrm{10,110,}10\right),\left(\mathrm{11,160,}20\right),\left(\mathrm{14,310,}120\right),\left(\mathrm{15,360,}120\right);\text{}\left(C\right){P}_{T}={Q}_{T}=\{\left(\mathrm{10,110,}10\right),\left(\mathrm{11,160,}20\right),\left(\mathrm{14,310,}120\right),\left(\mathrm{15,360,}120\right)$.
Step 2. $O=\varnothing $. Stop. Current solution {(10, 110, −10), (11, 160, −20), (14, 310, −120), (15, 360, −120)}. There are four shortest paths from the source node S to the sink node T respecting the resource constraints: (1) Path: S → D0 → L 0 → T with cost −10; (2) Path: S → D1 → L 3 → T with cost −20; (3) Path: S → D1 → L 1 → L 0 → T with cost −120; (4) Path: S → D0 → L 2 → L 3 → T with cost −120. Further examination would suggest that Path 4 is dominated by Path 3. As a result, Paths 1, 2 and 3 are nondominated optimal and are eligible to be added into the path/route pools in the master problem.
Note that in this example, labels in P and Q are always having the same sets of labels and paths because the resource constraints never get violated. In the case where certain paths exceed the resource limit, those parts are deemed as infeasible and therefore excluded at the current stage. The lower bound Q, however, would keep this infeasible label for future treatments. This technique is extremely useful when the resource on a directed link is negative. In other words, the previously infeasible path may become feasible again by adding consecutive links to it. This means that before all the paths are examined, one simply cannot determine the best or even most feasible paths.
6 THE SCHEME OF BRANCH AND BOUND (B&B)
Commercial software CPLEX is used to solve the master problem (linear programming relaxation). The optimal solution is generally noninteger (fractional). The B&B scheme is thus invoked and embedded into the column generation process to obtain an integer solution. The following are the details for the implementation (refer to the flow chart). If there is any existing feasible solution can be extracted from the current truck operating plan, then it should be utilized as the initial solution to speed up the search process.
6.1 The Branching Strategy
It is easy to observe that if the solution is noninteger, there must exist at least a pair of consecutive load nodes t _{1} , t _{2} such that
where T R(t _{1}, t _{2}) is the set of routes in which t _{2} is executed immediately after t _{1}. Based on the set TR(t _{1} , t _{2}), the original problem is partitioned (branched) into two subproblems:
The testing shows that this strategy gives a more balanced search tree than default variable branching and generally finds an acceptable integer solution more quickly. Conceivably, this method is to branch on the relationship between two consecutive loads. It forces the link to be selected or to be eliminated in the subproblem.
6.2 How B&B Scheme is Embedded into the Column Generation Process?
When the B&B scheme is embedded into column generation process, a search tree is created. The nodes in the search tree (stored in a sorted queue) are treated one by one. If the queue is empty, the algorithm is terminated, and the optimal integer solution obtained so far is the solution we are seeking. By treating a node we mean that two branched nodes will be created from this node:

one with $\underset{j\in TR\left({t}_{1},{t}_{2}\right)}{{\displaystyle \sum {x}_{j}}}=0$ (left node),

and the other with $\underset{j\in TR\left({t}_{1},{t}_{2}\right)}{{\displaystyle \sum {x}_{j}}}=1$ (right node).
By creating a node we mean that the objective value and the solution at the node are found. Only the node corresponding to fractional solution with objective value smaller than that of the current optimal integer solution is inserted into the queue. Once a node corresponding to an integer solution is created and its objective value is smaller than that of the current optimal integer solution, then all untreated nodes in the queue with objective value not smaller than that of the created node will be pruned off. Whenever a node a created, a set of columns (for solving the master problem) and a network (for solving the subproblem) should be updated accordingly.
6.3 Information Storage and Retrieval
For solving the master problem and the subproblem at different nodes of the search tree, the information on columns and network must conform to the node to be created. There are two ways to get the information: one is from the root node; the other is from the node to be treated (parent node). With the first one, much more memory can be saved, but more computing time is needed (repeated computing), while with the second, the situation is reversed. There is a tradeoff between the memory and the computing time. In the first case, we only need to store the original network and the columns generated at root node. The search tree is a binary tree and the root node is at 0 level. We use a label consisting of (j + 1) identifiers to represent the location of a node at jth level: (j, j _{1} , j _{2} , ..., j _{j} ), where the first one j is a digit (or digits) which represents the level number of the node location; the second one j _{1} represents the location (left or right) of its ancestor at the first level; the third one j _{2} represents the location of its ancestor at second level, ..., and j _{j} represents its location at jth level. j _{1} , j _{2} , ..., j _{j} = L or R, L stands for.left branch (0branch); R stands for right branch (1branch). For example, a node with label (3, R, R, L) means that the node is at the third level of the search tree, and its ancestors are at the first and second levels on the right location. The thirdlevel node is at left location. In the second case, we only need to indicate where the treated node is located and then modify the columns and network of its parent node respectively to process.
6.4 The Modification of the Columns for the Master Problem
In the 0branch, we delete all columns i for which ${a}_{{t}_{1},i}=1$ and ${a}_{{t}_{2},i}=1$. In the 1branch, we delete all columns i for which ${a}_{{t}_{1},i}=1$ and ${a}_{{t}_{2},i}=0$, or ${a}_{{t}_{1},i}=0$ and ${a}_{{t}_{2},i}=1$, and the row corresponding to constraint for node t _{2} due to the redundancy. Figure 3 illustrates this process. The underlying assumption is that each load can be covered only once at most. The 2nd coverage for the same load will not bring additional revenue.
In 1branch, the dimension of the basis matrix is thus reduced by one. The modification of the columns depends on the location of the created node in the search tree and also depends on from where the information is obtained. For example, if we use the information on columns at root node and the created nodes are (3, R, R, L), then all columns at root node having ${a}_{{t}_{1},i}=1$ and ${a}_{{t}_{2},i}=0$, ${a}_{{t}_{1},i}=0$ and ${a}_{{t}_{2},i}=1$, ${a}_{{t}_{3},i}=1$ and ${a}_{{t}_{4},i}=0$, ${a}_{{t}_{3},i}=0$ and ${a}_{{t}_{4},i}=1$, ${a}_{{t}_{5},i}=1$ and ${a}_{{t}_{6},i}=1$ are deleted. The row corresponding to constraint for node t _{6} is also deleted due to the redundancy. We assume that the branching at first, second, and third level is based on T R(t _{1} , t _{2}), T R(t _{3} , t _{4}), T R(t _{5} , t _{6}) respectively.
6.5 Modification of the Network for the Subproblem
We must restrict the columns to be generated by the subproblem to those that are compatible with the current created node in the search tree. The structure of the network used to generate the feasible columns should be modified accordingly. In the 0branch, the columns covering consecutively nodes t _{1} and t _{2} are forbidden. As a reminder, there are the columns having t _{2} executed immediately after t _{1}. In the network we split the middle node M, where the link corresponding to t _{1} terminates and the link corresponding to t _{2} starts, into two nodes M _{1} and M _{2}; all links originally terminated at M are now moved to M _{1}, and all links originally started from M are moved to M _{2}. The resource constraints on M _{1} and M _{2} remain the same as M. In the 1branch, we force any route covering t _{1} to also cover t _{2} immediately. In the network, the two links corresponding to t _{1} and t _{2} are condensed into one link, and all the links originally terminated at the middle node of the previous two links are deleted. Figure 4 shows this modification. The cost and resource consumption for the condensed link are: (a) the cost of the newly created link equals to cost(t _{1}) + cost(t _{2}), where cost(t _{1}) and cost(t _{2}) are the cost of t _{1} and t _{2}, respectively; (b) the number of pieces of work still equals zero ; (c) the spread equals t _{1} + t _{2}; (d) the work time equals t _{1} + t _{2}.
With the same example, for the created node (3, R, R, L) the modified network at this node is obtained from the original network. After splitting the two nodes (one between the links corresponding to t _{1} and t _{2}, and the other between t _{3} and t _{4}) into four nodes, and condensing two links corresponding to t _{5} and t _{6} into one, all the links originally terminated at the middle node of the links corresponding to t _{5} and t _{6} are deleted.
In the iteration process, a set of the dual values are passed from master problem. The cost (t _{i} ) of the link i corresponding to load node t _{i} (ith constraint) is subtracted by dual value π before the subproblem algorithm is executed. Once the desirable columns are generated, the cost of each link i is set to cost (t _{i} ) (the original cost).
6.6 The Order of Choosing the Next Node and the Queue
The sum of fractional variables for a specified set T R(t _{1} , t _{2}) and the objective values at each node should determine the next node to branch on. The node with smallest sum of fractions is chosen as a priority. If there is a tie, the one with the smaller objective value is chosen. If there is still a tie, then choose any arbitrary one. One should always keep a sorting queue (nondecreasing sequence of the sum of fractional variables) for the nodes to be treated and take the first one from the queue. When a fractional solution occurs, and its objective value is smaller than that of current optimal integer solution value, then the node corresponding to this fractional solution is inserted into the queue. Before inserting, we need to know the sum of fractions. From the fractional solution, one must find the first column x _{i} < 1 and a set TR(t _{1} , t _{2}) consisting of two consecutive nodes t _{1} and t _{2} in column x _{i} such that
otherwise, the process is repeated until such set is found.
7 TESTING RESULTS AND REMARKS
To prove the capability of the prototyped routing optimizer, we select top 10 US cities with the largest populations and set Dallas as home based depot (as shown in Table 2 and Figure 5).
From these 10 cities we can have 45 nondirectional city pairs sorted in an alphabet order on the origin city and destination city (as shown in Table 3).
For each pair, we call a random number between 0 and 1, if this number is less than 0.5, set a direction from origin to destination (ex. Houston to San Diego); otherwise set an opposite direction from destination to origin (ex. San Diego to Houston), where Dallas is a home base depot that every driver must return back to Dallas within 7000 miles (approximately equivalent to2 weeks).
One of the performance measurements is Load Factor (LF). For each truck driver route, LF is the total loaded miles divide by total miles (loaded miles + empty miles) traveled starting from and returning to the depot: Dallas.
100 random generated data sets have been tested between current existing algorithm and our new columngeneration based algorithm to compare both number of drivers used and Load Factor (LF) to complete each set of 45 loads crossing those 10 cities. Table 4 gives the details on each test case. The average run time is very stable and is within a few minutes for the tested problem size. Table 5 shows the overall performance.
Our sample testing indicates that our new columngeneration based solution method will reduce about 16% of drivers and improve the Load Factor (LF) by about 9%.
The company also conducted an onsite experiment with real production data. About 3,500 externally contracted drivers, who were mainly independent business operators, and 1,500 corporate salaried drivers were considered in the problem for dispatching. Again, corporate drivers normally cost less than external contracts. About 15,000 loads are required to be covered. The decision variable is defined as a route, which is the combination of drivers, the loads and sequence of covering the loads. Without any of the preprocessing to eliminate routes, the number of decision variables is 2.5 × 10^{20} if the maximum of four loads are allowed in a single route. This number increases to 3.8 × 10^{24} if five loads are allowed in a single route. The prototype successfully considered all the major business criteria and constraints to have problem solved within 20 minutes on the company’s mainframe machine. After comparing with the company’s thencurrent operating plan, over 10% cost saving was achieved, which amounted to millions of dollars annually.
The proposed decomposition algorithm can be easily generalized to many other industries that share similar characteristics. For example, railway, airlines and pipeline share the similar network characteristics. Resource constraints can be modeled in terms of time, materials, vehicle or line capacity, distances, volumes or even headcount. Sometimes it is difficult to consider certain requirements such as recurring maintenance constraints and sophisticated union rules, solving the problem without these rules as in the case of this study at least can serve as a benchmark for the real problem. We leave as a future potential research effort to incorporate those additional constraints into routes developed.
REFERENCES

^{1}ASTIAN C & RINNOOY KAN AHG. 1992. The Stochastic Vehicle Routing Problem Revisited. European Journal of Operational Research, 56: 407412.

^{2}BERMAN O & SIMCHILEVI D. 1989. The Traveling Salesman Location Problem on Stochastic Networks. Transportation Science, 23: 5457.

^{3}BERTSIMAS DJ & HOWELL LH. 1993. Further Results on the Probabilistic Traveling Salesman Problem. European Journal of Operational Research , 65: 6895.

^{4}BRAEKERS K, RAMAEKERS K & VAN NIEUWENHUYSE I. 2016. The vehicle routing problem: State of the art classification and review. Computers & Industrial Engineering, 99: 300313.

^{5}BRAMEL J, COFFMAN JR EG, SHOR P & SIMCHILEVI D. 1992. Probabilistic Analysis of Algorithms for the Capacitated Vehicle Routing Problem with Unsplit Demands. Operations Research, 40: 10951106.

^{6}BRAMEL J, LI CL & SIMCHILEVI D. 1994. Probabilistic Analysis of the Vehicle Routing Problem with Time Windows. American Journal of Mathematical and Management Science, 13: 267322.

^{7}COOK TM & RUSSELL RA. 1978. A Simulation and Statistical Analysis of Stochastic Vehicle Routing with Timing Constraints. Decision Science, 9: 673687.

^{8}DERIGS U, PULLMANN M & VOGEL U. 2013. Truck and Trailer RoutingProblems, Heuristics and Computational Experience. Computers and Operations Research , 40: 536546.

^{9}DESROCHERS M, LENSTRA JK & SAVELSBERGH MWP. 1990. A Classification Scheme for Vehicle Routing and Scheduling Problems. European Journal of Operational Research , 46: 322332.

^{10}GOLDEN BL & STEWART WR. 1978. Vehicle Routing with Probabilistic Demands. HOGBEN D & FIFE D. (eds.). Computer Science and Statistics: Tenth Annual Symposium on the Interface. NBSS Special Publication, National Book Service, Toronto, Canada, pp. 252259.

^{11}LABADIE N & PRINS C. 2012. Vehicle Routing Nowadays: Compact Review and Emerging Problems. In: Production System and Supply Chain Management in Emerging Countries: Best Practices, Selected Papers From the International Conference on Production Research (ICPR). Springer, Berlin, pp. 141166.

^{12}LAPORTE G, LOUVEAUX FV & MERCURE H. 1994. A Priori Optimization of the Probabilistic Traveling Salesman Problem. Operations Research , 42: 543549.

^{13}LAPORTE G, TOTH P & VIGO D. 2013. Vehicle Routing: Historical Perspective and Recent Contributions. European Journal of Operational Research , 1: 14.

^{14}LI Y, MIAO Q & WANG X. 2014. HighSpeed Train Network Routing with Column Generation. Transportation Research Record: Journal of the Transportation Research Board, 2466: 5867.

^{15}LI Y, MIAO Q & WANG X. 2010. A Column Generation Method for the U.S. Army Logistics Air Fleet Scheduling. In: Transportation Research Record: Journal of the Transportation Research Board , 2197: 3642.

^{16}LI Y, MIAO Q & WANG X. 2012. U.S. Postal Airmail Routing Optimization. In: Transportation Research Record: Journal of the Transportation Research Board , 2234: 2128.

^{17}POWELL WB, JAILLET P & ODONI A. 1995. Stochastic and Dynamic Networks and Routing. BALL MO, MAGNANTI TL, MONMA CL & NEMHAUSER GL. (eds.). Handbooks in Operations Research and Management Science, 8. Network Routing. Elsevier (NorthHolland), Amsterdam, pp. 141296.

^{18}POWELL WB, SNOW W & CHEUNG RK. 2000. Adaptive Labeling Algorithms for the Dynamic Assignment Problem. Transportation Science , 34: 5066.

^{19}SOLOMON MM. 1987. Algorithms for the Vehicle Routing and Scheduling Problem with Time Window Constraints. Operations Research , 35: 254265.

^{20}STEWART WR & GOLDEN BL. 1983. Stochastic Vehicle Routing: A Comprehensive Approach. European Journal of Operational Research , 14: 371385.

^{21}YANG J, JAILLET P & MAHMASSANI H. 2004. RealTime Multivehicle Truckload Pickup and Delivery Problems. Transportation Science , 38(2): 135148.

^{22}YANG J, JAILLET P & MAHMASSANI HS. 1998. Online Algorithms for Truck Fleet Assignment and Scheduling Under Realtime Information. Transportation Research Record, 1667: 107113.
Publication Dates

Publication in this collection
MayAug 2018
History

Received
29 Apr 2017 
Accepted
17 Mar 2018