Flexible truckload pickup and delivery problem considering reserved orders and fuel consumption

Journal Pre-proofs Flexible truckload pickup and delivery problem considering reserved orders and fuel consumption Min Huang, Hui Zhang, Hanbin Kuang, Yang Yu, Loo Hay Lee, Xingwei Wang PII: DOI: Reference:

S0360-8352(19)30586-8 https://doi.org/10.1016/j.cie.2019.106117 CAIE 106117

To appear in:

Computers & Industrial Engineering

Received Date: Revised Date: Accepted Date:

21 November 2018 6 September 2019 6 October 2019

Please cite this article as: Huang, M., Zhang, H., Kuang, H., Yu, Y., Hay Lee, L., Wang, X., Flexible truckload pickup and delivery problem considering reserved orders and fuel consumption, Computers & Industrial Engineering (2019), doi: https://doi.org/10.1016/j.cie.2019.106117

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Ltd.

Flexible truckload pickup and delivery problem considering reserved orders and fuel consumption Min Huang a, Hui Zhang a, Hanbin Kuang* a, Yang Yu a, Loo Hay Lee b, Xingwei Wang c

a. College of Information Science and Engineering ； State Key Laboratory of Synthetical Automation for Process Industries, Northeastern University, Shenyang, Liaoning, 110819, China b. Department of Industrial and Systems Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore c. College of Computer Science and Engineering, Northeastern University, Shenyang, Liaoning, 110169, China Corresponding author: [email protected]

Acknowledgements This work is supported by the NSFC Major International (Regional) Joint Research Project Grant No. 71620107003; the Program for Liaoning Innovative Research Team in University Grant No. LT2016007; the Fundamental Research Funds for State Key Laboratory of Synthetical Automation for Process Industries Grant No. 2013ZCX11.

Abstract: The orders of a transportation enterprise can be fulfilled by the private fleet or outsourced to a common carrier for a flexible transportation planning. In reality, some of the orders have to be reserved for the private fleet due to some constraints (such as trade secret), which can reduce the flexibility of transportation, thereby increasing the operational cost. To handle this problem, a flexible truckload pickup and delivery problem is proposed, where time widows, outsourcing, reserved orders, and fuel consumption of the private fleet are considered, and truck speed is taken as a decision variable instead of a parameter. Then, the proposed problem is described as a directed graph, and its attributes are analyzed to reduce the solution space; consequently, a mixed integer nonlinear programming model based on truck visiting time is established. After that, an updated time window partition based (UTWPB) method is developed, which is described by a sparse matrix to further reduce the solution space, and it can obtain a lower bound of the proposed problem by relaxing the constraints of the partitioned time windows. Numerical experiments demonstrate that the UTWPB method performs well compared with existing approximate methods, and it can offer 1

reasonably accurate solutions for real-life scale instances. Additionally, sensitivity analysis is also conducted to confirm the effectiveness of the proposed model. Keywords ： Truckload pickup and delivery problem; flexible; reserved order, fuel consumption, time window partition, optimizing speed

1 Introduction The vehicle routing problem with private fleet and common carrier (VRPPC) (Chu, 2005) is an extension of the classical VRP (Dantzig and Ramser, 1959), in which, a transportation enterprise can outsource some of its orders to a common carrier, to make the transportation planning feasible or efficient when its capacity is insufficient, or some orders are unprofitable. In particular, to meet customers’ service level or to protect trade secret, the transportation enterprise may have to reserve some orders for its private fleet, which can increase the complexity of solving the problem and the operational cost (Ziebuhr and Kopfer, 2016). Generally speaking, the operational cost of the VRPPC includes the charges by using a common carrier and the cost of the private fleet. Specifically, the former depends on outsourced orders, while the latter is mainly determined by fuel consumption, which is affected by transportation distance, truck load, and truck speed (Bektaş and Laporte, 2011). Although the VRPPC has been well studied in many papers, rare research considers the factor of truck speed in the operational cost. As an important factor influencing transportation capacity and fuel consumption, truck speed can be determined to reduce the operational cost of the VRPPC. Firstly, truck speed can affect the transportation capacity of the private fleet when orders must be served within given time windows; and secondly, truck speed can affect marginal profits of orders. Specifically, a truck with a higher speed may avoid the violations of time windows of orders, and increase the feasibility of trucks’ routes (Kramer et al., 2015), thereby fulfilling more orders within a given time. Whereas, if truck speed is set to be too high, it may result in too much fuel consumption (Bektaş and Laporte, 2011). Consequently, the transportation enterprise needs to design a proper truck speed for its private fleet, and it is meaningful to treat the truck speed as a decision variable to reduce the operational cost. However, when truck speed is treated as a decision variable for fuel consumption cost, the problem is difficult to solve because of its nonlinear objective function (Demir et al., 2012). The VRPPCs can be mainly divided into less truckload VRPPCs and truckload VRPPCs, but only a little research focuses on truckload VRPPCs. Because of the differences between orders, the models and effective algorithms for less truckload VRPPCs cannot be used directly into truckload VRPPCs. Considering that many road transportation problems are truckload transportation, such as container drayage problem (Zhang et al., 2010) and timber transportation problem (Derigs et al., 2012), in this paper, we study a flexible truckload pickup and delivery problem with reserved orders and fuel consumption (FTPDP-RF), where both time widows and outsourcing are considered, and truck speed is taken as a decision variable instead of a parameter. The main contributions of this paper are threefold: (1) We propose the FTPDP-RF, which is an extension of the traditional truckload pickup and delivery problem (TPDP); (2) We establish a mixed integer nonlinear programing (MINLP) model for the FTPDP-RF 2

through identifying trucks’ visiting time; (3) We propose an updated time window partition based (UTWPB) method to cope with the characteristics of the FTPDP-RF and to linearize the MINLP model, as well as offering a lower bound through relaxing the constraints of partitioned time windows. The research is organized as follows: Section 2 presents a brief literature review and justifies our contributions; Section 3 describes the FTPDP-RF by using a directed graph, analyzes the attributes of truck speed, truck travel time, time windows of orders, and then establishes a MINLP model based on truck visiting time; Section 4 offers a linear approximation and a lower bound for the MINLP model, and then introduces a UTWPB method based on an iterative framework; Section 5 designs numerical experiments to evaluate the proposed model and the solution method; moreover, Section 6 provides conclusions and research suggestions for future work.

2 Literature review Considering that the proposed FTPDP-RF is relevant to the VRPPC and the green vehicle routing problem (GVRP), a brief literature review is conducted to provide an investigation on previous work, and to justify our contributions.

2.1 VRPPC The VRPPC is proposed by Chu (2005). In this problem, a transportation enterprise needs to select a sub set of its orders for a common carrier to reduce its operational cost. According to types of orders, the literature on VRPPC can be divided into two categories: less truckload VRPPC and truckload VRPPC. The research on less truckload VRPPC has two main branches: one is to develop an efficient algorithm to solve the classical less truckload VRPPC, such as saving algorithm (Bolduc et al., 2007; Chu, 2005), tabu search algorithm (Côté and Potvin, 2009; Potvin and Naud, 2011), genetic algorithm (Kratica et al., 2012), evolutionary algorithm (Euchi et al., 2011), hybrid meta-heuristic algorithm (Euchi, 2017); the other is to update the classical VRPPC model by introducing extra features, such as heterogeneous vehicles (Bolduc et al., 2008), less truckload pickup and delivery orders with time windows (Krajewska and Kopfer, 2009; Ziebuhr and Kopfer, 2016), and many other outsourcing modes(Gahm et al., 2017; Krajewska and Kopfer, 2009; Ziebuhr and Kopfer, 2016). For outsourcing modes, Krajewska and Kopfer (2009) proposed three main types of outsourcing modes that were paid on distance basis, daily basis, and a function of distance and freight weight basis, respectively; Gahm et al. (2017) proposed an additional outsourcing mode where subcontractors can offer volume discounts; Ziebuhr and Kopfer (2016) considered both long-term subcontractors and common carriers in their problem, and they introduced three types of orders including reserved orders, when trade secret is considered. Compared with less truckload VRPPC, there is less research focusing on truckload VRPPC. Among them, Liu et al. (2010) considered both the outsourcing cost and contracting earning in their objective function, and they proposed a memetic algorithm to solve their

3

problem. However, they did not consider reserved orders, and they just took the truck speed as a parameter. As can be seen, rare research about VRPPC discusses the effect of truck speed on the operational cost, and a little research considers the reserved orders. To fill the gap, we present the FTPDP-RF, which has three features: (1) the objective is a nonlinear function; (2) the transportation enterprise needs to determine average truck speed of the private fleet on each route segment; (3) some orders are reserved. These features lead the FTPDP-RF to be difficult to solve.

2.2 GVRP GVRP is a variant of the VRPs, compared with which, GVRP considers the factor of carbon emissions (or fuel consumption) in the objective function. According to the setting of vehicle speed, GVRP can be classified into two categories: one takes vehicle speeds as parameters (GVRP-p), and the other takes vehicle speeds as decision variables (GVRP-v). The objective of the GVRP-p can be modeled as a linear function (Kopfer et al., 2014). The objective of the GVRP-v, however, is a nonlinear function (Bektaş and Laporte, 2011). The proposed FTPDP-RF is a variant of GVRP-v. There are three main groups of literature on solving variants of GVRP-v: approximate solutions (Bektaş and Laporte, 2011; Franceschetti et al., 2013), meta-heuristic solutions (Demir et al., 2012, 2014; Franceschetti et al., 2017; Jabali et al., 2012; Koç et al., 2014; Kramer et al., 2015), and exact solutions (Dabia et al., 2017). Approximate solutions are usually used to solve small size instances. Bektaş and Laporte (2011) proposed the well-known pollution routing problem (PRP) using a comprehensive fuel consumption model, and adopted a speed discretization based (SDB) method to handle the small scale instances with up to 20 customers. Franceschetti et al. (2013) employed the SDB method to solve and analyze the PRP considering traffic congestion with up to 20 customers, too. Meta-heuristic solutions can be used to handle large size instances. Demir et al. (2012), Demir et al. (2014), and Koç et al. (2014) presented adaptive large neighbor search meta-heuristics to solve different variants of PRP, namely, PRP considering low truck speed, bi-objective PRP, and fleet size and mix PRP; the sizes of instances are up to 200 customers. Jabali et al. (2012) analyzed the travel time and carbon emissions of trucks in traffic congestion for determining optimal truck routes, but they did not consider time windows of customers; they provided a tabu search to solve the instances with up to 80 customers. Kramer et al. (2015) proposed a departure time and speed optimization algorithm and embedded it into an iterated local search to solve the PRP with up to 200 customers. Exact solutions are often used for special problems. Dabia et al. (2017) proposed a branch and price algorithm to solve a special case of the original PRP, where the speed was a decision variable, but it is set to be the same for any route; the instances are up to 100 customers. In this paper, the setting of truck speed is similar as that in Demir et al. (2012). The truck speed can be determined during a given range on any arc. Instead of determining the truck speed directly as that in the existing literature, we treat the visiting time for customers as

4

decision variables to determine each truck’s travel time, thereby determining the truck’s speed.

3 Problem description and formulation This section describes and formulates the FTPDP-RF. The FTPDP-RF is formulated by using a directed graph, and interpretations are given by using an example in Subsection 3.1. A fuel consumption model and its transformation are presented in Subsection 3.2. Subsection 3.3 investigates the attributes of the FTPDP-RF to reduce solution space. Subsection 3.4 develops a mixed integer nonlinear programming model.

3.1 Description of the FTPDP-RF In the FTPDP-RF, a transportation enterprise has some orders and some trucks that are parked in its depot. It is assumed that each order must be picked up at its origin and directly delivered to its destination with a full truckload request during the corresponding time windows. The orders have two types: reserved orders and selective orders. For reserved orders, they must be fulfilled by the private fleet; for selective orders, they can either be fulfilled by the private fleet or be outsourced to a common carrier. In addition, if an order is fulfilled by the private fleet, the enterprise needs to pay for the fuel consumption cost; otherwise, service fee needs to be paid to the common carrier. The operational cost of the enterprise is the sum of the fuel consumption cost and the outsourcing charges. The objective of the FTPDP-RF is to minimize the operational cost by selecting orders to fulfill and optimizing the truck speed as well as routes for the private fleet. Inspired by Zhang et al. (2009) and Gendreau et al. (2015), the FTPDP-RF can be described with a directed graph 𝐺 = (𝑁0, 𝐴), where 𝑁0 is the set of nodes, and 𝐴 represents the set of arcs. 𝑁0 = {0} ∪ 𝑁, and {0} ∩ 𝑁 = ∅, where 0 is the depot and 𝑁 is the set of order nodes. 𝐴 = {(𝑖,𝑗)|𝑖 ∈ 𝑁0， 𝑗 ∈ 𝑁0, 𝑖 ≠ 𝑗}. The depot 0 is a start/return node for parking 𝐾 homogenous trucks. The working time for each truck is within [0, 𝐻]. An order node 𝑖 ∈ 𝑁 represents the fulfillment of order 𝑖 by loading from origin point 𝑂𝑖, traveling to destination point 𝐷𝑖, and unloading. The time windows of 𝑂𝑖 and 𝐷𝑖 are [𝑎𝑂𝑖 , 𝑏𝑂𝑖](𝑎𝑂𝑖 ≤ 𝑏𝑂𝑖) and [𝑎𝐷𝑖, 𝑏𝐷𝑖] (𝑎𝐷𝑖 ≤ 𝑏𝐷𝑖), respectively. A truck has to wait if it arrives at the origin or the destination ahead of corresponding time windows. The loading time at 𝑂𝑖 is 𝑠𝑂𝑖, and the unloading time at 𝐷𝑖 is 𝑠𝐷𝑖. The payload distance is 𝑑𝑖 = 𝐷𝑖𝑠(𝑂𝑖,𝐷𝑖), where 𝐷𝑖𝑠 (𝑂𝑖,𝐷𝑖) represents the distance from 𝑂𝑖 to 𝐷𝑖. An arc (𝑖,𝑗) ∈ 𝐴, where 𝑖,𝑗 ∈ 𝑁0, represents the traveling from 𝐷𝑖 to 𝑂𝑗, and the travel distance is 𝑑𝑖𝑗 = 𝐷𝑖𝑠(𝐷𝑖,𝑂𝑗). For the depot, we define 𝐷0=𝑂0=0. Since the FTPDP-RF considers fuel consumption and outsourcing, some extra attributes are needed compared with the traditional TPDP. Specifically, for a reserved order 𝑖 ∈ 𝑁, 𝑅𝑒𝑠𝑖 = 1, and for a selective order 𝑖 ∈ 𝑁, 𝑅𝑒𝑠𝑖 = 0. An outsourced order 𝑖 ∈ 𝑁 will be charged 𝜃𝑓𝑖𝑑𝑖 by the common carrier, where 𝜃 is the unit price of outsourcing, and 𝑓𝑖 is the payload of order 𝑖; the cost of order 𝑖 ∈ 𝑁 conducted by the private fleet is determined by the fuel consumption. The unite price of fuel consumption per liter is denoted by 𝛿; the total weight of a truck travelling from 𝑂𝑖 to 𝐷𝑖 is (𝑤 + 𝑓𝑖), where 𝑤 is the weight of an 5

empty truck, and the total weight of the truck travelling from 𝐷𝑖 to 𝑂𝑗 is 𝑤 (the truck is unloaded), and the range of truck speed on each road segment is [𝑣𝐿, 𝑣𝑈](0 < 𝑣𝐿 ≤ 𝑣𝑈). Note that a road segment indicates a link connecting any two separate locations among the origins, the destinations, and the depot. Considering the time windows for transportation of orders, an order 𝑖 should satisfy 𝑎𝑂𝑖 + 𝑠𝑂𝑖 + 𝑑𝑖/𝑣𝑈 ≤ 𝑏𝐷𝑖, or the order cannot be fulfilled. From the perspective of the introduced graph, the objective of the FTPDP-RF is to identify a subset of order nodes conducted by the private fleet (the rest of the order nodes are conducted by a common carrier), a subset of arcs connecting the depot and these order nodes (including reserved order nodes), and the truck speed on these nodes and arcs, thereby minimizing the operational cost in a time horizon (one day). To further interpret the FTPDP-RF, an illustrative example is shown in Figure 1, where, the square is the depot, denoted as 0, with a private truck. The enterprise has 3 orders and order 𝑖 (𝑖 = 1, 2, 3) has origin 𝑂𝑖 and destination 𝐷𝑖, and the freight for each order must be transported from its origin to its destination directly. Among the orders, Order 1 is a reserved order (𝑅𝑒𝑠1 = 1), and the others are selective orders (𝑅𝑒𝑠2 = 𝑅𝑒𝑠3 = 0). A solid line with an origin and a destination represents an order node, and a broken line indicates an arc. Then, a feasible solution for this example is that Order 3 is selected to be outsourced, and Order 1 and Order 2 are conducted by the private truck following Route 1. Note that if there is no reserved order, it is better for the private truck to conduct Order 2 and Order 3 considering shorter travelling distance.

D1 O2

Order 1 Order 2

Route 1

O1

D3

D2

Order 3

O3

An order node

0 An arc

Figure 1. An illustrative example

3.2 Calculation of fuel consumption According to the description of Demir et al. (2012), we can estimate the fuel consumption rate as 𝐹(𝑣) = ∆1𝑑/𝑣 + ∆2(𝑤 + 𝑓)𝑑 + ∆3𝑑𝑣2, (1) where the unites of 𝑑, 𝑣, 𝑤 and 𝑓 are kilometer, kilometer/minute, ton, and ton, respectively. According to the given typical values in Demir et al. (2012), we can obtain that 6

∆1=0.061058345 liter/minute, ∆2=0.008403232 liter/(ton ∙ kilometer), and ∆3= 0.039108399 liter ∙ minute2/kilometer3. Since travel time 𝑡 = 𝑑/𝑣 (minute), formula (1) can be rewritten as 𝐹(𝑡) = ∆1𝑡 + ∆2(𝑤 + 𝑓)𝑑 + ∆3𝑑3/(𝑡 + 𝜖)2.

(2)

Note that if 𝑑 = 0, then 𝑡 = 0. In this case, the fraction of formula (2) no longer has any meaning without ϵ. Hence, a very small value ϵ(ϵ > 0) is introduced in formula (2).

3.3 Analysis of the FTPDP-RF The solution space will increase when the truck speed on each road segment is treated as a decision variable. To properly reduce the solution space, the attributes of the FTPDP-RF are investigated. Specifically, the range of truck speed is reduced; the properties of truck travel time are analyzed; the ranges of time windows of orders are reduced; as well as, infeasible arcs are eliminated.

3.3.1 Range of truck speed Demir et al. (2012) proposed a proposition for the unconstrained optimal speed in the process of minimizing the truck fuel consumption on a road segment, and the optimal truck ∗ = 3 ∆1/(2∆3). This proposition holds on the assumption that the speed was defined as 𝑣𝑢𝑐 ∗ ∗ value of 𝑣𝑢𝑐 falls in the range of truck speed. In some cases, however, 𝑣𝑢𝑐 may be out of the range because 𝐹(𝑣) is a convex function of truck speed 𝑣 (kilometer/minute) (𝑣 > 0) (Demir et al., 2012) according to formula (1). For example, in China, the maximum of truck ∗ speed is limited to a low level (may be lower than 𝑣𝑢𝑐 ) on rural and urban roads for traffic safety, while the minimum of truck speed is limited to a relatively high level (may be higher ∗ than 𝑣𝑢𝑐 ) in highway for transportation efficiency. Hence, we extend the proposition as follows. Proposition 3.1. If the truck speed 𝑣 is constrained to fall between 𝑣𝐿 and 𝑣𝑈, the optimal speed 𝑣𝑐∗ for minimizing truck fuel consumption has the following properties: if 𝑣𝑈 ∗ ∗ ∗ ∗ ≤ 𝑣𝑢𝑐 < 𝑣𝑈, then 𝑣𝑐∗ = 𝑣𝑢𝑐 ≤ 𝑣𝐿, then 𝑣𝑐∗ = 𝑣𝐿. , then 𝑣𝑐∗ = 𝑣𝑈; if 𝑣𝐿 < 𝑣𝑢𝑐 ; and if 𝑣𝑢𝑐 Proof of Proposition 3.1. See the appendix. Based on Proposition 3.1, Proposition 3.2 is introduced to reduce the range of truck speed. Proposition 3.2. Given the description as shown in Subsection 3.1 and Subsection 3.2, the truck speed on each road segment of the optimal solution of the FTPDP-RF must fall in the range of [𝑣𝑐∗ , 𝑣𝑈]. Proof of Proposition 3.2. See the appendix. Consequently, the range of truck speed [𝑣𝐿,𝑣𝑈 ] can be rewritten as [𝑣𝐿𝑟, 𝑣𝑈𝑟], where,

7

∗ 𝑣𝑈, if 𝑣𝑈 ≤ 𝑣𝑢𝑐 ∗ ∗ , if 𝑣𝐿 ≤ 𝑣𝑢𝑐 < 𝑣𝑈, 𝑣𝐿𝑟 = 𝑣𝑐∗ = 𝑣𝑢𝑐 𝐿 ∗ 𝑣, if 𝑣𝑢𝑐 < 𝑣𝐿

{

(3)

𝑣𝑈𝑟 = 𝑣𝑈.

(4)

3.3.2 Range of truck travel time Similar as the truck speed, the travel time also has some properties. Proposition 3.3. Given the description as shown in Subsection 3.2, a truck travels 𝑑 (kilometer). Then, 𝐹(𝑡) is a convex function of travel time 𝑡 (𝑡 > 0) (minute) and the ∗ ∗ = 𝑑/𝑣𝑢𝑐 optimal travel time for minimizing truck fuel consumption is 𝑡𝑢𝑐 . Proof of Proposition 3.3. See the appendix. Following Proposition 3.3, Corollary 3.1 is introduced as follows. ∗ Corollary 3.1. If 𝑡 ≥ 𝑡𝑢𝑐 , then 𝐹(𝑡) is a monotone increasing function, whereas if 𝑡 ∗ ≤ 𝑡𝑢𝑐 , then 𝐹(𝑡) is a monotone decreasing function. Proposition 3.3 holds on the assumption that 𝑡 > 0. When 𝑡 ∈ [𝑑/𝑣𝑈𝑟, 𝑑/𝑣𝐿𝑟] , the following propositions are introduced. Proposition 3.4. Given the description as shown in Subsection 3.2, a truck travels 𝑑 (kilometer). Then, 𝐹(𝑡) is a monotone non-increasing function of 𝑡 ∈ [𝑑/𝑣𝑈𝑟, 𝑑/𝑣𝐿𝑟] (minute), and the optimal travel time 𝑡𝑐∗ for minimizing truck fuel consumption is 𝑡𝑐∗ = 𝑑/ 𝑣𝐿𝑟. Proof of Proposition 3.4. See the appendix. Proposition 3.5. When the travel time falls in the range of [0, 𝑇], the optimal travel time ∗ 𝑡𝑐 follows: if 𝑇 ≥ 𝑑/𝑣𝐿𝑟, then 𝑡𝑐∗ = 𝑑/𝑣𝐿𝑟; if 𝑑/𝑣𝑈𝑟 ≤ 𝑇 < 𝑑/𝑣𝐿𝑟, then 𝑡𝑐∗ = 𝑇; if 𝑇 < 𝑑/𝑣𝑈𝑟, 𝑡𝑐∗ is nonexistent. Proof of Proposition 3.5. See the appendix.

3.3.3 Time windows of orders Time windows of truckload orders can be reduced by proper methods when truck speeds are treated as parameters (Jula et al., 2005; Zhang et al., 2010). However，when truck speeds on road segments are considered as decision variables, these methods can no longer be directly used. Hence, a modified method is developed in this research. Let 𝑡𝑂𝐷 be the travel time from origin 𝑂𝑖 to destination 𝐷𝑖 of order node 𝑖, 𝑖 ∈ 𝑁. 𝑖 Note that 𝑠𝑂𝑖 represents the loading time at 𝑂𝑖. If the time window of 𝐷𝑖 is shifted back in 𝑂 𝐷 𝑂 𝑂𝐷 𝐷 𝑂 𝑂𝐷 time by 𝑡𝑂𝐷 𝑖 + 𝑠𝑖 , the range of time window becomes [𝑎𝑖 ― 𝑠𝑖 ― 𝑡𝑖 ，𝑏𝑖 ― 𝑠𝑖 ― 𝑡𝑖 ]. As shown in Sub-subsection 3.3.1, the range of truck speed is [𝑣𝐿𝑟, 𝑣𝑈𝑟]. Then, 𝑑𝑖/𝑣𝑈𝑟 ≤ 𝑡𝑂𝐷 𝑖 ≤ 𝐿𝑟 𝑑𝑖/𝑣 , and the range of the shifted time window of 𝐷𝑖 becomes 𝐷 𝑂 𝐿𝑟 𝐷 𝑂 𝑈𝑟 𝐷 [𝑎𝑖 ― 𝑠𝑖 ― 𝑑𝑖/𝑣 , 𝑏𝑖 ― 𝑠𝑖 ― 𝑑𝑖/𝑣 ]. If a truck starts at 𝑂𝑖 before the time of 𝑎𝑖 ― 𝑠𝑂𝑖 ― 𝑑𝑖/𝑣𝐿𝑟, the truck has to wait at 𝐷𝑖, and the fuel consumption of the truck will remain unchanged; if a truck starts at 𝑂𝑖 after the time of 𝑏𝐷𝑖 ― 𝑠𝑂𝑖 ― 𝑑𝑖/𝑣𝑈𝑟, the time window of 𝐷𝑖 will be violated. To avoid the unnecessary waiting time and the violations of time windows, the overlap of the shifted time window of 𝐷𝑖 and the time window of 𝑂𝑖 is taken as the new time window of 𝑂𝑖. After the transformation, the time window of 𝑂𝑖, [𝑎𝑂𝑖,𝑏𝑂𝑖], can be rewritten as 𝑂𝑟 [𝑎𝑂𝑟 𝑖 ,𝑏𝑖 ], where,

8

{ {

𝐷 𝑂 𝑎𝑂𝑟 𝑖 = min max 𝑎𝑖 ― 𝑠𝑖 ―

{

𝐷 𝑂 𝑏𝑂𝑟 𝑖 = min 𝑏𝑖 ― 𝑠𝑖 ―

,𝑎𝑂𝑖 ,𝑏𝑂𝑖 , 𝑣𝐿𝑟

} }

(5)

}

(6)

𝑑𝑖

𝑑𝑖

,𝑏𝑂𝑖 . 𝑣𝑈𝑟

𝐷𝑟 Similarly, [𝑎𝐷𝑖,𝑏𝐷𝑖] of order 𝑖 ∈ 𝑁 can be rewritten as [𝑎𝐷𝑟 𝑖 ,𝑏𝑖 ], where,

{

𝑂𝑟 𝑂 𝑎𝐷𝑟 𝑖 = max 𝑎𝑖 + 𝑠𝑖 +

{

}

𝑑𝑖

,𝑎𝐷𝑖 , 𝑣𝑈𝑟

{

𝐷 𝑂𝑟 𝑂 𝑏𝐷𝑟 𝑖 = max 𝑎𝑖 ,min 𝑏𝑖 + 𝑠𝑖 +

𝑑𝑖

,𝑏𝐷𝑖 𝑣𝐿𝑟

(7)

}}

.

(8)

3.3.4 Graph arcs 𝐷 𝑂𝑟 𝑑 𝑈𝑟 Based on the description above, if 𝑎𝐷𝑟 𝑖 + 𝑠𝑖 + 𝑖𝑗 𝑣 > 𝑏𝑗 , arc(𝑖,𝑗) ∈ 𝐴 is infeasible. 𝐴 Hence, arc set can be rewritten as 𝐷 𝑈𝑟 𝑂𝑟} {(𝑖,𝑗)|𝑖 = 0, 𝑗 ∈ 𝑁, or 𝑖 ∈ 𝑁, 𝑗 = 0, or 𝑖 ∈ 𝑁, 𝑗 ∈ 𝑁, 𝑖 ≠ 𝑗, 𝑎𝐷𝑟 + 𝑠 + 𝑑 /𝑣 ≤ 𝑏 . After 𝑖 𝑖 𝑖𝑗 𝑗 that, arc set 𝐴 can be described as a sparse matrix to reduce the decision variables.

3.4 Mixed integer nonlinear programming model According to the formulation above, truck speed and truck travel time can be converted to each other. Therefore, we can determine truck speed by determining truck travel time. In this subsection, a mixed integer nonlinear programming model is established for the FTPDP-RF. Decision variables are listed as follows: 1,𝑖𝑓 𝑎𝑟𝑐 （𝑖,𝑗） ∈ 𝐴 is selected in a feasible solution, 𝑥𝑖𝑗 = 0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒,

{

𝑢𝑂𝑖 ∈ 𝑹, the start time of loading at 𝑂𝑖 of order 𝑖 ∈ 𝑁, 𝑢𝐷𝑖 ∈ 𝑹, the start time of unloading at 𝐷𝑖 of order 𝑖 ∈ 𝑁. Moreover, three auxiliary variables are introduced to provide a clear description: 1, 𝑖𝑓 𝑜𝑟𝑑𝑒𝑟 𝑖 ∈ 𝑁 is selected to be outsourced, 𝑦𝑖 = 0, 𝑖𝑓 𝑜𝑟𝑑𝑒𝑟 𝑖 ∈ 𝑁 is selected to be fulfilled by a private fleet,

{

𝑡𝑖𝑗 ∈ 𝑹, the travel time of a private truck on arc （𝑖,𝑗） ∈ 𝐴, 𝑡𝑖 ∈ 𝑹, the travel time of a private truck fulfilling order 𝑖 ∈ 𝑁. Based on formula (2) and the directed graph as shown in Subsection 3.1, Model-T is established as: Minimize

∑(

𝛿 ∆1𝑡𝑖 + ∆2(𝑤 + 𝑓)𝑑𝑖 + ∆3

∑ （𝑖,𝑗） ∈ 𝐴

9

) )

(𝑡𝑖 + 𝜖)2

𝑖∈𝑁

+

𝑑3𝑖

(

𝛿 ∆1𝑡𝑖𝑗 + ∆2𝑤𝑑𝑖𝑗 + ∆3

𝑑3𝑖𝑗

(𝑡𝑖𝑗 + 𝜖)2

(1 ― 𝑦𝑖)

𝑥𝑖𝑗

(9)

(10)

+

∑ 𝜃𝑓 𝑑 𝑦

(11)

𝑖 𝑖 𝑖

𝑖∈𝑁

Subject to

∑ 𝑥 ≤ 𝐾, ∑ 𝑥 = ∑ 𝑥 , ∀𝑖 ∈ 𝑁 , 𝑦 = 1 ― ∑ 𝑥 , ∀𝑖 ∈ 𝑁, ∑ 𝑥 = 1, ∀𝑖 ∈ 𝑁 and 𝑅𝑒𝑠 = 1, 0𝑗

(12)

𝑗∈𝑁

𝑖𝑗

(𝑖,𝑗) ∈ 𝐴

0

𝑗𝑖

(13)

(𝑗,𝑖) ∈ 𝐴

𝑖

𝑖𝑗

(14)

(𝑖,𝑗) ∈ 𝐴

𝑖𝑗

𝑖

(15)

(𝑖,𝑗) ∈ 𝐴 𝑂 𝑂𝑟 𝑎𝑂𝑟 (16) 𝑖 ≤ 𝑢𝑖 ≤ 𝑏𝑖 , ∀𝑖 ∈ 𝑁, 𝐷𝑟 𝐷 𝐷𝑟 𝑎𝑖 ≤ 𝑢𝑖 ≤ 𝑏𝑖 , ∀𝑖 ∈ 𝑁, (17) 𝐷 𝑂 𝑢0 = 0, and 𝑢0 = 𝐻 (18) 𝐷 𝑂 𝑂 𝑡𝑖 ≤ 𝑢𝑖 ― 𝑢𝑖 ― 𝑠𝑖 , ∀𝑖 ∈ 𝑁 (19) 𝑂 𝐷 𝐷 (𝑡𝑖𝑗 ― (𝑢𝑗 ― 𝑢𝑖 ― 𝑠𝑖 ))𝑥𝑖𝑗 ≤ 0,∀(𝑖, 𝑗) ∈ 𝐴 (20) 𝑑𝑖 𝑑𝑖 ≤ 𝑡𝑖 ≤ 𝐿𝑟, ∀𝑖 ∈ 𝑁 (21) 𝑣𝑈𝑟 𝑣 𝑑𝑖𝑗 𝑑𝑖𝑗 ≤ 𝑡 ≤ , ∀(𝑖, 𝑗) ∈ 𝐴 (22) 𝑖𝑗 𝑣𝑈𝑟 𝑣𝐿𝑟 𝑥𝑖𝑗 ∈ {0，1}, 𝑡𝑖𝑗 is a positive real variable, ∀(𝑖,𝑗) ∈ 𝐴, (23) 𝑂 𝐷 𝑦𝑖 ∈ {0，1},𝑡𝑖,𝑢𝑖 , 𝑢𝑖 are positive real variables, ∀𝑖 ∈ 𝑁. (24) The objective function is to minimize the operational cost, where terms (9) and (10) represent the fuel consumption cost of the private fleet for payload and deadhead kilometers, respectively, and term (11) represents the outsourcing charges. Constraint (12) limits the number of trucks of the private fleet. Constraint (13) balances the truck flow and ensures that each vehicle starts and ends at the depot, and the number of trucks entering an order node is equal to that exiting the node. Constraint (14) ensures that each order must be served exactly once, either by the private fleet or by the common carrier. Constraint (15) makes sure that reserved orders must be served by the private fleet. Constraints (16) and (17) guarantee that the visiting time of the origin and the destination must fall in the corresponding time windows, respectively. For depot 𝑖 = 0, we define 𝑢𝐷0 = 0 and 𝑢𝑂0 = 𝐻. Constraint (18) implies the working hours of a truck cannot exceed 𝐻. Constraints (19) and (20) indicate that the visiting time of any two locations must satisfy the continuity requirement, and the vehicle cannot serve any order more than once, thereby eliminating the sub tours. Constraints (21) and (22) limit the truck speed. Constraints (23) and (24) are the ranges of decision variables. It is clear that Model-T is a mixed integer nonlinear programming model. Due to constraints (19) and (21), for a feasible solution, we have 𝑑𝑖/𝑣𝑈𝑟 ≤ 𝑡𝑖 ≤ min {𝑑𝑖/𝑣𝐿𝑟,𝑢𝐷𝑖 ― 𝑢𝑂𝑖 ― 𝑠𝑂𝑖 }, ∀𝑖 ∈ 𝑁 Similarly, 𝑑𝑖𝑗/𝑣𝑈𝑟 ≤ 𝑡𝑖𝑗 ≤ min {𝑑𝑖𝑗/𝑣𝐿𝑟,𝑢𝑂𝑗 ― 𝑢𝐷𝑖 ― 𝑠𝐷𝑖 }, ∀(𝑖, 𝑗) ∈ 𝐴 According to proposition 3.5, for a feasible solution, the optimal truck travel time on the nodes and arcs should satisfy formulas (25) and (26), respectively

10

{ {

𝑡𝑖 = max min

{ {

𝑡𝑖𝑗 = max min

𝑑𝑖

𝑑𝑖𝑗

} } 𝑑𝑖

,𝑢𝐷𝑖 ― 𝑢𝑂𝑖 ― 𝑠𝑂𝑖 ,𝑣𝑈𝑟 , ∀𝑖 ∈ 𝑁

𝑣𝐿𝑟

} } 𝑑𝑖𝑗

,𝑢𝑂𝑗 ― 𝑢𝐷𝑖 ― 𝑠𝐷𝑖 , 𝑣𝑈𝑟 , ∀(𝑖, 𝑗) ∈ 𝐴

𝑣𝐿𝑟

(25) (26)

Remark 3.1: By substituting formula (14) into term (11), variable 𝑦𝑖 can be eliminated; by substituting formulas (25) and (26) into Model-T, variables 𝑡𝑖 and 𝑡𝑖𝑗, constraints (21) and (22) can be eliminated. Thus, we only need to determine the optimal arcs and the truck visiting time at nodes. Note that the truck speed on node 𝑖 ∈ 𝑁, 𝑣𝑖, can be determined by 𝑣𝑖 = 𝑑𝑖 𝑡𝑖, and similarly, the truck speed on arc (𝑖, 𝑗) ∈ 𝐴 can be determined by 𝑣𝑖𝑗 = 𝑑𝑖𝑗/𝑡𝑖𝑗.

4 Time window partition based method Because Model-T is a mixed integer nonlinear model, it is difficult to be solved directly by commercial optimization software. Hence, we adopt an updated time window partition based (UTWPB) method to find the near-optimal arcs and the start working time at nodes. Simultaneously, a lower bound is provided by relaxing the sub time windows constraints. The time window partition based (TWPB) method is initially proposed by Wang and Regan (2002) to linearize the mathematical model of a truckload pickup and delviery problem with time windows. The TWPB method is an iterative framework and consists of an upper bound model and a lower bound model. At each iteration, the time windows of orders are re-partitioned by a predefined width, and this procedure repeats until a stop condition is satisfied. Afterwards, Wang and Regan (2009) proved the covergence of the TWPB method: when the width of patitioned time windows trends to zero, the optimal solution of the lower bound model or the upper bound model trends to be the optimal solution of the original problem. Moreover, Zhang et al. (2010) pointed out that each iteration of the TWPB method can be implemented individually due to the independency of the iterations. In general, the smaller the width of sub time windows is, the more accurate solution can be obtained; accordingly, the more computer memory is required, and the more computational time will be consumed. To reduce the computational time, Zhang et al. (2010) estimated the minimum width of sub time windows in advance, and then solved the upper bound model and lower bound model only once. The numerical experiments for the static and dynamic container drayage problems(Zhang et al., 2014; Zhang et al., 2010) confirm the effectiveness of the modified strategy. Because the opitmal visiting time must fall in a sub time window, it is not necessary to partition all parts of the time widows if we can find out the most potential range where the optimal visiting time falls. Following this idea, we proposed the UTWPB method which is also an iterative framework. Additionally, we use a sparse matrix to describe the UTWPB method to further reduce the solution space. The structure of this section is orgized as follows. Subsection 4.1 presents the methods for partitioning time windows, and offers a equlivent model of Model-T, namely, Model-S. Subsectin 4.2 provides an upper bound for Model-S, namely, Model-UB. Subsection 4.3 estimates a lower bound for Model-S, namely, Model-LB. Subsection 4.4 gives an example

11

for illustrating Model-UB and Model-LB. An iterative framework of the UTWPB method is introduced in Subsection 4.5.

4.1 Processes of time window partition Two partition methods are introduced in this subsection. One is a constant width method, similar as the existing literature (Wang and Regan, 2002; Zhang et al., 2014; Zhang et al., 2010), and the other is a bisection method. The constant width method is described as follows. Given a length Λ, we can partition the time windows for origins and destinations of all the orders. Then, two groups of sub time windows can be generated: one is for the origins, and the other is for the destinations. After that, an additional sub time window, whose width is zero, is added to each group. Specifically, 𝑂𝑟 𝑂𝑟 𝑂𝑟 𝑂𝑟 𝑂𝑟 for the origin of order 𝑖, 𝑖 ∈ 𝑁, [𝑎𝑂𝑟 𝑖 , 𝑏𝑖 ] can be partitioned to be [𝑎𝑖 ，𝑎𝑖 ], [𝑎𝑖 ，𝑎𝑖 + 𝑂𝑟 𝑂𝑟 𝑂𝑟 𝑂𝑟 𝑂𝑟 Λ],[𝑎𝑂𝑟 𝑖 + Λ,𝑎𝑖 +2Λ], … , [𝑎𝑖 + (⌈(𝑏𝑖 ― 𝑎𝑖 ) Λ⌉ ― 1)Λ, 𝑏𝑖 ]. In total,

𝑂𝑟 ⌈(𝑏𝑂𝑟 𝑖 ― 𝑎𝑖 ) Λ⌉+1

of sub time windows can be obtained for the origin of order 𝑖. When each sub time window is marked with a unique identification number, the set of sub nodes of origins, 𝑂𝑠, is formed, 𝑂 𝑂𝑠 and node 𝑗 ∈ 𝑂𝑠 with sub time window [𝑎𝑂𝑠 𝑗 , 𝑏𝑗 ] is a sub origin node of order 𝑖 = 𝛩 (𝑗) ∈ 𝑁. For example, sub node 3 ∈ 𝑂𝑠 and order 𝛩𝑂(3) = 1 ∈ 𝑁 indicate that the time window of 𝑠 𝑂𝑠 sub node 3 ∈ 𝑂𝑠 is [𝑎𝑂𝑠 3 , 𝑏3 ], and the location of sub node 3 ∈ 𝑂 is the origin of order 1; similarly, 𝐷𝑠 is the set of sub nodes of destinations, and node 𝑘 ∈ 𝐷𝑠 with sub time window 𝑂 𝐷𝑠 [𝑎𝐷𝑠 𝑘 , 𝑏𝑘 ] is a sub destination node of order 𝑖 = 𝛩 (𝑘) ∈ 𝑁. The bisection method is described as follows: if the width of a selected sub time window is greater than the given minimum width ∇, the selected sub time window is partitioned into two equal sections; otherwise, the sub time window keeps unchanged. After that, once the selected sub time window is partitioned, it is removed from the set of sub time windows, and the sub-sub time windows generated by the bisection method are added into the set of sub time windows. Note that a sub-sub time window can generate a sub origin node or a sub destination node. After the time windows are partitioned, let 𝐺𝑆 be a directed graph, where 𝐺𝑆 = 𝑈𝑟 𝑂 𝐷𝑠 ({0}, 𝑁𝑆,𝐴𝑆), 𝑁𝑆 = {(𝑖, 𝑗)|𝑖 ∈ 𝑂𝑠，𝑗 ∈ 𝐷𝑠,and 𝛩𝑂(𝑖) = 𝛩𝐷(𝑗); 𝑎𝑂𝑠 𝑖 + 𝑠𝛩𝑂(𝑖) + 𝑑𝛩𝑂(𝑖)/𝑣 ≤ 𝑏𝑗 } and 𝐴𝑆 = {(𝑖,𝑗)|𝑖 = 0, 𝑗 ∈ 𝑂𝑠;or 𝑖 ∈ 𝐷𝑠, 𝑗 = 0;or 𝑖 ∈ 𝐷𝑠,𝑗 ∈ 𝑂𝑠and 𝛩𝐷(𝑖) ≠ 𝛩𝑂(𝑗); 𝑎𝐷𝑠 𝑖 + 𝑈𝑟 𝐷𝐷 𝑂𝑠 𝑠𝛩 (𝑖) + 𝑑𝛩𝐷(𝑖)𝛩𝑂(𝑗)/𝑣 ≤ 𝑏𝑗 }. Arc (𝑖, 𝑗) ∈ 𝑁𝑆 means that a truck starts at origin 𝑂𝛩𝑂(𝑖) during sub time window [𝑎𝑂𝑠 𝑖 , 𝑆 𝑂𝑠 𝐷𝑠 𝐷𝑠 𝑏𝑖 ], and arrives at destination 𝐷𝛩𝐷(𝑗) during sub time window [𝑎𝑗 , 𝑏𝑗 ]. 𝑁 is equivalent to 𝑁. Arc (𝑖, 𝑗) ∈ 𝐴𝑆 means that a truck starts at destination 𝐷𝛩𝐷(𝑖) during sub time window 𝐷𝑠 𝑂𝑠 𝑂𝑠 [𝑎𝐷𝑠 𝑖 , 𝑏𝑖 ], and arrives at origin 𝑂𝛩𝑂(𝑖) during sub time window [𝑎𝑗 , 𝑏𝑗 ]. For the depot, we 𝑆 𝐷𝑠 𝐷𝑠 𝑂𝑠 𝑂𝑠 define 𝑎0 = 𝑏0 = 0, 𝑎0 = 𝑏0 = 𝐻. 𝐴 is equivalent to 𝐴. Based on the description above, decision variables for Model-S are introduced as follows: 𝜔𝑖𝑗 =

solution , {1, 𝑖𝑓 (𝑖,𝑗) ∈ 𝑁 is selected in0,the𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

12

𝑠

𝑥𝑖𝑗 =

solution , {1, 𝑖𝑓 (𝑖,𝑗) ∈ 𝐴 is selected in0,the𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 𝑠

𝑠 𝑢𝑂𝑠 𝑖 , the start time of loading at sub node 𝑖 ∈ 𝑂 , 𝑠 𝑢𝐷𝑠 𝑖 , the start time of unloading at sub node 𝑖 ∈ 𝐷 . Moreover, two auxiliary variables are introduced: 𝑠 𝑡𝑂𝐷 𝑖𝑗 , truck travel time on (𝑖,𝑗) ∈ 𝑁 , 𝑠 𝑡𝐷𝑂 𝑖𝑗 , truck travel time on (𝑖,𝑗) ∈ 𝐴 . Then, Model-S can be established as:

∑

Minimize

𝛿

(𝑖, 𝑗) ∈ 𝑁

∑

+

(

∆1𝑡𝑂𝐷 𝑖𝑗

𝑠

(

∆1𝑡𝐷𝑂 𝑖𝑗

𝛿

（𝑖,𝑗） ∈ 𝐴𝑠

+ ∆2(𝑤 + 𝑓)𝑑𝛩𝑂(𝑖) + ∆3

𝑑3𝛩𝑂(𝑖)

)

2 (𝑡𝑂𝐷 𝑖𝑗 + 𝜖)

+ ∆2𝑤𝑑𝛩 +

𝐷

𝑂

(𝑖)𝛩 (𝑗)

+ ∆3

𝑑3𝛩𝐷(𝑖)𝛩𝑂(𝑗)

)

2 (𝑡𝐷𝑂 𝑖𝑗 + 𝜖)

𝑤𝑖𝑗

𝑥𝑖𝑗

(27)

(28)

∑ 𝜃𝑓 𝑑 𝑦

(29)

𝑖 𝑖 𝑖

𝑖∈𝑁

Subject to

∑𝑥 𝑗∈𝑂

∑

0𝑗

≤ 𝐾,

∑

𝜔𝑗𝑖 =

(30)

𝑠

𝑥𝑖𝑗, ∀𝑖 ∈ 𝐷𝑠,

(31)

（𝑖,𝑗） ∈ 𝐴𝑠

(𝑗,𝑖) ∈ 𝑁𝑠

∑

𝑥𝑗𝑖 =

（𝑗,𝑖） ∈ 𝐴𝑠

∑

∑

𝜔𝑖𝑗, ∀𝑖 ∈ 𝑂𝑠,

(32)

(𝑖,𝑗) ∈ 𝑁𝑠

𝜔𝑖𝑗 + 𝑦𝑘 = 1, ∀𝑘 ∈ 𝑁,

(33)

𝑠

(𝑖,𝑗) ∈ 𝑁 𝛩𝑂(𝑖) = 𝛩𝐷(𝑗) = 𝑘

∑

𝜔𝑖𝑗 = 1, ∀𝑘 ∈ 𝑁, and 𝑅𝑒𝑠𝑘 = 1, 𝑠

(34)

𝑂𝑠 𝑂𝑠 𝑠 𝑎𝑂𝑠 𝑖 ≤ 𝑢𝑖 ≤ 𝑏𝑖 , ∀𝑖 ∈ 𝑂 ,

(35)

𝐷𝑠 𝐷𝑠 𝑠 𝑎𝐷𝑠 𝑖 ≤ 𝑢𝑖 ≤ 𝑏𝑖 ,∀𝑖 ∈ 𝐷 ,

(36)

𝐷𝑠 𝑂𝑠 𝑂 𝑠 (𝑡𝑂𝐷 𝑖𝑗 ― (𝑢𝑗 ― 𝑢𝑖 ― 𝑠𝛩 (𝑖)))𝑤𝑖𝑗 ≤ 0, ∀(𝑖, 𝑗) ∈ 𝑁 ,

(37)

𝑂𝑠 𝐷𝑠 𝐷 𝑠 (𝑡𝐷𝑂 𝑖𝑗 ― (𝑢𝑗 ― 𝑢𝑖 ― 𝑠𝛩 (𝑖)))𝑥𝑖𝑗 ≤ 0,∀(𝑖, 𝑗) ∈ 𝐴 ,

(38)

(𝑖,𝑗) ∈ 𝑁 𝛩 (𝑖) = 𝛩𝐷(𝑗) = 𝑘 𝑂

𝑂

𝐷

{ {

𝑡𝑂𝐷 𝑖𝑗 = max min

{ {

𝑡𝐷𝑂 𝑖𝑗 = max min 13

𝑑𝛩𝑂(𝑖) 𝑣𝐿𝑟

𝑑𝛩𝐷(𝑖)𝛩𝑂(𝑗) 𝑈𝑟

𝑣

} }, ∀(𝑖, 𝑗) ∈ 𝑁 , }, }, ∀(𝑖, 𝑗) ∈ 𝐴 .

𝑂𝑠 𝑂 ,𝑢𝐷𝑠 𝑗 ― 𝑢𝑖 ― 𝑠𝛩𝑂(𝑖) , 𝐷𝑠 𝐷 ,𝑢𝑂𝑠 𝑗 ― 𝑢𝑖 ― 𝑠𝛩𝐷(𝑖)

𝑑𝛩𝑂(𝑖)

𝑠

𝑣𝑈𝑟

𝑑𝛩𝐷(𝑖)𝛩𝑂(𝑗) 𝑣𝑈𝑟

𝑠

(39) (40)

Terms (27)-(29) are corresponding to terms (9)-(11). Constraint (30) is corresponding to constraint (12). Constraints (31) and (32) are corresponding to constraint (13). Constraints (33) - (38) are corresponding to constraints (14)-(20), respectively. Formulas (39) and (40) are corresponding to formulas (25) and (26) which can eliminate constraints (21) and (22). Because the objective function, the solution space, and the constraints of Model-S are equivalent to Model-T, Model-S is equivalent to Model-T. After the partition of time windows, Model-S is still a mixed integer nonlinear programming model. Further linearization method is required.

4.2 Upper bound estimation for Model-S If the decision variables of visiting time at sub nodes are replaced with discrete parameters, the upper bound for Model-S can be determined because the decision space for the upper bound is a sub set of that for Model-S (Fagerholt et al., 2010; Wang and Regan, 2002; Zhang et al., 2010). 𝑠 𝑠 𝑂𝑠 𝐷𝑠 𝐷𝑠 For constraints (35) and (36), let 𝑢𝑂𝑠 𝑖 = 𝑏𝑖 , ∀𝑖 ∈ 𝑂 ; let 𝑢𝑗 = 𝑏𝑗 , ∀𝑗 ∈ 𝐷 . According to formula (39), we can determine 𝑡𝑂𝐷 𝑖𝑗 as

{ {

𝑡𝑂𝐷 𝑖𝑗 = max min

} }

𝑑𝑝

𝑑𝑝 𝑂𝐷 ,𝑇 , , ∀ (𝑖, 𝑗) ∈ 𝑁𝑠 𝑖𝑗 𝑣𝑈𝑟 𝑣𝐿𝑟

(41)

𝑠 𝑂 𝐷𝑠 𝑂𝑠 𝑂 where 𝑇𝑂𝐷 𝑖𝑗 = 𝑏𝑗 ― 𝑏𝑖 ― 𝑠𝑝 , ∀ (𝑖, 𝑗) ∈ 𝑁 , 𝛩 (𝑖) = 𝑝. According to formula (40), we can determine 𝑡𝐷𝑂 𝑗𝑖 as

{ {

𝑡𝐷𝑂 𝑗𝑖 = max min

} }

𝑑𝑞𝑝

𝑑𝑞𝑝 𝐷𝑂 ,𝑇 , ,∀(𝑗, 𝑖) ∈ 𝐴𝑠 𝑗𝑖 𝑈𝑟 𝑈𝑟 𝑣 𝑣

(42)

𝑠 𝑂 𝐷 𝑂𝑠 𝐷𝑠 𝐷 where 𝑇𝐷𝑂 𝑗𝑖 = 𝑏𝑖 ― 𝑏𝑗 ― 𝑠𝑞 , ∀(𝑗,𝑖) ∈ 𝐴 , 𝛩 (𝑖) = 𝑝, 𝛩 (𝑗) = 𝑞. Substitute formula (41) into formula (27), and then we can calculate the minimum fuel consumption cost on (𝑖,𝑗) ∈ 𝑁𝑠 as

{(

(

)

𝑑𝑝 2 𝛿 ∆1 𝐿𝑟 + ∆2(𝑤 + 𝑓𝑝)𝑑𝑝 + ∆3𝑑𝑝(𝑣𝐿𝑟) ， 𝑣

𝑄𝑖𝑗 =

𝛿

∆1𝑇𝑂𝐷 𝑖𝑗

+ ∆2(𝑤 + 𝑓𝑝)𝑑𝑝 + ∆3

if 𝑇𝑂𝐷 𝑖𝑗 ≥

𝑑𝑝

. 𝑣𝐿𝑟

)

(𝑑𝑝)3

2 (𝑇𝑂𝐷 𝑖𝑗 + 𝜖)

𝑑𝑝 𝑑𝑝 ,if 𝑈𝑟 ≤ 𝑇𝑂𝐷 𝑖𝑗 < 𝐿𝑟. 𝑣 𝑣

(43)

Note that the truck is loaded when it travels from the origin to its destination of an order. Substitute formula (42) into formula (26), and then we can calculate the minimum fuel consumption cost on (𝑗,𝑖) ∈ 𝐴𝑠 as

(

{(

)

𝑑𝑞𝑝 2 𝛿 ∆1 𝐿𝑟 + ∆2𝑤𝑑𝑞𝑝 + ∆3𝑑𝑞𝑝(𝑣𝐿𝑟) ， 𝑣

𝑄𝑗𝑖 =

𝛿

14

∆1𝑇𝐷𝑂 𝑗𝑖

+ ∆2𝑤𝑑𝑞𝑝 + ∆3

(𝑑𝑞𝑝)3

)

2 (𝑇𝐷𝑂 𝑗𝑖 + 𝜖)

if 𝑇𝐷𝑂 𝑗𝑖 ≥

𝑑𝑞𝑝 . 𝑣𝐿𝑟

𝑑𝑞𝑝 𝑑𝑞𝑝 < . ,if 𝑈𝑟 ≤ 𝑇𝐷𝑂 𝑗𝑖 𝑣𝐿𝑟 𝑣

(44)

Note that the truck is unloaded when it travels from the destination of an order to the origin of another order. Based on the description above, the upper bound for Model-S can be determined by the following model, namely, Model-UB:

∑

Minimize

𝑄𝑖𝑗𝜔𝑖𝑗 +

∑

𝑄𝑗𝑖𝑥𝑗𝑖 +

（𝑗,𝑖） ∈ 𝐴𝑠

(𝑖,𝑗) ∈ 𝑁𝑠

∑ 𝜃𝑓 𝑑 𝑦

𝑖 𝑖 𝑖

𝑖∈𝑁

(45)

Subject to (30) - (34)

(

𝑑𝛩𝑂(𝑖)

(

𝑣𝑈𝑟

𝑑𝛩𝐷(𝑗)𝛩𝑂(𝑖) 𝑣𝑈𝑟

)

𝑠 ― 𝑇𝑂𝐷 𝑖𝑗 𝑤𝑖𝑗 ≤ 0, ∀(𝑖, 𝑗) ∈ 𝑁 ,

)

𝑠 ― 𝑇𝐷𝑂 𝑗𝑖 𝑥𝑗𝑖 ≤ 0,∀（𝑗,𝑖） ∈ 𝐴 .

(46)

(47)

The three terms of objective function (45) are corresponding to terms (27) - (29). Constraints (46) and (47) can ensure the feasibility of formulas (41) and (42). The reason is 𝑈𝑟 𝐷𝑂 that if 𝑑𝛩𝑂(𝑖)/𝑣𝑈𝑟 > 𝑇𝑂𝐷 𝑖𝑗 , then 𝑤𝑖𝑗 = 0; if 𝑑𝛩𝐷(𝑗)𝛩𝑂(𝑖)/𝑣 > 𝑇𝑗𝑖 , then 𝑥𝑗𝑖 = 0. Hence,

constraints (46) and (47) can be removed by redefining 𝑁𝑠 and 𝐴𝑠 as sparse matrices: 𝑂 𝑈𝑟 𝐷𝑠 𝑁𝑈𝐵 = {(𝑖,𝑗)│𝑖 ∈ 𝑂𝑠, 𝑗 ∈ 𝐷𝑠,𝛩𝑂(𝑖) = 𝛩𝐷(𝑗) and 𝑏𝑂𝑠 𝑖 + 𝑠𝛩𝑂(𝑖) + 𝑑𝛩𝑂(𝑖)/𝑣 ≤ 𝑏𝑗 }, 𝑈𝑟 𝑂𝑠 𝐷 𝐴𝑈𝐵 = {(𝑗,𝑖)|𝑗 ∈ 𝐷𝑠,𝑖 ∈ 𝑂𝑠;𝛩𝐷(𝑗) ≠ 𝛩𝑂(𝑖), and 𝑏𝐷𝑠 𝑗 + 𝑠𝛩𝐷(𝑗) + 𝑑𝛩𝐷(𝑗)𝛩𝑂(𝑖)/𝑣 ≤ 𝑏𝑖 }. 𝑈𝐵 𝑂𝑠 Arc (𝑖, 𝑗) ∈ 𝑁 means that a truck starts at origin 𝑂𝛩𝑂(𝑖) on time 𝑏𝑖 , and arrives at 𝑈𝐵 destination 𝐷𝛩𝐷(𝑗) on time 𝑏𝐷𝑠 means a truck starts at destination 𝐷𝛩𝐷(𝑗) 𝑗 . Arc (𝑗, 𝑖) ∈ 𝐴 𝐷𝑠 on time 𝑏𝑗 , and arrives at origin 𝑂𝛩𝑂(𝑖) on time 𝑏𝑂𝑠 𝑖 . After the transformation, Model-UB is a binary linear programming model, and the sparse matrix strategy can be used to improve the efficiency of the TWPB method (See Subsection 5.2.2). Remark 4.2: For any (𝑖,𝑗) ∈ 𝑁𝑈𝐵 or (𝑗,𝑖) ∈ 𝐴𝑈𝐵, it has two dimensions: time and space. If (𝑖,𝑗) ∈ 𝑁𝑈𝐵 or (𝑗,𝑖) ∈ 𝐴𝑈𝐵 is selected in a solution, the visiting time increases from 𝑖 to 𝑗 𝑈𝑟 𝑈𝑟 𝑂 𝐷𝑠 𝑂𝑠 𝐷 due to 𝑏𝑂𝑠 and 𝑏𝐷𝑠 𝑖 + 𝑠𝛩𝑂(𝑖) + 𝑑𝛩𝑂(𝑖)/𝑣 ≤ 𝑏𝑗 𝑗 + 𝑠𝛩𝐷(𝑗) + 𝑑𝛩𝐷(𝑗)𝛩𝑂(𝑖)/𝑣 ≤ 𝑏𝑖 . Thus, sub tour will not be generated. Constraints (31) and (32) ensure that: the visiting time of a route is continuous, which can be observed visually in Figure 2 (see Subsection 4.4). That’s why Model-UB is feasible, and it is the upper bound of Model-S. Accordingly, this optimal solution can be selected as an approximate optimal solution of Model-S.

4.3 Lower bound estimation for Model-S As shown in Model-S, the cost of the private fleet contains the fuel consumption cost on both selected orders and arcs. Therefore, if the lower bound of the fuel consumption cost for each order or each arc can be obtained, one can obtain the lower bound of the fuel consumption cost. The fuel consumption cost for each order is dependent on the truck travel time, so we can relax the time window constraints to obtain the lower bound. As shown in Model-S, 15

following constraints (35) and (37), we have 𝑠 𝐷 𝑂 𝑂 𝐷𝑠 𝑂𝑠 𝑂 𝑡𝑂𝐷 𝑖𝑗 ≤ 𝑢𝑗 ― 𝑢𝑖 ― 𝑠𝑝 ≤ 𝑏𝑗 ― 𝑎𝑖 ― 𝑠𝑝 , ∀(𝑖,𝑗) ∈ 𝑁 , where 𝑝 = 𝛩𝑂(𝑖). Then, 𝑡𝑂𝐷 𝑖𝑗 can be relaxed by 𝐷𝑠 𝑂𝑠 𝑂 𝑡𝑂𝐷 𝑖𝑗 ≤ 𝑏𝑗 ― 𝑎𝑖 ― 𝑠𝑝 .

(48)

According to Proposition 3.5, we have, 𝑠 𝐿𝑟 𝐷𝑠 𝑂𝑠 𝑂 𝑈𝑟 𝑡𝑂𝐷 𝑖𝑗 = max{min {𝑑𝑝/𝑣 ,𝑏𝑗 ― 𝑎𝑖 ― 𝑠𝑝 },𝑑𝑝/𝑣 }, ∀(𝑖, 𝑗) ∈ 𝑁 .

(49)

Similarly, 𝐿𝑟 𝑂𝑠 𝐷𝑠 𝐷 𝑈𝑟 𝑠 𝑡𝐷𝑂 𝑗𝑖 = max{min {𝑑𝑞𝑝/𝑣 ,𝑏𝑖 ― 𝑎𝑗 ― 𝑠𝑞 },𝑑𝑞𝑝/𝑣 },(𝑗,𝑖) ∈ 𝐴 . 𝐷

𝑡𝑂𝐷 𝑖𝑗

𝑡𝐷𝑂 𝑗𝑖

𝑄𝐿𝐵 𝑖𝑗 𝜔𝑖𝑗 +

∑

(50)

where 𝑞 = 𝛩 (𝑗). Substitute and into the objective function, and then we can obtain the lower bound of the fuel consumption cost on (𝑖,𝑗) ∈ 𝑁𝑠 and on (𝑗,𝑖) ∈ 𝐴𝑠, 𝐿𝐵 namely, 𝑄𝐿𝐵 𝑖𝑗 and 𝑄𝑗𝑖 , respectively. Then, similar as Model-UB, we can define two sparse matrices for 𝑁𝑠 and 𝐴𝑠 as 𝑂 𝑈𝑟 𝐷𝑠 𝑁𝐿𝐵 = {(𝑖,𝑗)│𝑖 ∈ 𝑂𝑠, 𝑗 ∈ 𝐷𝑠,𝛩𝑂(𝑖) = 𝛩𝐷(𝑗) and 𝑎𝑂𝑠 𝑖 + 𝑠𝛩𝑂(𝑖) + 𝑑𝛩𝑂(𝑖)/𝑣 ≤ 𝑏𝑗 }, 𝑈𝑟 𝑂𝑠 𝐷 𝐴𝐿𝐵 = {(𝑗,𝑖)|𝑗 ∈ 𝐷𝑠,𝑖 ∈ 𝑂𝑠;𝛩𝐷(𝑗) ≠ 𝛩𝑂(𝑖), and 𝑎𝐷𝑠 𝑗 + 𝑠𝛩𝐷(𝑗) + 𝑑𝛩𝐷(𝑗)𝛩𝑂(𝑖)/𝑣 ≤ 𝑏𝑖 }. After the relaxation and refinement, we can estimate the lower bound for Model-S by the following model, Model-LB: Minimize

∑ (𝑖,𝑗) ∈ 𝑁

𝐿𝐵

（𝑖,𝑗） ∈ 𝐴

𝑄𝐿𝐵 𝑖𝑗 𝑥𝑖𝑗 + 𝐿𝐵

∑ 𝜃𝑓 𝑑 𝑦

𝑖 𝑖 𝑖

𝑖∈𝑁

(51)

Subject to (30)-(34) As can be seen, Model-LB is a binary linear programming model.

4.4 An illustrative example An illustrative example as shown in Figure 2 is presented for the UTWPB method. Figure 2 shows a solution in terms of visiting time, distance, and truck speed for the path 𝑂1 ― 𝐷1 ― 𝑂2 in Figure 1. We can see that Order 3 is outsourced due to the violation of time window. The black line represents the path of the optimal solution; the green line and the red line are the upper bound and lower bound of the optimal solution, respectively. The slope of a line indicates the truck speed. 𝑂𝑟 𝐷𝑟 𝑂𝑟 𝑂𝑟 [𝑎𝐷𝑟 In this example, each of time windows [𝑎𝑂𝑟 1 , 𝑏1 ] , 1 , 𝑏1 ], and [𝑎2 , 𝑏2 ] is partitioned into four sub time windows (the width of one sub time window is zero), 𝐷𝑠 respectively. A truck travels from 𝐷1 during sub time window [𝑎𝐷𝑠 3 , 𝑏3 ] to 𝑂2 during sub 𝐷𝑠 time window [𝑎𝐷𝑠 7 , 𝑏7 ]. In this road segment, the truck is unloaded. In the optimal solution, 𝑂𝑠 it should start at 𝐷1 on time 𝑢𝐷𝑠 3 , and arrives at 𝑂2 on time 𝑢7 ; in the upper bound 𝑂𝑠 solution, it has to start at 𝐷1 on time 𝑏𝐷𝑠 3 , and arrives at 𝑂2 on time 𝑏7 ; and in the lower 𝐷𝑠 bound solution, it has to starts at 𝐷1 on time 𝑎3 , and arrives at 𝑂2 on time 𝑏𝑂𝑠 7 . Without loss generality, since the service time at 𝐷1 is the same, we set it to be 0. Then, we can see 𝐷𝑠 𝑂𝑠 𝐷𝑠 𝑂𝑠 𝐷𝑠 that 𝑏𝑂𝑠 7 ― 𝑏3 ≤ 𝑢7 ― 𝑢3 ≤ 𝑏7 ― 𝑎3 , and the fuel consumption cost can be calculated as 𝑈𝐵 𝐷𝑠 𝑄𝐿𝐵 37 ≤ 𝑄37 ≤ 𝑄37 . For the upper bound solution, the truck arrives at 𝐷1 on time 𝑏3 from 𝑂1, and starts from 𝐷1 on 𝑏𝐷𝑠 3 to 𝑂2, the visiting time is continuous and increased. That is why the upper bound solution is feasible and no sub tour exists. For the lower bound solution, 16

𝐷𝑠 𝐷𝑠 𝐷𝑠 the truck arrives at 𝐷1 on time 𝑏𝐷𝑠 3 from 𝑂1, but starts from 𝐷1 on 𝑎3 (𝑎3 ≤ 𝑏3 ). As can be seen, the lower bound solution may be an infeasible solution. Distance

Optimal solution for load truck

O3 a7Os

b7Ds

for unload truck Upper bound solution for load truck

O2

for unload truck

a3Ds

b3Ds Lower bound solution for load truck

D1

for unload truck Time window and sub time widow

O1

u3Ds

u7Os Time

a1Or

b1Or

Figure 2. An illustrative example for TWPB method

4.5 Framework of the UTWPB method The framework of the UTWPB method is as shown in Algorithm 4.1. Lines 1-3 generate an initial solution, and then lines 4 – 8 improve the initial solution. To be specific, the method proposed by Zhang et al. (2010) is adapted to generate the initial solution, and the constant width partition method is used to reduce the number of iterations. Then, the selected time windows by Model-LB are repartitioned by the bisection method to update the set of partitioned time windows, based on which, Model-LB and Model-UB are resolved. This process repeats until one of the following three conditions is satisfied: (1) LB/UB is greater than a given value, where LB and UB are the optimal solutions obtained by solving Model-LB and Model-UB, respectively; (2) the iterations reach a given number, MaxIter; (3) the computational time of the iterations is greater than a given time, MaxTim. Algorithm 4.1 The framework of the UTWPB method 1 Partition the time windows with the constant width method (see Subsection 4.1); 2 Solve Model-UB; 3 Solve Model-LB; 4 Repeat 5 Partition the selected time windows with the bisection method (see Subsection 4.1); 6 Solve Model-UB; 7 Solve Model-LB; 8 Until a stop condition is satisfied In the process of Algorithm 4.1, the solution with the largest objective value obtained by Model-LB, and the solution with the minimum objective value obtained by Model-UB are 17

taken as the lower bound and upper bound of the optimal solution of Model-T, respectively. For Algorithm 4.1, we have Proposition 4.1, Proposition 4.1 In Algorithm 4.1, according to the iterations progress, the upper bound of the optimal solution obtained by Model-UB is monotonically non-increasing, and the lower bound of the optimal solution obtained by Model-LB is monotonically non-decreasing. Proof of Proposition 4.1. See the appendix. We can see that the time complexity of Algorithm 4.1 is dependent on Model-UB and Model-LB, and the time complexity of Model-UB and Model-LB are dependent on the size of sub nodes. If an instance of the FTPDP-RF consists of 𝑛 orders, then, ∑𝑛

(⌈(𝑏𝑂𝑟 ― 𝑎𝑂𝑟) Λ⌉ + 1)

𝑖=1

𝑛 sub origins and ∑𝑖 = 1(⌈(𝑏𝐷𝑟 ― 𝑎𝐷𝑟) Λ⌉ + 1) sub destinations

can be generated by the constant width method according to Subsection 4.1. Moreover, for each iteration, no more than 𝑛 sub origins and 𝑛 sub destinations can be generated by the bisection method. To reduce the computational time, we present Algorithm 4.2, where Model-UB is not necessary to be solved at each iteration. According to Proposition 4.1, the optimal solution obtained by Model-UB is monotonically non-increasing. Then, we need only to record the results of Model-UB at the final iteration. Accordingly, based on the property of Model-LB in iteration process, the stop condition (1) can be modified by that LBIter - 1/LBIter is greater than a given value ε, where LBIter - 1 and LBIter are the optimal objective values obtained by solving Model-LB at the previous iteration and the current iteration, respectively. Algorithm 4.2 The framework of the improved UTWPB method 1 Partition the time windows with the constant width method; 2 Solve Model-LB; 3 Repeat` 4 Partition the selected time windows with the bisection method. 5 Solve Model-LB 6 Until a stop condition is satisfied 7 Solve Model-UB

5 Numerical experiments The numerical experiments are designed to analyze the performance of the UTWPB method, and to compare it with an existing TWPB method and the SDB method (Subsection 5.2); then, a sensitivity analysis of the FTPDP-RF is conducted, and management insights are provided (Subsection 5.3). Most of the existing literature on traditional full truckload pickup and delivery problem adopt randomly generated instances to evaluate their models and algorithms (Gendreau et al., 2015; Jula et al., 2005; Lai et al., 2017; Li et al., 2015; Liu et al., 2010; Zhang et al., 2014). In this research, we also employ the similar method for our experiments (Subsection 5.1).

5.1 Generation of instances The parameters for calculating fuel consumption come from existing literature and the reality. Specifically, the parameters, ∆𝑖(𝑖 = 1, 2, 3), 𝑤 = 6.35 ton, 𝑣𝐿 = 20 kilometer/hour, 18

are inferred from the expressions and typical values in Demir et al. (2012). According to traffic laws about trucks in China, 𝑣𝑈= 90 kilometer/hour, and 𝜃 = 1.5 RMB/(kilometer * ton), and 𝛿 = 6.5 RMB/Liter are designed as local prices for outsourcing and diesel, respectively. The method for generating truckload orders is similar as Zhang et al. (2010) and Zhang et al. (2014). The parameters are shown in Table 1. The location of the depot is set to be (90, 90) since the depot usually locates in the central region. The locations of the origin and the destination for each order are uniformly distributed between 0 and 180. All the orders must be served during a working day (within eight hours, from 60 to 540 minutes). The working time for each truck is limited within ten hours (during 600 minutes) according to relevant labor laws. Since the start time of loading must be ahead of the start time of unloading, the upper bound of the time window for an origin is designed to be 420 minutes. For order 𝑖 ∈ 𝑁, 𝑎𝑂𝑖 is a random value falling within [60, OTWaMax] and 𝑏𝑂𝑖 = 𝑎𝑂𝑖 + TWWidth. Then, 𝑎𝐷𝑖 = 𝑎𝑂𝑖 + 𝑠𝑂𝑖 + 𝑑𝑖/𝑣𝑈 and 𝑏𝐷𝑖 = 𝑎𝐷𝑖 + TWWidth. If 𝑏𝐷𝑖 > 𝐷𝑇𝑊𝑏𝑀𝑎𝑥, the order cannot be fulfilled on time. Hence, the time windows for order 𝑖 ∈ 𝑁 need to be regenerated. The algorithm is coded in C++. Model-UB and Model-LB are solved by calling a MIP solver of CPLEX 12.5, namely, branch and cut solver. All the experiments are implemented on a 64bit Windows 10 personal computer with the processor, Inter(R) Core(TM) i7-2600, @3.4 GHz and 8GB memory. Table 1. Parameters of instances Name Description Value Xmax maximum length of the Euclidean plain 180 kilometer Ymax maximum width of the Euclidean plain 180 kilometer rateR/O rate of reserved orders [0.1, 0.3] rateNT/NDT satisfaction rate of the demand for trucks [0.5, 1.0] 𝐻

maximum of working time of trucks

600 minutes

OTWaMax DTWbMax TWWidth 𝑠𝑂𝑖 𝑠𝐷𝑖 𝑓𝑖

upper bound of the time window for an origin upper bound of the time window for a destination width of time window service time at the origin for order 𝑖 ∈ 𝑁 service time at the destination for order 𝑖 ∈ 𝑁 weight of freight for order 𝑖 ∈ 𝑁

[60, 420] minutes 540 minutes [30, 60] minutes [10, 30] minutes 𝑠𝐷𝑖 = 𝑠𝑂𝑖 [4, 12] tons

5.2 Analysis and performance of the solution methods To analyze and show the performance of the UTWPB method, we present the following experiment. We generate four sets of instances with different sizes at first. Then, we provide an illustrative example to demonstrate the validity of Proposition 4.1 and 4.2 (Subsection 5.2.1). After that, we show the effectiveness of using spare matrix to describe Model-UB (Subsection 5.2.2). Finally, we compare the performance of the UTWPB method with that of other two existing methods (Subsection 5.2.3).

19

5.2.1 Analysis of the UTWPB method To show the relation between the width of sub time window 𝛬 and the accuracy of solutions, we solve a 25-order instance using Algorithm 4.1, the setting of which is as follows: ∇ = 1 minute, LB/UB ≥ 0.9999, MaxIter = 10, MaxTim = 3600s. The results are shown in Figure 3, where 𝐿𝐵𝑖 and 𝑈𝐵𝑖, 𝑖 = 2, 3,4,5, represent the optimal solution at corresponding iteration obtained by solving Model-LB and Model-UB, in which 𝛬 = 𝑖. Figure 3 confirms the validity of Proposition 4.1. Meanwhile, we noticed that the smaller 𝛬 is, the less of the gap between LB and UB is. The values of UB decrease slowly and that of LB increase slowly after 5 iterations. Analogously, we can acquire the relation between the width of sub time window 𝛬 and the computational time. As shown in Figure 4, 𝐼𝑇𝑖 and 𝑇𝑇𝑖, 𝑖 = 2, 3,4,5, represent the computational time for the corresponding individual iteration and total iterations, respectively. As can be seen, with the iterations increase, the computational time of an individual iteration increases gently, but the total computation time increases rapidly. Specifically, the smaller 𝛬 is, the more rapidly the time increases. The reasons for these results can be interpreted as follows: the number of sub nodes increases gently at each iteration; the total computational time is the accumulation of the computational time of individual iteration; the smaller 𝛬 is, the more sub nodes are generated. 6030

Cost(RMB)

6010 LB5

5990

LB4

5970

LB3

5950

LB2 UB5

5930

UB4

5910

UB3

5890 0

1

2

3

4 5 Iteration

6

7

8

9

10

Figure 3. Relation among width of initial sub time windows, iterations, and solutions

20

UB2

Computional time (s)

600 500

IT5

400

IT4 IT3

300

IT2

200

TT5 TT4

100

TT3

0

TT2 0

1

2

3

4 5 Iteration

6

7

8

9

10

Figure 4. Relation among width of initial sub time windows, iterations, and computational time In addition, a comparison between Algorithm 4.1 and Algorithm 4.2 is conducted and the results are shown in Table 2. Algorithm 4.2 requires less computational time compared with Algorithm 4.1 without increasing the UB nor distinctly decreasing the LB. There are two explanations for this result: Model-UB in Algorithm 4.2 is only solved once; the iteration process for Algorithm 4.2 may stop earlier considering its stop condition (1). Table 2. Comparison between Algorithm 4.1 and Algorithm 4.2 　 Algorithm 4.1 Algorithm 4.2 Imp% Λ

UB

2 3 4 5

5996.1 5996.2 5996.2 6021.1

LB Iter 5924.5 5925.8 5923.5 5923.0

10 10 10 10

CT

UB

569 270 183 157

5996.1 5996.2 5996.2 6021.1

LB Iter 5923.3 5925.8 5919.6 5923.0

CT

UB

LB

CT

6 182 10 150 6 55 10 73

0% 0% 0% 0%

-0.02% 0.00% -0.07% 0.00%

-68% -44% -70% -53%

Iter represents the number of iterations. CT represents the computational time. Imp% represents the improvement of calculation results by using Algorithm 4.2 over Algorithm 4.1.

5.2.2 Effectiveness of the spare matrix A comparison is conducted between the TWPB method described by a spare matrix (TWPBM-S) and the TWPB method (TWPBM-N) employed by Zhang et al. (2010) and Zhang et al. (2014). In this case, twenty instances are generated, and each of them contains ten orders. Setting Λ=1, these instances are solved by Model-UB. To easily record the corresponding data, we adopt Lingo 11.0 as the solver. The average results are shown in Table 3: the TWPBM-S, in average, can reduce the number of decision variables by 82.36%, the number of constraints by 99.56%, the number of nonzero coefficients by 77.14%, the

21

memory usage by 68.18%, and the computational time by 36%. The results confirm the effectiveness of TWPBM-S. Table 3. TWPBM-N VS. TWPBM-S TWPBM-N

TWPBM-S

GAP

NumVar

232330.95

40979.8

-82.36%

NumCom

211657.25

932.25

-99.56%

NumCoe

944677.15

215977.3

-77.14%

64078

20392.55

-68.18%

22.9

14.65

-36.03%

MomUse (KB) CT (s)

TWPBM-N represents the existing TWPB method not using a spare matrix. TWPBM-S represents the TWPB method using a spare matrix NumVar is the average number of decision variables of Model-UB. NumCom is the average number of constraints of Model-UB. NumCoe is the average number of nonzero coefficient of Model-UB. MomUse is the average usage of memory in Model-UB. CT is average computational time when solving Model-UB. GAP = (TWPBM-S – TWPBM-N) / TWPBM-N *100%.

5.2.3 Comparison among TWPB methods and SDB method In this subsubsection, we compare the performance among TWPB methods (TWPBM-S and UTWPB method) and SDB method (Bektaş and Laporte, 2011) since all of them are approximate methods. Four sets of instances are generated: 25-order instances, 50-order instances, 100-order instance, and 200-order instance. Each set contains twenty instances. These four sets of instances are solved as follows: (1) The 25-order instances are solved by the SDB method with setting the width of sub speeds to be 0.01 kilometer/minute, and by the TWPBM-S with setting the width of sub time windows, Λ =1 minute, respectively. Because the size of orders and Λ are very small, a good initial solution can be obtained by the TWPB method in a short time. Therefore, we don’t need to employ the UTWPB method to process the iterations. The results are shown in Table 4. (2) The 50-order instances, 100-order instances, and 200-order instances are solved by the SDB method, the TWPBM-S, and the UTWPB method, respectively. Accordingly, the widths of sub speeds are set to be 0.01 kilometer/minute, 0.01 kilometer/minute, and 0.05 kilometer/minute, respectively; the widths of sub time windows for the constant width method are set to be 2 minutes, 4 minutes, and 10 minutes, respectively; the minimum width for the bisection method is set to be 1 minutes; moreover, ε ≥ 0.999, MaxIter = 5, MaxTim = 3600s, 3600s, and 7200s, respectively. The results are shown in Table 5 – 7. For the 25-order instances and 50-order instances, both of the SDB method and the TWPB methods can obtain near-optimal solutions. Meanwhile, the average performance of the SDB method is better than the TWPB methods in terms of solutions and computational time, but some of the solutions obtained by the TWPB methods are better than that obtained

22

by the SDB method. Moreover, the average gaps between the computational results of them are very small (less than 0.12%). For the 100-order instances and 200-order instances, the average performance of the TWPB methods is better than the SDB method. The UTWPB method outperforms the TWPBM-S since the iteration process can decrease the upper bound and increase the lower bound. Moreover, with the increasing of widths of sub time windows of the UTWPB method, the gap between the upper bound and the lower bound increases, but it is acceptable even when the size of orders adds up to 200 orders. These results indicate that both the SDB and the UTWPB method are valid, but their performance varies in different circumstances. Specifically, for instances the size of which is less than 50 orders, the SDB method is the better choice; whereas, for instances the size of which is between 100 and 200 orders, the proposed UTWPB method is the better choice. Table 4. Result of 25-order instances Ins.

Cost (RMB) DV

UB0

GAP

LB0

Computational time(s)

DV\UB0

UB0\LB0

DV

UB0

LB0

1

5469.6

5468.6

5460.9

0.02%

0.14%

11.2

215.6

211.7

2

9464.4

9464.0

9315.1

0.00%

1.57%

14.6

173.9

112.6

3

6126.7

6126.5

6108.0

0.00%

0.30%

15.3

161.0

151.8

4

5381.5

5380.9

5364.6

0.01%

0.30%

10.8

112.5

108.9

5

5244.8

5244.3

5238.3

0.01%

0.11%

10.4

161.5

128.2

6

5970.6

5969.6

5959.0

0.02%

0.18%

10.5

104.5

128.6

7

5504.5

5504.0

5491.3

0.01%

0.23%

15.7

148.1

136.5

8

10440.7

10439.7

10420.4

0.01%

0.18%

13.5

146.6

129.9

9

7990.2

7989.2

7974.8

0.01%

0.18%

17.6

132.8

117.6

10

7052.9

7052.0

7041.8

0.01%

0.15%

14.5

126.5

100.2

11

10970.2

10969.9

10951.1

0.00%

0.17%

15.6

135.3

131.5

12

6459.5

6459.0

6448.5

0.01%

0.16%

16.6

79.0

71.7

13

8329.0

8328.2

8317.7

0.01%

0.13%

18.0

96.1

113.1

14

10163.4

10162.6

10144.5

0.01%

0.18%

21.5

100.7

111.6

15

5458.1

5457.7

5442.8

0.01%

0.27%

21.1

184.0

170.0

16

5825.3

5824.4

5810.8

0.02%

0.23%

19.1

112.4

98.6

17

7613.6

7612.8

7591.7

0.01%

0.28%

28.5

107.4

81.2

18

9221.6

9221.0

9208.3

0.01%

0.14%

24.4

118.8

111.9

19

5954.1

5953.2

5941.7

0.02%

0.19%

24.2

124.7

116.0

20

5939.1

5996.1

5927.4

-0.95%

1.14%

33.1

131.2

86.2

average

7229.0

7231.2

7207.9

-0.03%

0.32%

17.8

133.6

120.9

Ins. represents the serial number of instance, and DV represents the results obtained by the SDB method. UB0, LB0 represent the results by solving Model-UB and Model-LB of the TWPBM-S, respectively. Note that DV\UB0 = (DV ― UB0 )/ UB0 * 100%, and UB0\LB0 = (UB0 ― UB0 )/ UB0 * 100%.

23

Table 5. Result of 50-order instances Ins.

Cost (RMB)

GAP

DV

UB0

UB1

LB0

1

12955.7

12975.5

12974.9

12552.5

2

12852.6

12854.7

12854.3

3

9906.2

10050.3

4

16437.9

5

LB1

Computational time(s)

DV\UB1

UB0\UB1

UB1\LB1

LB0\LB1

DV

UB0

LB0

12872.7

-0.15%

0.00%

0.79%

12751.5

12761.5

-0.01%

0.00%

10049.6

9838.6

9846.3

-1.43%

16454.5

16452.7

16301.5

16369.6

11211.4

11210.7

11210.2

11157.8

6

11371.8

11373.1

11371.9

7

12208.2

12216.2

8

14973.4

UB1-LB1

-2.49%

232.64

168.74

153.4

631.352

0.72%

-0.08%

150.48

161.21

190.5

631.352

0.01%

2.02%

-0.08%

149.67

142.71

119.1

482.042

-0.09%

0.01%

0.51%

-0.42%

74.422

136.22

148.4

658.078

11163.8

0.01%

0.00%

0.41%

-0.05%

59.479

93.631

98.8

352.948

11293.5

11295.8

0.00%

0.01%

0.67%

-0.02%

54.892

77.488

77.41

251.388

12215.6

12126.8

12134.1

-0.06%

0.00%

0.67%

-0.06%

78.647

100.72

104.7

337.695

14974.2

14972.6

14877.1

14885.5

0.01%

0.01%

0.58%

-0.06%

108.55

104.26

128.1

427.419

9

11740.6

11787.0

11786.9

11680.5

11682.6

-0.39%

0.00%

0.88%

-0.02%

71.496

156.37

107.8

399.131

10

11107.8

11107.0

11106.2

11061.9

11065.0

0.01%

0.01%

0.37%

-0.03%

113.78

126.09

235.4

515.823

11

24175.8

24176.3

24175.1

24107.9

24113.8

0.00%

0.00%

0.25%

-0.02%

100.79

103.3

98.8

326.259

12

21054.5

21055.3

21054.5

20659.1

20739.1

0.00%

0.00%

1.50%

-0.39%

89.013

109.66

119.6

624.178

13

14598.4

14598.2

14597.8

14544.7

14549.9

0.00%

0.00%

0.33%

-0.04%

96.235

79.005

77.68

293.037

14

14137.5

14136.6

14135.5

14059.8

14062.4

0.01%

0.01%

0.52%

-0.02%

81.415

86.49

87.95

295.613

15

21695.6

21695.2

21694.5

21608.2

21611.9

0.01%

0.00%

0.38%

-0.02%

159.51

133.36

89.52

342.103

16

15173.4

15187.6

15187.4

14977.9

14986.2

-0.09%

0.00%

1.32%

-0.06%

167.67

216.34

156.9

557.845

17

15751.5

15751.1

15750.9

15635.8

15641.1

0.00%

0.00%

0.70%

-0.03%

83.845

158.16

99.38

328.122

18

12722.0

12721.9

12721.2

12647.4

12653.9

0.01%

0.01%

0.53%

-0.05%

134.55

95.194

100.5

334.441

19

15406.5

15481.9

15480.6

14969.0

15004.5

-0.48%

0.01%

3.08%

-0.24%

149.45

146.39

116.2

587.946

20

12658.7

12659.6

12658.8

12536.0

12547.9

0.00%

0.01%

0.88%

-0.09%

160.06

114.31

82.34

277.202

average

14607.0

14623.3

14622.6

14469.4

14499.4

-0.11%

0.01%

0.84%

-0.21%

115.83

125.48

119.6

432.699

UB1, LB1 represent the results obtained by Algorithm 4.2. Note that DV\UB1 = (DV ― UB1 )/ UB1 * 100%, UB0\UB1 = (UB0 ― UB1 )/ UB1 * 100%, UB1\LB1 = (UB1 ― LB1 )/ LB1 * 100%, and LB0\LB1 = (LB0 ― LB1 )/ LB1 * 100%. UB1-LB1 indicates the computational time by using Algorithm 4.2. 24

Table 6. Result of 100-order instances Ins.

Cost (RMB)

GAP

DV

UB0

UB1

LB0

LB1

1

27810.6

25707.1

25338.4

24390.6

2

19218.6

19270.1

19207.6

3

28408.3

28015.4

4

31799.0

5

Computational time(s)

DV\UB1

UB0\UB1

UB1\LB1

LB0\LB1

DV

UB0

LB0

UB1-LB1

24703.1

9.76%

1.46%

2.51%

-1.27%

3600

256.347

178.25

1311.33

18898.3

19047.6

0.06%

0.33%

0.83%

-0.78%

3600

245.552

302.84

2706.06

27990.8

27173.4

27579.2

1.49%

0.09%

1.47%

-1.47%

3600

226.474

296.01

2399.88

30706.9

30353.8

29213.6

29845.3

4.76%

1.16%

1.68%

-2.12%

3600

96.429

178.32

1389.77

34496.0

32164.1

31838.7

30389.4

31433.7

8.35%

1.02%

1.27%

-3.32%

3600

228.56

315.55

2551.61

6

25167.4

23384.5

23251.1

22310.5

22719.1

8.24%

0.57%

2.29%

-1.80%

3600

153.919

241.34

3167.87

7

22526.3

22647.8

22599.8

22291.8

22397.7

-0.33%

0.21%

0.89%

-0.47%

3600

89.48

191.9

1546.94

8

22255.7

22238.0

22227.8

21787.8

21968.1

0.13%

0.05%

1.17%

-0.82%

3600

218.381

157.61

1020.07

9

23224.2

23304.5

23294.5

22724.2

22961.4

-0.30%

0.04%

1.43%

-1.03%

3600

259.233

149.61

1569.19

10

24991.2

24699.6

24569.3

23557.2

23819.7

1.72%

0.53%

3.05%

-1.10%

3600

118.079

136.55

1754.51

11

29937.5

29690.1

29678.7

27795.9

28980.8

0.87%

0.04%

2.35%

-4.09%

3600

214.005

220.04

1407.11

12

26022.4

26063.8

26048.0

25548.0

25782.2

-0.10%

0.06%

1.02%

-0.91%

3600

181.405

268.98

2075.82

13

21997.2

22008.5

21992.7

21540.2

21813.0

0.02%

0.07%

0.82%

-1.25%

3600

181.771

251.85

1987.79

14

21787.1

21787.5

21779.8

21602.9

21674.6

0.03%

0.04%

0.48%

-0.33%

3600

257.909

262.38

1147.08

15

21360.9

21441.4

21376.4

21091.6

21191.0

-0.07%

0.30%

0.87%

-0.47%

3600

110.905

128.46

1077.93

16

22930.0

22960.3

22904.2

22596.6

22755.5

0.11%

0.24%

0.65%

-0.70%

3600

126.084

247.9

2772.32

17

19829.0

19855.1

19836.9

19603.2

19704.2

-0.04%

0.09%

0.67%

-0.51%

3600

120.758

191.94

2074.56

18

40698.7

41053.3

40636.3

39845.9

40285.8

0.15%

1.03%

0.86%

-1.09%

3600

120.426

135.34

1506.58

19

29409.8

28554.2

28456.1

27438.2

27979.2

3.35%

0.34%

1.68%

-1.93%

3600

279.164

324.75

4315.88

20

21413.4

21393.5

21378.4

21126.9

21237.6

0.16%

0.07%

0.66%

-0.52%

3600

144.765

133.86

1315.77

average

25764.2

25347.3

25238.0

24546.3

24893.9

2.08%

0.43%

1.36%

-1.40%

3600

181.482

215.67

1954.904

Table 7. Result of 200-order instances 25

Ins.

Cost (RMB)

GAP

Computational time(s)

DV

UB0

UB1

LB0

LB1

DV\UB1

UB0\UB1

UB1\LB1

LB0\LB1

DV

UB0

LB0

UB1-LB1

1

88015.4

68253.9

67059.5

60644.8

64135.9

31.25%

1.78%

4.36%

-5.44%

7200

125.73

288.443

2753.07

2

65075.7

45711.8

44766.9

40745.8

43179.2

45.37%

2.11%

3.55%

-5.64%

7200

162.87

160.723

4296.89

3

87669.2

63519.7

62688.7

56353.4

61052.2

39.85%

1.33%

2.61%

-7.70%

7200

185.48

231.534

4323.34

4

46044.1

40493.0

40287.5

39050.6

39695.3

14.29%

0.51%

1.47%

-1.62%

7200

96.165

112.981

2508.99

5

76617.0

55105.6

54282.2

49548.2

52421.1

41.15%

1.52%

3.43%

-5.48%

7200

217.17

151.034

2880.96

6

67634.7

54179.9

53423.9

48873.4

51656.7

26.60%

1.42%

3.31%

-5.39%

7200

91.008

143.968

2958.26

7

65017.8

45074.6

44518.1

41167.4

43018.8

46.05%

1.25%

3.37%

-4.30%

7200

195.46

172.583

3916.19

8

42474.5

39720.2

39536.7

38439.5

38439.5

7.43%

0.46%

2.78%

0.00%

7200

101.64

136.164

1998.29

9

76176.0

66709.2

65547.1

60660.8

63708.4

16.22%

1.77%

2.81%

-4.78%

7200

218.36

225.935

3185.52

10

79221.2

53975.4

52342.2

47348.0

50241.8

51.35%

3.12%

4.01%

-5.76%

7200

138.5

209.369

4072.78

11

99790.4

82905.9

81199.6

73828.1

79190.0

22.90%

2.10%

2.47%

-6.77%

7200

143.12

176.738

3027.01

12

50967.7

40617.8

40404.3

39187.3

39811.0

26.14%

0.53%

1.47%

-1.57%

7200

121.54

230.704

2828.15

13

69777.8

46711.4

45549.0

41826.5

43940.7

53.19%

2.55%

3.53%

-4.81%

7200

200.36

250.657

5470.95

14

116593.0

97994.5

95047.7

85417.9

92435.8

22.67%

3.10%

2.75%

-7.59%

7200

97.171

146.131

1870.85

15

111859.0

89941.1

88748.7

80611.9

86148.5

26.04%

1.34%

2.93%

-6.43%

7200

112.19

355.107

2437.55

16

39460.6

39367.1

39240.2

38431.0

38903.8

0.56%

0.32%

0.86%

-1.22%

7200

109.42

153.351

1917.31

17

41657.0

41243.4

41075.9

40073.9

40641.0

1.41%

0.41%

1.06%

-1.40%

7200

157.36

156.222

2200.14

18

82249.6

63077.3

62037.7

56230.1

60005.1

32.58%

1.68%

3.28%

-6.29%

7200

198.69

267.291

3272.64

19

76886.7

44866.4

44149.7

41138.0

42781.8

74.15%

1.62%

3.10%

-3.84%

7200

179.69

172.407

3927.38

20

91479.0

71019.6

70228.9

63039.3

67575.7

30.26%

1.13%

3.78%

-6.71%

7200

149.81

233.028

3805.7

73733.3

57524.4

56606.7

52130.8

54949.1

30.26%

1.62%

2.93%

-5.13%

7200

150.09

198.719

3182.6

average

26

5.3 Analysis of the model 5.3.1 Effect of truck speed on operational cost This sub-subsection analyzes the effect of truck speed on the total operational cost of the FTPDP-RF. We generate twenty random instances, and each instance contains ten orders. The truck speeds are set to be 60, 70, 80, 90 and variable, respectively. Note that ‘variable’ represents that truck speeds on road segments are set to be decision variables. Meanwhile, two types of outsourcing charges are investigated: high price of outsourcing and low price of outsourcing. The computational results are shown in Figure 5. From Figure 5(a), when the outsourcing price is set to be high in the market, along with the increasing of truck speeds from 60 to 90, the outsourcing charges decrease, while the fuel cost increases. Moreover, the least operational cost can be reached only when the truck speeds are considered as decision variables. The reason for this result is that the private fleet with a high truck speed can serve more orders in the given market. When truck speeds are considered as decision variables, the outsourcing cost can be reduced by increasing the truck speeds on some road segments, and the fuel cost can be reduced by decreasing the truck speeds on some other road segments. From Figure 5(b), we can see that: when the outsourcing price is set to be low in the market, along with the increasing of truck speeds from 60 to 90, the outsourcing charges increase, and the fuel cost grows insignificantly. The reason for this result is that when the speed of private trucks increases from 60-90, more orders will be outsourced to decrease the fuel consumption in the given market. Same as Figure 5(a), the least operational cost can be reached only when the truck speeds are considered as decision variables.

5.3.2 Effect of truck speed on flexibility In this sub-subsection, we investigate the relation between truck speed and the flexibility of the FTPDP-RF. Two experiments, 𝐸1 and 𝐸2, are designed. The purposes of 𝐸1 and 𝐸2 are to investigate how the changes of truck number and the rate of reserved orders can affect the flexibility of the FTPDP-RF, respectively. To implement the experiments, we compute the average operational cost of twenty 10-order instances in two scenarios: one is to set the truck speed as an economic speed, 70 kilometer/hour, and the other needs to determine the truck speed on each road segment. The results are shown in Figure 6 and Figure 7, respectively. In Figure 6, y1 and y2 are the trend lines of the operational cost when the number of trucks changes from enough to insufficient. Specifically, y1 is corresponding to the scenario when truck speed is set to be a given parameter, and y2 is corresponding to the scenario when the truck speed is treated as a decision variable. In Figure 7, y1 and y2 are the trend lines of the operational cost when reserved orders are included in the instances and excluded from the instances, respectively. As can be seen, the lower operational cost can be reached when the truck speeds are considered as decision variables. Meanwhile, in both of the scenarios, the rate of y1 is greater than that of y2. The 27

phenomenon indicates that when constraints are changed, the FTPDP-RF model that treating truck speed as a decision variable is more stable.

Average Cost ( RMB )

4000 3500 3000 2500 2000

Outsourcing Charges

1500

Fuel Cost

1000 500 0 60 70 80 90 variable Truck Speed ( kilomiters per hour )

(a) High price of outsourcing

Average Cost (RMB)

4000 3500 3000 2500 2000

Outsourcing Charges

1500

Fuel Cost

1000 500 0 60 70 80 90 variable Truck Speed (kilometer per hour)

(b) Low price of outsourcing Figure 5. Effect of truck speed on the cost

28

3600 70 km/h

3400

Cost(RMB)

3200

varible speed

3000

70 km/hytrend

2800

line

2600

varible speed y trend line

2400 2200 2000 Enough Trucks

Insufficient Trucks

Figure 6. Relation between truck number and flexibility 3600

70 km/h

Cost (RMB)

3400 3200

varible speed

3000

70 km/hytrend line

2800

variable speed y trendline

2600 2400 2200 2000 without reserved orders with reserved orders

Figure 7. Relation between reserved order and flexibility

6 Conclusions and future suggestions This paper proposes a flexible truckload pickup and delivery problem considering fuel consumption. In this problem, some of the orders can be outsourced, while some of the orders are reserved; the operational cost is considered as a nonlinear function of the decision variable of truck speed, which causes the challenge of solving the proposed problem. To address the challenge, we describe the proposed problem as a directed graph and then analyze its attributes to reduce the solution space. Subsequently, we establish a mixed integer nonlinear programming model which is based on truck visiting time. Furthermore, we present a UTWPB method and provide a lower bound of the proposed problem by relaxing partitioned time windows. Finally, the UTWPB method is evaluated on randomly instances, and 29

sensitivity analysis is conducted. The results of the numerical experiments are as follows. (1) The proposed model and solution method are valid and efficient. For small scale instances, each of which contain less than 50 orders, both the SDB method and the UTWPB method can obtain near-optimal solutions (average gap less than 1%) within a short computational time, and the SDB method performs better; for real-life scale instances, each of which contains 100-200 orders, the proposed UTWPB method outperforms the existing TWPB method and the SDB method. (2) Compared with that considering truck speed as a parameter, the strategy that takes truck speed as a decision variable can significantly reduce the operational cost and upgrade the flexibility of the FTPDP-RF. Because the application of the UTWPB is limited by the scale of instances, meta-heuristic or column generation algorithms will be designed in our future work. Other extensions of the FTPDP-RF, considering more complex situations, such as dynamic situation and traffic congestion, will also be conducted.

Appendix. Proof of the main results ∗ Proof of Proposition 3.1. According to Demir et al. (2012), when 𝑣 ≤ 𝑣𝑢𝑐 , 𝐹(𝑣) is a ∗ monotone decreasing function, and when 𝑣 ≥ 𝑣𝑢𝑐, 𝐹(𝑣) is a monotone increasing function; and then, Proposition 3.1 is proved. ∎ ∗ ≤ 𝑣𝐿, 𝑣𝑐∗ = 𝑣𝐿, and [𝑣𝐿, Proof of Proposition 3.2. According to proposition 3.1, if 𝑣𝑢𝑐 ∗ ∗ ∗ < 𝑣𝑈, then 𝑣𝐿 < 𝑣𝑐∗ ≤ 𝑣𝑢𝑐 𝑣𝑈] can be rewritten as [𝑣𝑐∗ , 𝑣𝑈]; and if 𝑣𝐿 < 𝑣𝑈 ≤ 𝑣𝑢𝑐 , or 𝑣𝐿 < 𝑣𝑢𝑐 . 𝐿 ∗ Hence, Proposition 3.2 is true if only if 𝑣 < 𝑣𝑐 . Define 𝑅𝑆(𝑆,𝑣𝑐) as any road segment of the feasible solution 𝑆 for the FTPDP-RF where the truck speed on the road segment 𝑅𝑆 is 𝑣𝑐. The objective value of 𝑆 is 𝐶(𝑆) = 𝐹𝐶(𝑅𝑆(𝑆,𝑣𝑐)) + 𝐶(𝑆\𝑅𝑆(𝑆,𝑣𝑐)) where 𝐶(𝑆) represents the operational cost of the feasible solution 𝑆, 𝐹𝐶(𝑅𝑆(𝑆,𝑣𝑐)) is the cost of fuel consumption on 𝑅𝑆 and 𝐶 (𝑆 \𝑅𝑆(𝑆,𝑣𝑐)) is the operational cost except for the fuel consumption cost on 𝑅𝑆. Define 𝑆 ∗ as another solution which is same as 𝑆 except for the truck speed on 𝑅𝑆. Then, 𝑅𝑆(𝑆 ∗ ,𝑣𝑐∗ ) represents RS of 𝑆 ∗ where the truck speed on 𝑅𝑆 is 𝑣𝑐∗ . Hence, the only difference between 𝑅𝑆(𝑆,𝑣𝑐) and 𝑅𝑆(𝑆 ∗ ,𝑣𝑐∗ ) is the corresponding truck speed on RS. If 𝑣𝑐 < 𝑣𝑐∗ , then 𝑆 ∗ must be also a feasible solution, because a truck at a higher speed will spend less travel time, thereby not violating time windows of corresponding orders. Since 𝐶 (𝑆 ∗ \𝑅𝑆(𝑆 ∗ ,𝑣𝑐∗ )) = 𝐶(𝑆\𝑅𝑆(𝑆,𝑣𝑐)), the objective value of 𝑆 ∗ is 𝐶(𝑆 ∗ ) = 𝐹𝐶(𝑅𝑆(𝑣𝑐∗ )) + 𝐶 (𝑆 ∗ \𝑅𝑆(𝑆 ∗ ,𝑣𝑐∗ )) = 𝐹𝐶(𝑅𝑆(𝑣𝑐∗ )) + 𝐶(𝑆\𝑅𝑆(𝑆,𝑣𝑐)). Arguing by contradiction, assume that S is the optimal solution, while there exists 𝑣𝑐∗ and 𝑣𝑐 < 𝑣𝑐∗ on 𝑅𝑆. Since S and 𝑆 ∗ are same except for the truck speed, then 𝐹𝐶 (𝑅𝑆 (𝑆,𝑣𝑐)) = 𝛿𝐹(𝑣𝑐) and 𝐹𝐶(𝑅𝑆 (𝑆 ∗ ,𝑣𝑐∗ )) = 𝛿𝐹(𝑣𝑐∗ ). Considering 𝐹(𝑣) is a monotone decreasing function, when 𝑣 < 𝑣𝑐∗ , 𝐹(𝑣𝑐∗ ) < 𝐹(𝑣𝑐) and then 𝐹𝐶(𝑅𝑆 (𝑆 ∗ ,𝑣𝑐∗ )) < 𝐹𝐶 (𝑅𝑆 (𝑆,𝑣𝑐)). Then, for the feasible solution 𝑆 ∗ , 𝐹𝐶(𝑅𝑆 (𝑆 ∗ ,𝑣𝑐∗ )) + 𝐶(𝑆\𝑅𝑆(𝑆,𝑣𝑐))=𝐶(𝑆 ∗ ) < 𝐶(𝑆) = 𝐹𝐶(𝑅 𝑆(𝑆,𝑣𝑐)) + 𝐶(𝑆\𝑅𝑆(𝑆,𝑣𝑐)). Thus, 𝑆 is not the optimal solution. This contradicts the assumption; therefore, Proposition 3.2 follows the proof.∎ Proof of Proposition 3.3. Since 𝐹(𝑡) = ∆1𝑡 + ∆2(𝑤 + 𝑓)𝑑 + ∆3𝑑3/𝑡2, take the first derivative of 𝐹(𝑡), 𝐹(𝑡)‘ = ∆1 ―2∆3𝑑3/𝑡3; and take the second derivative of 𝐹(𝑡), 𝐹(𝑡)

30

" = 6∆3𝑑3/𝑡4. Because 𝑡>0, ∆3 > 0, 𝑑>0, we have 𝐹(𝑡)‘’ > 0; and then 𝐹(𝑡) is a convex 3

∗) = 0. ∎ function. formula (14) follows from 𝐹′(𝑡) = ∆1 ―2∆3(𝑑/𝑡𝑢𝑐

Proof of Proposition 3.4. According to formulas (3) and (4), if 𝑡 ∈ [𝑑/𝑣𝑈𝑟, 𝑑/𝑣𝐿𝑟], then ∗ we have: (i) if 𝑣𝑈 ≤ 𝑣𝑢𝑐 , then 𝑣𝐿𝑟 = 𝑣𝑈𝑟 = 𝑣𝑈, thus 𝑡 = 𝑑/𝑣𝑈𝑟 = 𝑑/𝑣𝐿𝑟, and 𝐹(𝑡) is a ∗ ∗] ∗ ∗ ≤ 𝑣𝑈, then 𝑡 ∈ [𝑑/𝑣𝑈𝑟, 𝑑/𝑣𝑢𝑐 = 𝑡𝑢𝑐 constant value; (ii) if 𝑣𝐿 ≤ 𝑣𝑢𝑐 , thus 𝑡 ≤ 𝑑/𝑣𝑢𝑐 , and then ∗ 𝐹(𝑡) is a monotone decreasing function of 𝑡 following Corollary 3.1; and (iii) if 𝑣𝑢𝑐 ≤ 𝑣𝐿, ∗ ∗ = 𝑡𝑢𝑐 then 𝑡 ∈ [𝑑/𝑣𝑈𝑟, 𝑑/𝑣𝐿𝑟],where 𝑑/𝑣𝐿𝑟 = 𝑑/𝑣𝐿. Since 𝑡 ≤ 𝑑/𝑣𝐿 ≤ 𝑑/𝑣𝑢𝑐 , 𝐹(𝑡) is a monotone decreasing function of 𝑡 following Corollary 3.1. To sum up, when 𝑡 ∈ [𝑑/𝑣𝑈𝑟 , 𝑑/𝑣𝐿𝑟], 𝐹(𝑡) is a monotone non-increasing function and then 𝑡𝑐∗ is the optimal travel time for minimizing fuel consumption. ∎ Proof of Proposition 3.5. Considering the constraint of truck speed, travel time 𝑡 can be neither more than 𝑑/𝑣𝐿𝑟 nor less than 𝑑/𝑣𝑈𝑟, 𝑡 ∈ [𝑑/𝑣𝑈𝑟, 𝑑/𝑣𝐿𝑟]. If 𝑇 ≥ 𝑑/𝑣𝐿𝑟, then 𝑡𝑐∗ = 𝑑/𝑣𝐿𝑟 following Proposition 3.4; if 𝑑/𝑣𝑈𝑟 ≤ 𝑇 < 𝑑/𝑣𝐿𝑟, then 𝑡 ∈ [𝑑/𝑣𝑈𝑟, 𝑇], and 𝑈𝑟 ∗ thereby 𝑡𝑐 = 𝑇 following Proposition 3.4; and 𝑇 < 𝑑/𝑣 violates the constraint of truck speed; thus, 𝑡𝑐∗ is nonexistent when 𝑇 < 𝑑/𝑣𝑈𝑟.∎ Proof of Proposition 4.1. Three types of arcs (including node arcs) will be generated after an iteration: (1) arcs that the time windows of two nodes are not partitioned, (2) arcs that the time window of one node is partitioned but the other is not, (3) arcs that the time windows of both nodes are partitioned. After the time windows of arcs are partitions, the distance and truckload are not changed; therefore, the fuel consumption of an arc is dependent on the time windows. For each arc of type (1), either the upper bound or the lower bound of the fuel consumption remains unchanged. First, we calculate the upper bound of the fuel consumption of each arc of type (3). According to the description of bisection method in Subsection 4.1, at each iteration, a time window may be divided into two sub time windows. For example, [𝑎𝑟𝑖, 𝑏𝑟𝑖] is divided into [𝑎𝑟𝑖, 𝑐𝑟𝑖] and [𝑐𝑟𝑖, 𝑏𝑟𝑖], where 𝑐𝑟𝑖 = 𝑎𝑟𝑖 +(𝑏𝑟𝑖 ― 𝑎𝑟𝑖)/2 . In the same way, another time window [𝑎𝑟𝑗, 𝑏𝑟𝑗] is divided into [𝑎𝑟𝑗, 𝑐𝑟𝑗] and [𝑐𝑟𝑗, 𝑏𝑟𝑗]. According to Subsection 4.2, the fuel consumption of an arc (𝑖, 𝑗) of Model-UB is dependent on the upper bounds of the sub time windows, 𝑏𝑟𝑖 and 𝑏𝑟𝑗; therefore, at the next iteration, there exists an arc (𝑖′, 𝑗′) the fuel consumption of which is equal to that of the arc (𝑖, 𝑗) because 𝑏𝑟𝑖 and 𝑏𝑟𝑗 still exist at the next iteration. Meanwhile, some new arcs are generated by the new generated sub nodes; thus, the decision space of Model-UB will increase. For each arc of type (2), we have the similar result. Consequently, the optimal solution obtained by solving Model-UB is monotonically non-increasing. Second, we calculate the lower bound of the fuel consumption of each arc of type (3). According to Subsection 4.3, the fuel consumption of an arc (𝑖, 𝑗) of Model-LB is dependent on 𝑎𝑟𝑖 and 𝑏𝑟𝑗; and then, at the next iteration, when the time windows [𝑎𝑟𝑖, 𝑏𝑟𝑖] and [𝑎𝑟𝑗, 𝑏𝑟𝑗] are partitioned, we can get four sub time windows, [𝑎𝑟𝑖, 𝑐𝑟𝑖] and [𝑐𝑟𝑖, 𝑏𝑟𝑖], [𝑎𝑟𝑗, 𝑐𝑟𝑗] and [𝑐𝑟𝑗, 𝑏𝑟𝑗], which can generate four sub arcs by the following combinations: 𝑎𝑟𝑖 and 𝑐𝑟𝑗, 𝑎𝑟𝑖 and 𝑏𝑟𝑗, 𝑐𝑟𝑖 and 𝑐𝑟𝑗, and 𝑐𝑟𝑖 and 𝑏𝑟𝑗. Exclude the arc (𝑖, 𝑗) which starts at 𝑎𝑟𝑖 and end at 𝑏𝑟𝑗, the fuel consumption of the rest three sub arcs are greater than that of arc (i, j) according to Subsection 4.3 because 𝑐𝑟𝑖 ≥ 𝑎𝑟𝑖 and 𝑐𝑟𝑗 ≤ 𝑏𝑟𝑗 and then less travel time can be used. Therefore, 31

none of the lower bounds of the four sub arcs is less than that of the origin arc (𝑖, 𝑗). This may result in a non-decreasing optimal solution. Moreover, a truck must enter and leave the same sub node; compared with the time windows before being partitioned, the end time at the node may be earlier and the start time at the node may be later. For example, according the description of the relaxed constraints in Subsection 4.3, before the time windows are partitioned, a truck starts at 𝑎𝑟𝑖 at node 𝑖, and ends at 𝑏𝑟𝑗 at node j, and then starts at node 𝑗 at 𝑎𝑟𝑗; after the time windows are partitioned, a truck should start at node 𝑗 at 𝑐𝑟𝑗, where 𝑐𝑟𝑗 ≥ 𝑎𝑟𝑗 according to constraints (31) and (32). These constraints may result in an infeasible route according to the definitions of 𝑁𝐿𝐵 and 𝐴𝐿𝐵; thereby, an arc with a higher fuel consumption may be selected. For each arc of type (2), we have the similar result. In conclusion, the optimal solution obtained by solving Model-LB is monotonically non-decreasing.∎

References Bektaş, T., Laporte, G., 2011. The Pollution-Routing Problem. Transportation Research Part B: Methodological 45(8), 1232-1250. Bolduc, M.-C., Renaud, J., Boctor, F., 2007. A heuristic for the routing and carrier selection problem. European Journal of Operational Research 183(2), 926-932. Bolduc, M.C., Renaud, J., Boctor, F., Laporte, G., 2008. A perturbation metaheuristic for the vehicle routing problem with private fleet and common carriers. Journal of the Operational Research Society 59(6), 776-787. Chu, C.-W., 2005. A heuristic algorithm for the truckload and less-than-truckload problem. European Journal of Operational Research 165(3), 657-667. Côté, J.-F., Potvin, J.-Y., 2009. A tabu search heuristic for the vehicle routing problem with private fleet and common carrier. European Journal of Operational Research 198(2), 464-469. Dabia, S., Demir, E., Woensel, T.V., 2017. An Exact Approach for a Variant of the Pollution-Routing Problem. Transportation Science 51(2), 607-628. Dantzig, G.B., Ramser, J.H., 1959. The Truck Dispatching Problem. Management Science 6(1). Demir, E., Bektaş, T., Laporte, G., 2012. An adaptive large neighborhood search heuristic for the Pollution-Routing Problem. European Journal of Operational Research 223(2), 346-359. Demir, E., Bektaş, T., Laporte, G., 2014. The bi-objective Pollution-Routing Problem. European Journal of Operational Research 232(3), 464-478. Derigs, U., Pullmann, M., Vogel, U., Oberscheider, M., Gronalt, M., Hirsch, P., 2012. Multilevel neighborhood search for solving full truckload routing problems arising in timber transportation. Electronic Notes in Discrete Mathematics 39, 281-288. 32

Euchi, J., 2017. The vehicle routing problem with private fleet and multiple common carriers: Solution with hybrid metaheuristic algorithm. Vehicular Communications 9, 97-108. Euchi, J., Chabchoub, H., Yassine, A., 2011. New Evolutionary Algorithm Based on 2-Opt Local Search to Solve the Vehicle Routing Problem with Private Fleet and Common Carrier. International Journal of Applied Metaheuristic Computing 2(1), 58-82. Fagerholt, K., Laporte, G., Norstad, I., 2010. Reducing fuel emissions by optimizing speed on shipping routes. Journal of the Operational Research Society 61(3), 523-529. Franceschetti, A., Demir, E., Honhon, D., Van Woensel, T., Laporte, G., Stobbe, M., 2017. A metaheuristic for the time-dependent pollution-routing problem. European Journal of Operational Research 259(3), 972-991. Franceschetti, A., Honhon, D., Van Woensel, T., Bektaş, T., Laporte, G., 2013. The time-dependent pollution-routing problem. Transportation Research Part B: Methodological 56, 265-293. Gahm, C., Brabänder, C., Tuma, A., 2017. Vehicle routing with private fleet, multiple common carriers offering volume discounts, and rental options. Transportation Research Part E: Logistics and Transportation Review 97, 192-216. Gendreau, M., Nossack, J., Pesch, E., 2015. Mathematical formulations for a 1-full-truckload pickup-and-delivery problem. European Journal of Operational Research 242(3), 1008-1016. Jabali, O., Van Woensel, T., de Kok, A.G., 2012. Analysis of Travel Times and CO2Emissions in Time-Dependent Vehicle Routing. Production and Operations Management 21(6), 1060-1074. Jula, H., Dessouky, M., Ioannou, P., Chassiakos, A., 2005. Container movement by trucks in metropolitan networks: modeling and optimization. Transportation Research Part E: Logistics and Transportation Review 41(3), 235-259. Koç, Ç., Bektaş, T., Jabali, O., Laporte, G., 2014. The fleet size and mix pollution-routing problem. Transportation Research Part B: Methodological 70, 239-254. Kopfer, H.W., Schönberger, J., Kopfer, H., 2014. Reducing greenhouse gas emissions of a heterogeneous vehicle fleet. Flexible Services and Manufacturing Journal 26(1-2), 221 - 248. Krajewska, M.A., Kopfer, H., 2009. Transportation planning in freight forwarding companies. European Journal of Operational Research 197(2), 741-751. Kramer, R., Maculan, N., Subramanian, A., Vidal, T., 2015. A speed and departure time optimization algorithm for the pollution-routing problem. European Journal of Operational Research 247(3), 782-787. Kratica, J., Kostic, T., Tosic, D., Dugosija, D., Filipovic, V., 2012. A genetic algorithm for the routing and carrier selection problem. Computer Science and Information Systems 9(1), 49-62. 33

Lai, M., Cai, X., Hu, Q., 2017. An iterative auction for carrier collaboration in truckload pickup and delivery. Transportation Research Part E: Logistics and Transportation Review 107, 60-80. Li, J., Rong, G., Feng, Y., 2015. Request selection and exchange approach for carrier collaboration based on auction of a single request. Transportation Research Part E: Logistics and Transportation Review 84, 23-39. Liu, R., Jiang, Z., Liu, X., Chen, F., 2010. Task selection and routing problems in collaborative truckload transportation. Transportation Research Part E: Logistics and Transportation Review 46(6), 1071-1085. Potvin, J.Y., Naud, M.A., 2011. Tabu search with ejection chains for the vehicle routing problem with private fleet and common carrier. Journal of the Operational Research Society 62(2), 326-336. Wang, W.F., Regan, A.C., 2002. Local truckload pickup and delivery with hard time window constraints. Transportation Research Part B: Methodological 36(2), 97-112. Zhang, R., Lu, J.-C., Wang, D., 2014. Container drayage problem with flexible orders and its near real-time solution strategies. Transportation Research Part E: Logistics and Transportation Review 61, 235-251. Zhang, R., Yun, W.Y., Kopfer, H., 2010. Heuristic-based truck scheduling for inland container transportation. OR Spectrum 32(3), 787-808. Zhang, R., Yun, W.Y., Moon, I., 2009. A reactive tabu search algorithm for the multi-depot container truck transportation problem. Transportation Research Part E: Logistics and Transportation Review 45(6), 904-914. Ziebuhr, M., Kopfer, H., 2016. Solving an integrated operational transportation planning problem with forwarding limitations. Transportation Research Part E: Logistics and Transportation Review 87, 149-166.

34

Highlights 

A mixed integer nonlinear programming model is established



Several strategies are used to reduce the solution space



An efficient updated time window partition based method is designed



A lower bound of the problem is proposed



Determining truck speed can enhance the flexibility of transportation planning

35