The heterogeneous green vehicle routing and scheduling problem with time-varying traffic congestion

Transportation Research Part E 88 (2016) 146–166

Contents lists available at ScienceDirect

Transportation Research Part E journal homepage: www.elsevier.com/locate/tre

The heterogeneous green vehicle routing and scheduling problem with time-varying traffic congestion Yiyong Xiao a, Abdullah Konak b,⇑ a b

School of Reliability and System Engineering, Beihang University, Beijing 100191, China Information Sciences and Technology, Penn State Berks, Tulpehocken Road, P.O. Box 7009, Reading, PA 19610-6009, United States

a r t i c l e

i n f o

Article history: Received 20 July 2015 Received in revised form 31 December 2015 Accepted 26 January 2016

Keywords: CO2 emissions Vehicle routing Vehicle scheduling Green logistics Mixed integer programming Hybrid optimization Matheuristics

a b s t r a c t The green vehicle routing and scheduling problem (GVRSP) aims to minimize green-house gas emissions in logistics systems through better planning of deliveries/pickups made by a fleet of vehicles. We define a new mixed integer liner programming (MIP) model which considers heterogeneous vehicles, time-varying traffic congestion, customer/vehicle time window constraints, the impact of vehicle loads on emissions, and vehicle capacity/range constraints in the GVRSP. The proposed model allows vehicles to stop on arcs, which is shown to reduce emissions up to additional 8% on simulated data. A hybrid algorithm of MIP and iterated neighborhood search is proposed to solve the problem. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Climate scientists undoubtedly point out that the concentration of greenhouse gases in the atmosphere, particularly CO2 emitted by human activities, is the primary cause of global warming. The Keeling Curve (http://keelingcurve.ucsd.edu/) shows that the concentration of CO2 in the atmosphere has been growing in the past half-century and even faster in the last decade. According to the US Environment Protection Agency (US EPA, 2014), transportation is one of the major sectors emitting CO2, counting for 27% of the total US emissions in 2013, and almost three-quarters of these emissions from transportation were due to road transportation. Today, all large and growing urban areas experience high levels of road traffic congestion. Road traffic congestion, often companied with frequent acceleration and deceleration, contributes to CO2 emissions significantly (Barth and Boriboonsomsin, 2008; Franceschetti et al., 2013). According to the International Road Transport Union (IRTU), road traffic congestion could increase CO2 emissions by 300% and was responsible for 100 billion liters of wasted fuel, or 250 billion tons of CO2 emissions in the U.S. (IRTU, 2012). In China, cars and trucks are also responsible for 40% of the PM2.5 air pollution, which is caused by particles less than 2.5 micrometers in diameter, observed in many central Chinese cities such as Beijing where 175 days in 2014 were reported to have dangerous levels of PM2.5 air pollution (Chai, 2015). Traffic congestion may be caused by various factors such as road capacities, rush-hour, school zones, road work, tidal roads, accidents, and traffic controls. Whatever the cause is, traffic congestion is almost always time dependent, which we refer to as time-varying traffic conditions, indicating that the average travel speed on a road may change from one period to another with predictable patterns. Considering time-varying traffic conditions in vehicle routes and schedules provides an ⇑ Corresponding author. E-mail addresses: [email protected] (Y. Xiao), [email protected] (A. Konak). http://dx.doi.org/10.1016/j.tre.2016.01.011 1366-5545/Ó 2016 Elsevier Ltd. All rights reserved.

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

147

opportunity for reducing vehicle emissions in logistic systems by implementing Operations Research techniques to avoid traffic congestion and meet delivery requirements at the same time. Several recent papers related to this issue can be found in the literature, including Figliozzi (2010, 2012), Kuo (2010), Franceschetti et al. (2013), Kwon et al. (2013), Gaur et al. (2013), Demir et al. (2014a, 2014b), Soysal et al. (2015), and Xiao and Konak (2015a). In recent years, green logistics, which aims at improving environmental sustainability by reducing fossil fuel consumption in freight transportation systems, has been a subject undergoing intense study and attracted the attention of Operations Research professionals (Demir et al., 2014a,b; Rodrigue et al., 2013; Eglese and Bektasß, 2014; Bektasß et al., 2016). In the context of the vehicle routing problem (VRP), traditional logistics activities are often optimized in a way that customer requirements are met at minimum cost in terms of solely monetary measures. On the other hand, the green vehicle routing problem (GVRP) is characterized by the objective of balancing the environmental and economic costs by implementing effective vehicle routes and schedules. Considering the fuel consumption cost (or the CO2 emissions cost) in the GVRP may lead to quite different logistic plans for customer assignments, route selections, and schedules than the traditional VRP (Lin et al., 2014; Xiao and Konak, 2015a). In this paper, we present a comprehensive mixed integer linear programming (MILP) model to formulate the GVRP for a fleet of heterogeneous vehicles operating under time-varying traffic conditions. The formulated problem is referred to as Heterogeneous Green Vehicle Routing and Scheduling Problem (HGVRSP). In the MILP model, we use total CO2 emissions as the objective function to be minimized by various operational decisions including customer-vehicle assignment, route selection, and travel time scheduling. The main advantages of our MILP model lie on three aspects: (1) idle/waiting times can be scheduled at any point during tours (i.e., the non-stopping assumption is relaxed), (2) heterogeneous vehicles are considered with individualized features including vehicle types, CO2 emissions rates/models, load capacities, fuel tank capacities, and time availabilities, and (3) vehicle emissions are calculated by considering their dynamic payload weights and travel speeds. In addition, the model also considers general traffic congestion patterns with several periods and time windows dictated by customer requirements or the availability of vehicles. Therefore, the proposed model is very comprehensive and general. For the solution approach, we propose a hybrid algorithm of partial MIP optimization and iterative neighborhood search (P-MIP-INS) based on the concepts from variable neighborhood search. Comprehensive computational experiments with different problem sizes are carried out to study the effectiveness and efficiency of the P-MIP-INS algorithm with respect to several well-known construction heuristics from the literature. The rest of the paper is organized as follows. In Section 2, we review the previous related work and point out the drawbacks of previous models with the non-stop assumption. In Section 3, we remodel the time-dependent travel and formulate a comprehensive model for the HGVRSP. In Section 4, we carry out computational experiments on 20 small-sized problems to compare the solutions found with three different objective functions—(1) the total travel distance, the total travel time, and the total CO2 emissions. In Section 5, we develop a hybrid solution approach of partial MIP optimization and variable neighborhood search (P-MIP-INS) for the proposed HGVRSP model in this paper and introduce our modifications on the construction heuristics from the literature. In Section 6, computational experiments are carried out on 140 problem instances (small, medium, and large). Finally, we conclude the paper in Section 7. 2. Literature review Since its introduction by Dantzig and Ramser (1959), the VRP has been extensively studied in the literature. Various versions of the VRP have been developed for different applications, such as pickup and delivery VRP, capacitated VRP, multiple depot VRP, VRP with time windows, split delivery VRP, and time-dependent VRP. Surveys on various VRP formulations and algorithms can be found in Laporte (1992), Laporte et al. (2000), Toth and Vigo (2002, 2014), Golden et al. (2008), and Eksioglu et al. (2009). In the following, we review the previous VRP work considering time constraints. 2.1. Vehicle routing and scheduling problems The vehicle routing and scheduling problem (VRSP), first defined by Bodin and Golden (1981), refers to the case when customers have specific service time requirements. Therefore, a solution to the VRSP includes both vehicle routes and schedules. The route of a vehicle is defined by the sequence of pickup and/or delivery points which the vehicle visits, and the schedule of the vehicle indicates the associated arrival and departure times to and from these points. Solomon (1983) first presented an MIP formulation for the VRSP with time window constraints (VRSPTW) where some customers impose deadlines or earliest time constraints on their pickups/deliveries. Since then, the VRSPTW have been a subject of intensive research efforts for both construction heuristics (Solomon, 1987; Christofides et al., 1979; Van Landeghem, 1988; Malandraki, 1989; Ioannou et al., 2001) and meta-heuristics (Garcia et al., 1994; Thangiah, 1995; Vidal et al., 2013). Comprehensive surveys on the VRPTW can be found in Olli and Michel (2005a,b) and Desaulniers et al. (2014). The time-dependent VRPTW (TD-VRPTW) is an extension of VRSPTW such that the travel time of an arc is considered as a function of the departure time. Malandraki (1989) and Malandraki and Daskin (1992) first formulated the TD-VRPTW as a MIP model and proposed a greedy nearest-neighbor heuristic as the solution approach. Malandraki and Dial (1996) proposed a restricted dynamic programming heuristic for constructing the route for a time-dependent traveling salesperson problem. Chen et al. (2006) studied the real-time TD-VRPTW. Maden et al. (2010) proposed a heuristic algorithm for the single-depot

148

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

VRPTW utilizing the Road TimeTableTM and reported 7% reduction in CO2 emissions in a case study. Figliozzi (2012) proposed a heuristic construction/improvement algorithm for the TD-VRP with hard or soft time window constraints. In recent years, another raising issue related to the TD-VRP is how to minimize CO2 emissions in logistic systems. Kara et al. (2007) proposed a cost function to minimize energy consumption in the context of the VRP. To the best of our knowledge, the earliest paper addressing CO2 emissions optimization in the context of the VRSP under time-varying traffic conditions is the work by Figliozzi (2010). Kuo (2010) also proposed a model for minimizing fuel consumption in the TD-VRP and used a simulated annealing algorithm as the solution approach. Bektasß and Laporte (2011) presented the pollution routing problem (PRP) where the optimal travel speed of each arc (i; j) is determined to minimize a comprehensive objective function including fuel consumption, emissions, and driver costs. Demir et al. (2012) developed a two-stage heuristic algorithm to solve the PRP for large-size problem instances with up to 200 nodes. In their approach, an adaptive large neighborhood search heuristic is used to solve the VRPTW in the first stage, and the optimal travel speed for each selected arc is determined in the second stage. Jabali et al. (2012) developed a model for an emission-based TD-VRP where a uniform speed limit is determined for all arcs. Franceschetti et al. (2013) extended the PRP to Time-Dependent PRP (TD-PRP) by considering traffic congestion variations in two periods (a congested period and a pre-congested period). Lin et al. (2014) presented a review of the GVRSP and described its relationships with other versions of the TD-VRP. Kramer et al. (2015) proposed a matheuristic approach combining a local search-based metaheuristic with MIP for the PRP and reported better performance than the algorithm proposed by Maden et al. (2010) on the same benchmark problems. Soysal et al. (2015) proposed a MILP model for a two-echelon time-dependent capacitated VRP with the objective of minimizing fuel consumption such that freights are delivered to customers via intermediate depots rather than direct shipments. Several other researchers also focus on minimizing CO2 emissions or fuel consumption in the context of the VRP, but without considering time-varying traffic congestion (e.g., Xiao et al. (2012), Kwon et al. (2013), Gaur et al. (2013), and Demir et al. (2014a,b)). Xiao et al. (2012) incorporated fuel consumption due to a vehicle’s load into the capacitated VRP. Gaur et al. (2013) studied the cumulative VRP with the objective of minimizing fuel consumption and proposed an approximation algorithm for the cumulative VRP. Kwon et al. (2013) developed a heterogeneous-fleet VRP model with a fixed carbon emissions rate (per unit distance) and considered the carbon emissions trade cost in the objective function. Demir et al. (2014a,b)proposed an adaptive large neighborhood search algorithm (ALNS) to minimize the fuel consumption and the driving time with Pareto optimality for the PRP. Xiao and Konak (2015a) presented a MIP model for the TD-PRP with hierarchical objectives (including the total CO2 emissions, the arrival time, and the total travel time on road) and a weighted tardiness penalty under a carbon-trade system and developed a simulated annealing algorithm as the solution approach. Xiao and Konak (2015b) also proposed a hybrid approach of genetic algorithms and dynamic programming to the TD-PRP by considering discrete time points for scheduling vehicles.

Fig. 1. Existing methods of modeling time-dependent travel times.

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

149

2.2. Models with time-varying travel times In most VRSPTW cases, the travel time t ij of each arc (i; j) is assumed to be constant throughout the planning horizon. In time-varying traffic conditions, however, travel times are time-dependent, i.e., each arc (i; j) has a dynamic travel time t kij depending on departure time k from node i. Malandraki (1989) and Malandraki and Daskin (1992) first treated the travel time as a step function of the departure time. The step function was also used in the time-dependent traveling salesperson problem by Malandraki and Dial (1996) and in the real-time TD-VRPTW by Chen et al. (2006). Hill and Benton (1992) assumed that each customer node has a constant travel speed around its vicinity and used the average speed between node pairs to estimate time-dependent travel times. The earlier models for the TD-VRSP did not include the non-passing property (Ahn and Shin, 1991) or the First-in-First-out (FIFO) property (Ichoua et al., 2003), which is generally accepted by the literature for modeling time-dependent travel times in the latter. The non-passing property states that among two identical vehicles that travel arc (i; j), the one that departs from node i first will always arrive at node j earlier than the other vehicle. To guarantee the non-passing property, Ichoua et al. (2003) modeled travel speed (instead of travel time) as a step function of time and used a simple method to obtain continuous travel times. Fleischmann et al. (2004) proposed a smoothed version of the time-dependent travel time function and showed that if the slope of the function is less than 45°, then the non-passing property will hold. Eglese et al. (2006) proposed Road TimeTableTM, which is updated in real-time, to calculate the time-dependent travel time between two nodes at any time of the day. Figliozzi (2012) assumed constant speed intervals to guarantee the FIFO property for the TD-VRP with hard or soft time windows. Jabali et al. (2012) considered a uniform speed limit for all arcs, and used a continuous piecewise linear function over the starting time to estimate the time-dependent travel time. Franceschetti et al. (2013) used the same method with Jabali et al. (2012) to estimate time-dependent travel times. In the following, we summarize the existing methods of modeling the time-dependent travel time under time-varying traffic congestion.

(i) A step function of the departure time as shown in Fig. 1(A). This model was first proposed by Malandraki (1989) and used by Malandraki and Daskin (1992), Hill and Benton (1992), Malandraki and Dial (1996), Chen et al. (2006). The time-dependent travel times depend on the departure time. (ii) A smoothed step function of the departure time as shown in Fig. 1(B). This model was first proposed by Ahn and Shin (1991) and improved in Fleischmann et al. (2004). The angle of the slopes are set lower than 45° to ensure the non-passing property. (iii) A function of the departure time and arrival times in multiple intervals (periods), proposed by Hill and Benton (1992). The travel speed is the mean value of the speed around the vicinity of customer i in the period of departure time b and the speed around the vicinity of customer j in the period of arrival time e, i.e., rijb ¼ ðri;TðbÞ þ ri;TðeÞ Þ=2.The speed in the vicinity of each customer is a step function of travel time of arc (i; j) is then calculated as b þ dij =r ijb . (iv) A step function of the time-dependent travel speed as shown in Fig. 1(C). This model was first proposed by Ichoua et al. (2003) and used by Eglese et al. (2006), Maden et al. (2010), Kuo (2010), Figliozzi (2012), and Soysal et al. (2015). The model assumes that the travel speed is constant within a period and changes at the boundary of two consecutive periods, and then the model uses a simple cumulative method to obtain the travel time. This method has been widely used in the recent work related to time-varying traffic congestion. (v) The model in Fig. 1(D) was proposed first by Jabali et al. (2012) and used by Franceschetti et al. (2013). The planning horizon is divided into three consecutive periods: (1) a congestion period with a steady speed v c (a constant), (2) a transition period, and (3) a free-flow period with a steady speed v f (a variable). The travel time of arc (i; j) with distance d and departing time w is calculated as follows:

tðwÞ ¼

8 > < vd=vv c f

c

if w 6 ða  d=v c Þ

þ

d

v f ða  wÞ þ v f > : d=v f

if ða  d=v c Þ < w < a : þ

if w P a

2.3. Models with/without the non-stopping rule All of the models summarized above, including Malandraki (1989), Ahn and Shin (1991), Malandraki and Daskin (1992), Hill and Benton (1992), Malandraki and Dial (1996), Ichoua et al. (2003), Chen et al. (2006), Eglese et al. (2006), Maden et al. (2010), Figliozzi (2012, 2010), Kuo (2010), Jabali et al. (2012), Franceschetti et al. (2013), and Soysal et al. (2015), do not allow a vehicle to stop during the travel of an arc. Once a vehicle starts traveling arc (i; j), the vehicle must travel the whole arc distance without a waiting/idle time during its travel. This property refers to as the non-stopping assumption as

150

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

discussed before. The pollution routing models by Ichoua et al. (2003), Figliozzi (2010, 2012), and Franceschetti et al. (2013) allow idle/waiting times only at customer premises while the models by Kuo (2010), Bektasß and Laporte (2011), and Jabali et al. (2012) do not allow idle/waiting times even at customer premises. The validity of the non-passing property of Ahn and Shin (1991) and FIFO property of Ichoua et al. (2003) is also based on the non-stopping assumption. The non-stopping assumption is reasonable in the TD-VRSP models where travel time or travel distance related costs are usually used as the objective function to be minimized. However, it can be a drawback when the objective of minimizing CO2 emissions (or fuel consumption) is considered. Bellow, we explain this drawback using the example given in Fig. 2 where (A) illustrates a typical rush-hour travel speed pattern with three periods such that each period has a different maximum travel speed that is below the most economical speed. A vehicle is going to travel arc (i; j) starting at time t in period 1. Suppose the vehicle cannot travel the whole arc in the remaining time of period 1, which is e1  t, and has to travel it in the following periods. Fig. 2(B) shows the non-stopping travel time path, which is a continuous tour starting at time t in period 1, going through period 2, and ending at time t 00 in period 3. While Fig. 2(C) presents another travel time path where the vehicle starts at the same time t in period 1, stops at the beginning of period 2, then restarts to travel the remaining distance in period 3 (after the rush-hour), and finally arrives at the destination at time t0 . Obviously, Fig. 2(B) will lead to a shorter travel time than that of Fig. 2(C), i.e., t00 < t 0 , which is preferable when the objective is minimization of the total travel time. However, if the objective is emissions reduction, Fig. 2(C) will result in lower emissions than that of Fig. 2(B) because the vehicle travels the same distance with a lower emissions rate per unit distance. In Fig. 2(C), note that the vehicle is able to travel with travel speeds closer to the economical level. Our computational experiments in Section 4.3 show that CO2 emissions can be reduced by 8.0% on the average (up to 41.1%) if idle times are allowed during tours. This relaxation of the non-stopping property is particularly important for arcs representing long-haul roads. 2.4. Models with heterogeneous fleets Another extension of the VRP related to the problem defined in this paper is known as Heterogeneous Vehicle Routing Problem (HVRP), first presented by Golden et al. (1984) and lately reviewed by Baldacci et al. (2008) and Koç et al. (in press), in which one must additionally decide on the fleet composition and the customer/route selections for each vehicle. The current state-of-art solution approach for the HVRP with time window is a hybrid evolutionary algorithm (HEA) with an ALNS procedure developed by Koç et al. (2015). The heterogeneity of the vehicle fleet is particularly important for the Green VRPs since different vehicles may have different emission rates, and logistic companies may have a heterogeneous fleet of various vehicles, including new and old ones, different engine types (fuel, electric, or hybrid), and even different sources. Recent works related to the HVRP with CO2 emissions reduction can be found in Kwon et al. (2013) and Koç et al. (2014). However, these approaches do not consider time-varying traffic congestion. In the HGVRSP, we consider a general description of vehicle heterogeneity, including various load capacities, fuel tank (or electric power) capacities, dynamic emissions rates (the fuel consumption rates and fuel-to-CO2 conversion rates) with respect to vehicle loads and speeds, and availabilities (with individual time-window). To the best of our knowledge, a general heterogeneity of fleet vehicles are considered for the first time in the context of Green VRSP with time-varying traffic congestions. 3. Problem description and formulation The HGVRSP is defined on a complete directed graph G = (N; A) where N = {0, 1, . . . , n} is the set of nodes, and A = {ði; jÞ : i – j 2 N} is the set of arcs. Node 0 denotes the depot, and the other nodes represent customers. If a vehicle h is used, it starts its tour from the depot and returns back to it after visiting all the assigned customers. Each customer i is associated with a non-negative demand ri in weight, a service time si , and a service time window [Si ; Ei ]. The shortest path for each pair of two customers, i.e., arc ði; jÞ 2 A, is a fixed route that has a fixed distance Dij . The road network has

(A)

(B) Fig. 2. A comparison of arc travel with and without stop in the middle.

(C)

151

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

time-varying traffic conditions. The planning horizon is partitioned into a set of periods M = {[b1 ; e1 ], [b2 ; e2 ], . . ., [bm ; em ]} where bk and ek represent the starting and ending time of period k, respectively. The travel speed (v ijk ) on each arc (i; jÞ 2 A in each period k 2 M is assumed to be known and constant during the period. Each customer must be served by only a single vehicle. The set of fleet vehicle types (U) includes heterogeneous vehicles. Each vehicle type u 2 U has a maximum load capacity C u , a fuel tank capacity limit Lu , a given function of fuel consumption rate FCRu ðv ; f ) depending on travel speed v, load weight f, and a known fuel-to-CO2 conversion rate cu . In addition, each vehicle h is associated with a time window [F h ; T h ] indicating its availability to serve customers. The objective function is to minimize the total CO2 emissions of all used vehicles. 3.1. Modeling time-dependent travel of arcs Under time-varying traffic conditions, the travel time, fuel consumption, and CO2 emissions of traveling arc (i; j) depend on when the vehicle departs from node i, which is referred to as time-dependent travel. As discussed in the introduction section, the widely accepted non-stopping assumption in the literature does not guarantee the optimal paths and schedules in terms of fuel consumption or CO2 emissions. In the following, we introduce a mathematical model for the HGVRSP by relaxing the non-stopping assumption. In order to determine the travel schedule of arc (i; j), we introduce non-negative continuous variables li which denotes the starting time of travelling arc (i; j) and dijk which indicates the distance traveled on arc (i; j) in period k. If arc (i; j) is selected in P a tour, k2M dijk ¼ Dij must be satisfied. The travel time in each period can be calculated by sijk ¼ dijk =v ijk where v ijk is the P travel speed of arc (i; j) in period k. Then, the total travel time (tij ) of arc (i; j) is calculated as t ij ¼ k2M sijk . Note that sijk may not fully utilize all available time in period k. There may be multiple breaks while traveling arc (i; j), which are waiting/idle times scheduled to avoid traffic congestion, as shown in Fig. 3. We use binary variable xijk to indicate whether arc (i; j) is traveled in period k (xijk ¼ 1) or not (xijk ¼ 0). Variable xijk is bound to dijk as follows:

dijk 6 Dij xijk

8k2M

Note that if arc (i; j) is not selected by any vehicle, then xijk ¼ 0 8 k 2 M. The relationship between sijk and departure time li and arrival time aj can be expressed as follows:

li 6 ek  sijk þ em ð1  xijk Þ 8 k 2 M and aj P bk þ sijk  em ð1  xijk Þ 8 k 2 M; where bk and ek are the beginning and ending time of period k, respectively, and em denotes the ending time of the last period (i.e., period m). There may be many pairs of feasible (li ; aj ) that will result in minimum CO2 emissions. This can be explained as the time flexibility that the driver can utilize to schedule waiting/idle time in the tour as long as the arrival time aj satisfies Sj 6 aj 6 Ej where [Sj ; Ej ] is the time window of customer j. The approach summarized in this section leads to a linear model because v ijk is a parameter that can be calculated based on the historical data before solving the problem. Thereby, it is possible to express the non-linear relationship between CO2 emissions and travel speed as a linear constraint without requiring auxiliary binary variables for linearization. To achieve a linear formulation, it is assumed that the travel speed of an arc is constant within a period. This assumption is well justified in the case of real-life operations of commercial fleets that must strictly follow speed limits. In addition, the travel speed is

(A)

(B)

Fig. 3. An example of a time path without the non-stopping assumption.

152

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

mainly dictated by traffic and road conditions, not a decision variable in many real-life road systems. The values of v ijk can be derived from historical records and geographical information systems including traffic data. For the sake of presentation clarity, the planning horizon is divided into the same set of periods (i.e., M) for all arcs in the MIP model of the HGVRSP presented in this section. However, we need to point out that period set M can be tractably replaced with individual period set M ij for each arc (i; j) in the model, and arcs can have different numbers of periods with unequal lengths. In addition, the number of periods can be increased to better approximate actual traffic congestion patterns. These changes will make the formulation more realistic in an area with various traffic congestion patterns, but reduce the computational tractability of the problem significantly. Our solution approach given in Section 5 can also work with multiple set of periods. It is also needed to point out that in the time-dependent travel modeled above each arc (i; j) is treated as a unique and non-dividable segment like many traditional arc-level formulation of the VRP. Therefore, we simply assume arc (i; j) represents the shortest path from origin i to destination j, and the v ijk is the average travel speed on the arc in period k. Thereby, the travel time can be quickly estimated by dijk =v ijk . In the case of urban areas, a street-level model that allows the arc (i; j) to be represented as many different physical street paths with different traffic patterns may be more realistic but also significantly increases the computational complexity as well. For this case, a pre-calculated function table will be needed for each arc (i; j) for accessing the non-linear travel time of any given sub-distance (i.e., dijk ) started at any given time point, instead of calculating it directly by dijk =v ijk . This will obviously increase the computational cost for solving the MIP model formulated in the next sub-section. 3.2. The description of the mixed integer linear programming model Below we formulate the HGRVSP as a MIP. The parameters of the model are as follows: /⁄- - - - - - - - - - - - - -notations for customer locations and requirements- - - - - - - - - - - - - - - - - -⁄/ i; j index of nodes, i = 0, 1, 2, . . . , n (the depot is represented by 0) n total number of nodes (including the depot) N set of nodes including the depot N0 set of nodes excluding the depot, N 0 ¼ N n f0g A set of arcs formed by all pairs of nodes, A ¼ fði; jÞ : 8 i 2 N; j 2 N; i – jg Dij distance of arc (i; j) ri demand of customer i si service time of customer i ½Si ; Ei  time-window for starting to serve customer i /⁄- - - - - - - - - - - - - -notations for time-varying traffic conditions- - - - - - - - - - - - - - - - - - - - - - - - - - - -⁄/ k index of period, k = 1, 2, . . . , m m number of time periods, m = card(M) M set of periods, k 2 M ½bk ; ek  beginning and ending time of period k v ijk travel speed on arc (i; j) in period k /⁄- - - - - - - - - - - - - -notations for heterogeneous vehicles - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -⁄/ h index of vehicle, h = 1, 2, . . . , q H set of vehicles, h 2 H q total number of vehicles, q = card(H) ½F h ; T h  The available time window for vehicle h /⁄- - - - - - - - - - - - - -notations for heterogeneous vehicle types- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -⁄/ u index of vehicle type, u = 1, 2, . . . , card(U) U set of vehicle types Cu maximum payload of vehicle type u Lu Fuel tank capacity of vehicle type u whu {0, 1} value indicating whether vehicle h belongs to type u /⁄- - - - - - - - - - - - - -notations for fuel consumption and CO2 emissions - - - - - - - - - - - - - - - - - - - -⁄/ g ijku fuel consumption rate (FCR) of vehicle type u for traveling arc (i; j) in period k with zero payload /u additional fuel consumption rate of vehicle type u when traveling with one more unit of load (kg) cu fuel-to-CO2 conversion rate of vehicle type u (kg/kg) The binary and continuous decision variables of the model are defined as follows: X ij {0, 1} variable indicating whether arc (i; j) is traveled (1) or not (0) yijh {0, 1} variable indicating whether arc (i; j) is traveled by vehicle h (in any period) (1) or not (0) xijkh {0, 1} variable indicating whether arc (i; j) is traveled in period k by vehicle h (1) or not (0)

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

dijkh

continuous continuous continuous continuous continuous

sijkh li ai f ijh

variable variable variable variable variable

indicating indicating indicating indicating indicating

the the the the the

153

traveled distance of arc (i; j) in period k by vehicle h time spent on arc (i; j) in period k by vehicle h departure time from node i arrival time (or service starting time) at node i load of vehicle h on arc (i; j)

We use the fuel consumption function developed by Demir et al. (2012) and Franceschetti et al. (2013), which is based on the comprehensive module emissions model (CMEM) by Barth et al. (2005) and Barth and Boriboonsomsin (2008), to estimate vehicles’ fuel consumption rates (FCR) as a function of travel speed v and travel load f. For a heavy-duty vehicle type u traveling on a road with no slope, the FCR function is given as follows:

FCRu ðv ; f Þ ¼ au v 1 þ bu v 2 þ cu þ uu f ;

ð1Þ

where au , bu , cu ; and uu ; are the coefficients related to the vehicle type u. We can calculate the FCR value, denoted as g ijku , for an unloaded vehicle of type u traveling arc (i; j) in period k using the first three terms of Eq. (1). In addition, we can calculate the FCR due to the vehicle’s travel load using the last term of Eq. (1) as uu f ijh where f ijh is the travel load vehicle h when it is traveling arc (i; j). Therefore, the objective function that minimizes the total CO2 emissions can be expressed as follows:

Min CO2 ¼

X X XX

cu g ijku whu dijkh þ

h2H ði;jÞ2A k2M u2U

X X X cu uu whu f ijh Dij

ð2Þ

h2H ði;jÞ2A u2U

The first term represents the emissions induced by the vehicles, and the second term is the emissions induced by the vehicles’ loads. Note that the objective function is linear because only dijkh and f ijh are variables, and all other parameters can be pre-calculated before solving the problem. We use continuous variable dijkh as the primary decision variable to divide arc (i; j) into several segments that can be traveled in different periods. An advantage of this approach is that the whole arc does not need to be continuously traveled; thereby, traffic congestion in some periods can be scheduled off to save fuel or reduce CO2 emissions. In real-life scenarios, a vehicle can experience different levels of traffic congestion while traveling a road, particularly a long haul one. Since CO2 emissions are pre-calculated based on periods, such scenarios can be considered in the proposed model. This is another advantage of the modeling approach in this paper. In most of the existing PRP models in the literature, vehicles are not often allowed to make stops once they have started their tours (Kuo, 2010; Bektasß and Laporte, 2011; Jabali et al., 2012; Soysal et al., 2015) or are allowed to have a waiting time at customer nodes, but not during traveling arcs (Ichoua et al., 2003; Figliozzi, 2010; Franceschetti et al., 2013). Note that the proposed model can also be easily modified to prevent stops on arcs. Variable f ijh is a dynamic load of vehicle h along its tour and determined by the demands of the customers that the vehicle is going to serve. Binary variable xijkh indicates whether variable dijkh has a positive value (indicated by xijkh ¼ 1) or zero (by xijkh ¼ 0) as follows:

ð1Þ dijkh 6 Dij  xijkh

8 ði; jÞ 2 A; k 2 M; h 2 H

Constraint (1) requires that variable xijkh must be 1 when variable dijkh is a positive number. Note that dijkh cannot exceed the distance of arc (i; j). To express route constraints, two other binary variables, yijh and X ij , are introduced. Variable yijh is tightly bound to xijkh by Constraints (2) and (3), and variable X ij is tightly bound to yijh by Constraint (4).

ð2Þ yijh P xijkh ð3Þ yijh 6

m X xijkh

8 ði; jÞ 2 A; k 2 M; h 2 H 8 ði; jÞ 2 A; h 2 H

k¼1

ð4Þ X ij ¼

q X yijh

8 ði; jÞ 2 A

h¼1

The following Constraints (5)–(8) are used to construct the routes of heterogonous vehicles. Constraint (5) ensures that each vehicle departs the depot at most once. Constraint (6) and (7) are typical node balance constraints guaranteeing that each customer is visited only once by a single vehicle. Constraint (8) ensures the same vehicle enters to and departs from a customer node. Constraint (9) ensures that the total distance of a selected arc (i; j) is fully traversed. Constraint (10) states that if arc (i1 ; i2 ) is traveled by vehicle h in time period k2 , all arcs departing from node i2 cannot be traveled in a time period k1 such that k1 < k2 . Constraint (10) also eliminates any sub-tours in the solution.

154

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

ð5Þ

n X y0jh 6 1

8h2H

j¼1

ð6Þ

n X X ij ¼ 1

8 i 2 N0

j¼0

ð7Þ

n X X ij ¼ 1

8 j 2 N0

i¼0

ð8Þ

n X

0

8 h 2 H; j 2 N0

yji0 h

0

i 2N;i –j

i2N:i–j

ð9Þ

n X

yijh ¼

q m X X

8 ði; jÞ 2 A

dijkh ¼ X ij Dij

k¼1 h¼1

ð10Þ xi2 ;i3 ;k1 ;h 6 2  xi1 ;i2 ;k2 ;h  yi1 ;i2 ;h

8 h 2 H; ði1 ; i2 Þ; ði2 ; i3 Þ 2 A; k1 ; k2 2 K; k1 < k2

Travel time sijkh in minutes is calculated by Constraint (11). Note that Constraint (11) is linear because travel speed v ijk is a parameter (in km/h) in the model. Constraint (12) restricts the total travel time of vehicle h in period k to the length of the period.

ð11Þ ð12Þ

sijkh ¼ 60  dijkh =v ijk X

8 ði; jÞ 2 A; k 2 M; h 2 H

sijkh 6 ek  bk

8 k 2 M; h 2 H

ði;jÞ2A

Constraints (13)–(17) given below are used to compute departure time li and arrival time ai of node i.

ð13Þ li 6 ek  sijkh þ em  ð1  xijkh Þ

8 ði; jÞ 2 A; k 2 M; h 2 H

ð14Þ aj P bk þ sijkh  em  ð1  xijkh Þ ð15Þ aj P li þ

q m X X

8 ði; jÞ 2 A; k 2 M; h 2 H

sijkh  em  ð1  X ij Þ

8 ði; jÞ 2 A

k¼1 h¼1

ð16Þ ai þ si 6 li ð17Þ Si 6 ai 6 Ei

8 i 2 N0 8 i 2 N0

Constraint (13) ensures that if arc (i; j) is traveled in period k by vehicle h (i.e., xijkh ¼ 1), the departure time of node i must be before the end time of period k minus the travel time of arc (i; j) in that period. Similarly, Constraint (14) states that the arriving time of node j must be greater than or equal to the start time of period k plus the travel time of arc (i; j) in that period. Note that these two constraints are enforced only for xijkh ¼ 1. For xijkh ¼ 0, Constraints (13) and (14) still hold because em is a large number that ensures the validity of the inequalities. Constraint (15) is a disjunctive constraint to calculate the arrival times at the nodes. Note that Constraint (15) is only active for the arcs that are selected (i.e., X ij ¼ 1). Constraint (16) ensures that the vehicle’s departure time from a node is after its arriving time plus the service time. Constraint (17) imposes the upper and lower bounds on arrival time ai to serve customer i. So far, Constraints (1)–(17) and the objective function in Eq. (2) define a basic MIP model for the HGVRSP under timevarying traffic conditions. In the following, we introduce the constraints to consider real-life scenarios such as the effect of vehicle loads on fuel consumption, vehicles’ capacity limits, vehicles’ fuel tank capacity (or travel length) limits, customers’ time-windows, and vehicles’ individual availability time windows. Constraint (18) is used to compute the vehicle’s payload along the tour as well as to eliminate any sub-tours. When this constraint is used, Constraint (10) can be removed from the model. Constraint (19) limits the vehicle’s load by its maximum load capacity C h and also forces f ijh to be zero if arc (i; j) is not traveled by vehicle h (i.e., yijh ¼ 0).

ð18Þ

X

f ijh 

h2H;i2N;i–j

ð19Þ f ijh 6 C h  yijh

X

f jih ¼ r j

h2H;i2N;i–j

8 h 2 H; ði; jÞ 2 A

8 j 2 N0

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

155

Constraint (20) enforces vehicles’ maximum travel distance limits in terms of the fuel tank capacity. The right-hand side of Constraint (20) is the total amount of fuel consumed by vehicle h and the left-hand side is its fuel tank capacity.

ð20Þ

X X XX X X whu  Lu P whu  g ijku  dijkh þ uu  whu  f ijh  Dij u2U

ði;jÞ2A k2K u2U

8h2H

ði;jÞ2A u2U

Finally, considering that not all vehicles are always available for the entire planning horizon, each vehicle h has a time window ½F h ; T h  indicating its availability. This constraint is quite realistic from a practical point of view because logistics companies may face vehicle shortages in busy seasons, rent vehicles with different schedules, or outsource some services to drivers. The following Constraints (21) to (24) enforce the availability of vehicle h.

ð21Þ F h 6 ek  s0jkh þ em  ð1  x0jkh Þ

8 j 2 N0 ; k 2 M; h 2 H

ð22Þ F h 6 aj  s0jkh þ em  ð1  x0jkh Þ

8 j 2 N0 ; k 2 M; h 2 H

ð23Þ T h P bk þ si0kh  em  ð1  xi0kh Þ ð24Þ T h P li þ si0kh  em  ð1  xi0kh Þ

8 i 2 N0 ; k 2 M; h 2 H 8 i 2 N 0 ; k 2 M; h 2 H

Finally, Constraint (25) defines the ranges of the decision variables in the model.

(

ð25Þ

X ij 2 f0; 1g; yijh 2 f0; 1g; xijkh 2 f0; 1g;

dijkh P 0; sijkh P 0; ai P 0; li P 0; f ijh P 0

8 i; j; k; h

Objective function Eq. (2) and Constraints (1)–(25) form a MIP model for the time-dependent, heterogeneous, and capacitated vehicle routing and scheduling problem with time-windows for CO2 emissions minimization. Because the model is linear, small-sized problems can be solved optimally by most existing MIP solvers such as CPLEX and Lingo. The presented model can be easily extended to a case where arcs have different number of time periods and/or time periods with different lengths. Such an extension would be useful in a large area where traffic congestion might follow different patterns depending on locations. In addition to the constraints above, the following inequalities are useful to cut off infeasible solutions and develop more efficient heuristics. X ij þ X ji 6 1 X ij P yijh X P xijkh Pijn i¼1 yi0h 6 1 P Pn whu C u P whu ni¼0 j¼1 r j  yijh Pn 1 6 i¼1 X i0 6 q Pn 1 6 j¼1 X 0j 6 q

8 8 8 8 8

ði; jÞ 2 A ði; jÞ 2 A; h 2 H ði; jÞ 2 A; k 2 M; h 2 H h2H h 2 H; u 2 U

//One arc be traveled by at most once in one direction. //If an arc (i; j) is not selected, all yijh should be zero. //If an arc (i; j) is not selected, all xijhk should be zero. //Each vehicle returns to depot only once //The load capacity of a vehicle should not be violated //The total number of returns to the deport is in [1, q] //The total number of departures from the deport is in [1, q]

4. Comparisons of the CO2 emissions and the traditional VRP objectives The HGRVSP model described in the previous section aims to minimize CO2 emissions. However, traditional VRP objectives such as the total travel distance or the total travel time can also be used as follows:

Min Trav el TimeðTTÞ ¼

X X X

tijkh

ð3Þ

h2H ði;jÞ2A k2M

Min Trav el DistanceðTDÞ ¼

X

Dij  X ij

ð4Þ

ði;jÞ2A

In this section, we compare the solutions found with three different objective functions, namely the total CO2 emissions, total travel distance, and total travel time in order to demonstrate to what extent CO2 emissions can be reduced at the expense of longer travel distances or travel times. We used AMPL/CPLEX (version 12.4.0.1) to solve the HGVRSP model with the three objective functions defined above. 4.1. Description of the test problems We adopted the PRP data generated by Demir et al. (2012), which includes real geographical distances between randomly selected 200 cities in UK, random time windows, demands, and service times for customers (http://www.apollo.management.soton.ac.uk/prplib.htm). In this test problem set, there are 9 problem groups of different sizes ranging from 10 customers to 200 customers, and each problem group contains 20 instances. The average distance

156

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

between the nodes is about 80 km, and the planning horizon is 9 h. We divided the 9-h delivery window into five equallength periods and associated each period with a travel speed range (km/h) as follows: [70, 100], [30, 60], [10, 30], [30, 60], and [70, 100] to create time-varying traffic conditions based on a V-shape traffic congestion pattern. The traffic speed for each arc in each period was randomly and uniformly selected from the corresponding range. For each instance, the number of available vehicles is set to dn=5e. We used two types of heavy-duty vehicles with different emissions models. Type I is a 6.35-ton (curb-weight) diesel truck whose fuel consumption function was developed by Franceschetti et al. (2013) based on the Comprehensive Modal of Emissions Model (CMEM) of Barth et al. (2005) and Barth and Boriboonsomsin (2008). Type II is a 9.58-ton (curb-weight) diesel truck whose fully-loaded (26 tons) fuel consumption rate was determined as 7.96 MPG through real-life empirical studies by Coyle (2007), and its curb-weight CO2 emissions model is calculated based on the MEET model, which is developed by Hickman (1999) and widely used in the literature. The maximum capacities of Type I and II vehicles were set as 10 tons and 15 tons, respectively, about 1.5 times of their curb-weight (heavy-duty trucks can have a maximum payload up to 2.5–3.0 times of their curb-weights). In order to make the problems as close to real-life as possible, customers’ demands were trebled, and the fuel tank limits were set as 150 kg and 200 kg. Thereby, the vehicles may be subject to both load capacity and fuel capacity limits. Respectively, the emission rates of Type I and II vehicles under different speeds and payloads are numerically expressed as follows:

ERI ðv ; f Þ ¼ 9:75862  v 1 þ 2:90264  105  v 2 þ 0:142139 þ 2:23840  105  f

ð5Þ

ERII ðv ; f Þ ¼ 871  16:0  v þ 0:143  v 2 þ 32031  v 2 þ 1:79961  105  f

ð6Þ

where ERI and ERII are emissions rates in kg/km, v is speed in km/h, and f is the load in kg. In the fleet, we set odd numbered vehicles as Type I and the even numbered ones as Type II. We assigned time windows [0, 540] and [30, 540] for Type I and Type II vehicles, respectively. For the same travel speed, Type II vehicles have a higher emissions rate than Type I vehicles when they are empty or lightly-loaded. Therefore, Type I vehicles may have lower emissions than Type II vehicles with lighter loads. However, since Type II vehicles have a lower marginal emissions per additional unit of payload, i.e., uI = 2.23840  105 kg/km/kg versus uII = 1.79961  105 kg/km/kg, they would be superior to Type I vehicles with heavy loads. Fig. 4 illustrates the emissions functions of Type I and II vehicles with respect to different travel speeds when they are empty (Fig. 4A) and with respect to different payload weights when they travel at a constant speed of 70 km/h (Fig. 4B).

4.2. Comparison of different objectives We optimally solved 20 instances of problem group 10  5 (n  m) using the two vehicle types for the three objectives independently, the total CO2 emissions (CO2), total travel time (TT), and total travel distance (TD). All instances (excluding No. 11 which had no feasible solution) were optimally solved and the results are summarized in Table 1. The optimal results found for the objectives CO2, TT, and TD are given under columns CO2 , TT⁄, and TD⁄, respectively. For the optimal solutions found using a particular objective, the corresponding values of the other objectives are also provided in Table 1. For example in Problem No 1, the optimal solution found with the CO2 objective has CO2 = 477.1, and the values of the other two objectives are TT = 562.8 and TD = 421.6 for this optimal solution. Columns D% indicate the percent deviations of the objective values from their corresponding optimal values. In Problem No 1 for example, the optimal solution found using the CO2 objective has 53.5% higher travel time than its optimal travel time which is TT⁄ = 366.7. When the TD objective is used, the schedule of a vehicle is irrelevant. For these solutions, we determined the schedules that yielded the best and worst possible CO2 emissions and the best and worst possible total travel times. These values are also given in subcolumns Min. and Max. of CO2 and TT under the Total travel distance objective in Table 1. Finally, the t:ðmÞ column is the CPU time in minutes. If the CO2 values are used as benchmarks, the solutions that minimized the TT objective had 44.7% higher CO2 emissions on the average. Similarly, the solutions that minimized the TD objective yielded 19.4% higher CO2 emissions on the average, ranging between 5.1% and 33.8%. The comparisons in Table 1 justify the merits of using emissions-based objective functions under time varying traffic conditions. In particular, solutions even with the same travel distance may have very different CO2 emissions due to scheduling of the vehicles in different periods. The significant gap between the best and worst CO2 emissions of the solutions found by the TD objective supports the importance of proper scheduling decisions in addition to routing decisions in order to minimize CO2 emissions under time-varying traffic conditions. Note that solutions found by the CO2 objective lead to an average increase of 2.1% in travel distance and an average increase of 45.4% in travel time. The computational time varied significantly depending on instances with the overall average of 12.4 CPU minutes. We were not able to solve larger-sized problems even within a 24-h of CPU lime. CPU time requirement increased significantly even with the addition of a single node. For larger-sized problems, therefore, it is practical to obtain near-optimal solutions using heuristics.

157

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166 1.5

0.8

1.4 1.3 1.1 1 0.9 0.8 Vehicle II (9.58T)

0.7 0.6 0.5 0.4

Vehicle I (6.35T)

CO2 emission rate (Kg/Km)

CO2 emission rate (kg/km)

1.2

0.7

0.6 Vehicle II

0.5 Vehicle I

0.3 0.2

0.4

0.1 0

0

20

40

60

80

100

120

Vehicle's speed (km/hour)

0

2

4

6

8

10

12

14

16

Payload (Tone)

(B) With constant speed v=70km/h

(A) With empty load

Fig. 4. Vehicle types’ CO2 emission rates with different speeds and payloads.

4.3. Comparison with the non-stopping assumption Next, we consider the non-stopping assumption by adding the following constraint to the model.

ð26Þ aj  li 6

XX

sijkh þ em  ð1  X ij Þ

8 ði; jÞ 2 A

k2K h2H

The constraint above enforces the time difference between arrival time aj and departure time li to be equal to the traveling time of arc (i; j) so that stops are not allowed while traveling arc (i; j). We solved the same 20 instances used in the previous section. In Table 2, we compare the solutions found with and without the non-stopping assumption. It can be observed that the solutions with the non-stopping assumption are always worse or equal to (only instance No. 20) the solutions without the non-stopping assumption. The deviation from optimality is 8.0% on the average. For the worst case (instance No. 19), the deviation is even as high as 41.1%, i.e. 41.1% more CO2 will be emitted compared to the optimal solution without the non-stopping assumption. The model presented in this paper allows vehicles to stop while traveling arcs. This is particularly important in real-life cases where some arcs may represent long-haul roads. 4.4. A restricted-MIP method The HGVRSP model defined in this paper can be solved to optimality only for small-sized problems. The computational efficiency of CPLEX degrades dramatically when the problem size grows slightly larger. In our experiments, it took 12 min on the average to solve small problems with size of 10  5, but larger problems with size of 15  5 could not be solved optimally within a 24-h CPU time limit. Therefore, finding near-optimal solutions within acceptable computational time is an alternative approach for the HGVRSP studied in this paper. In this section, we introduce a restricted-MIP method to be used as a benchmark for the heuristic algorithms described in Section 5. Our objective is to demonstrate that good solutions can be found by considering only a subset of arcs. We also use a variation of this approach to design a hybrid algorithm to solve large-sized problems. In the restricted-MIP method, the arc selection variables, i.e., X ij , are restricted to only a small set of promising arcs, eliminating most of the other unpromising, long-distance arcs from the consideration. Each node reserves only a fixed number of its promising arcs. The underlying justification of this approach is that nodes are much more likely to connect to their close neighbors rather than farther ones in an optimal solution. In Section 4.2, the distance-optimized solutions have 5.1% (for the best case) higher CO2 emissions than those of CO2-optimized solutions on the average. This observation also implies that a shorter travel distance is likely to result in lower CO2 emissions as well. Therefore, we use the following two inequities to restrict the arcs that can be selected for each node to only its k nearest neighbors.

158

Table 1 Experimental results of three different objectives for 20 instances of 10  5 group (two vehicles). CO2 emissions

Total travel time

Total travel distance

No.

CO2

TT

D%

TD

D%

t (m)

TT⁄

CO2

D%

TD

D%

t (m)

TD⁄

Min.

D%

Max.

D%

Min.

D%

Max.

D%

t (m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Avg

477.1 786.7 630.5 584.1 517.1 781.1 591.4 738.6 508.0 580.3 – 512.3 588.6 540.4 364.6 540.3 534.5 446.2 443.3 434.8

562.8 736.5 655.6 633.6 640.1 706.3 674.7 687.1 610.7 680.9 – 632.0 679.9 655.2 494.2 587.3 666.3 544.6 599.9 534.1

53.5 31.0 39.6 37.2 64.9 16.5 34.9 19.8 56.4 48.4 – 50.4 38.0 38.2 65.0 61.5 67.1 44.9 34.6 60.2 45.4

421.6 602.9 516.2 486.4 463.1 590.5 513.1 567.8 461.5 485.2 – 466.2 510.7 507.8 354.4 447.7 479.5 429.9 433.2 387.0

0.5 2.1 1.5 0.0 3.6 1.0 0.0 0.1 1.2 0.0 – 1.3 0.0 6.0 12.0 3.5 2.3 3.5 0.0 1.8 2.1

11.6 12.6 21.4 7.4 31.2 14.0 4.8 2.1 3.0 3.0 – 2.8 6.2 11.0 47.1 44.2 2.5 11.4 0.4 7.5 12.9

366.7 562.1 469.7 461.8 388.1 606.2 500.0 573.6 390.5 458.7 – 420.1 492.6 474.0 299.6 363.7 398.8 375.9 445.7 333.4

770.0 888.6 841.3 877.9 812.3 880.5 726.4 833.3 829.3 867.4 – 743.5 932.3 860.1 542.9 839.2 787.2 693.4 621.8 708.1

61.4 13.0 33.4 50.3 57.1 12.7 22.8 12.8 63.2 49.5 – 45.1 58.4 59.2 48.9 55.3 47.3 55.4 40.3 62.9 44.7

432.6 599.8 537.6 504.1 460.6 590.5 513.1 567.8 472.9 533.2 – 465.8 526.6 514.3 348.2 465.6 490.5 437.0 434.0 404.9

3.1 1.5 5.7 3.7 3.1 1.0 0.0 0.1 3.7 9.9 – 1.2 3.1 7.3 10.0 7.6 4.6 5.2 0.2 6.5 4.1

1.0 5.2 1.8 1.3 1.6 3.4 2.2 0.4 0.4 0.8 – 0.3 2.5 10.5 8.1 1.1 0.4 1.2 1.3 1.0 2.3

419.5 590.7 508.6 486.3 446.9 584.5 513.1 567.5 456.1 485.2 – 460.3 510.7 479.1 316.4 432.7 468.8 415.2 433.2 380.2

495.5 829.7 641.4 652.9 571.8 810.2 591.4 783.1 516.9 580.3 – 584.1 604.9 552.4 377.9 632.5 539.6 472.5 443.3 456.3

3.9 5.5 1.7 11.8 10.6 3.7 0.0 6.0 1.8 0.0 – 14.0 2.8 2.2 3.6 17.1 1.0 5.9 0.0 4.9 5.1

733.7 927.8 696.6 812.9 708.0 900.7 729.3 867.2 776.6 720.3 – 627.7 762.0 813.8 506.1 726.1 760.6 670.9 593.6 643.1

53.8 17.9 10.5 39.2 36.9 15.3 23.3 17.4 52.9 24.1 – 22.5 29.5 50.6 38.8 34.4 42.3 50.4 33.9 47.9 33.8

390.0 662.1 577.7 488.4 489.7 630.8 500.0 593.4 456.2 516.3 – 551.7 566.1 491.5 335.5 421.5 422.8 440.0 500.3 381.0

6.4 17.8 23.0 5.8 26.2 4.1 0.0 3.4 16.8 12.6 – 31.3 14.9 3.7 12.0 15.9 6.0 17.1 12.3 14.3

597.9 723.7 680.6 650.2 667.7 711.6 680.0 669.3 655.6 726.5 – 600.7 694.3 655.8 505.5 592.4 656.1 579.8 614.6 563.6 12.8

63.1 28.8 44.9 40.8 72.1 17.4 36.0 16.7 67.9 58.4 – 43.0 41.0 38.4 68.7 62.9 64.5 54.3 37.9 69.0 48.7

22.0 57.6 33.8 11.2 41.3 24.7 11.3 3.5 6.3 9.1 – 4.5 10.9 37.7 67.2 49.7 8.4 21.3 2.2 7.9 22.7

CO2

Note: bold face indicates optimal values.

TT

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

Obj

159

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166 Table 2 Comparison between solutions with and without non-stopping assumption. 10  5

CO2

With non-stopping

No. 1 2 3 4 5 6 7 8 9 10

477.1 786.7 630.5 584.1 517.1 781.1 591.4 738.6 508.0 580.3

10  5

CO2

D%

t (m)

No.

512.6 805.4 649.7 631.0 534.9 840.8 663.5 774.3 515.7 599.1

7.4 2.4 3.0 8.0 3.5 7.6 12.2 4.8 1.5 3.2

42 32 23 37 55 43 27 4 3 8

11 12 13 14 15 16 17 18 19 20 Avg

CO2

With non-stopping CO2

D%

t (m)

– 512.3 588.6 540.4 364.6 540.3 534.5 446.2 443.3 434.8

– 533.2 633.0 593.3 421.3 543.6 573.5 497.4 625.4 434.8

– 4.1 7.5 9.8 15.5 0.6 7.3 11.5 41.1 0.0

8 17 85 895 123 11 24 27 24

557.9

599.1

8.0

78

Xij 6 k þ n  ð1  X ij Þ 8 ði; jÞ 2 A; i – 0; j – 0; and Xij 6 k0 þ n  ð1  X ij Þ 8 ði; jÞ 2 A; ði ¼ 0 or j ¼ 0Þ: where

Xij

the distance ranking of arc (i; j) with respect to node i, e.g., Xij =1 means that arc (i; j) is the shortest distance arc emanating from node i k parameter that restricts the number of the adjacent arcs to customer nodes to its k nearest customers. k0 parameter that restricts the depot (i = 0) being connected to only its k0 nearest customers.

The constraints above reduce the solution space by fixing variables X ij to zero for all arcs that have Xij > k. Clearly, this approach can be directly implemented on the problem data by eliminating such arcs from arc set A. We solved the 20 small-sized problems using CPLEX with k0 = n and k = 2, . . . , 7 to investigate the tradeoff between computational efficiency and solution quality. The detailed results are given in Table 3, and the summary of the results is illustrated in Fig. 4. It can be observed that as k decreases from its original value 10 to 5, the computational time improves dramatically by 87% while the solution quality degrades only by 1.5%. For k = 4, the computational time even improves by 96% at the expense of 3.8% degradation in the objective function. This experiment indicates that the restricted-MIP method is able to find good solutions in relatively short times (see Fig. 5).

Table 3 Results of restricted-MIP on problem group 10  5 (two vehicles). 10  5

Opt.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Avg

CO2 k=2

No. 477.1 786.7 630.5 584.1 517.1 781.1 591.4 738.6 508.0 580.3 – 512.3 588.6 540.4 364.6 540.3 534.5 446.2 443.3 434.8

– – – – – – – – 572.8 625.4 – – – 580.4 410.5 – – – – – 10.1%

t (s)

k=3

t (s)

k=4

t (s)

k=5

t (t)

k=6

t (s)

k=7

t (s)

– – – – – – – 1 1 – – – 1 1 – – – – – 1

507.4 – 832.1 584.1 551.9 – 620.9 738.6 508.0 610.4 – 627.7 637.0 545.1 371.5 575.1 640.6 503.6 693.3 502.0 11.7%

4 – 9 3 2 – 6 4 6 7 – 3 6 4 3 5 12 2 9 1 7

477.1 801.9 678.1 584.1 548.0 798.7 591.4 738.6 508.0 580.3 – 627.7 591.8 540.4 371.5 551.7 552.3 503.6 461.3 464.6 3.8%

27 21 50 76 22 89 21 12 11 16 12 27 7 16 46 12 45 12 7 7 28

477.1 801.9 641.4 584.1 535.5 781.1 591.4 738.6 508.0 580.3 – 512.3 591.8 540.4 364.6 547.1 534.5 503.6 443.3 464.6 1.5%

135 164 174 100 182 150 71 19 86 73 40 8 22 88 141 229 88 83 7 163 101

477.1 801.9 641.4 584.1 535.5 781.1 591.4 738.6 508.0 580.3 – 512.3 591.8 540.4 364.6 547.1 534.5 446.2 443.3 464.6 0.8%

286 318 173 153 685 316 113 34 117 109 70 107 127 181 479 630 124 31 7 291 217

477.1 786.7 630.5 584.1 525.4 781.1 591.4 738.6 508.0 580.3 – 512.3 591.8 540.4 364.6 547.1 534.5 446.2 443.3 434.8 0.2%

346 863 386 268 1082 318 100 45 146 127 88 115 147 212 1045 820 150 41 31 235 328

Note: Bold face indicates optimal solutions. The notations ‘–’ indicates infeasible solution.

160

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166 900

Deviation from optimal (%)

Time (second)

800 700 600 500 400 300 200 100 0 1

2

3

4

5

6

7

(λ) 8

9

10

11

14 12 10 8 6 4 2 0 1

2

Numer of neighbors

3

4

5

6

7 (λ) 8

9

10

11

Numer of neighbors

Fig. 5. Computational time and the deviations (%) from optimality with respect to various k values.

5. A partial-MIP method with iterated neighborhood search (P-MIP-INS) In this section, we present a hybrid algorithm using partial-MIP and iterated neighborhood search (INS) to solve largesized HGVRSP instances. We use term ‘‘partial-MIP” to indicate that the problem is solved only for a sub-set of the binary variables during the search process. The proposed algorithm will be referred to as P-MIP-INS hereinafter. The basic idea of the P-MIP-INS is to fix a set of decision variables and to use partial-MIP to optimize the model for the remaining variables. Note that once the vehicle types and routes are known, the scheduling decision variables can be determined rather efficiently. We use this observation in partial-MIP optimization to generate new solutions from an incumbent solution based on yijh variables. Let S represent a solution to the HGVRSP and yijh (S) be the value of binary variable yijh for solution S. Given solution S, MIP is used to find a new neighbor S 0 as follows: (1) select a sub-set F of yijh decision variables; (2) fix yijh

yijh (S)

for all (i; j) and h; (3) unfix yijh for all yijh 2 F; and (4) Call the MIP solver to optimize the problem and return S 0 . It is important to select yijh decision variables in set F in a logical way that the resulting new solution after solving partialMIP can be a different and likely a better one than the incumbent solution. Therefore, set F is selected using three operators such that it can be possible to form new routes and/or to change vehicle types in the incumbent solution. As described in the following section, these operators do not make perturbations on the incumbent solution, but indicate which part of the incumbent solution can be improved by partial-MIP. We use Variable Neighborhood Search (VNS) to alternate the operators so that different solution structures can be investigated during the research. VNS, first proposed by Mladenovic and Hansen (1997), is a top-level meta-heuristic that guides local search to explore globally optimized solutions through a systematical change of searching in predefined neighborhoods. Many VNS-based approaches can be found in the literature for various VRP related optimization problems (Hansen et al., 2008; Mladenovic et al., 2012; Labadie et al., 2012; Xiao et al., 2014). One of the contribution of the P-MIP-INS is the integration of VNS and MIP in order to solve large size problem instances. 5.1. Description of the neighborhood search operators In the P-MIP-INS, we use three operators to select set F for partial-MIP optimization, and these operators also define the corresponding neighborhood structures for implementing INS. The first two operators are based on a string representation of solutions. The routes of the selected vehicles is represented as string S = {0, s11 , . . . , 0, . . . , 0, sh1 , . . . , 0, . . . , 0, sq1 , . . . , 0} where 0 represents the depot, shp is the index of the pth customer node visited by vehicle h, and q is the number of available vehicles. For example, S = {0, 3, 1, 4, 0, 0, 2, 5, 0} with three vehicles (q = 3) is decoded as follows: the first vehicle starts from the depot, visits customers 3, 1, and 4 consecutively, and returns back to the depot; the second vehicle is not used; and the third vehicle starts from the depot and visits customers 2 and 5, and returns back to the depot. The first operator selects a random set of W consecutive nodes on string S as shown in Fig. 6A. Therefore, it is called singlerange-selection (SRS). Similarly, the second operator, which is called two-range-selection (TRS), selects two disjoints sets of W/2 consecutive nodes on string S as shown in Fig. 6 (B). The third operator, which is called nearest-neighbor-selection (NNS), selects a set of W nodes within a randomly placed circle on the network map as shown in Fig. 6 (C). In the case where customer coordinates are not available (i.e., only the distances among customers are provided in the problem data), a random node and its W  1 closest neighbors are selected. To find a new solution from the current solution, all binary decision variables related with the selected nodes are unfixed and optimized by the MIP solver while the other binary variables are fixed to their respective values in the incumbent solution. After a solution is found and accepted as the new incumbent, new string S is constructed for this new incumbent. These three operators are used in a framework of INS which is introduced in the following section. We also compare the effectiveness of these three operators through computational experiments. We control the size of set F by parameter W. If W = 8, for example, the SRS operator selects an interval with 8 nodes on string S, the DRS operator selects two non-overlapped intervals with 4 nodes in each, and the NNS operator randomly selects a random customer node and its 7 closest neighbors. Generally, a larger value of parameter W leads to a steeper improve-

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

161

Fig. 6. Three operators for conducting the partial-MIP optimization.

ment over the incumbent solution, but also significantly increases the computational time. Based on our preliminary computational experiments, parameter W should be set to a value smaller than 8 independently from the problem size. Given set F and incumbent solution S, finding a new solution involves unfixing decision variables yijh among all pairs of the nodes in set F and between the nodes in set F and the other nodes including the depot for all vehicles used in set F. For example, the SRS operator in Fig. 6 (A) will unfix decision variables y611 ; y612 ; y601 ; y602 ; y641 ; y642 ; y691 ; y692 ; y161 ; y162 ; y101 ; y 102; y141 ; y142 ; y191 ; y192 ; y061 ; y062 ; y011 ; y012 ; y041 ; y042 ; y091 ; y092 ; y461 ; y462 ; y411 ; y412 ; y401 ; y 402; y491 ; y492 ; y961 ; y962 ; y911 ; y912 ; y901 ; y902 ; y941 ; y942 ; y261 ; y211 ; y201 ; y241 ; y291 ; y632 ; y132 ; y032 ; y432 , and y932 . This enables the MIP solver to freely swap or relocate nodes in the selected range in order to achieve a local optimum. 5.2. Iterated neighborhood search (INS) The pseudo code of the P-MIP-INS algorithm is given in Fig. 7. Based on the above three operators (O = 1, 2, and 3), we define three neighborhoods of the incumbent solution. The P-MIP-INS algorithm starts from an initial solution and systematically searches for better solutions in these three pre-defined neighborhoods in turn. If a solution better than the incumbent solution is found, then it is accepted as the new incumbent solution, which resets the counter for non-improvements (N 0); otherwise, the search is moved to the next neighborhood (i.e., O ð ðO mod 3Þ þ 1 Þ. The stopping condition can be either a maximum time limit (T max ) or a maximum number of iterations (N max ) without improvement on the best solution so far (Sbest ). In Fig. 7, step 3A acts as the shaking phase, and step 3D is the local search using the MIP solver. Parameter T max and N max can be used to control the stopping condition, and parameter W is for controlling the number of selected nodes for partial MIP-optimization. 5.3. Greedy construction heuristics Constructing a feasible initial solution for the HGVRSP is quite difficult due to the several hard constraints such as customer and vehicle time windows, vehicle capacity limits, and vehicle range limits. The time-varying traffic conditions make it even more difficult to design a problem-specific construction heuristic that can guarantee feasible solutions. In the literature, several construction heuristics are proposed for the VRPTW such as the sweep algorithm (Christofides et al., 1979), the nearest-neighbor heuristic (Solomon, 1987), the savings heuristic (Clarke and Wright, 1964; Van Landeghem, 1988), and the insertion heuristic (Solomon, 1987; Ioannou et al., 2001). Olli and Michel (2005a) report that the insertion heuristic improved by Ioannou et al. (2001) based on the generic insertion framework of Solomon (1987) provide the best results for the VRPTW. Using CO2 emissions as the evaluation criterion, we can also utilize these construction heuristics to find initial solutions for the P-MIP-INS algorithm as described below. (1) The Insertion Heuristic: The insertion heuristic first initializes a route with a ‘seed’ customer, and then repeatedly uses some criteria to select the best customer to be inserted into the best location in the current partial route. If no more customers can be inserted into the current route, a new route is started with another ‘seed’ customer. In this paper, we modify this heuristic by considering multiple partial routes. The criterion for evaluating the best customer and the best location among multiple partial routes is the minimum increase in CO2 emissions after an insertion. Our insertion heuristic starts with an empty string S and then repeatedly inserts the best customer into the string at the best location until all customers are inserted. To determine the best customer and the best location, we use the HGVRSP model defined in this paper. In the case of very constrained problems, the algorithm may not be able to insert several last

162

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

Fig. 7. Framework of iterative neighborhood search.

customers into the existing routes due to constraint violations. In this case, we add additional virtual vehicles to the fleet so that all customers can be served. These virtual vehicles are penalized by setting their CO2 emissions rates very high. If the final solution includes any virtual vehicles, then it is considered as infeasible. (2) The Savings Heuristic: The savings heuristic is a well-known construction method for the VRP (Olli and Michel, 2005a). It starts with a solution in which each customer is served by a separate vehicle and then combines two routes {0, . . . , i, 0} and {0, j, . . . , 0} into a single route {0, . . . , i; j, . . . , 0} that yields the maximum distance saving of di0 þ d0j  dij among all route pairs. In this paper, we redefine the cost saving as the amount of CO2 emissions reduction resulted from combining of any two routes. Similarly, virtual vehicles with higher CO2 emissions rates are added to the fleet at the beginning of the savings heuristic, and the final solution is considered as infeasible if a virtual vehicle is used. (3) The Nearest-neighbor Heuristic: The nearest-neighbor heuristic starts a route by finding the unserved customer that is the closest to the depot. In the subsequent iterations, the heuristic searches for the closest customer to the last customer added and adds it to the route. If a customer cannot be added to the current route due to a capacity constraint, a new route is started. Similarly, we redefine the measurement of the ‘‘closest” in terms of CO2 emissions, i.e., adding the ‘‘closest” customer to the end of the route should result in the least increase in CO2 emissions. 6. Computational experiments Since the HGVRSP model defined in this paper has not been studied in the literature previously, we cannot directly compare the performance of the P-MIP-INS algorithm using results from the literature. Therefore, the P-MIP-INS algorithm is compared to other methods introduced in this paper using the PRP data of Demir et al. (2012), which is also used in Section 3.1. The test problem set includes seven groups of problems (10  5, 15  5, 20  5, 25  5, 50  5, 75  5 and 100  5), and each problem group has 20 problem instances. The solutions found by P-MIP-INS algorithm are compared to the solutions found by five other methods: (1) MIP/CPLEX, (2) Restricted-MIP/CPLEX, (3) the nearest-neighbor heuristic, (4) the savings heuristic, and (5) the insertion heuristic. The MIP/CPLEX method is using CPLEX to solve the MIP model. The Restricted-MIP/CPLEX method is described in Section 4.4, and it was used to solve the problem groups 10  5, 15  5, 20  5, and 25  5 with k = 5 and for the CPU time limit of 7200 s. All methods were implemented in AMPL, and the resulting subproblems of the three construction heuristics were solved by CPLEX (version 12.4.0.1). In the construction heuristics, subproblems have the same MIP formulation given in Section 3.2, except that decision variables yijh are all fixed as parameters whose values are pre-determined by the construction heuristics before calling the solver. Particularly, the P-MIP-INS algorithm used the insertion heuristic to find an initial solution. For the problem groups of 10  5, 15  5, 20  5, and 25  5, the P-MIP-INS algorithm was run for five times and obtained five solutions for each instance. The average and best objective function values are presented in Table 4 for comparisons. For the other large-sized instances, including 50  5, 75  5 and 100  5, each instance was solved only once. We used parameters N max ¼ n and T max = 7200 s as the stopping conditions.

163

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166 Table 4 Comparison of results for 140 problem instances. Feasible/instances

Vehicles usedb

Avg. obj.

Avg. D. (%)

Best/opt. rate (%)

Avg. time

19/20 19/20 19/20 19/20 19/20 19/20 19/20 19/20 5/20 7/20 8/20

– – – – 2/2.0/2 – 2/2.0/2 2/2.0/2 2/2.0/2 2/2.0/2 2/2.0/2

605.7 562.5 561.2 557.9 558.6 557.9 557.9 565.3 635.9 637.5 767.1

7.2 0.8 0.6 0.0 0.1 0 0 1.5 30.8 22.8 61.8

16 89 68 100 79 100 100 63 0 0 0

7s

20/20 20/20 20/20 20/20 20/20 20/20 17/20 19/20 12/20 10/20 10/20

– – – – 3/3.0/2 – 3/2.9/2 3/2.9/2 3/2.9/2 3/3.0/3 3/2.9/2

949.6 834.7 827.0 753.0 770.9 746.1 733.6 735.2 855.3 719.9 790.2

21.7 9.0 8.6 1.1 3.1 0.3 4.7 1.6 39.3 22.4 38.1

0 30 10 85 25 95 15 40 0 0 0

20/20 20/20 20/20 20/20 20/20 20/20 17/20 9/20 9/20 8/20

– – – – 4/3.9/3 – 4/3.9/3 4/3.9/3 4/4.0/4 4/3.9/3

1072.3 981.5 985.2 962.5 971.2 958.7 1099.5 1173.8 1119.0 1084.2

10.4 2.2 2.7 0.4 1.3 0.0 17.6 45.2 32.6 35.2

0 25 5 70 5 100 0 0 0 0

R-MIP/CPLEX (k = 5) Nearest-Neighbor Savings heuristic Insertion heuristic

20/20 20/20 20/20 20/20 20/20 20/20 17/20 17/20 16/20 16/20

– – – – 5/4.7/3 – 5/4.8/4 5/4.4/3 5/4.9/4 5/4.8/3

1008.0 962.7 971.7 944.9 952.3 936.0 1093.9 1355.2 1066.3 1187.0

7.8 2.9 3.8 0.9 1.7 0.0 21.7 56.8 21.9 38.0

0 5 0 35 0 100 0 0 0 0

120 m⁄ <1 m 1m <1 m

50  5

P-MIP-INS (W = 4) P-MIP-INS (W = 6) P-MIP-INS (W = 8) Nearest-Neighbor Savings heuristic Insertion heuristic

20/20 20/20 20/20 20/20 19/20 18/20

– – 10/8.6/8 10/8.6/7 10/9.6/9 10/8.9/8

1832.3 1741.1 1710.2 2584.6 2080.3 2321.2

7.3 1.9 0.0 52.0 21.8 37.2

0 35 100 0 0 0

4m 10 m 59 m 1m 2.8 m 3m

75  5

P-MIP-INS (W = 4) P-MIP-INS (W = 6) P-MIP-INS (W = 8) Nearest-Neighbor Savings heuristic Insertion heuristic

20/20 20/20 20/20 20/20 20/20 20/20

– – 15/12.1/11 15/13.0/11 15/14.3/13 15/12.7/11

2613.4 2520.9 2483.5 3806.5 2977.4 3308.9

5.4 1.6 0.0 53.7 20.1 33.9

0 0 100 0 0 0

12 m 21 m 89 m 4m 7m 10 m

100  5

P-MIP-INS (W = 4) P-MIP-INS (W = 6) P-MIP-INS (W = 8) Nearest-Neighbor Savings heuristic Insertion heuristic

20/20 20/20 20/20 20/20 18/20 20/20

– – 20/15.5/13 20/15.8/12 20/18.1/15 20/15.9/14

3337.0 3190.4 3149.8 4794.3 3766.3 4311.8

6.0 1.3 0.0 52.6 20.1 37.3

0 10 95 0 0 0

35 m 50 m 98 m 8m 13 m 31 m

Problem group

Methods

10  5

P-MIP-INS (W = 4) P-MIP-INS (W = 6) P-MIP-INS (W = 8)

Avg. of 5 runs Best. of 5 runs Avg. of 5 runs Best. of 5 runs Avg. of 5 runs Best. of 5 runs

MIP/CPLEX R-MIP/CPLEX (k = 5) Nearest-Neighbor Savings heuristic Insertion heuristic 15  5

P-MIP-INS (W = 4) P-MIP-INS (W = 6) P-MIP-INS (W = 8)

Avg. of 5 runs Best. of 5 runs Avg. of 5 runs Best. of 5 runs Avg. of 5 runs Best. of 5 runs

MIP/CPLEX R-MIP/CPLEX (k = 5) Nearest-Neighbor Savings heuristic Insertion heuristic 20  5

P-MIP-INS (W = 4) P-MIP-INS (W = 6) P-MIP-INS (W = 8)

Avg. of 5 runs Best. of 5 runs Avg. of 5 runs Best. of 5 runs Avg. of 5 runs Best. of 5 runs

R-MIP/CPLEX (k = 5) Nearest-Neighbor Savings heuristic Insertion heuristic 25  5

P-MIP-INS (W = 4) P-MIP-INS (W = 6) P-MIP-INS (W = 8)

Avg. of 5 runs Best. of 5 runs Avg. of 5 runs Best. of 5 runs Avg. of 5 runs Best. of 5 runs

Notes: bold face indicates best value. a Indicates limited computational time without optimality guarantee. b x/x/x, represents the number of available vehicles/the average number of vehicles used/the minimum number of vehicle used.

41s 464s 774s 3m <1 m <1 m <1 m 20s 86s 775s 120 ma 120 ma <1 m <1 m <1 m 32s 124s 1079s 120 ma <1 m <1 m <1 m 53s 239s 1864s

164

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

We also tested W ¼ f4; 6; 8g for all problem instances. The experiments were carried out on a Mac computer with 2.50 GHz IntelÒ CoreTM i5-2400S CPU. The results are summarized in Table 4. In the table, the column Feasible/Instances indicates the number of the instances for which the algorithm was able to find feasible solutions and the number of instances in the problem group. Column Vehicles used displays three numbers formatted as x/x/x, representing the number of available vehicles/average number of vehicles used/minimum number of vehicle used, respectively. Column Avg. obj. indicates the average amount of CO2 emissions (kg) found for the feasible solutions. Column Avg. D. (%) indicates the average percent difference from the best solution (or the optimal solution for group 10  5) found by all methods. Column Best/Opt. Rate (%) indicates the ratio of the number of the best/optimal solutions found by the method to the total number of feasible solutions. Column Avg. time indicates the average CPU time (wall clock time). It can be seen that P-MIP-INS (W = 8) provided the best solution quality in all problem groups. The solutions found by P-MIP-INS (W = 6) were slightly inferior to the ones by P-MIP-INS (W = 8), but they required significantly less computational time. P-MIP-INS (W = 4) was very fast, but performed poorly. From the computational cost/solution quality point of view, the method P-MIP-INS (W = 6) was the most balanced method. It found good quality solutions in relatively short CPU times. Both P-MIP-INS (W = 6) and P-MIP-INS(W = 8) were able to find 100% of optimal solutions for the 20 small instances of 10  5 group. The MIP/CPLEX method could also find optimal solutions for group 10  5, but it could not solve group 15  5 or larger ones within reasonable CPU times. The R-MIP/CPLEX (k = 5) method also performed well for small problems up to size of 25  5, but its performance deteriorated in the case of the large-sized problems. The construction heuristics were also efficient for small-sized problems. However, their computational times increased dramatically for the large-sized problems. Generally, the savings heuristic was the best construction heuristic, outperforming the other two with about 15% average difference from best solutions found. Next, we study the performances of the three search operators used in the P-MIP-INS algorithm. In the experiments, we recorded the number of improvements on the incumbent solution and the number of times that the operators are used during the search. Table 5 presents the average improvement ratios (the ratio of the number of times that an operator improved the incumbent solution and the number of times that the operator was used) of the three operators for different problem groups. Fig. 8 illustrates the overall trends of the ratios for all problem groups. It can be observed that the three operators contributed to the search process equally in the case of small-sized problems, i.e., group 10  5, 15  5, and 20  5. However, for large-sized problem groups, i.e., group 50  5, 75  5, and 100  5, the third operator, i.e., the NNS operator, outperformed the SRS and TRS operators. For example, the NNS operator had an improvement ratio of 56% for problem group 100  5.

Table 5 Comparison of the local search operators. Problem group

Feasible instances

Average improvements

Ratio % SRS (O = 1) (%)

TRS (O = 2) (%)

NSS (O = 3) (%)

10  5 15  5 20  5 25  5 50  5 75  5 100  5

20 19 19 20 20 20 20

6.6 12.3 14.9 20.7 43.2 63.2 91.2

31 28 26 24 25 23 20

34 33 35 29 27 25 24

35 39 39 46 48 52 56

Fig. 8. Improvement ratio of the three different operators.

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

165

7. Conclusions In this paper, we present a comprehensive mixed integer linear model for time-dependent, heterogeneous, capacitated vehicle routing and scheduling problem with time windows and time-varying traffic congestion in order to minimize CO2 emissions. The new proposed model has some new features as follows: 1. A general case for time-varying traffic congestion in multiple periods is considered. 2. The travel times of arcs depend on the actual periods that they are traveled, not on their travel starting time. This enables to schedule idle/waiting times both at customer premises and on the roadside (in long-haul roads) in order to avoid traffic congestion. 3. The non-stopping assumption accepted by most of the TD-VRP models in literature is relaxed. We have showed that the non-stopping assumption might reduce the solution quality 8.0% on the average when a CO2 emissions-based objective is used. Relaxing this assumption is important in the case of arcs that represent long-haul roads. 4. The fuel tank capacity is modeled as a constraint, which is practical for alternative fuel vehicles under time varying traffic conditions. 5. A fleet of vehicles with general heterogeneities is considered. Vehicles can have different fuel consumption rates, fuel-toCO2 conversion rates (due to using different types of fuel), load capacities, fuel tank capacities, and even different time windows. We developed six approaches to solve the HGVRSP, which are MIP/CPLEX, Restricted-MIP/CPLEX, Insertion heuristic, Nearest-neighbor heuristic, Savings heuristic, and P-MIP-INS, and studied their performances using 140 benchmark problem instances with up to 100 nodes. The experimental results showed that the P-MIP-INS algorithm provided the best solutions for all groups of problem instances. The P-MIP-INS algorithm is a combination of MIP and a framework of neighborhood search, which is a promising approach for solving MIP optimization problems with difficult constraints. Studying the computational aspects of the proposed model, formulating the problem in a multi-objective framework, and considering a streetlevel model are interesting further research venues. Acknowledgement This work is partly supported by the National Natural Science Foundation of China under Grants No. 71271009. References Ahn, B.-H., Shin, J.-Y., 1991. Vehicle-routing with time windows and time-varying congestion. J. Oper. Res. Soc. 42 (5), 393–400. Barth, M.J., Boriboonsomsin, K., 2008. Real-world CO2 impacts of traffic congestion. Transp. Res. Rec.: J. Transp. Res. Board 2058, 163–171. Barth, M.J., Younglove, T., Scora, G., 2005. Development of a Heavy-Duty Diesel Modal Emissions and Fuel Consumption Model. Technical Report. UCBITSPRR-2005-1, California PATH Program, Institute of Transportation Studies, University of California at Berkeley. Baldacci, R., Battarra, M., Vigo, D., 2008. Routing a heterogeneous fleet of vehicles. In: The Vehicle Routing Problem: Latest Advances and New Challenges. Springer, pp. 3–27. Bektasß, T., Laporte, G., 2011. The pollution-routing problem. Transp. Res. Part B 45, 1232–1250. Bektasß, T., Demir, E., Laporte, G., 2016. Green vehicle routing. In: Green Transportation Logistics. Springer International Publishing, pp. 243–265. Bodin, L., Golden, B., 1981. Classification in vehicle routing and scheduling. Networks 11 (2), 97–108. Chai, J., 2015. Under The Dome: Report of PM2.5 air pollution in China. (access September 9, 2015). Chen, H.-K., Hsueh, C.-F., Chang, M.-S., 2006. The real-time time-dependent vehicle routing problem. Transp. Res. Part E 42, 383–408. Clarke, G., Wright, J.W., 1964. Scheduling of vehicles from a central depot to a number of delivery points. Oper. Res. 12, 568–581. Christofides, N., Mingozzi, A., Toth, P., 1979. The vehicle routing problem. In: Christofides, N., Mingozzi, A., Toth, P., Sandi, C. (Eds.), Combinatorial Optimization. Wiley, Chichester, pp. 315–338. Coyle, M., 2007. Effects of payload on the fuel consumption of trucks. UK Department of Transport Report. Dantzig, G.B., Ramser, J.H., 1959. The truck dispatching problem. Manage. Sci. 6 (1), 80–91. Demir, E., Bektasß, T., Laporte, G., 2012. An adaptive large neighborhood search heuristic for the pollution-routing Problem. Eur. J. Oper. Res. 223, 346–359. Demir, E., Bektasß, T., Laporte, G., 2014a. A review of recent research on green road freight transportation. Eur. J. Oper. Res. 237 (3), 775–793. Demir, E., Bektasß, T., Laporte, G., 2014b. The bi-objective pollution-routing problem. Eur. J. Oper. Res. 232, 464–478. Desaulniers, G., Madsen, O.B., Røpke, S., 2014. The vehicle routing problem with time windows. In: Toth, P., Vigo, D. (Eds.), Vehicle Routing: Problems, Methods, and Applications, . second ed., MOS-SIAM Series on Optimization second ed. SIAM, Philadelphia, pp. 119–159. Eksioglu, B., Vural, A.V., Reisman, A., 2009. The vehicle routing problem: a taxonomic review. Comput. Ind. Eng. 57 (4), 1472–1483. Eglese, R., Madena, W., Slaterb, A., 2006. A Road TimetableTM to aid vehicle routing and scheduling. Comput. Oper. Res. 33, 3508–3519. Eglese, R., Bektasß, T., 2014. Chapter 15: green vehicle routing. In: Toth, P., Vigo, D. (Eds.), Vehicle Routing: Problems, Methods, and Applications, . second ed., MOS-SIAM Series on Optimization second ed. SIAM, Philadelphia, pp. 437–458. Figliozzi, M.A., 2010. Vehicle routing problem for emissions minimization. Transp. Res. Rec.: J. Transp. Res. Board 2197, 1–7. Figliozzi, M.A., 2012. The time dependent vehicle routing problem with time windows: benchmark problems, an efficient solution algorithm, and solution characteristics. Transp. Res. Part E 48, 616–636. Fleischmann, B., Gietz, M., Gnutzmann, S., 2004. Time-varying travel times in vehicle routing. Transp. Sci. 38, 160–173. Franceschetti, A., Honhon, D., Woensel, T.V., Bektasß, T., Laporte, G., 2013. The time-dependent pollution-routing problem. Transp. Res. Part B 56, 265–293. Garcia, B.L., Potvin, J.Y., Rousseau, J.M., 1994. A parallel implementation of the tabu search heuristic for vehicle routing problems with time window constraints. Comput. Oper. Res. 21, 1025–1033. Gaur, D.R., Mudgal, A., Singh, R.R., 2013. Routing vehicles to minimize fuel consumption. Oper. Res. Lett. 41, 576–580. Golden, B.L., Assad, A.A., Levy, L., Gheysens, F., 1984. The fleet size and mix vehicle routing problem. Comput. Oper. Res. 11, 49–66. Golden, B.L., Raghavan, S., Wasil, E., 2008. The Vehicle Routing Problem: Latest Advances and New Challenges. Springer, Berlin. Hansen, P., Mladenovic, N., Pérez, J.A.M., 2008. Variable neighborhood search. Eur. J. Oper. Res. 191, 593–595.

166

Y. Xiao, A. Konak / Transportation Research Part E 88 (2016) 146–166

Hickman, A.J., 1999. Methodology for calculating transport emissions and energy consumption, Transportation Research Laboratory Report, No: SE/491/98. Hill, A.V., Benton, W.C., 1992. Modeling intra-city timedependent travel speeds for vehicle scheduling problems. J. Oper. Res. Soc. 43 (4), 343–351. Ichoua, S., Gendreau, M., Potvin, J.-Y., 2003. Vehicle dispatching with time-dependent travel times. Eur. J. Oper. Res. 144, 379–396. International Road Transport Union, 2012. CO emission. (accessed 16.02.14). Ioannou, G., Kritikos, M., Prastacos, G., 2001. A greedy look-ahead heuristic for the vehicle routing problem with time windows. J. Oper. Res. Soc. 52, 523– 537. Jabali, O., Van Woensel, T., de Kok, A.G., 2012. Analysis of travel times and CO2 emissions in time-dependent vehicle routing. Prod. Oper. Manage. 21 (6), 1060–1074. Kara, I., Kara, B.Y., Yetis, M.K., 2007. Energy minimizing vehicle routing problem. In: Dress, A., Xu, Y., Zhu, B. (Eds.), Combinatorial Optimization and Applications: FIRST International Conference, LNCS 4616. Springer, Berlin/Heidelberg, pp. 62–71. Koç, Ç., Bektasß, T., Jabali, O., Laporte, G., 2014. The fleet size and mix pollution-routing problem. Transp. Res. Part B: Methodol. 70, 239–254. Koç, Ç., Bektasß, T., Jabali, O., Laporte, G., 2015. A hybrid evolutionary algorithm for heterogeneous fleet vehicle routing problems with time windows. Comput. Oper. Res. 64, 11–27. Koç, Ç., Bektasß, T., Jabali, O., Laporte, G., 2015. Thirty years of heterogeneous vehicle routing. Eur. J. Oper. Res. (in press), http://dx.doi.org/10.1016/j.ejor. 2015.07.020. Kramer, R., Subramanian, A., Vidal, T., Lucídio dos Anjos, F.C., 2015. A matheuristic approach for the pollution-routing problem. Eur. J. Oper. Res. 243 (2), 523–539. Kuo, Y., 2010. Using simulated annealing to minimize fuel consumption for the time-dependent vehicle routing problem. Comput. Ind. Eng. 59, 157–165. Kwon, Y.J., Choi, Y.J., Lee, D.H., 2013. Heterogeneous fixed fleet vehicle routing considering carbon emission. Transp. Res. Part D 23, 81–89. Labadie, N., Mansini, R., Melechovsky´, J., Calvo, R.W., 2012. The team orienteering problem with time windows: an LP-based granular variable neighborhood search. Eur. J. Oper. Res. 220 (1), 15–27. Laporte, G., 1992. The vehicle routing problem: an overview of exact and approximate algorithms. Eur. J. Oper. Res. 59, 345–358. Laporte, G., Gendreau, M., Potvin, J.Y., Semet, F., 2000. Classical and modern heuristics for the vehicle routing problem. Int. Trans. Oper. Res. 7 (4–5), 285– 300. Lin, C., Choy, K.L., Ho, G.T., Chung, S.H., Lam, H.Y., 2014. Survey of green vehicle routing problem: past and future trends. Expert Syst. Appl. 41 (4), 1118– 1138. Maden, W., Eglese, R., Black, D., 2010. Vehicle routing and scheduling with time-varying data: a case study. J. Oper. Res. Soc. 61, 515–522. Malandraki, C., 1989. Time Dependent Vehicle Routing Problems: Formulations, Solution Algorithms and Computational Experiments. Ph.D. Dissertation, Northwestern University, Evanston, Illinois. Malandraki, C., Daskin, M.S., 1992. Time-dependent vehicle-routing problems – formulations, properties and heuristic algorithms. Transp. Sci. 26 (3), 185– 200. Malandraki, C., Dial, R.B., 1996. A restricted dynamic programming heuristic algorithm for the time dependent traveling salesman problem. Eur. J. Oper. Res. 90 (1), 45–55. Mladenovic, N., Hansen, P., 1997. Variable neighborhood search. Comput. Oper. Res. 24, 1097–1100. Mladenovic, N., Urosevic, D., Hanafi, S., Ilic, A., 2012. A general variable neighborhood search for the one-commodity pickup-and-delivery travelling salesman problem. Eur. J. Oper. Res. 220, 270–285. Olli, B., Michel, G., 2005a. Vehicle routing problem with time windows, part I: route construction and local search algorithms. Transp. Sci. 39 (1), 104–118. Olli, B., Michel, G., 2005b. Vehicle routing problem with time windows, part II: metaheuristics. Transp. Sci. 39 (1), 119–139. Rodrigue, J.P., Comtois, C., Slack, B., 2013. The geography of transport systems. Routledge, London. Solomon, M., 1983. Vehicle Routing and Scheduling with Time Window Constraints: Models and Algorithms. Ph.D. Dissertation, Dept. of Decision Sciences, University of Pennsylvania. Solomon, M., 1987. Algorithm for the vehicle routing and scheduling problems with time window constraints. Oper. Res. 35 (2), 254–265. Soysal, M., Bloemhof-Ruwaard, J.M., Bektasß, T., 2015. The time-dependent two-echelon capacitated vehicle routing problem with environmental considerations. Int. J. Prod. Econ. 164, 366–378. Thangiah, S.R., 1995. Vehicle routing with time windows using genetic algorithms. In: Chambers, L. (Ed.), Application Handbook of Genetic Algorithms: New Frontiers, vol. II. CRC Press, BocaRaton, FL, pp. 253–277. Toth, P., Vigo, D., 2002. An overview of vehicle routing problems. In: Toth, P., Vigo, D. (Eds.), The Vehicle Routing Problem, SIAM Monographs on Discrete Mathematics and Applications, vol. 9. SIAM, pp. 1–26. US EPA, 2014. Inventory of U.S. Greenhouse Gas Emissions and Sinks: 1990–2012. USEPA #EPA 430-R- 14-003. . Toth, P., Vigo, D., 2014. Vehicle Routing Problems, Methods, and Application, MOS-SIAM Series on Optimization, second ed. SIAM, Philadelphia, . Van Landeghem, H., 1988. A Bi-criteria heuristic for the vehicle-routing problem with time windows. Eur. J. Oper. Res. 36 (2), 217–226. Vidal, T., Crainic, T.G., Gendreau, M., Prins, C., 2013. A hybrid genetic algorithm with adaptive diversity management for a large class of vehicle routing problems with time-windows. Comput. Oper. Res. 40, 475–489. Xiao, Y., Zhao, Q., Kaku, I., Xu, Y., 2012. Development of a fuel consumption optimization model for the capacitated vehicle routing problem. Comput. Oper. Res. 39 (7), 1419–1431. Xiao, Y., Zhang, R., Zhao, Q., Kaku, I., Xu, Y., 2014. A variable neighborhood search with an effective local search for uncapacitated multilevel lot-sizing problems. Eur. J. Oper. Res. 235 (1), 102–114. Xiao, Y., Konak, A., 2015a. A simulating annealing algorithm to solve the green vehicle routing & scheduling problem with hierarchical objectives and weighted tardiness. Appl. Soft Comput. 34, 372–388. Xiao, Y., Konak, A., (2015b). Green vehicle routing problem with time-varying traffic congestion. In: Proceedings of 14th INFORMS Computing Society Conference, Richmond, Virginia, January 11–13, 2015, pp. 134–148.