University of Southern California Dissertations and Theses

Routing problems for fuel efficient vehicles

(USC Thesis Other)

Routing problems for fuel efficient vehicles

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Routing Problems for Fuel Efficient Vehicles
By
Yihuan (Ethan) Shao
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(INDUSTRIAL SYSTEMS AND ENGINEERING)
December 2017
Committee:
Dr. Maged M. Dessouky, Chair
Dr. Fernando Ordonez
Dr. Ketan Savla
Contents
1 Introduction 4
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Research gap and contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Organization of the dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Literature Review 10
2.1 Green Vehicle Routing Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Location Routing Problem for Fuel Stations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Pickup and Delivery Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 CNG Truck Routing Problem with Fuel Stations 18
3.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Mixed Integer Programming Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 Heuristic Solution Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4 CNG Truck Location Routing Problem 45
4.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2 Mixed Integer Programming Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3 Heuristic Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.4 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5 Hybrid Transit System 61
5.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.2 Mixed Integer Programming Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3 Heuristic Solution Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.4 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6 Conclusion and Future Directions 76
References 78
1
List of Figures
1 An example of a truck’s multiple visits to a fuel station . . . . . . . . . . . . . . . . . . . . . 24
2 An example of route k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3 Steps of constructing an initial feasible solution . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4 Steps of the Adaptive Large Neighborhood Search (ALNS) . . . . . . . . . . . . . . . . . . . . 30
5 Local search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
6 Average total working time with different number of CNG fueling stations . . . . . . . . . . . 40
7 Average total working time with difference tank capacities . . . . . . . . . . . . . . . . . . . . 41
8 Average total working time with different tank capacities for light duty trucks . . . . . . . . . 43
9 Average total working time with different sizes of light duty truck fleets . . . . . . . . . . . . 44
10 Steps of the heuristic method for the CTLRP . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
11 Convergence of the ALNS method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
12 Locations of fuel stations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
13 Sensitivity analysis of the CTLRP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
14 Bus map for small instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
15 Bus map for large instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
16 Average results with different ratios of
1
=
2
. . . . . . . . . . . . . . . . . . . . . . . . . . . 74
17 Average results with different headway . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
2
List of Tables
1 Results on small CTRPFS instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2 Results comparison between different solution methods . . . . . . . . . . . . . . . . . . . . . . 38
3 Results for different number of CNG fueling stations . . . . . . . . . . . . . . . . . . . . . . . 40
4 Results for different fuel tank capacities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5 Results on small CTLRP instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6 Results comparison between different solution methods . . . . . . . . . . . . . . . . . . . . . . 57
7 Comparison of different branch-and-cut algorithms . . . . . . . . . . . . . . . . . . . . . . . . 58
8 Notations for the exact formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
9 Results comparison between different solution methods . . . . . . . . . . . . . . . . . . . . . . 72
10 Results comparison with different ratios of
1
=
2
. . . . . . . . . . . . . . . . . . . . . . . . . 73
11 Results comparison with different headway . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3
1 Introduction
In this chapter, we first introduce the background of the Green Vehicle Routing Problem (GVRP), which is
related to our research. Then we highlight the research gaps and our contributions. In the last part, we list
the organization of this dissertation.
1.1 Background
The climate of our planet keeps changing. We are facing the risks due to these changes, such as temperature
rising, rain pattern shifting. All these changes are believed to be due to the increasing level of Greenhouse
Gases (GHG) in the atmosphere, which is mainly caused by human activities. Within these activities,
transportation accounted for a large part in GHG emissions. In 2014, transportation contributed about
26% of the total U.S. GHG emissions, which makes it the second largest among the five major economic
sectors (Others are electricity, industry, commercial & residential, and agriculture) (U.S. Environmental
Protection Agency, 2016). Globally speaking, transportation emissions are also an important driver of
increasing GHG emissions (Intergovernmental Panel on Climate Change, 2015). It accounted for 14% of
global GHG emissions in 2010 (U.S. Environmental Protection Agency, 2017). All these factors call for
better planning on transportation. So far, a significant amount of effort has been made on transportation
to reduce its negative impact on the environment. In this dissertation, we mainly focus on efforts in two
directions, which are exploring alternative fuels and changing people’s behavior. The first direction focuses
on reducing emissions per traveled mile. It introduces clean energy technologies like biofuel, electric vehicles
etc. The second direction aims at reducing the number of vehicle miles traveled. It covers topics such as
ridesharing and bus scheduling.
For the direction of reducing emissions per traveled mile, various kinds of green energies like biodiesel,
electricity, ethanol, hydrogen, methanol, natural gas and propane have been used for alternative fuel ve-
hicles (Erdoğan and Miller-Hooks, 2012). Within these alternative fuel energies, Compressed Natural Gas
(CNG) is considered as a possible solution for fossil fuel substitution because of its wide availability, engine
compatibility, and low operations cost (Khan et al., 2015). In this dissertation, we build models based on
the background of CNG trucks. However, the models can also be used for other alternative fuel vehicles. In
2015, CNG vehicles accounted for 14% of the total Alternative Fuel Vehicles (AFVs) inventory in the U.S.
and contributed 55% of the petroleum savings due to AFVs (Johnson and Singer, 2016). More importantly,
life cycle analysis/assessment (LCA) of CNG vehicles shows they can reduce GHG emissions (Shahraeeni
4
et al., 2015; Tong et al., 2015). CNG vehicles accounted for 15% of the GHG tons reduced by AFVs in 2015
in the U.S. (Johnson and Singer, 2016) The first natural gas powered vehicle can be dated back to the 1930s.
Since then renewed attention is being given to CNG vehicles especially when the gasoline price rises (Khan
et al., 2015). Today there are more than 22 million CNG vehicles distributed in over 80 countries (Altfuels
Communications Group, 2017). With the increased emphasis on sustainability, some government systems
and companies have already purchased CNG vehicles to support their daily operations (Dessouky and Wang,
2009). Successful stories can be found and more CNG vehicles and fueling facilities are deploying now (U.S.
Department of Energy, 2017d).
CNG vehicles are expected to contribute significantly in reducing GHG emissions; however, they still face
several challenges on the way to replacing traditional diesel or gasoline fuel. One of the biggest challenges
is the sparse, uneven allocation of CNG fueling stations, which may lead to a long detour for refueling. It is
reported by the Los Angeles Clean Cities Coalition for the U.S. Department of Energy, there are 98 natural
gas fueling stations in the City of Los Angeles, which is not even comparable to electric fueling stations
who have 3,625 spots around the city (U.S. Department of Energy, 2017c). In fact, the detour to a CNG
fueling station could increase overall energy consumption and finally offset the low emission benefit from
the CNG fuel. The long waiting time in CNG fueling stations is also another important problem. The
waiting time may due to the long queues for the pumps (Ahmed et al., 2013; Khan et al., 2015) or the refill
time for the CNG storage vessels in the stations (Dessouky and Wang, 2009; U.S. Department of Energy,
2017b). Additional challenges exist like limited tank capacity (Erdoğan and Miller-Hooks, 2012; Salimifard
and Raeesi, 2014). For instance, the Freightliner M2 tractor is estimated to have a range of 275 to 325 miles
per fill (Laughlin and Burnham, 2016) and it is hard to extend further because of the large gas tank volume.
SomebiginstitutionsbuildtheirprivateCNGfuelingstations. However, manyofthelocalfreightdelivery
operations are performed by independent truckers who are too small to manage their own fuel station and
the public fuel stations are the only option for them. In Southern California, even for drivers who have
CNG fueling stations in the base depots, they are also reliant on public stations for half of their refueling
needs (Kelley and Kuby, 2015). Empirical research shows that drivers do consider station locations and
convenient location is the most important reason for choosing a fuel station (Kelley and Kuby, 2015). Thus,
it becomes imperative to develop efficient deployment strategies for refueling related decisions that will have
minimal impact on their routing operations. First, we introduce the CNG Truck Routing Problem with
Fueling Stations (CTRPFS) to model decisions to be made with regards to the vehicle routes including
the choice of fuel stations. Second, we propose a CNG Truck Location and Routing Problem (CTLRP) to
5
simultaneously solve the location problem for CNG fueling stations together with the routing problem for
CNG trucks. The CTRPFS is a special case of the well-known Vehicle Routing Problem (VRP). It combines
classical restrictions such as capacity, daily traveling distance as well as new constraints for CNG trucks,
such as fuel tank capacity, and waiting time for refueling. Compared with other VRPs that include fuel
stations like the routing problem for electric vehicles, CNG vehicles have a special cost due to the fixed
fueling time, which occurs every time a truck visits a fuel station. The fixed fueling time mainly comes
from the long waiting time at the fuel stations (Ahmed et al., 2013; Khan et al., 2015; U.S. Department of
Energy, 2017b), and is independent of the gas volume a truck refuels. This fixed fueling time calls for special
attention, because the detour is no longer the only factor considered in making fueling decisions. The model
needs to keep both the fueling frequency and detour at low levels. CTRPFS is designed for small fleets
of CNG trucks who have to use the public fuel stations for their daily operations, while CTLRP is a new
business model targeting a group of small fleets who work together to build some shared private CNG fueling
stations in their delivery regions, and independently deliver products to customers. The idea of cooperation
is not new. In a report published by the National Renewable Energy Laboratory for the U.S. Department
of Energy, a model called Vehicle Infrastructure and Cash-Flow Evaluation (VICE) is developed to evaluate
the financial soundness for fleets of CNG vehicles to build a CNG fueling infrastructure, and the location of
CNG fueling station is simplified to sit at the center of the region (Mitchell, 2015). Some organizations exist
to unite a diverse group of stakeholders. For instance, the Natural Gas Vehicle Technology Forum (NGVTF)
was established in 2008. One purpose of the forum is to enable fleets and large purchasers to aggregate
demand for natural gas vehicles and equipment (U.S. Department of Energy, 2017e). Though the CTLRP
and CTRPFS are designed for CNG trucks, they can be used to solve problems for other alternative fuel
vehicles, who also face the challenges from small tank capacity and sparsely allocated fuel stations.
For the other direction, changing people’s behavior, ridesharing is one of the hottest methods. Rideshar-
ing, or more specifically real-time ridesharing, has enjoyed a boost in recent years because of the development
of navigation and location sharing technologies on mobile devices (Furuhata et al., 2013). Many companies
have started businesses on ridesharing in the last five years and people even created a new word "transporta-
tion network companies" for this category of service. Reported by some of these companies, the ridesharing
companies are not only new occupants of the traditional taxi market, but also enlarge the market by at-
tracting more people at low traveling cost. In fact, ridesharing is changing people’s behavior into a more
environmentally friendly way by reducing the number of private cars needed for people’s daily travel. A
case study in Dublin confirms the contribution of ridesharing in reducing GHG emissions and the vehicle
6
kilometers traveled (Caulfield, 2009). Also, the reduction in private cars in turn can solve other problems like
traffic congestion, limited parking spaces, caused by traditional private cars (Stiglic et al., 2015). With higher
occupancy rates the energy efficiency of shared vehicles is usually higher than private cars, and compared
with vehicles in the public transit system, like buses or subways, the energy efficiency of shared vehicles is
still low. However, fixed-route vehicles are troubled with the famous last-mile problem, which lowers the
attraction of this transportation method. In order to combine these two systems together, we propose a
hybrid transit system in which the routes of shared vehicles are scheduled together with the public transit
system. Passengers are allowed to use a combination of shared vehicles and buses or subways so that the total
traveling cost for passengers is reduced and the fuel efficiency of the whole system can be further improved.
Moreover, the bus stops or subway stations can serve as meeting points for the ridesharing systems and may
result in additional feasible matches between drivers and riders (Stiglic et al., 2015).
1.2 Research gap and contribution
Though the Vehicle Routing Problem has been a hot topic for decades, the GVRP still calls for further
research. Some of the features of AFVs may have large impacts on the daily operations in reality and those
features may change the models fundamentally. For instance, in the classic VRP, the tank capacity is usually
not considered into the model due to the large mileage range and easy access to gas stations. However, for
most of the AFVs, this assumption no longer holds. The models have to include this factor. As to the pickup
and delivery problem, which is a classical topic in VRP, the hybrid transit system is a new background
setting. This hybrid scheduling on one hand, is meaningful to society. It should improve the fuel efficiency
of our transit system. Also based on current navigation and location sharing technologies on mobile devices,
the implementation of the coordination is not hard. On the other hand, the given solution methods, both
exact algorithms and metaheuristics, for the pickup and delivery problem are no longer suitable for the new
hybrid system. In the new system, the vertices representing bus stops may or may not be visited by shared
vehicles. This modification brings in a fundamental change in the solution method.
In this dissertation, we make contributions to the research about GVRP by proposing three models,
which are the CTRPFS, the CTLRP and the hybrid transit system.
WefirstintroducetheCTRPFStomodeldecisionstobemadewithregardstothevehicleroutesincluding
the choice of fuel stations. A fixed refueling cost, which may represent the waiting time at the fuel station,
is considered for the fueling procedure. This fixed refueling cost is unique for CNG trucks and few papers
add it into their models. Although only few researchers touch this topic, the motivation for fixed refueling
7
cost comes from the real world. Users may wait for hours in long queues to get their trucks refueled at
CNG fueling stations (Khan et al., 2015). Meanwhile the fixed refueling cost can have a large influence
on the solution method. A fuel station with a shorter detour may not lead to a better solution. We
need to balance the detour and the frequency of visiting fuel stations. In order to find exact solutions for
the problem, we provide a MIP model. Our contribution is to introduce a new decision variable setting
method with a variable elimination method and valid inequalities. Small instances are solved via CPLEX
and serve as benchmarks for the heuristic method. To deal with instances in realistic size, we propose a
new hybrid heuristic method which combines an ALNS with a local search and MIP model. Compared
with other methods, our unique embedded MIP model can solve a problem for a more general situation.
In computational experiments, we modify benchmark instances from Augerat et al. (1995) with constraints
from CNG trucks. The computational experiments explore some insights. Based on the estimation of daily
goods delivery data from the Ports of Los Angeles and Long Beach, we conduct experiments comparing the
total traveling distance using CNG trucks with the current use of diesel trucks for local delivery. The purpose
of these experiments is to show how much the losses of route efficiency for CNG trucks can be reduced with
an efficient deployment.
Based on the CTRPFS, we propose the CTLRP, in which multiple CNG truck fleets work together to
build private CNG fueling stations in an area to support their daily operations. The CTLRP includes two
problems which are the location problem for the CNG fueling stations and the routing problem for the trucks.
Considering the inter-dependence of the location and routing problem, we simultaneously solve these two
problems. We propose a MIP formulation with valid inequalities and an adaptive branch-and-cut algorithm
to solve the problem exactly, and a hybrid heuristic algorithm to solve large scale problems. Experiments
show that our valid inequalities and adaptive branch-and-cut algorithm are efficient and can improve the
solving speed of the MIP. The hybrid heuristic algorithm can find better solutions than well-known methods
such as the heuristic methods for the p-median problem.
The last model, our hybrid transit system, combines a ridesharing system with a fixed-route public
transit system. The hybrid transit system is designed to solve the routing problem for the shared vehicles
who receive delivery request from riders. The model includes time windows from customers and shared
vehicles, as well as the schedule from the public transit system. We propose a MIP with valid inequalities
and a branch-and-cut algorithm to solve the problem exactly. Since the problem is computationally hard
to solve, we also develop a heuristic method. The computational results show our valid inequalities and
the branch-and-cut algorithm can save CPU time. Large scale problems solved by the proposed heuristic
8
method show some insights of the hybrid transit system.
1.3 Organization of the dissertation
This dissertation is organized as follows. In Chapter 2 a review of the related literature is presented.
Chapter 3 describes the CTRPFS model including a general MIP formulation and a heuristic solution
method to solve the problem. Similarly, Chapter 4 gives the CTLRP model as well as a MIP formulation
and a heuristic solution method. The hybrid transit system is proposed in Chapter 5. In Chapter 6, we give
the conclusion and future directions.
9
2 Literature Review
In this chapter, we review papers for the three models, including CTRPFS, CTLRP and the hybrid transit
system, mentioned in the previous chapter. CTRPFS can be categorized as a special case of the Green
Vehicle Routing Problem (GVRP), which recently is a hot topic in the well-known VRP. Our review starts
from papers about GVRP with different factors considered in the routing problem, and focuses more on the
research which is similar to our CTRPFS. Our second model, CTLRP, is related to several kinds of research
streams. It is a VRP when considering the routing part of the model and it is also a refueling location
problem when we focus on the facility location part. In the second part of the literature review chapter,
papers about the Location Routing Problem (LRP) are presented. The hybrid transit system is in fact a
Pickup and Delivery Problem (PDP). Papers about PDP cover both heuristic and exact methods. We also
list papers whose backgrounds are quite similar to the hybrid transit system.
2.1 Green Vehicle Routing Problem
VRP is a central problem in transportation (Bektaş et al., 2016). It first appeared in a paper by Dantzig
and Ramser (1959). It solves a routing problem for trucks to meet the demands from customers located in
differentplaceswhileminimizingthetotalmileagetraveledbythefleet. Sincethen, VRPhasbeenahottopic
in Operations Research and extended to a number of variants with respect to different real-world constraints.
Basic variants of VRP includes Capacitated VRP (CVRP), VRP with Time Windows (VRPTW), Pickup
and Delivery Problems (PDP), etc. Nowadays, with increasing concerns on the sustainability of the logistics
system, current research on VRP has received more interest on environmental issues, often referred as the
Green Vehicle Routing Problem (GVRP). The standard objective function for traditional VRP is to minimize
the total traveling distance, while many of the GVRP papers consider the reduction of pollutants, such as
nitrogen oxides (N
2
O), particulate matter (PM), as well as greenhouse gases (GHG). Some papers directly
add those items into the objective function, while many more consider fuel consumption related terms (Demir
et al., 2014), since fuel consumption can sometimes be used as a surrogate measure for the emissions of air
pollutants.
GVRP covers a wide range of topics, and factors which influence fuel consumption are considered in
GVRP papers. Based on a survey paper, these factors can be divided into five groups, including vehicle,
environment, traffic, driver and operations related categories (Demir et al., 2014). Speed and congestion in
the traffic related category are common factors in many GVRP related papers (Qian and Eglese, 2016; Xiao
10
and Konak, 2016). Some studies (Bektaş and Laporte, 2011; Xiao et al., 2012) also seek to add a payload
factor, which belongs to the operations related category, into their models. Various kinds of emissions models
have been developed to integrate these factors together for fuel consumption or emissions calculations. For
instance, one of the most widely used emissions models is called MEET (Demir et al., 2014), which was first
proposed by Hickman et al. (1999). Xiao and Konak (2016) adopt a MEET model to calculate the emissions
associated with the arcs between cities based on factors such as speed and payload. Instead of traveling
distances, emissions are the new costs for the arcs and a comprehensive MIP formulation is proposed to
solve the GVRP. They use CPLEX to solve the MIP and also develop some heuristic methods to deal with
large cases. Bektaş and Laporte (2011) present another emissions model to incorporate factors like speed and
payload into the cost of emissions, and finally into the routing problem. A MIP formulation with strengthen
constraints are given by the authors. Computational results shed light on the tradeoffs between various
performance measures.
Speed and congestion as well as payload are very important factors for fuel consumption. However,
for alternative fuel vehicles like CNG trucks, some other factors such as fuel stations, may also have a
large influence on fuel consumption. Compared with the traditional diesel or gasoline vehicles, the sparse
allocation of fuel stations and small fuel tank capacity of alternative fuel vehicles call for an extra attention
on the refueling process. Traditional vehicle routing problems assume that the fuel is adequate for covering
the whole tour and the vehicles can be refueled any time. However, such an assumption may not hold for
CNG trucks. Surveys (Ahmed et al., 2013; Kuby et al., 2013; Ma et al., 2013) show that CNG vehicle drivers
suffer from issues such as long waiting time at fuel stations and long detour off their least-time routes. Since
GVRP for alternative fuel vehicles is a relatively new area, the existing research regarding fuel stations is
limited.
So far only few papers have considered fuel stations in their routing problems. In those papers, different
models and solution methods are proposed based on the unique features for their specific background.
Erdoğan and Miller-Hooks (2012) extend CVRP to GVRP where the routing problem for alternative fuel-
powered vehicle fleets is solved. Limited by its tank capacity, an alternative fuel vehicle has to visit one or
more times the fuel stations during its tour to serve the customers. Two construction heuristic methods and
a customized improvement technique are developed to minimize the total distance traveled by the vehicles.
Schneider et al. (2014) present a hybrid heuristic that combines a Variable Neighborhood Search (VNS)
algorithm with a Tabu Search (TS) to solve an electric vehicle routing problem. They construct an initial
solution without guarantee of feasibility and set penalties on the violations of feasibility so that the following
11
search algorithm tends to focus on feasible solutions. The hybrid heuristic is tested on newly designed
instances based on the traditional cases, as well as special benchmark cases from other papers like Erdoğan
and Miller-Hooks (2012). These tests show that hybridization has a positive effect on the solution quality.
Hiermann et al. (2016) propose a more complex model to solve a routing problem for a heterogeneous fleet.
They develop a hybrid heuristic which combines an ALNS with an embedded local search and labeling
procedure for intensification. The choice of recharging stations in each electric vehicle route is explicitly
handled with a labeling procedure. Their method takes advantage of tight time windows, whereas our
proposed method does not depend on time windows to reduce the search space.
2.2 Location Routing Problem for Fuel Stations
There is one research stream which is closely related to our model, which is called Location Routing Problem
(LRP). The LRP has received extensive attention and it inspired a growing stream of research in the last
decade. Several papers (Drexl and Schneider, 2015; Min et al., 1998; Nagy and Salhi, 2007; Prodhon and
Prins, 2014) have conducted comprehensive surveys and present detailed categorizations of this research
topic. A classical LRP problem solves two inter-dependent problems, the depot location and vehicle routing
problems, simultanousely. It solves the location problem for the depots where trucks depart from and return
to. Meanwhile, the customers allocated around are assigned to a depot nearby and their demand is fulfilled
by trucks from that depot. Starting from uncapacitated depots (Tuzun and Burke, 1999), the background of
the problem gets more complicated. In some papers, more contraints, such as depot and/or vehicle capacity
limitation (Barreto et al., 2007; Doulabi and Seifi, 2013; Prins et al., 2007), are added into the models so that
the solutions are closer to realistic problems. Some papers modify the structure of the problems from the
regular two-layer (depots and customers) problems to the multi-layer (factories, warehouses and customers)
problems (Baldacci et al., 2013; Gendron and Semet, 2009). As to the methodology, both exact and heuristic
algorithms are proposed. The first exact algorithm for the general LRP is by Laporte and Nobert (1981).
They use a branch-and-bound algorithm to solve a single depot location and vehicle routing problem. Other
exact algorithms like branch-and-cut, column generation are developed (Baldacci et al., 2013; Belenguer
et al., 2011; Contardo et al., 2013). Since LRP is a combination of two already NP-hard problems (Prodhon
and Prins, 2014), the size of the problem which can be solved by exact algorithm is limited. Thus, there
is a significant amount of papers on heuristic methods that have been published. These heuristic methods
are categorized based on different dimensions in survey papers (Nagy and Salhi, 2007; Prodhon and Prins,
2014).
12
Though the LRP has been studied for decades, research on the LRP usually does not include intermediate
depots or refill points because the location decisions are implicitly determined by the routing (Drexl and
Schneider, 2015). The intermediate depots or refill points related problems are special compared with regular
depots in the classical LPR. In the classical LRP a truck is required to visit a depot at the beginning and
end of the day, while in those problems a truck can visit a fuel station at any time based on its gas level
and it may visit a fuel station for multiple times in a single day. Since our CTLRP are about the location
problems for fuel stations, we focus more on the LRP papers for fuel stations.
The research on the LRP for fuel stations can be divided into two categories, vertex based and flow
based (Schiffer and Walther, 2017). The vertex based approaches are mainly developed from location models,
including maximum covering, p-center, p-median and so on. Details about the location models can be found
in the book by Daskin (2011). Within these location models, the p-median model is perhaps the most
popular one and the implementations of the p-median model on the location decision for fuel stations is
supported by its appeal of locating stations convenient to where people live (Upchurch and Kuby, 2010). A
classical p-median model solves a location problem of the p facilities in a network so that the total cost is
minimized. The total cost is given by the summation of the product of the demand at each node and the
distance from that node to the nearest facility Daskin (2011). It is believed that the p-median model is first
introduced by Goodchild and Noronha (1987) to solve the location problem for fuel stations. Since then,
p-median model has been widely used in solving the fuel station location problems. For instance, Nicholas
et al. (2004) and Nicholas and Ogden (2006) use p-median model to evaluate the locations of hydrogen
stations. Wang and Wang (2010) modify the p-median model to solve a dual-objective problem, in which
they include the construction cost of the fuel stations as well as the operation cost from the product of the
population of a node and the distance to its corresponding fuel station. In Upchurch and Kuby (2010)’s
work, the solutions from the p-median model served as a benchmark for the location problem. Though
the distance factor is considered in the vertex based approaches, most of the studies assume the demands
are fulfilled by vehicles directly from the fuel stations to demands, namely in the format of full truckload
(FTL). The routing problem is not simultaneously solved with the location problem. In some papers, the
calculation of the routing part takes place at a preprocessing step, which is not solved with the location
problem simultaneously (Schiffer and Walther, 2017).
The other research stream is called flow based model. Hodgson (1990) is one of the earliest researchers
who propose the concept of Flow Capturing Location Allocation Model (FCLM). Different from the vertex
based models which assume that demands are located at the node of the network, the FCLM assumes that
13
the demands are allocated by the origin-destination (O-D) paths and are captured by the facilities on those
paths. In Hodgson (1990)’s work, the flow-capturing problem is modified as a maximum covering problem,
in which the locations of the p facilities are chosen so that the captured flows are maximized. Based on
the FCLM, Kuby and Lim (2005) propose a new model called Flow Refueling Location Model (FRLM) to
solve the location problem for alternative fueling stations. In their model, multiple facilities are allowed on
a path and the problem is solved by exogenously determining the potential facility combination on a path.
Kim and Kuby (2012) extend FRLM by allowing drivers to deviate from their shortest paths to get refueled.
The demand flows captured by the stations are assumed to decrease as the deviation increases. Capar et al.
(2013) develop a new solution method, which focuses on the covering of arcs that comprise each path, to
solve the FRLM more efficiently. Flow based models are different from our CTLRP because in most of
the flow based models, demands, or called flows in many of those papers, are given in the format of O-D
pairs. The models are designed to maximize the flows covered by the fuel stations, or minimize the detours
of vehicles from the directed O-D pairs. However, in our CTLRP model, demands are allocated at nodes
within an area and they are fulfilled by LTL shipment. O-D pairs are unknown and implicitly determined
by the routing problem. Thus, flow based models can not be directly used to solve our problem.
As far as we know, only a few papers are similar to our CTLRP model. Yang and Sun (2015) consider
a LRP for electric vehicles battery swap stations. They propose a MIP formulation to optimally solve the
problem, as well as a four-phase heuristic to solve problems of large size. Schiffer and Walther (2017) present
another LRP under the background of charging stations with the allowance of partial recharging. They
develop a MIP formulation to solve the problem and some strengthen constraints to improve the speed.
Different objective functions are discussed in their paper. To date, no paper has solved the LRP with fixed
fueling time, which is not negligible for CNG fueling stations. Moreover, the exact and heuristic algorithms
of the LRP for fuel stations can be further improved.
2.3 Pickup and Delivery Problem
2.3.1 Exact Algorithms for Pickup and Delivery Problem
There are two primary approaches to solve the pickup and delivery problem optimally, which are branch-and-
price algorithm and branch-and-cut algorithm (Ropke and Cordeau, 2009). In general, a branch-and-price
algorithm has a better performance when dealing with a routing problem with a small search region, for
instance, a routing problem with tight time windows (Dumas et al., 1991). A branch-and-cut algorithm has
advantages in solving a routing problem with a large feasible region. In the rest of this section, we review
14
papers on these two algorithms.
Branch-and-price algorithm, or called column generation in some papers, was first proposed by Dumas
et al. (1991) in solving a pickup and delivery problem for a fleet of heterogeneous trucks with known demands
and time windows from the customers. In the paper, Dumas et al. (1991) separate the problem into a
master problem and a subproblem. For the master problem, also called set partitioning formulation, the
variables correspond to admissible routes for vehicles. The subproblem, solved by a dynamic programming
algorithm, generates the admissible routes for the master problem. Many efforts have been made to improve
the performance of this method. Instead of a dynamic programming algorithm, Xu et al. (2003) use a
fast heuristic method to solve the subproblem. Ropke and Cordeau (2009) propose several kinds of valid
inequalities and dynamically add them into the set partitioning formulation. They name their method as
branch-and-cut-and-price algorithm. Baldacci et al. (2011) reduce the search area of the set partitioning
formulation by computing a near-optimal dual solution of its linear relaxation. Stålhane et al. (2012) use
the branch-and-price algorithm to solve the pickup and delivery problem with time windows and split loads.
They modify the subproblem to be a multi-dimensional knapsack problem. Recently, Mahmoudi and Zhou
(2015) propose a relatively new algorithm, which is different from the traditional branch-and-price algorithm.
They use a Lagrangian relaxation approach to decompose the multi-vehicle routing problem to be a single
vehicle routing subproblem and the subproblem is solved by a dynamic programming algorithm with a new
state labeling method.
The other exact algorithm, called branch-and-cut algorithm, solves the problem by adding valid inequali-
ties to the MIP formulation. These inequalities are redundant for the formulation, but can reduce the search
space. Lu and Dessouky (2004) formulate the pickup and delivery problem by introducing extra binary vari-
ables to represent the sequence and precedence constraints and solved it by a branch-and-cut algorithm with
four kinds of valid inequalities. Based on previous work on other similar problems, Cordeau (2006) develops
a branch-and-cut algorithm to solve a pickup and delivery problem. In his branch-and-cut algorithm, valid
inequalities from time, capacity, subtour elimination, precedence and other constraints are added into the
MIP formulation. More cuts from subtour elimination, generalized order and other constraints are selected
by a heuristic method and dynamically added into the brand-and-bound procedure. Ropke et al. (2007)
introduce a new formulation for the pickup and delivery problem with time windows and solve it by adding
families of valid inequalities. Dumitrescu et al. (2010) presents new polyhedral results, valid inequalities and
separation procedures for the traveling salesman problem with pickup and delivery.
15
2.3.2 Metaheuristics for Pickup and Delivery Problem
Parragh et al. (2007)’s survey shows Tabu Search and Adaptive Large Neighborhood Search (ALNS)/Large
Neighborhood Search (LNS) are two popular metaheuristic solution methods proposed for the pickup and
deliveryproblem. Anumberofpapersshowtheefficiencyofthesetwomethods. Inrecentyears,metaheuristic
solution methods are getting more and more sophisticated. Researchers combine different methods together
to take advantage of the structures of the problems. In the rest of this section, a review of Tabu Search and
ALNS/LNS is given.
Tabu Search was developed by Glover in his papers (Glover, 1989, 1990). It is a strengthen form of local
search, which allows worsening moves and prohibits revisiting a neighbor state in a short time. Tabu search
is first used by Landrieu et al. (2001) to solve a pickup and delivery problem with time windows. They are
able to solve a case with a single vehicle and up to 15 requests. Nanry and Barnes (2000) present a reactive
Tabu Search approach with three dedicated move neighborhoods that capitalize on the dominance of the
precedence and coupling constraints. They solve problem instances for multiple vehicles with 100 requests.
Li and Lim (2003) propose a Tabu-embedded Simulated Annealing method with a restart trigger. Many
other papers use Tabu Search to solve the pickup and delivery problem under different backgrounds (Cordeau
and Laporte, 2003; Wang et al., 2016).
Shaw (1998) is the first person who uses the LNS method to solve vehicle routing problems. Compared
with the neighborhood search, LNS enlarges the search space by allowing destroy and repair operators.
Based on Shaw’s work, Ropke and Pisinger (2006) propose ALNS for the pickup and delivery problem with
time windows. ALNS enriches the destroy and repair operators. Every time when needed, ALNS selects
one destroy operator and one insert operator from its multiple candidates based on the operators’ previous
contributions. Bent and Van Hentenryck (2006) create a two-stage hybrid algorithm for the routing problem.
The first stage is a Simulated Annealing method and the second stage uses LNS. Recently ALNS/LNS has
been enriched by solving routing problems with a specific objective function or constraints (Hiermann et al.,
2016; Masson et al., 2013).
2.3.3 Related Alternative Pickup and Delivery Problems
Thepickupanddeliveryproblemhasbeenintensivelystudiedinrecentyears, whilethehybridtransportation
system has been approached by only few researchers. Aldaihani and Dessouky (2003) extend the general
vehicle routing problem to an integrated problem dealing with combining a fixed route service with on-
demand vehicles. The solution method includes three phases, which are candidate fixed route selection,
16
request insertion and solution improvement. Quadrifoglio et al. (2007) develop a similar system called
Mobility Allowance Shuttle Transit (MAST) system in which vehicles may deviate from a fixed path so
that customers within an area can be picked up or dropped off with flexibility. The routing problem for
the MAST system is solved by an insertion based heuristic method. Later on, Quadrifoglio et al. (2008)
develop an exact algorithm, which is a MIP model, to solve the routing problem for the MAST system.
They propose sets of logic cuts by using reasonable assumptions on passengers’ behavior. Masson et al.
(2014) solve a pickup and delivery problem with shuttle routes, in which a customer’s trip may be divided
into two trips, first trip for vehicles and second trip for fixed-route shuttles. They solve the problem with
a branch-and-cut-and-price algorithm with five families of valid inequalities for the master problem and a
new labeling algorithm for the subproblem. Ghilas et al. (2016a) propose a slightly different model, in which
freight flows, instead of passengers, are integrated with scheduled public transportation services. The model
is called the pickup and delivery problem with time windows and scheduled lines (PDPTW-SL) and solved
by a MIP formulation with strengthen constraints. Ghilas et al. (2016b) extend the PDPTW-SL further with
stochastic demand. They develop an adaptive large neighborhood search embedded with a sample average
approximation method to solve the problem in a scenario-based way.
Although only a few papers solve exactly the same model as ours, the bus stops in the hybrid transit
system are mathematically similar to regular transfer points. There are some papers about the pickup and
delivery problem with transfer points. Mitrovic-Minic and Laporte (2006) solve this specific kind of problem
by an insertion heuristic algorithm. The solution quality is improved by Masson et al. (2013) via using an
ALNS embedded with a new destruction and repair heuristic dedicated to the use of transfer points. Cortés
et al. (2010) propose a branch-and-cut algorithm which is solved by Benders’ decomposition with a master
problem for network flow balance and a feasibility problem for logical constraints. Unique valid inequalities
are also added to speed up the search for the optimal solution.
17
3 CNG Truck Routing Problem with Fuel Stations
This chapter first describes the CTRPFS problem, and then a MIP model is presented to optimally solve
the problem. In order to improve the solving speed, variable elimination and valid inequalities are provided
and an adaptive branch-and-cut algorithm is proposed. We also develop a heuristic method modified from
ALNS to solve a problem of large size. In the last part of this chapter, we show the computational results
of the proposed methods and explore some insights from the results.
3.1 Problem Statement
CTRPFS is designed to solve the vehicle routing problem for a distributor that employs a fleet of CNG
trucks to daily deliver products to customers. Customers are widely distributed in a region, and demands
from customers are known ahead of time. We assume the distributor’s depot does not own CNG fueling
stations by itself, which means its trucks have to use the sparsely distributed public CNG fueling stations
during their delivery routes. Thus, we need to develop route plans for the fleet which minimizes the total
working time while meeting constraints such as fuel tank capacity and load capacity limitations. A typical
route for the fleet may start from a departure depot, visit some customers and fuel stations if necessary, and
return back to the departure depot at the end of the day.
The main objective for this problem is to minimize the total working time for the trucks, while satisfying
three kinds of constraints which are load capacity, daily working time and fuel tank capacity limitations.
We define the working time as the time when the trucks are away from the depot. It includes traveling time
in the routes, fixed waiting time in the stations and variable fueling time which depends on the volume of
the fuel that needs to be refilled. To simplify the model, we assume all of the trucks travel at a constant
speed. Furthermore, we assume fuel is also consumed at a constant rate. However, our solution methods
are arc based, so we can easily bring in nonlinearity by calculating the time and gas needed for each arc
individually. As previously mentioned, the fixed fueling time could be the waiting time for the availability of
the pump, which can be as long as hours in the case of CNG vehicles (Kuby et al., 2013). In this model,
we use the expected waiting time and ignore the congestion effect in the fuel stations. For
simplification, we assume the fixed fueling time is identical for all fuel stations. The variable
fueling time is set to be proportional to the volume of fuel that needs to be refilled. Another assumption
in our model is a truck can visit at most one fuel station between two adjacent customers in its route. In
some very special cases, the optimal route for a truck may contain two successive visits of fuel stations. For
18
instance, two customers are so far away that a truck has to visit multiple fuel stations in order to travel from
one customer to another. This may happen when two customers are located in two different cities. Since
our problem solves the daily routing problem for a fleet within a city and we assume the tank capacity is
high enough to travel from any fuel station in the region to any demand point in the region, we ignore such
a situation. Initially, all trucks are located at the depot and with a given level of gas left in their tank and
all trucks are required to return back to the depot with a minimal level of gas so that it can at least go to
the nearest fuel station on the next day.
3.2 Mixed Integer Programming Model
An instance network of CTRPFS consists of a set of vertices V and a set of arcs A. Vertices in set V are
marked from 0 to n. Vertices v
0
and v
n
are the departure and arrival depots respectively. Though v
0
and
v
n
are two different vertices in our network, in reality they could be located at the same place. The rest of
the vertices are divided into two parts, customers’ vertices C and fuel stations vertices F. Since a truck is
allowed to visit a station multiple times in its single route, we replicate each station enough times in the
network, so that we can keep track of the status (load, gas level, etc) when a truck visits a fuel station. An
arc a
k
i;j
in set A represents a directed route from vertex v
i
to vertex v
j
for truck k in fleet M. The arcs
in the network include arcs from departure vertex v
0
to any vertex in C[F, from any vertex in C[F to
arrival vertex v
n
, arcs between any two different vertices inC, and arcs between a vertex inC and a vertex
in F. Note that we do not include arcs between vertices in F, namely we do not allow successive visits of
fuel stations. We introduce more notation for the MIP formulation. Let V
+
i
andV

i
be the sets of vertices
which have arcs from and to vertex v
i
, respectively. Overall, the fleet has m trucks and all of them need to
depart from vertex v
0
and arrive at vertex v
n
at the end of the day.
In the MIP model, binary variable x
k
i;j
is equal to 1 if arc a
k
i;j
is traveled by truck k in the solution and
0 otherwise. Continuous variables t
k
i
, g
k
i
and c
k
i
keep track of the working time from the departure depot,
remaining gas level and current delivered load of truck k at vertex v
i
, respectively.
Let constant R
1
represent the fixed fueling time for a truck’s single visit to a fuel station. Constants T
and G are the limitations for daily working time (in minutes) and fuel tank capacity (in gallons). Constant
L is the limitation of a truck’s cargo load. We assume all trucks have identical T, G and L. Let
1
,
2
and
3
be the speed of the trucks (miles/minute), fuel consumption rate (miles/gallon) and fueling speed
(gallons/minute) respectively.
i;j
is the distance from vertex v
i
to vertex v
j
.
i
is the load for vertex v
i
. It
is the customer’s demand if the vertex is a customer vertex and 0 otherwise. All of the customers’ demands
19
are known in our model. Constants
k
are the initial fuel level for the trucks.
k
should be set to be a
number between min
j2F

j;n

2

andG since within this range it guarantees a truck has enough fuel to go to
the nearest fuel station at the beginning of the period.
3.2.1 Basic Formulation
The basic MIP model for CTRPFS is as follows.
min
X
k2M
t
k
n
(1)
s.t.
X
i2V

j
x
k
i;j

X
i2V
+
j
x
k
j;i
= 0 8j2C[F;k2M (2)
X
i2V

j
X
k2M
x
k
i;j
= 1 8j2C (3)
X
j2V
+
0
x
k
0;j
1; 8k2M (4)
X
i2V

j
x
k
i;j
1; 8j2F;k2M (5)
t
k
i
+M

x
k
i;j
1

+

i;j

1
t
k
j
; 8i2f0g[C;j2V
+
i
;k2M (6)
t
k
i
+M

x
k
i;j
1

+

i;j

1
+
Gg
i

3
+R
1
t
k
j
; 8i2F;j2V

i
;k2M (7)
0t
k
i
T; 8i2V;k2M (8)
g
k
i
+M

1x
k
i;j

i;j

2
g
k
j
; 8i2f0g[C;j2V
+
i
;k2M (9)
M

1x
k
i;j

+G

i;j

2
g
k
j
; 8i2F;j2V
+
i
;k2M (10)
0g
k
i
G; 8i2V;k2M (11)
g
k
0
=
k
; 8k2M (12)
g
k
n
min
j2F

j;n

2

; 8k2M (13)
c
k
i
+M

1x
k
i;j

i
c
k
j
; 8i2f0g[C[F;j2V
+
i
;k2M (14)
0c
k
i
L; 8i2V;k2M (15)
x
k
i;j
2f0; 1g; 8i2V;j2V
+
i
;k2M (16)
Objective function (1) is to minimize the total working time, which includes traveling time, fixed waiting
20
time and variable fueling time. Constraints (2) and (3) are the flow conservation equations that ensure that
the solution defines a route plan for the fleet to visit every customer. Constraint (4) restricts that each truck
can be used at most once. Constraint (5) guarantees each replicate of the fuel station can be visited by the
same truck at most once. Constraints (6) and (8) are working time related constraints. The daily working
time increases due to the traveling time only if the truck travels through the regular arcs, which is represented
inconstraint(6). Itincreasesduetothetravelingtime, variablefuelingtimeandfixedfuelingtimeifittravels
through arcs adjacent to the fuel stations, which is shown in constraints (7). Constraint (8) bounds the daily
working time between 0 and T. The fuel tank capacity constraint is covered by constraints (9) to (13).
Constraint (9) decreases the remaining gas level if a truck travels through a regular arc. Constraint (10)
means if a truck travels through an arc starting from a fuel station, the remaining gas level will be set to
the full gas level minus the gas used for that arc, namely the truck is fully refueled at the fuel station. The
remaining gas level boundary, initial gas level and final gas level are set in constraints (11), (12) and (13)
respectively. Constraints (14) and (15) guarantee the demand for each customer can be fulfilled by the
assigned truck. Constraint (16) is the domain constraint.
3.2.2 Preprocessing and Valid Inequalities
Before solving the MIP formulation, some variables are eliminated based on the feasibility of the arcs,
and parameters like big M are appropriately selected; and valid inequalities which are redundant but can
strengthen the formulation are added to the model.
Variable elimination. Variables are eliminated based on different feasibility rules. The first elimination
rule is based on the tank capacity limitation. Variable x
k
i;j
is eliminated if the corresponding arc distance

i;j
satisfies one of the conditions (17) to (20).

i;j
>
2
G (17)
min
i
0
2V

i

i
0
;i
+
i;j
+ min
j
0
2V
+
j

j;j
0 >
2
G; i2C;j2V
+
i
[C (18)

i;j
+ min
j
0
2V
+
j

j;j
0 >
2
G; i2F;j2V
+
i
(19)
min
i
0
2V

i

i
0
;i
+
i;j
>
2
G; i2C;j2V
+
i
[F (20)
Condition (17) is intuitive. A truck can not travel on an arc if it is more than the tank capacity limitation
G. Condition (18) is a strengthened version of condition (17), which means an arc between two customers is
21
not reachable if a truck can not travel from any other place to the arc and go to any other place. Similarly, in
conditions (19) and (20), an arc away from a fuel station, or to a fuel station is not reachable if the minimal
distance of an adjacent arc combined with this arc is larger than the tank capacity limitation.
The second rule is for symmetry elimination. Note that this rule works only if the fleet is homogeneous,
namely all of the trucks are identical, including the initial gas level. If the fleet is homogeneous, two solutions
are identical if we can transfer from one to the other by changing the numbering of the trucks. The following
symmetry elimination conditions (21) to (23) are modified from Cortés et al. (2010). Condition (21) implies
that customer vertices inC are visited in order by each one of the trucks. For instance, for customer 2, the
demand can only be fulfilled by truck 1 or 2. Similarly, conditions (22) and (23) imply that replicates of
the fuel stations are visited in order. In these two equations, we introduce notation (i) to represent that
the corresponding vertex in the network is the (i)th replicate of a fuel station. An instance of these two
constraints is a truck may only go to the first and second copies of the fuel stations immediately after or
before visiting customer 2’s vertex.
x
k
i;j
= 0; 8i2C;j2V
+
i
;k2M;k>i (21)
x
k
i;j
= 0; 8i2C;j2V
+
i
[F; (j)>i;k2M (22)
x
k
i;j
= 0; 8j2C;i2V

j
[F; (i)>j;k2M (23)
The third elimination rule is induced by the domination of the fuel stations. We say a fuel station vertex
v
f1
dominates another v
f2
between a pair of customer vertices v
i
and v
j
, if
i;f1

i;f2
and
f1;j

f2;j
are
satisfied. A fuel station vertex v
f2
can never be an option between v
i
and v
j
in the optimal solution if it is
dominated by any other fuel station. Variables x
k
i;f2
is eliminated if the fuel station vertex v
f2
is dominated
by other fuel stations in any pair between v
i
and a vertex successive to v
i
, namely condition (24) holds.

i;f1

i;f2
;
f1;j

f2;j
; i2C;f
1
;f
2
2F;f
1
6=f
2
;8j2V
+
f1
[V
+
f2
(24)
The fourth elimination rule is induced by the remaining gas level after visiting a fuel station. Variablex
k
i;j
representing an arc from customer vertexv
i
to fuel station vertexv
j
by truckk is eliminated if condition (25)
holds.

i;j

i;j
0

2
>G

j;j
0

2
; 8j
0
2V
+
j
(25)
22
A truck can travel from vertex v
i
to v
j
only if the remaining gas can cover distance
i;j
. The left-hand side
of condition (25) is the minimal remaining gas at vertex v
j
0 if the truck goes directly from vertex v
i
to v
j
0,
while the right-hand side is the remaining gas if the truck goes to the fuel station. A visit from customer
vertex v
i
to a fuel station vertex v
j
is not necessary if the remaining gas level after the visit is always lower
no matter where the truck goes.
The fifth elimination rule is also related to the visiting of a fuel station. A visit from customer vertex v
i
to a fuel station vertexv
j
is impossible if the minimal tour including this arc is larger than the daily working
time limitation. The right-hand side of condition (26) shows the minimal tour.

0;i
+
i;j
+
j;0
>
1
(TR
1
) (26)
The last elimination rule concerns the load capacity limitation. x
k
i;j
is eliminated if condition (27) is
satisfied, namely the demands
i
and
j
are more than a truck’s load capacity.

i
+
j
>L (27)
Parameter setting Parameters like bigM and the number of replicates of each fuel station are explicitly
handled. In theory the big M in constraints (6), (7), (9), (10) and (14) can be any number as long as it
is sufficiently large. However, an arbitrary big M may bring in some numerical problems and weaken the
bound from the linear relaxation of the MIP. We need to choose M as small as possible to accelerate the
computation speed. For constraints (6) and (7) the minimal M should satisfy that even if t
k
i
reaches its
maximal possible value, the minimal possible value oft
k
j
will not be influenced whenx
k
i;j
is 0. Equations (28)
and (29) give the detailed representation ofM for constraints (6) and (7). Similarly, bigM in constraints (9)
23
and (10) are given in equations (30) and (31). For constraint (14) big M is simply set to be 1.
M =

T

i;n

1

0;j

1
+

i;j

1
(28)
M =

T

i;n

1

0;j

1
+

i;j

1
+
G

3
+R
1
(29)
M = maxfg
1
;g
2
g min
i
0
2V
+
i
[F

i;i
0

2
+

i;j

2
;
where g
1
=

0;j

2
; g
2
=G min
j
0
2V

j
[F

j
0
;j

2
(30)
M = maxfg
1
;g
2
gG +

i;j

2
; (31)
where g
1
=

0;j

2
; g
2
=G min
j
0
2V

j
[F

j
0
;j

2
(32)
As we stated at the beginning of Section 3.2, we replicate the fuel stations in case a truck needs to visit
a fuel station multiple times. Too many replicates will increase the number of variables dramatically. The
upper bound on the number of times a truck visits a fuel station is given byd
2T
1

2G
e + 1.
Figure 1: An example of a truck’s multiple visits to a fuel station
Figure 1 shows an example of a truck’s multiple visits to a fuel station. A truck starts from a depot,
goes through a fuel station, visits some customers in the subtours and goes back to the depot. We claim
that only one subtour’s distance can be less than

2G
2
. Otherwise, the subtours can combine together with
less distance because of the triangle inequality with one less visit to a fuel station. Thus the total number
of times a truck visits a fuel station is bounded byd
2T
1

2G
e + 1.
Valid inequalities Some valid inequalities are added to the model to strengthen the linear relaxation
of the problem. The valid inequalities we add are from the tank capacity limitation. If the total distance of
two successive arcs is larger than
2
G, a truck can not travel through both of the arcs. Suppose the vertices
24
are v
i1
;v
i2
and v
i3
, Inequality (33) is added to the formulation if
i1;i2
+
i2;i3
>
2
G.
x
k
i1;i2
+x
k
i2;i3
1; 8k2M (33)
3.2.3 Adaptive Cut Generation
The branch-and-cut algorithm is a popular algorithm for solving combinatorial optimization problems. We
generate the initial pool of cuts based on the load capacity limitation, and the dominance relationship in the
fueling process as well as the regular subtour elimination from Cordeau (2006).
The first group of cuts is from the load capacity limitation. If a subsetC
0
of the customers’ cargo load is
more than a truck’s load capacity, those subsets of customers can not be served by a truck. We select those
subsets whose cardinality is 3, namelyjC
0
j = 3, and add inequality (34) into the initial pool of cuts for all
suchC
0
.
X
i2C
0
X
j2C
0
;i6=j
x
k
i;j
jC
0
j 2 8k2M (34)
The second group of cuts is with the remaining gas level after visiting a fuel station, which is similar
to inequality (25). Let v
i1
;v
i2
represent two customer vertices, and v
j
represents a fuel station vertex.
Inequality (35) is added into the initial pool of cuts if
2
(
i1;j

i1;i2
)G
2

j;i2
.
x
k
i1;j
+x
k
j;i2
1; 8k2M (35)
We generate the last group of cuts based on the subtour elimination inequality (36), whereC
0
is a subset
of the customers.
X
i;j2C
0
;i6=j
x
k
i;j
jC
0
j 1 (36)
Cuts are added to the pool if the cardinality ofC
0
is 2, and the corresponding two vertices do not appear in
any set from (34).
Preliminary experiments show that a significant amount of time is spent in searching for violated in-
equalities in the pool of cuts. We propose an adaptive cut generation to improve the searching efficiency.
The idea of adaptive cut generation is a group of cuts will be checked with higher probability during the
branch-and-bound process if this group of cuts are likely to be violated. We set p
i;t
to be the probability of
25
checking groupi’s cuts in iterationt. In our model,i = 1; 2; 3, refer to cut the groups from the load capacity
limitation, remaining gas level and subtour elimination respectively. Initially p
i;0
= 1;i = 1; 2; 3. We define
lhs
i
and rhs
i
;i = 1; 2; 3 as the left-hand and right-hand side of inequalities (34), (35) and (36), respectively.
Take one cut j in inequality (34) where i = 1 as an example,
lhs
j
1
=
X
i2C
0
X
j2C
0
;i6=j
x
k
i;j
(37)
rhs
j
1
=jC
0
j 2 (38)
Let lhs
j
i;t
be the value of lhs
j
i;t
in iteration t by substituting x
k
i;j
with its value in the iteration. Since the
right-hand side of each group of cuts has the same constant value, we simply use rhs
i
;i = 1; 2; 3 to denote
the right-hand side values. We update p
i;t
based on equation (39)
p
i;t+1
= (1r)p
i;t
+r
Y
j
lhs
j
i;t
rhs
i
(39)
wherer is a constant number and we set it to 0:01 for our problem. From equation (39) we know the score of
a cut is defined by the
lhs
j
i;t
rhsi
part. It is bigger than 1 if the cut is violated. Thus the corresponding group of
cuts has a higher probability to be checked in the next iteration if the cut is violated in the current iteration.
3.3 Heuristic Solution Method
ALNS is an efficient heuristic method for solving the classical truck routing problem. It is proposed by Ropke
and Pisinger (2006) from Large Neighborhood Search (LNS), which is first introduced by Shaw (1998) for
the routing problem. Hiermann et al. (2016) extend it by embedding a local search and labeling procedure
to solve a truck routing problem with time windows and recharging stations. In this section, we propose a
hybrid heuristic algorithm that combines an ALNS with a MIP model, which is different from the previous
research, to solve the routing problem with fuel stations. Our new method has better performance when the
search space is large. In our method, the search is conducted within the feasible region. This procedure has
an advantage that whenever the method is stopped, we will always get a feasible solution for the fleet.
Based on the settings in the problem statement, we introduce some more notations. Let r
k
be the route
fortruckk. r
k
(i)
istheithnodeinroutek. Figure2showsanexampleofroutek. Fromtheproblemstatement
in the Section 3.1, the first node r
k
(0)
and the last node r
k
(j)
are the depot vertices v
0
and v
n
. Within the
route, some nodes, for instance r
k
(i+1)
, may be fuel stations. Let Time

r
k

be the total working time for
26
Figure 2: An example of route k
route k and Load

r
k

be the total load for route k. Let Time

r
k
(i)

and Load

r
k
(i)

be the working time
and cargo load from the start depot until the current node r
k
(i)
. Let Fuel

r
k
(i)

be the remaining fuel level
at node r
k
(i)
. Note Time

r
k

and Load

r
k

are bounded by T and L respectively. Fuel

r
k
(i)

is bounded
byG andFuel

r
k
(0)

=,Fuel

r
k
(j)

min
f2F

f;n

2

. All of these boundaries are identical with the ones
in the MIP formulation in Section 3.2. In the heuristic solution method, between any two successive nodes
i and i + 1, only the fuel station with the smallest total detour are considered during the insertion. Let

+
r
k
(i)
;r
k
(i+1)
denote the distance from a customer node r
k
(i)
to the fuel station and

r
k
(i)
;r
k
(i+1)
denote distance
from the fuel station to a customer node r
k
(i+1)
. Letc
k
be the cost of truckk andC be the fleet’s total cost,
then equation (40) is the objective function of the hybrid heuristic algorithm. Thus the objective is identical
with the MIP model.
C =
X
k2M
c
k
=
X
k2M
Time

r
k

(40)
3.3.1 Initial Feasible Solution Construction
The construction of an initial feasible solution is modified from Ropke and Pisinger’s parallel insertion
heuristic (Ropke and Pisinger, 2006). Figure 3 shows the general steps for the construction process. We
start from an empty fleet. In each loop, we pick a customer who is not inserted into the fleet yet, go through
all positions in the current fleet for the customer to check if the route plan is still valid after the insertion,
select the position with the smallest additional traveling distance or add one more truck if there is not a
valid insertion position for the current customer, insert the customer into the fleet and repair the influenced
route if necessary. We repeat the procedures until all customers are inserted. The following section gives
a more detailed explanation on valid insertion, additional traveling distance, customer insertion and route
repair.
Suppose customer vertex v
i
is going to be inserted into truck k’s route, after the jth node. A valid
insertion means it meets four constraints, which are working time, fuel tank capacity, load capacity and
location constraints. An insertion is also valid if it only violates the fuel tank capacity limitation, which
can be simply repaired by adding a fuel station in front of the inserted customer while keeping all the other
27
Figure 3: Steps of constructing an initial feasible solution
constraints valid. Details about the four constraints are listed in the following equations (41) to (46).
Working time constraint. Equation (41) shows the criterion for the working time constraint. The daily
working time can not exceed T.
Time

r
k

+

r
k
(j)
;i
+
i;r
k
(j+1)

r
k
(j)
;r
k
(j+1)

1
T (41)
Fuel tank capacity constraint. For the fuel tank capacity constraint, equations (42) to (44) show its
criteria. Equation (42) is easy to understand. It ensures the truck has enough fuel to go to the fuel station.
In equation (43), Fuel

r
k
(j+1)

is the gas level at node r
k
(j+1)

where node r
k
(j+1)

is the first fuel station
after node r
k
(j+1)
, or the last node if there is no fuel station after it. The quantity in the left-hand side of
equation (43) is the gas needed to cover the range between the newly inserted fuel station and the nearest
following fuel station, or the last node r
k
(j+1)

. This quantity should be no bigger than the tank capacity
G. It eliminates a special case in which the distance from the fuel station back to vertex v
i
is so large that
inserting a fuel station into the route still can not make it feasible. Equation (44) ensures the truck still
meets the daily working time constraint if it needs to refuel. The working time changes due to the extra
traveling time and extra fueling time, and the extra fueling time includes changes in the variable fueling time
28
and the fixed fueling time.
Fuel

r
k
(j)

+
r
k
(j)
;i

2
(42)

r
k
(j)
;i

2
+Fuel

r
k
(j+1)

Fuel

r
k
(j+1)

G (43)
Time

r
k

+
t
+
f
T
where
t
=

+
r
k
(j)
;i
+

r
k
(j)
;i
+
i;r
k
(j+1)

r
k
(j)
;r
k
(j+1)

1

f
=
GFuel

r
k
(j)

1
2

+
r
k
(j)
;i

3

G
1
2

r
k
(j)
;i
+
i;r
k
(j+1)

Fuel

r
k
(j+1)

3
+R
1
(44)
Load capacity constraint. The load capacity constraint is explained in Equation (45). The load of truck
k can not be more than L.
Load

r
k

+
i
L (45)
Location constraint. The last constraint for a valid insertion is the location constraints. Given the
assumption that fuel stations can not be visited successively, node r
k
(j)
can not be a fuel station vertex.
r
k
(j)
= 2F (46)
Let d
k
i
denote the change in distance incurred by a valid insertion of customer vertex v
i
into route k at
the position that increases the distance for truck k the least. The change, or increase in distance, is simply
measured by
r
k
(j)
;i
+
i;r
k
(j+1)

r
k
(j)
;r
k
(j+1)
. Initially d
k
i
is set to be +1 and d
k
i
is updated by equation (47).
d
k
i
= min
j2J
k
i

r
k
(j)
;i
+
i;r
k
(j+1)

r
k
(j)
;r
k
(j+1)

(47)
whereJ
k
i
is the set of valid insertion positions for truck k for customer vertex v
i
. Finally we choose truck k
that minimizes d
k
i
and insert customer vertexv
i
at the minimum cost position in that truck. If d
k
i
= +1
for all k, there is no valid insertion for customer vertex v
i
and we need to use another truck route.
After inserting customer vertexv
i
into truckk, we need to process an inspection on the truck’s fuel tank
capacity constraint. If there exists a node r
k
(j)
whose Fuel

r
k
(j)

is less than 0, an insertion of the fuel
station is required. Though equations (42) to (44) ensure that inserting a fuel station with minimum detour
29
in front ofv
i
will make the route feasible again, we still want to find a feasible position j so thatFuel

r
k
(j)

is bigger than 0 while as small as possible. In reality, this strategy means that we try to delay the fueling as
much as possible. The truck will keep working until its fuel tank is almost empty. Intuitively, this delay will
minimize the number of visits to the fuel stations. The search of position j starting from the first node r
k
(m)
in router
k
whoseFuel

r
k
(m)

is less than 0. We setj =m and if two adjacent nodesr
k
(j1)
andr
k
(j)
satisfy
equations (42) to (44), we find feasible position j. Otherwise we setj =j 1 and check again, until we find
a feasible position. In the worst case, we will stop at the position where we just inserted the customer vertex
v
i
.
The insertion will stop until all customers are inserted into the fleet.
3.3.2 Adaptive Large Neighborhood Search
The ALNS is adapted from Hiermann et al. (2016). However, our search is conducted only within the feasible
region, which is different from the paper by Hiermann et al. (2016). Figure 4 shows the steps of one iteration
of the ALNS. We first do local search to intensify the search in each iteration. The next step is the destroy
and repair procedures from the general methodology of the ALNS. We also try to make the used fleet size as
small as possible in the following step. The last one is using a MIP model to find the optimal fueling plan
for each route.
Figure 4: Steps of the Adaptive Large Neighborhood Search (ALNS)
Local Search
Weintroduce some well-known localsearchmethods like shift, swap, 2-opt and2-opt* toimprovethesolution
quality in each search iteration. These methods are widely used in various kinds of heuristic methods for
truck routing problems (Chen and Wu, 2006; Hiermann et al., 2016; Montané and Galvao, 2006). In our
30
(a) Shift (b) Switch
(c) 2-Opt (d) 2-Opt*
Figure 5: Local search
paper, graphs from Hiermann et al. (2016) are used for illustration purposes.
Shift. Figure 5a shows the shift movement. We randomly select a node, for instance the ith node from
route k, and insert it into another randomly selected route k
0
.
Switch. In the switch movement, two nodes i and j are randomly picked from two routes k and k
0
, and
switch positions, as shown in Figure 5b.
2-Opt. In this movement, we first randomly select two positionsi andj, reverse the order of nodes between
i and j, and insert into the original position. Figure 5c shows this movement.
2-Opt*. 2-Opt* movement, also called 2-Exchange, is a movement within two routes. Figure 5d is an
example of a 2-Opt* movement. Two routes k and k
0
exchange the rest of the nodes at positions i and j.
To select the local search methods to use, we adapt the weighted roulette wheel selection method and
adaptive weight adjustment from Ropke and Pisinger (2006). Similar to their work, we define a segment as a
number of iterations of the local search; here we define a segment as 10,000 iterations. Let w
s
i;t
represent the
weight for local search method i at segment t. We select method i with probability given by equation (48)
at segment t. The weights are updated via equation (49) and initial weights w
s
i;0
;i = 1; 2; 3; 4 for the four
methods are equal.
w
s
i;t
P
4
j=1
w
s
j;t
i = 1; 2; 3; 4 (48)
w
s
i;t+1
=w
s
i;t
(1r
s
) +r
s

i;t

i;t
i = 1; 2; 3; 4 (49)
In equation (49), r
s
is a constant number,
i;t
is the score of method i obtained during segment t; here we
set it to be the decrease in the objective function, and
i;t
is the number of times we have attempted to use
31
method i during segment t. In our experiments, we select w
s
i;0
= 100 and r
s
= 0:5.
Note that the search is conducted within the feasible region. The first acceptance criterion is that the
new routes after the local search should still be feasible. The feasibility check procedures are similar to
equations (41) to (44) in Section 3.3.1. The second criterion is from the Simulated Annealing method, which
is also used by Ropke and Pisinger (2006). To avoid getting trapped in a local minimum, an inferior solution
may also be accepted with some probability. A solution will be accepted with probability
e
(C
0
C)=T
s
(50)
where C
0
is the cost of the new routes and C is the original one. T
s
is the temperature and starts with an
initial valueT
s
0
. In each iterationT
s
is updated viaT
s
=
s
T
s
, where 0<
s
< 1. We selectT
s
= 10 and

s
= 0:9999 for our experiments.
Since local search can be conducted very fast, we run
1
iterations in every ALNS iteration. If after a
local search iteration a route is empty, or has fuel stations only, it will be deleted. Based on our test, we set
= 100; 000 in the experiments.
Destroy and Repair Operators
The destroy and repair operators we introduce here are mainly from Ropke and Pisinger (2006). We get a
new solution by removing some nodes using the destroy operator, followed by the repair operator. These
operators will extend the search space by bringing a larger search neighborhood. We implement three destroy
operators and two repair operators. Notice only the nodes representing customer vertices will be considered
in the destroy and the following repair operators. We want to generate a feasible solution whenever the
algorithm stops. A destroy of a fuel station might make it hard to generate a feasible solution during the
repair procedure.
Random destroy operator. A random destroy operator is the simplest operator. Every node representing
a customer vertex is chosen with equal probability.
Related destroy operator. The idea of the related destroy operator is to remove similar nodes together,
and reinsert them into routes again. Intuitively if the removed nodes are far away from each other, the
solution after the repair operators may still be similar to the original one, or even worse. The relatedness
Related (i;j) for two vertices v
i
and v
j
is from Shaw (1998).
Related (i;j) =
1

0
i;j
+V
i;j
(51)
32
In equation (51),
0
i;j
is the normalized distance from vertex v
i
to vertex v
j
and
0
i;j
=
i;j
maxi;j
, where
i;j
is
the distance between vertex v
i
and vertex v
j
and max
i;j
is the maximal one for all of the distances. V
i;j
is
an indicator variable. It is 1 if vertex v
i
and vertex v
j
are in the same route and 0 otherwise. When using
the related destroy operator, the first node is randomly selected with the random destroy operator, followed
by 4 nodes, selected by the relatedness of the first node. The higher value the relatedness is, the higher
probability the node will be chosen.
Worst destroy operator The idea of the worst destroy operator is if the node leads to a high detour in
its route, it is likely to be misplaced and should be removed. This operator is used by Ropke and Pisinger
(2006). First we calculate the detour Detour(i) of the ith node in route k via equation (52).
Detour(i) =
r
k
(i1)
;r
k
(i)
+
r
k
(i)
;r
k
(i+1)

r
k
(i1)
;r
k
(i+1)
(52)
Then we choose the node to be removed based onDetour(i). The higher the detour is, the higher probability
the node will be chosen.
Cheapest route repair operator. This repair operator is exactly the same as the insertion method in
Section 3.3.1.
Cheapest customer repair operator. In this operator, we first calculate the cheapest insertion cost for
every customer removed by the destroy operators. The procedure is the same as the cheapest route repair
operator and is given in Section 3.3.1. Next we choose the customer with the lowest overall cost and insert
it into the corresponding route using the same weighted roulette selection method as in Section 3.3.1.. The
parameters in equations (48) and (49) for the destroy operator arew
d
i;t
;i = 1; 2; 3 and
d
. The parameters for
the repair operators arew
r
i;t
;i = 1; 2 and
r
. The segment is redefined as 100 iterations for both the destroy
and repair operators. After tuning the parameters, the initial values for w
d
i;0
;i = 1; 2; 3 andw
r
i;0
;i = 1; 2 are
all 100, and
d
=
r
= 0:8.
Every time before using the destroy and repair operators, we first decide how many nodes are going to
be removed. In our case, we randomly generate a number within
1
jCj and
2
jCj, where
1
and
2
are
constants and 0 <
1
<
2
< 1 andjCj is the total number of customers. Next we choose one destroy
and one repair operator for this iteration. The selection method is similar to the weighted roulette wheel
selection method and adaptive weight adjustment described in the Local Search part in this section. If a
route after the destroy operator is empty, or has fuel stations only, it will be deleted. In our experiments,

1
= 0:05 and
2
= 0:15.
33
Fleet Size Minimization
A lower bound for the fleet size can be calculated from the customers’ cargo load. LetS

be the lower bound.
S

=d
P
i2C

i
L
e (53)
If the current fleet size S is bigger than S

, we select SS

routes with the smallest number of customers,
delete these routes and put the customers into the request pool again. We insert those customers into the
fleet via repeating the construction procedure in Section 3.3.1. The new fleet will be accepted if its size S
0
is smaller thanS, otherwise we keep the original one. Though the new fleet size is usually smaller, the total
cost C is likely to be larger because of the relatively narrow view on the insertion procedure.
Preliminary experiments show the fleet size minimization procedure will dramatically slow down the
computing speed and increase the total cost. Thus, we set the frequency of this procedure to be very low.
We run it every
2
ALNS iterations. We set
2
to be 200 in our experiments so that the fleet size is usually
minimized, while the computing speed is not affected too much.
Fuel Station Optimization
Hiermann et al. (2016) also add the recharging stations into their model and implement a special labeling
algorithm to explicitly handle these stations. The idea of the labeling algorithm is to use a specific design
label to denote all of the possible recharging plans and pick the one with the lowest total cost. The success
of the labeling algorithm heavily depends on the elimination of the search space, namely the so-called
dominance checks in their paper. In their model, each customer has a relatively tight time window so the
search space is significantly reduced. In our setting, we do not have any time window constraints which
makes the dominance checks to be very weak. Without time window constraints, a detour to a fuel station
between two customers has a disadvantage on longer total traveling distance, but has likely an advantage on
having more gas left, which may save future detours to the fuel stations. This relationship holds for every
pair of successively visited customers, so the labels will extend exponentially even if we only consider one
possible fuel station between them. For instance, if a truck visits 35 stops a day, the search space may have
as many as 2
35
labels. In order to resolve this issue, we develop a MIP model to deal with the fuel stations.
Different from the MIP model in Section 3.2, the choice of fuel stations for each route is independent
now and each route can be solved separately. We take route r
k
for instance. First, delete all of the existing
fuel stations inr
k
and leave the nodes with the customer and depot vertices only. Let r
k
0
be the route after
deletion and suppose there are n
d
nodes left. In the heuristic method, for any two successive nodes in r
k
0
,
34
we just select the fuel station with the smallest detour as the only candidate. Let binary decision variable
z
k
(i)
denote the fueling decision for routek after theith node, which is customer noder
k
(i)
. Letz
k
0
(i)
be 1 if the
corresponding truck goes to fuel station after r
k
0
(i)
and be 0 otherwise. We follow the notation in Section 3.2.
For simplification, let z
i
represent z
k
0
(i)
, i represents customer node r
k
0
(i)
, t
i
represents t
k
0
i
, and g
i
represents
g
k
0
i
. The MIP is as follows.
min t
n
d
(54)
s.t. t
i
+

i;i+1

1
t
i+1
; 8i = 1;:::;n
d
1 (55)
t
i
+

+
i;i+1
+

i;i+1

1
+R
1
+
Gg
i
+
1
2

+
i;i+1

3
+M(z
i
1)t
i+1
; 8i = 1;:::;n
d
1 (56)
0t
i
T; 8i = 1;:::;n
d
(57)
g
i

i;i+1

2
+Mz
i
g
i+1
; 8i = 1;:::;n
d
1 (58)
M(1z
i
) +G

i;i+1

2
g
i+1
; 8i = 1;:::;n
d
1 (59)
g
i

+
i;i+1
z
i
; 8i = 1;:::;n
d
1 (60)
g
1
= (61)
g
dn
min
j2F

j;n

2

(62)
z
i
2f0; 1g 8i = 1;:::n
d
(63)
The objective function is to minimize the total working time. Constraints (55) and (56) keep track of the
working time from node to node. If the vehicle does not go to the fuel station, the working time increases
by the traveling time
i;i+1

1
, which is presented in constraint (55). Otherwise the working time increases by
the traveling time, fixed and variable fueling time, which is presented in constraint (56). Constraint (57)
gives the boundary of working time. Constraints (58) to (62) are the same as constraints (9) to (13), which
ensure the route meets the fuel tank capacity constraint.
In practice, we can easily construct an initial feasible solution by checking the original route plan before
the deletion of the fuel stations. Let z
i
be 1 if the truck goes into fuel station in the original route plan and
0 otherwise. Since the original route plan is feasible, the constructed solution must be a feasible solution for
the MIP. In our computational experiments, the MIP can be solved quickly for reasonable sized problems
(see Section 3.4). Similar to the fleet size minimization step, we run it every
3
ALNS iterations. In our
experiments, we set
3
to be 100.
35
Acceptance and Stopping Criteria
Similar to the rule used by Ropke and Pisinger (2006), in every ALNS iteration we accept the new solution
if it is better than the previous one, or with a probability calculated by equation (50) if it is worse. However,
the parameters T
a
0
and
a
are different from T
s
0
and
s
. In fact, these two probabilities are updated
independently.
The ALNS stops if the improvement of every iterations is less than a percentage, or the total number
of iterations reaches . We set T
a
0
= 1000,
a
= 0:999, = 100, = 0:1 and = 5; 000.
3.4 Computational Results
Data about CNG trucks are collected to help us set the parameters in our model. First we run the proposed
ALNS on small instances. These small instances are also solved by the MIP formulation using a standard
commercial optimization solver CPLEX. The small instances are used to show the efficiency of the pre-
processing and valid inequalities presented in Section 3.2.2 and our adaptive branch-and-cut algorithm in
Section 3.2.3 in reducing the computation time in finding an optimal solution. Results also serve as bench-
mark cases to show that our ALNS can get high quality solutions in a short amount of computation time.
Also we use our ALNS to solve moderate and large size instances to explore some interesting insights. The
moderate and large size instances include modified cases from Augerat et al. Augerat et al. (1995) and data
from freight demand originating from the Ports of Los Angeles and Long Beach in the Southern California
area.
All experiments are performed on a PC with Intel Core i7-4790K CPU and 32 GB RAM, running Linux
distribution Ubuntu 14.04 and CPLEX 12.5.
3.4.1 Results on Small Instances
We first create small instances by randomly generating 8 to 10 customers on a 100*100 map. A depot and 3
or 4 CNG fueling stations are also randomly picked within the map. The demands of customers are generated
from a uniform distribution between 0 to 0.6 unit. We solve these small instances with CPLEX 12.5, as well
as our proposed ALNS. In order to show the effectiveness of our valid inequalities and adaptive branch-and-
cut algorithm, three kinds of formulations, basic, improved and adaptive branch-and-cut, are tested. The
basic formulation includes the base MIP formulation and the big M selection. The improved formulation
includes everything presented in Section 3.2.2. The adaptive branch-and-cut cases include the improved
formulation and the adaptive cut generation method. Note we stop each run after 2 CPU hours if the exact
36
solution procedure has not yet found the optimal solution. The experiments are conducted based on the
background of CNG trucks and the parameters are modified from related materials (Laughlin and Burnham,
2016; U.S. Department of Energy, 2014). We set
1
,
2
and
3
to be 35 miles/hour, 3:75 miles/GGE (Gallon
of Gasoline Equivalent) and 3 GGE/minute. We set R
1
(fixed fueling time) to 20 minutes. Namely every
time when a truck goes to a CNG fueling station, its daily working time increases by 20 minutes. The daily
working time limitationT and tank capacityG are set to 8 hours (480 minutes) and 53.3 GGE, respectively.
The load capacity L for each truck is 1. The initial gas level
k
is a random number between min
j2F

j;n

2

and G.
Name
CPU Time (s) Solution Value
Opt Basic Opt Imp Opt B&C ALNS Opt # Opt Obj Gap
8C3F3T1 250:4 48:3 29:8 5:1 749:9 8 762:3 1:6%
8C3F3T2 115:0 53:3 25:2 3:2 1;033:4 10 1;033:4 0:0%
8C3F3T3 240:9 69:6 52:4 2:4 966:8 10 966:8 0:0%
8C3F3T4 96:6 18:6 11:5 1:7 792:9 6 795:1 0:3%
8C3F3T5 200:4 82:5 30:8 4:2 1;001:7 7 1;010:0 0:8%
Avg 180:6 54:4 29:9 3:3 908:9 8:2 913:5 0:5%
9C3F3T1 389:1 94:1 40:8 6:4 1;003:5 7 1;006:0 0:2%
9C3F3T2 1;135:4 436:1 339:2 5:9 1;115:2 10 1;115:2 0:0%
9C3F3T3 767:7 91:2 67:5 10:4 738:5 5 746:5 1:1%
9C3F3T4 517:2 38:4 28:3 11:5 1;031:0 7 1;040:7 0:9%
9C3F3T5 4;716:7 914:5 873:6 14:2 1;054:2 6 1;070:8 1:6%
Avg 1;505:2 314:9 269:9 9:7 988:5 7:0 995:9 0:7%
10C3F4T1 1;344:0 361:1 328:8 11:5 1;011:1 7 1;014:4 0:3%
10C3F4T2 2;452:1 2;435:7 1;512:6 10:4 1;049:2 5 1;052:1 0:3%
10C3F4T3 1;433:9 841:3 578:4 5:7 875:4 10 875:4 0:0%
10C3F4T4 > 7;200:0 3;436:9 3;192:7 17:1 1;184:9 6 1;196:4 1:0%
10C3F4T5 > 7;200:0 1;695:3 681:1 9:1 1;457:9 10 1;457:9 0:0%
Avg >3;926:0 1;754:0 1;258:7 10:8 1;115:7 7:6 1;119:2 0:3%
Table 1: Results on small CTRPFS instances
Table 1 presents the results for the small instances. In total we conduct experiments with 3 different
settings and each setting has 5 random instances. Name 8C3F3T1 means it is the first instance with a
37
setting of 8 customers, 3 CNG fueling stations and 3 trucks. The first part CPU Time has four columns.
Column Opt Basic, Opt Imp, Opt B&C are the CPU times (in seconds) used by the basic, improved and
adaptive branch-and-cut formulations. Column ALNS is the CPU time (in seconds) for our ALNS. Note that
our adaptive branch-and-cut algorithm and ALNS have randomness parts, so we run the adaptive branch-
and-cut 5 times and ALNS 10 times, and record the average of the CPU times. The three columns, Opt
Basic, Opt Imp and Opt B&C show that our valid inequalities and adaptive branch-and-cut algorithm can
save CPU time. On average, our adaptive branch-and-cut algorithm with the improved formulation finds
the optimal solution in the shortest amount of time. The second part Solution Value contains the solution
values from the MIP and the ALNS. The solution values represent the total working time (in minutes) for
the fleets in different instances. Column Opt is the optimal solution value. For the ALNS, we run 10 times
for each instance. Column # Opt is the number of times the ALNS obtained the optimal solution. Column
Obj is the average of the 10 solution values. Column Gap is the gap between the average value from the
ALNS, which is column Obj, and the optimal solution value, which is column Opt. The Solution Value part
indicates that for each instance, our ALNS can achieve the optimal solution within 10 runs, and the gap
between the average value and optimal solution is also very small. In these scenarios, the average gap is
below 1.0% and at least one run of the 10 can identify the optimal solution for each instance.
Name
Basic Imp B&C ALNS
Upper Lower Gap Upper Lower Gap Upper Lower Gap Best Avg Gap
15C4F5T1 1;495:7 829:8 44:5% 1;483:6 917:4 38:2% 1;429:4 904:2 36:7% 1;417:9 1;452:8 37:8%
15C4F5T2 1;447:5 823:6 43:1% 1;496:2 955:8 36:1% 1;390:0 990:2 28:8% 1;389:0 1;413:5 30:5%
15C4F5T3 1;510:3 748:9 50:4% 1;507:9 793:6 47:4% 1;452:9 772:9 46:8% 1;507:9 1;509:1 47:5%
15C4F5T4 1;079:8 734:1 32:0% 1;079:8 836:1 22:6% 1;079:8 850:5 21:2% 1;079:8 1;082:7 21:5%
15C4F5T5 1;593:3 809:3 49:2% 1;593:3 906:0 43:1% 1;593:3 933:6 41:4% 1;593:3 1;611:9 42:6%
Avg 1;425:3 789:1 43:9% 1;432:1 881:8 37:5% 1;389:1 890:3 35:0% 1;397:6 1;414:0 36:0%
Table 2: Results comparison between different solution methods
We continue to run more instances to show the efficiency of our improve formulation and adaptive
branch-and-cut algorithm, as well as the ALNS. Table 2 shows the results of larger instances which can not
be optimally solved within 2 hours of CPU time. That is, after 2 hours of CPU time, we record the lower and
upper bounds. In total we run experiments with a setting which has 15 customers, 4 CNG fueling stations
and 5 trucks. We run 5 random instances and each instance is solved by three methods which are basic,
improved, adaptive branch-and-cut and ALNS. The three group of columns (Basic, Imp, B&C) in Table 2
shows the upper and lower bounds as well as the gaps for these three methods. Note that since the adaptive
cut generation has randomness, the columns in the group B&C are the average of 5 replicates. The last
group of column (ALNS) records the best (Best) and average (Avg) solutions within 10 replicates. It also
38
includes the gap (Gap) between the average solution and the best lower bound from the previous three exact
methods. Comparing the three exact methods, our adaptive branch-and-cut algorithm always has the lowest
upper bound and the gap is also the smallest. The ALNS has a good performance in terms of the solution
quality. On average, there is only a 1:8% gap between the average solution of the 10 ALNS replicates and
the best upper bound from the exact methods. In some instances, like 15C4F5T1, the best ALNS solution
may even be better than the best upper bound of the exact methods.
3.4.2 Results on Moderate Instances
In order to explore the effectiveness of the CNG trucks, we conduct experiments on moderate instances. The
moderate instances are modified from the test case A-n33-k5 by Augerat et al. (1995). In the original test
case, there are 32 customers and 1 depot randomly allocated on a 100*100 map. Each customer is assigned
with a random number representing the demand volume. Each truck has identical capacity limitation and
delivers cargo to customers with a tour starting and ending at the depot. The optimal solution for this
problem is known in the literature and in the optimal solution, 5 trucks are used to deliver products to the
customers. We modify the test case by adding the tank capacity and daily traveling distance constraints to
the trucksso thatthey aremore likeCNG trucks. The truck speed
1
, fuel consumptionrate
2
, fueling speed

3
, and fixed fueling time R
1
are set as 35 miles/hour, 3.75 miles/GGE, 3 GGE/minute and 20 minutes,
which are identical with Section 3.4.1. Note that the previous known optimal solution may no longer be
feasible with these added constraints. We also randomly select several points within the map to be CNG
fueling stations. We assume that there is no CNG fueling station at the depot and all trucks in the fleet
need to use the public CNG fueling stations. If we develop a route plan for a single day, all trucks are
likely to return to the depot with nearly empty tanks, which will have an influence on the next day. Thus
we treat the customers’ demands as daily, rather than a single time, and extend the schedule to be 5 days
with the same repeated daily demand. We require that all trucks need to return to the depot at the end
of each day. The heuristic solution method proposed in Section 3.3 can be used with a slight modification.
For the construction of initial feasible solution in Section 3.3.1, an extra requirement of minimal fuel level
is added at the end of a day so that the truck has enough fuel to go to the nearest CNG fueling station in
the next day. Only nodes that belong to the same day can be manipulated by the ALNS steps proposed in
Section 3.3.2, except for the CNG fueling station optimization in the Fuel Station Optimization part, which
is optimized with the nodes for 5 days, since the fueling decision for a truck is made based on its fuel level
within the 5 days. For daily working time T, tank capacity G and load capacity L, we set to 8 hours (480
39
minutes), 53.3 GGE, and 1, respectively. For initial gas level
k
, we generate it by using the same method
as in Section 3.4.1.
jFj
Instance 1 Instance 2 Instance 3 Instance 4 Instance 5
Time Fuel CPU Time Fuel CPU Time Fuel CPU Time Fuel CPU Time Fuel CPU
3 6;742:5 18:7 291:0 6;885:4 19:2 215:5 6;810:2 19:0 207:7 6;691:1 18:4 206:0 6;810:4 18:9 174:4
4 6;676:9 18:4 275:0 6;805:0 19:1 218:5 6;685:5 18:5 202:6 6;577:6 17:8 205:7 6;717:6 17:7 176:5
5 6;642:3 18:3 241:2 6;711:8 18:6 214:8 6;571:8 17:5 200:5 6;541:5 17:6 206:6 6;644:4 17:0 173:3
7 6;553:0 17:6 210:7 6;628:3 17:6 209:8 6;542:2 17:3 197:2 6;526:4 17:4 206:4 6;586:6 16:6 182:2
10 6;464:3 17:0 201:4 6;564:3 17:2 204:0 6;517:4 17:4 202:2 6;467:6 16:9 214:1 6;460:9 16:5 174:1
Table 3: Results for different number of CNG fueling stations
Figure 6: Average total working time with different number of CNG fueling stations
Table 3 presents the detailed results on different number of CNG fueling stations. We run 5 instances and
each instance has CNG fueling stations increasing from 3 to 10, by gradually adding extra randomly located
CNG fueling stations into the map. The instances along the column in Table 3 differ only in the locations
of the CNG fueling stations. ColumnjFj in Table 3 shows the cardinality of the CNG fueling stations set
F, namely the number of CNG fueling stations. For each instance with a certain number of CNG fueling
stations, we run the ALNS 10 times and column Time is the average total working time of these 10 runs.
Since we make a route plan for 5 days, the total working time includes the 5 trucks’ working time in 5 days.
Column Fuel is the average number of times the 5 trucks use a CNG fueling station in 5 days. Note that
the values of this column vary from 16.5 to 19.2, namely on average a truck visits CNG fueling stations less
than one time per day. Though the frequency is low, the fueling decision does influence the total working
time. Column CPU is the average CPU time used to solve the instance. From Table 3 we observe that when
increasing the number of CNG fueling stations, the total working time decreases, as we expect. The number
of times a truck visits a CNG fueling station also decreases. One potential reason is that trucks have more
flexibility on the choice of CNG fueling stations when the number of CNG fueling stations increases and
drivers are more likely to refuel when their tanks are almost empty. Another reason for the reduction of
40
visiting times is that the total traveling distance decreases. Figure 6 presents the average total working time
of these 5 instances. It has a rapid decrease when the number of CNG fueling stations increases from 3 to
4, which means the number of CNG fueling stations is an important factor when the total number is low,
and it also follows the law of diminishing marginal utility. Column CPU is the average CPU time used in
seconds. Onaverageittakesonlyfewminutes. Wecomparethesesolutionswiththeoptimalsolutionwithout
considering the daily working time and tank capacity constraints. A solution without these constraints can
represent a route plan for a fleet of diesel trucks and is a standard solution to a capacitated VRP problem.
The optimal solution for the original instance A-33n-5k can be found in the work from Augerat et al. (1995)
and the corresponding total working time is 5; 665:7 minutes, if the truck speed is 35 miles/hour. Compared
with the optimal solution, the average total working time with 3 CNG fueling stations is 19.8% more and
with 10 CNG fueling stations is 14.6% more.
G
Instance 1 Instance 2 Instance 3 Instance 4 Instance 5
Time Fuel CPU Time Fuel CPU Time Fuel CPU Time Fuel CPU Time Fuel CPU
46:7 6;796:5 21:2 209:4 6;791:7 21:1 196:0 6;668:1 20:2 192:4 6;673:6 20:5 206:1 6;722:2 21:7 169:3
53:3 6;569:0 17:8 211:5 6;613:1 17:6 191:8 6;533:7 17:6 193:9 6;510:6 17:1 207:8 6;582:6 16:9 169:5
60:0 6;493:4 16:3 206:3 6;488:9 14:8 192:1 6;442:3 15:7 191:8 6;436:9 15:0 204:7 6;401:8 14:7 168:7
66:7 6;393:8 13:5 207:0 6;412:9 12:9 189:9 6;362:6 13:6 187:7 6;360:4 12:4 204:1 6;347:9 12:8 167:9
73:3 6;314:5 11:6 206:6 6;352:5 11:1 188:5 6;298:3 11:7 189:7 6;301:7 10:9 203:6 6;273:8 10:7 169:1
Table 4: Results for different fuel tank capacities
Figure 7: Average total working time with difference tank capacities
Next we do a sensitivity analysis on the tank capacity to explore its influence on the route plan. We
modify the instances by fixing the number of CNG fueling stations at 7 and change the tank capacity from
46.7 GGE to 73.3 GGE. All the other parameters are identical with the previous group instances. Again we
run 5 instances and the only difference between these instances is the locations of the CNG fueling stations.
Within each instance we run the ALNS with different tank capacities and for each tank capacity level, we
41
run 10 times to get the average results. Table 4 and Figure 7 present the influence of the tank capacity.
Table 4 is quite similar to Table 3. Figure 7 shows the tank capacity has a larger influence on the route plan
when the tank capacity is low. When the tank capacity is set to 73.3 GGE, the average total working time
of 6; 308:2 is close to the total working time in the optimal solution without the refueling process and the
tank capacity constraints, which is 5; 665:7. Indeed, the average total working time only increased by 11.3%.
Considering the fixed and variable fueling time, this gap is small. The small gap indicates that refueling
is no longer an important factor for CNG trucks under this specific setting, as well as demonstrating the
efficiency of our proposed ALNS.
3.4.3 Results for Data from the Ports of Los Angeles and Long Beach
We run experiments on the estimation of daily goods delivery data for the Ports of Los Angeles and Long
Beach. The data is from the report "SCAG Regional Travel Demand Model and 2008 Model Valida-
tion" LSA Associates (2012), which was conducted by the Southern California Association of Governments
(SCAG) in 2012. Within this model the Southern California region is divided into 4; 109 blocks and each
block is associated with a pair of longitude and latitude, which represents its location. There is an estimation
of light duty diesel truck flow between any pair of blocks. We make some modifications and transfer the
demand model to be a routing problem for CNG trucks. The direct distance between every pair of blocks is
calculated based on their longitudes and latitudes and use it as the distance between them. For the locations
of CNG fueling stations, we collect their locations from the Alternative Fuels Data Center under the U.S.
Department of Energy U.S. Department of Energy (2017a) and locate them within one of the 4; 109 blocks.
In total there are 85 blocks that have publicly accessible CNG fueling stations in the Southern California
region. We select the block where the Ports of Los Angeles and Long Beach are located as the depot for our
fleet and use the light duty diesel trucks’ truckload trip rates between the depot and the rest of the 4; 108
blocks as the daily demands to meet. Within these 4; 108 blocks, only the blocks which can be reached by a
fully fueled CNG truck without having to stop for refueling are considered. A block which is far way from
the port, or far away from the CNG fueling stations are unreachable in our route plan. For each block, we
assume the integral part of the demand is fulfilled by FTL (Full Truckload) and the fractional part is fulfilled
by LTL (Less Than Truckload). Overall, there are 28 truckloads of FTL demand and 73.7 truckloads of LTL
demand allocated in 3; 743 blocks and we only consider the 73.7 truckloads of LTL demand for routing.
42
Figure 8: Average total working time with different tank capacities for light duty trucks
Similar to the experimental settings in Section 3.4.2, we make the route plan for 5 days with the same
repeated daily demand for light duty trucks to simulate the plan for the fleet over a week. ParameterT is set
to 8 hours and all the other parameters are the same as in Sections 3.4.1 and 3.4.2. In total, there are 73.7
truckloads of LTL demand, so at least the fleet size is 74 trucks. In this experiment, we allow the fleet to use
77 trucks at most, which is the smallest fleet size we can get from the ALNS method. We run experiments
with the tank capacity G increasing from 46.7 GGE to 80.0 GGE, eventually to positive infinity to simulate
the situation of diesel trucks where the fueling process is no longer a constraint. For a certain level of tank
capacity, we run 5 independent replicates and take the average of the best solution in each replicate as the
value in the figure. Figure 8 presents the number of refueling stops and the average total working time (in
hours) of the fleet in 5 days when changing the tank capacity from 46.7 GGE to 80.0 GGE. On average,
these 77 trucks need to refuel 243.3 times in 5 days, namely each truck refuels more than 3 times, when the
tank capacity is only 46.7 GGE. When the tank capacity increases to 80.0 GGE, on average, these 77 trucks
need to refuel 106.3 times in 5 days, namely each truck only needs to refuel less than twice. Similar to the
results in Section 3.4.2, the tank capacity has a larger influence on the total working time when it is low.
When the tank capacity increases from 46.7 GGE to 53.3 GGE, the total working time decreases by 67.3
hours, in which the savings from the fixed fueling time, 16.2 hours, make up only 24.1% of the total saving.
In fact, the decrease is mainly from the improvement of routing due to the flexibility from a larger tank
capacity. In our case, the benefits of flexibility still exist even when the tank capacity is large. The average
daily total working time is 1; 396:4 when the tank capacity is as large as 80 GGE. This number is quite close
to the one without tank capacity limitation, which is 1; 351:6, with only a difference of 3.3%. The program
finishes between 90 and 150 CPU minutes in each run, which is fast considering the size of the problem.
43
Figure 9: Average total working time with different sizes of light duty truck fleets
Another interesting thing to explore is the effect of the fleet size. In fact, the objective of the conventional
VRP with shortest traveling distance implies the minimum number of truck, if there are no extra constraints
such as tank capacity limitation. We use the same setting as the previous group of experiments and fix the
tank capacity at 53.3 GGE. Figure 9 shows the effect of the fleet size on the number of refueling stops and
the average total working time (in hours). As the figure shows, the average total working time decreases
when the fleet size increases. However, this reduction stops after the fleet size reaches 87.
44
4 CNG Truck Location Routing Problem
In this chapter, we describe the CTLRP and present a basic MIP formulation. Next we propose some
techniques including variable elimination and valid inequalities, to strengthen the model formulation. We
develop an adaptive branch-and-cut algorithm to further improve the speed. In order to solve problems of
large size, we modify the ALNS in Section 3.3. Computational results, as well as some insights from the
CTLRP are presented in the last part.
4.1 Problem Statement
Different from large CNG truck fleets, it is not cost efficient for companies with small fleets to build a private
CNG fueling station in their own depot. One option is to collaborate to build some shared private CNG
fueling stations allocated within their delivery regions. The CTLRP is designed to solve the location and
routing problem for a group of small CNG truck fleets. In the CTLRP problem, each fleet has some trucks
that need to deliver products from their own depot to their customers. The demands from customers are
known ahead of time, and are independent for each fleet, which means a truck from one fleet can not deliver
products to a customer who belongs to another fleet. These small fleets work together to build some private
CNG fueling stations within some candidate locations so that these fleets can be refueled during their routes.
The problem is to make the location and routing decisions so that the total traveling distance of the
fleets are minimized, while keeping the size of these fleets as small as possible, and satisfying three kinds of
constraints which are load capacity, daily traveling distance and fuel tank capacity. These fleets can only
build a certain number of shared private stations among the given candidate locations. We integrate the
fixed fueling time into our model via reducing the daily traveling distance. This fixed fueling time could be
the waiting time at a station because the station tank has to be refilled if it is recently used, which reduces
the productivity of CNG trucks. Our model reflects this reduction by decreasing the daily traveling distance
limitation every time when a truck goes to a fuel station. In our problem, a typical valid route for a truck is
that a truck starts from its own depot with a given level of gas left in the tank, delivers products to customers
allocated around the depot, refuels at the shared stations if necessary and comes back to the origin depot
with more than a minimal gas level, which allows it to go to the nearest shared station on the next day.
45
4.2 Mixed Integer Programming Model
Suppose there areM CNG truck fleets in total. Each fleet hasN
i
;i = 1;:::;M customers and in total there
are N =
P
i=1;:::;M
N
i
customers. This group of fleets has F candidate locations for CNG fueling stations
and plans to build K stations within these candidates. A network of the problem has a set of vertices V
and a set of arcsA. Vertices setV includes three kinds of subsets, set T representing departure and arrival
depots for the truck fleets, set C representing locations for the customers and set F representing candidate
locations for the CNG fueling stations. Since one CNG fueling station may be visited multiple times by the
same truck, set F includes enough replications for each candidate CNG fueling station so that each vertex
in F will be visited at most one time. Arc set A includes all directed arcs between any two vertices in set
V, except for those arcs which are not needed. These arcs contain arcs from all the other vertices to the
departure depot vertices, from the arrival depot vertices to the other vertices, and arcs between different
fleets. The last kind of arcs are unnecessary because all the fleets are independent and it does not make
any sense for a truck to visit a customer that belongs to another fleet, or travel between the other fleet’s
customers.
We introduce more notation for the MIP formulation. Let T
+
and T

denote vertices that correspond
to the departure and arrival depots, respectively. Let V
+
and V

represent vertices with positive out and
positive in degrees. Namely V
+
=T
+
[C[F and V

=T

[C[F. More specifically, let V
+
i
and V

i
be
the sets of vertices which have arcs from and to vertex v
i
, respectively. SetM is the set of CNG trucks and
M
i
is the collection of CNG trucks which go through vertex v
i
. Set L is the set of candidate locations for
CNG fueling stations andL
f
is location f’s replicated vertices in the graph.
In the MIP model, the binary variable x
k
i;j
is equal to 1 if truck k travels from vertex v
i
to v
j
in the
solution and 0 otherwise. Another binary variable y
f
is an indicator for the CNG fueling station. It is 1
if there is a CNG fueling station at location f and 0 otherwise. Note that in the formulation, variable y
f
can be relaxed as a continuous variable because y
f
can be only 0 or 1 in the optimal solution. Continuous
variablesd
k
i
,t
k
i
andc
k
i
keep track of the traveling distance from the departure depot, traveling distance since
the last fuel station visit and the current delivered load of truck k at vertex v
i
, respectively.
Constant R
1
represents the fixed fueling time measured in distance (i.e., fueling time multiplied by the
average traveling speed of the vehicle) for a truck’s single visit to a fuel station. It means every time a truck
visits a fuel station, its daily traveling distance limitation will decrease by R
1
. M is a big constant number.
Constants D and T are the limitations for daily traveling distance and the distance a truck can travel for a
full CNG fuel tank. ConstantL is the limitation of a truck’s cargo load.
i;j
is the distance from vertexv
i
to
46
vertex v
j
.
i
is the demand for vertex v
i
. Let
k
be the initial fuel level (in miles) for the trucks.
k
should
be a number between max
j2F

i;j
and T, where i is the departure vertex for truck k. Within this range, it
guarantees that the truck has enough fuel to go to any fuel station at the beginning of the day. Let T
k
be
the minimal fuel level (in miles) at the end of the day. T
k
= max
j2F

i;j
by the definition.
47
4.2.1 Basic Formulation
The basic MIP model for the CTLRP is as follows.
min
X
i2V
+
X
j2V
+
i
X
k2Mi

i;j
x
k
i;j
(64)
s.t.
X
j2V

i
x
k
j;i

X
j2V
+
i
x
k
i;j
= 0; 8i2C[F;k2M
i
(65)
X
j2V
+
i
X
k2Mi
x
k
i;j
= 1; 8i2C (66)
X
j2V
+
i
x
k
i;j
1; 8i2T
+
;k2M
i
(67)
X
j2V
+
i
x
k
i;j
1; 8i2F;k2M
i
(68)
d
k
i
+M

x
k
i;j
1

+
i;j
d
k
j
; 8i2T
+
[C;j2V
+
i
;k2M
i
(69)
d
k
i
+M

x
k
i;j
1

+
i;j
+R
1
d
k
j
; 8i2F;j2V
+
i
;k2M
i
(70)
0d
k
i
D; 8i2V;k2M
i
(71)
t
k
i
+M

1x
k
i;j

i;j
t
k
j
; 8i2T
+
[C;j2V
+
i
;k2M
i
(72)
M

1x
k
i;j

+T
i;j
t
k
j
; 8i2F;j2V
+
i
;k2M
i
(73)
0t
k
i
T; 8i2V;k2M
i
(74)

k
=t
k
i
; 8i2T
+
;k2M
i
(75)
t
k
i
T
k
8i2T

;k2M
i
(76)
x
k
i;j
y
f
; 8f2L;i2L
f
;j2V
+
i
;k2M
i
(77)
X
f2L
y
f
K (78)
c
k
i
+M

x
k
i;j
1

+
i
c
k
j
; 8i2V
+
;j2V
+
i
;k2M
i
(79)
0c
k
i
L (80)
x
k
i;j
2f0; 1g; y
f
0 (81)
Objective function (64) is to minimize the total traveling distance. Constraints (65) to (68) are the
flow conservation equations that ensure that the solution defines a routing plan for the fleets to visit every
customer. Constraints(69) and (71) are distanceconstraints. Constraint (69)increases the traveling distance
48
d
k
j
by
i;j
if the truck k travels from v
i
tov
j
. An extra distance penalty R
1
will be added in constraint (70)
if the truck goes through a CNG fueling station vertex. The daily total traveling distance for each truck is
bounded by D in constraint (71). The fuel tank capacity constraint is covered by constraints (72) to (76).
Similar to constraint (69), constraint (72) calculates the remaining fuel (in miles). Constraint (73) resets
the level if it visits a fuel station, which means the truck is fully refueled. Constraint (74) ensures that a
truck does not run out of gas. The initial and final gas levels are set in constraint (75) and constraint (76).
In constraint (77), indicator variable y
f
will be 1 if any replication of candidate CNG fueling station f is
visited by a truck, namely a CNG fueling station will be built at f. By the triangle inequality, in a optimal
solution, a truck will visit a candidate CNG fueling station only if the station is built there. Constraint (78)
limits the total number of stations built to be no more than K. Constraints (79) and (80) guarantee the
demand for each customer can be fulfilled by the assigned truck. Constraint (81) is the domain constraint.
4.2.2 Preprocessing and Valid Inequalities
Before solving the MIP formulation, some variables are eliminated based on the feasibility of the arcs, and
parameters like big M are selected; and valid inequalities which are redundant but can strengthen the
formulation are added to the model.
Variable elimination. Variables are eliminated based on different feasibility rules. All the rules are
similar to the ones in Section 3.2.2, except for the third rule.
The third elimination rule is induced by the domination of the fuel stations. We say a CNG fueling
stationf
1
dominates anotherf
2
between customer verticesv
i
andv
j
, if equations (82) and (83) are satisfied.

i;f1

i;f2
(82)

f1;j

f2;j
(83)
A fuel station f
2
can never be an option between v
i
and v
j
in the optimal solution if it is dominated by
more than FK fuel stations. If fuel station f
2
is dominated by any vertex successive to v
i
, all variables
x
k
i;f2
= 0; 8k2M
i
.
Parameter setting Parameters like bigM and the number of replicates of each fuel station are explicitly
handled. In theory the big M in constraints (69), (70), (72), (73) and (79) can be any number as long as
it is sufficiently large. However, an arbitrary big M may bring in some numerical problems and weaken the
bound from the linear relaxation of the MIP. We need to choose M as small as possible to accelerate the
49
computation speed. The minimal M for constraints (69)and (70) is D
j;n
, for constraints (72) and (73) is
T and for constraint (79) is 1.
We replicate the fuel stations in case a truck needs to visit a fuel station multiple times. The upper
bound of times for a truck to visit a fuel station is given byd
2D
=
T
e + 1. The proof is given in Section 3.2.2.
Valid inequalities Some valid inequalities are added to the model to strengthen the linear relaxation
of the problem. The valid inequalities we add are from the tank capacity limitation, which is similar to
inequality (33) in Section 3.2.2.
4.2.3 Adaptive Cut Generation
The branch-and-cut algorithm is a popular algorithm for solving combinatorial optimization problems. We
generate the initial pool of cuts based on load capacity limitation, and the domination relationship in the
fueling process as well as the regular subtour elimination from Cordeau (2006). The cuts are identical
with the ones in Section 3.2.3. Similarly, we use the adaptive cut generation method to solve the time in
searching for the cuts. We let p
i;t
be the probability of checking group i’s cuts in iteration t. Initially
p
i;0
= 1; 8i = 1; 2; 3. We update p
i;t
based on equation (39) in Section 3.2.3.
50
4.3 Heuristic Solution
Inspired by the hybrid heuristic method used in Section 3.3, we use a similar method to solve the CTLRP.
Figure 10: Steps of the heuristic method for the CTLRP
Figure 10 shows the main steps of the heuristic method for the CTLRP. In the first step, we formulate a p-
median problem to solve the initial locations for the private CNG fueling stations. In the p-median problem,
the demand of a customer is treated as the weight for that customer, and the Euclidean distance of any two
locations is used as the distance of the two locations. The whole problem can be formulated as a standard
p-median problem, which is pick up p locations for the fuel stations out of some candidate locations. The p-
median problem is solved by the heuristic method from Daskin (2011). In the second step, we independently
construct the initial feasible route plan for each fleet based on the location decision of the fuel stations in
the first step. For each truck, the construction method is identical with the one used in Section 3.3.1. In
the third step, the initial feasible route plans are improved by an ALNS method, which is almost the same
as the one used in Section 3.3.2. In the ALNS method, we first assume the locations of the fuel stations are
given and each fleet’s route plan is independently improved by the three parts from Section 3.3.2, which are
local search, destroy and repair, and fleet size minimization. Different from Section 3.3.2, the fuel station
optimization part here does not only optimize the choice of fuel stations, but also optimizes the location of
the fuel stations. We use a new MIP formulation to solve the fuel station optimization problem.
Similar to the fuel station optimization in Section 3.3.2, we start from route r
k
. First, delete all of the
existing fuel stations in r
k
and leave the nodes with the customer and depot vertices only. Let r
k
0
be the
route after deletion and suppose there are n
k
0 nodes left. In the heuristic method, for any two successive
nodes in r
k
0
, we only consider the fuel station candidates with the smallest detour. Let binary decision
51
variablez
k
0
(i)
denote the fueling decision for router
k
0
after theith node, which is the customer noder
k
0
(i)
. Let
z
k
0
(i)
be 1 if the corresponding truck goes to fuel station after r
k
0
(i)
and be 0 otherwise. Let Z
f
be the set of
allz
k
0
(i)
whose corresponding fuel stations are candidate f. For simplification, we usez
k
i
to representr
k
0
(i)
, use
i to represent customer node r
k
0
(i)
, use d
k
i
and t
k
i
to represent the accumulative traveling distance and gas
level of route r
k
0
at the ith node. Based on the formulation in Section 4.2.1, we use the following notation.
Decision variabley
f
is a binary variable. It is 1 if a fuel station is built at candidate location f, otherwise it
is 0. With the same reason in Section 4.2.1, it can be relaxed as a continuous variable. Let L be the set of
candidate locations,M be the set of all trucks in all fleets.
i;j
is the distance between i and j.
+
i;j
and

i;j
are the distances from i to the fuel station, and from the fuel station to node j, respectively. D;T;K are
constants for the daily traveling distance limitation, tank capacity, and the number of fuel stations to build.
The following MIP model is for the fuel station optimization problem. Note that different from the MIP
model in Section 3.3.2 in which each truck’s route is independent and solved separately, the MIP model
here solves the problem for all trucks’ routes together, since the choice of the locations of the fuel stations
depends on all the truck routes.
min
X
k2M
n
k
1
X
i=1

+
i;i+1
+

i;i+1

i;i+1
+R
1

z
k
i
(84)
s.t. d
k
i
+
i;i+1
+

+
i;i+1
+

i;i+1

i;i+1
+R
1

z
k
i
d
k
i+1
; 8k2M;i = 1;:::;n
k
1 (85)
0d
k
i
D; 8k2M;i = 1;:::;n
k
(86)
t
k
i
+ (
i;i+1
+T

i;i+1
)z
k
i;i+1

i;i+1
t
k
i+1
; 8k2M;i = 1;:::n
k
1 (87)
t
k
i

+
i;i+1
z
k
i;i+1
0; 8k2M;i = 1;:::n
k
1 (88)
0t
k
i
T; 8k2M;i = 1;:::n
k
1 (89)

k
=t
k
0
; 8k2M (90)
t
k
i
T
k
8k2M (91)
z
k
i
y
f
; 8z
k
i
2Z
f
;f2L (92)
X
f2L
y
f
K (93)
z
k
i
2f0; 1g; y
f
0 (94)
The objective function (84) is to minimize the total detour for all trucks due to the refueling process.
52
Constraints (85) and (86) are distance constraints. In constraint (85), the traveling distance increases
by
i;i+1
if the truck only does not visit a fuel station, and increases by
+
i;i+1
+

i;i+1
+R
1
if the truck
makes a detour for a fuel station, where R
1
is the cost in terms of distance due to the fixed fueling time.
Constraint (86) sets the boundary for the daily traveling distance. Constraints (87) to (91) are the tank
capacity constraints. In constraint (87), the remaining gas level t
k
i
decreases by
k
i;i+1
if the truck does not
go to a fuel station, and increases to T

i;i+1
otherwise. Constraint (88) guarantees the truck has enough
gas if it goes to a fuel station. Constraints (89), (90) and (91) set the boundary, initial and final value for
gas level t
k
i
. Constraints (92) ensure that a fuel station can be used only if it is built. Constraint (93) sets
the upper bound of the fuel stations to be built.
In practice, we can easily construct an initial feasible solution by checking the original route plans before
the deletion of fuel stations. We also set a upper bound of the CPU time spent for each MIP formulation.
Overall, the MIP can be solved quickly for reasonable sized problems (See Section 4.4.2).
53
4.4 Computational Results
In this section, we randomly generate cases to test the performance of our algorithms, and explore the
influences of different factors like tank capacity, number and location of fuel stations, on our objective. All
of the test cases are generated on a 100*100 map. Locations for depots, customers and candidate sites for the
fuel stations are randomly picked within the map. The demands of customers are random numbers between
0 and 1.
All experiments are performed on a PC with Intel Core i7-4790K CPU and 32 GB RAM, running Linux
distribution Ubuntu 14.04 and CPLEX 12.5.
4.4.1 Results for Small Instances
For small instances, the demands of customers are generated from a uniform distribution between 0 to 1.
Small instances which can be solved within 1 hour of CPU time are tested. In order to show the effectiveness
of our valid inequalities and adaptive branch-and-cut algorithm, three kinds of formulations basic, improved
and adaptive branch-and-cut, are tested. The basic formulation includes the base MIP formulation and big
M selection. The improved formulation includes everything presented in Section 4.2.2. The branch-and-cut
cases include the improved formulation and the adaptive cut generation method. Since the cut generation
part in our adaptive branch-and-cut algorithm has randomness, we run the algorithm 5 times and record
the average results in all of the following tables. The experiments are conducted based on the background
of CNG trucks and parameters are modified from related materials (Laughlin and Burnham, 2016; U.S.
Department of Energy, 2014). We set the daily traveling distance limitation D to 250 and tank capacity T
to 225. The parameter for the fixed fueling time, R
1
, is set to 20. Namely every time when a truck goes to
a fuel station, its daily working distance decreases by 20. The load capacity L for each truck is 1. In our
small instances, we only consider cases with 2 fleets. The initial fuel level
k
is a random number between 0
and max
i2T
+ (T max
j2F

i;j
).
54
Name
CPU Time Solution Value
Opt Basic Opt Imp Opt B&C Opt # Opt Obj Gap
2T5C1F3L1 73:7 31:7 6:4 711:9 10 711:9 0:0%
2T5C1F3L2 13:4 12:8 7:4 645:9 8 646:0 0:0%
2T5C1F3L3 29:2 25:5 10:7 693:2 8 695:3 0:3%
2T5C1F3L4 54:7 30:0 9:9 442:2 9 443:3 0:2%
2T5C1F3L5 115:8 56:9 25:2 613:5 10 613:5 0:0%
2T5C1F3L6 210:5 66:8 29:9 688:3 10 688:3 0:0%
2T5C1F3L7 29:2 30:0 8:6 592:6 10 592:6 0:0%
2T5C1F3L8 31:6 17:6 6:9 656:3 10 656:3 0:0%
2T5C1F3L9 146:9 102:2 33:5 683:0 8 684:3 0:2%
2T5C1F3L10 33:3 24:7 6:7 597:7 10 597:7 0:0%
2T5C1F3L11 62:0 42:3 32:8 734:4 10 734:4 0:0%
2T5C1F3L12 74:9 79:7 40:0 605:5 10 605:5 0:0%
2T5C1F3L13 31:9 23:0 9:5 660:1 4 679:2 2:9%
2T5C1F3L14 22:8 17:8 9:9 516:7 10 516:7 0:0%
2T5C1F3L15 26:7 18:6 6:3 619:9 10 619:9 0:0%
2T5C1F3L16 65:8 27:4 17:5 724:7 6 727:1 0:3%
2T5C1F3L17 30:5 19:5 6:8 694:9 8 701:0 0:9%
2T5C1F3L18 190:4 27:7 29:9 565:4 10 565:4 0:0%
2T5C1F3L19 54:4 52:3 13:1 691:9 8 700:4 1:2%
2T5C1F3L20 99:2 39:6 40:8 618:9 10 618:9 0:0%
Avg 69:8 37:3 17:6 637:9 9:0 639:9 0:3%
2T5C1F5L1 23:2 19:6 7:1 652:3 10 652:3 0:0%
2T5C1F5L2 678:5 663:8 160:7 652:8 10 652:8 0:0%
2T5C1F5L3 31:8 62:3 6:2 643:6 7 653:6 1:6%
2T5C1F5L4 1;924:1 171:4 174:3 635:7 10 635:7 0:0%
2T5C1F5L5 19:6 17:8 6:7 666:1 10 666:1 0:0%
2T5C1F5L6 53:4 25:7 11:8 562:4 10 562:4 0:0%
2T5C1F5L7 100:5 19:7 7:7 644:2 10 644:2 0:0%
2T5C1F5L8 27:1 23:8 10:1 511:2 10 511:2 0:0%
2T5C1F5L9 65:5 47:9 28:4 473:6 10 473:6 0:0%
2T5C1F5L10 13:2 17:7 4:5 428:5 10 428:5 0:0%
2T5C1F5L11 52:1 41:3 18:7 592:4 10 592:4 0:0%
2T5C1F5L12 37:0 43:6 24:0 546:5 7 560:9 2:6%
2T5C1F5L13 23:9 22:2 11:5 662:8 10 662:8 0:0%
2T5C1F5L14 27:6 26:7 15:6 674:3 4 682:0 1:1%
2T5C1F5L15 39:3 31:7 14:4 584:9 8 584:9 0:0%
2T5C1F5L16 251:6 235:3 152:6 637:7 7 648:5 1:7%
2T5C1F5L17 13:3 11:8 3:5 608:2 7 625:8 2:9%
2T5C1F5L18 197:1 193:6 73:9 624:0 8 624:0 0:0%
2T5C1F5L19 106:1 103:9 33:9 657:1 8 657:1 0:0%
2T5C1F5L20 37:4 32:1 24:0 522:6 8 522:6 0:0%
Avg 186:1 90:6 39:5 599:0 8:7 602:3 0:5%
Table 5: Results on small CTLRP instances
Table 5 presents the CPU time used to solve the small instances. In total we conduct experiments with 2
different settings and each setting has 20 random instances. Name 2T5C1F3L1 means it is the first instance
with a setting of 2 fleets, including 2 trucks and 5 customers for each fleet, and 1 CNG fueling station to
be built within 3 candidate locations. The first group of columns, CPU Time records the CPU Time (in
seconds) used for the three formulations. Column Opt Basic, Opt Imp and Opt B&C are the CPU time (in
seconds) for the basic, strength and adaptive branch-and-cut formulations, respectively. Since the adaptive
branch-and-cut algorithm includes randomness in the cut generation part, we run 5 replicates and record
55
the average CPU time in the table. Within these three formulations, the branch-and-cut is the fastest one
in most of the instances, while the basic formulation is usually the slowest one. For example, in the setting
of 2T5C1F3L, the branch-and-cut formulation is the fastest method in 18 out of the total 20 instances. In
contrast, the basic formulation is the slowest one in 18 instances. One average, the strength formulation
needs 37.3 CPU seconds to solve an instance, which is 46.6% faster than the basic formulation in terms of
the CPU time. The branch-and-cut formulation only needs 17.6 CPU seconds to solve an instance, which
is 74.8% faster than the basic formulation. The difference is even bigger when the size of the problem gets
larger. In the setting of 2T5C2F5L, on average, the strength formulation is 51.3% faster than the basic
formulation. The branch-and-cut formulation is the fastest method in 19 out of the 20 instances and it is
78.8% faster than the basic formulation. These small instances show the effectiveness of our valid inequalities
and adaptive branch-and-cut algorithm. The performance is stable in solving small instances. The second
group of columns, Solution Value records the solution values from the MIP model and the ALNS method.
Column Opt is the optimal solution from the MIP model. For the ALNS, we run 10 times for each instance.
Column # Opt is the number of times the ALNS obtained the optimal solution. Column Obj is the average
of the 10 solution values. Column Gap is the gap between the average value from the ALNS, which is column
Obj, and the optimal solution value, which is column Opt. The comparison of the solution values shows the
ALNS method has high solution quality in solving small instances. It can identify the optimal solution in
most of the runs. The gap between the average value and optimal solution is also very small. On average
the gap is only 0.3% and 0.5% for these two settings.
56
Name
Basic Imp B&C ALNS
Time Upper Lower Gap Time Upper Lower Gap Time Upper Lower Gap Avg Gap
3T7C1F5L1 3;600:0 873:1 832:0 4:7% 3;600:0 873:1 814:6 6:7% 1;140:4 873:1 873:0 0:0% 922:4 5:4%
3T7C1F5L2 3;600:0 750:0 534:2 28:8% 3;600:0 750:0 598:2 20:2% 3;600:0 750:0 709:5 5:4% 750:0 5:4%
3T7C1F5L3 3;600:0 845:8 699:2 17:3% 3;122:3 845:8 845:7 0:0% 1;303:3 845:8 845:7 0:0% 845:8 0:0%
3T7C1F5L4 3;600:0 878:3 631:6 28:1% 2;813:2 878:3 878:2 0:0% 1;347:4 878:3 878:2 0:0% 906:7 3:1%
3T7C1F5L5 3;600:0 774:5 610:3 21:2% 3;600:0 774:5 572:7 26:1% 3;600:0 774:5 702:2 9:3% 777:6 9:7%
Avg 3;600:0 824:3 661:5 20:0% 3;347:1 824:3 741:9 10:6% 2;198:2 824:3 801:7 3:0% 840:5 4:7%
3T7C2F5L1 3;600:0 784:7 608:5 22:5% 3;600:0 784:7 632:3 19:4% 3;600:0 784:7 667:2 15:0% 786:4 15:2%
3T7C2F5L2 3;600:0 861:2 535:1 37:9% 3;600:0 861:2 715:6 16:9% 3;600:0 861:2 786:6 8:7% 861:2 8:7%
3T7C2F5L3 3;600:0 951:1 718:7 24:4% 3;600:0 951:1 763:2 19:8% 3;600:0 951:1 753:3 20:8% 952:8 19:9%
3T7C2F5L4 3;600:0 943:1 670:5 28:9% 3;600:0 943:1 741:5 21:4% 3;600:0 943:1 740:4 21:5% 955:5 22:4%
3T7C2F5L5 3;600:0 749:9 605:7 19:2% 3;600:0 749:9 650:9 13:2% 3;600:0 749:9 669:6 10:7% 749:9 10:7%
Avg 3;600:0 858:0 627:7 26:6% 3;600:0 858:0 700:7 18:1% 3;600:0 858:0 723:4 15:3% 861:2 15:4%
4T9C1F5L1 3;600:0 1;310:0 725:3 44:6% 3;600:0 1;310:0 866:1 33:9% 3;600:0 1;310:0 1;020:5 22:1% 1;320:5 22:7%
4T9C1F5L2 3;600:0 1;042:5 711:3 31:8% 3;600:0 1;042:5 836:8 19:7% 3;600:0 1;042:5 851:7 18:3% 1;048:0 18:7%
4T9C1F5L3 3;600:0 1;198:2 726:4 39:4% 3;600:0 1;203:3 795:3 33:9% 3;600:0 1;198:2 914:4 23:7% 1;210:7 24:5%
4T9C1F5L4 3;600:0 973:4 757:5 22:2% 3;600:0 973:4 810:8 16:7% 3;600:0 973:4 835:1 14:2% 979:4 14:7%
4T9C1F5L5 3;600:0 844:9 681:5 19:3% 3;600:0 855:5 648:7 24:2% 3;600:0 844:9 679:9 19:5% 853:7 20:2%
Avg 3;600:0 1;073:8 720:4 31:5% 3;600:0 1;076:9 791:5 25:7% 3;600:0 1;073:8 860:3 19:6% 1;082:5 20:2%
4T9C2F5L1 3;600:0 1;084:8 657:7 39:4% 3;600:0 1;084:8 643:7 40:7% 3;600:0 1;084:8 793:2 26:9% 1;153:7 31:2%
4T9C2F5L2 3;600:0 1;208:1 651:2 46:1% 3;600:0 1;180:0 832:2 29:5% 3;600:0 1;180:0 853:6 27:7% 1;195:2 28:6%
4T9C2F5L3 3;600:0 1;094:7 696:0 36:4% 3;600:0 1;094:7 677:1 38:2% 3;600:0 1;094:7 691:5 36:8% 1;098:5 36:6%
4T9C2F5L4 3;600:0 1;332:1 707:5 46:9% 3;600:0 1;332:1 809:3 39:2% 3;600:0 1;332:1 794:9 40:3% 1;332:4 39:3%
4T9C2F5L5 3;600:0 1;059:7 759:2 28:4% 3;600:0 1;059:7 800:7 24:4% 3;600:0 1;056:8 847:5 19:8% 1;082:9 21:7%
Avg 3;600:0 1;155:9 694:3 39:4% 3;600:0 1;150:2 752:6 34:4% 3;600:0 1;149:7 796:2 30:3% 1;172:5 31:5%
Table 6: Results comparison between different solution methods
Table 6 shows the results in solving larger instances. This group of experiments have 4 different settings
and each setting has 5 random instances. Each instance is solved by three formulations, which are identical
withthepreviousexperiments, andtheALNSmethod. Thecolumnsinthetablearedividedintofourgroups,
which are Basic, Imp, B&C, and ALNS. The Basic, Imp and B&C groups, including CPU time (in seconds),
upper bound, lower bound and gap, record the results from the basic, improved and adaptive branch-and-cut
formulations. The ALNS group records the average value of the 10 runs of the ALNS method, and the gap
between the average value and the best lower bound from the previous three formulations. The results show
that the strength and branch-and-cut formulations can solve some instances which the basic formulation can
not solve within one CPU hour. For the instances which none of these methods can solve within one CPU
hour, the average gaps for the strength and branch-and-cut formulations are smaller than those for the basic
formulation. These instances show that our valid inequalities and adaptive branch-and-cut algorithm have
better performance in solving larger instances. We also observe that the average solution values of the ALNS
method are close to the upper bounds from the three formulations, which means our ALNS method can still
obtain high quality solutions.
57
Name
Adp B&C Std B&C
Time Upper Lower Gap Cuts # Cuts Time Time Upper Lower Gap Cuts # Cuts Time
3T7C1F5L1 1;140:4 873:1 873:0 0:0% 10:4 1:5 1;411:2 873:1 873:0 0:0% 11:0 21:5
3T7C1F5L2 3;600:0 750:0 709:5 5:4% 18:8 20:5 3;600:0 750:0 669:8 10:7% 22:0 365:9
3T7C1F5L3 1;303:3 845:8 845:7 0:0% 8:2 2:0 1;650:9 845:8 845:7 0:0% 9:0 48:2
3T7C1F5L4 1;347:4 878:3 878:2 0:0% 14:4 3:8 1;552:0 878:3 878:2 0:0% 18:0 57:4
3T7C1F5L5 3;600:0 774:5 702:2 9:3% 11:8 11:7 3;600:0 774:5 671:3 13:3% 14:0 210:9
Avg 2;198:2 824:3 801:7 3:0% 12:7 7:9 2;362:8 824:3 787:6 4:8% 14:8 140:8
3T7C2F5L1 3;600:0 784:7 667:2 15:0% 12:8 25:2 3;600:2 784:7 651:1 17:0% 14:0 390:4
3T7C2F5L2 3;600:0 861:2 786:6 8:7% 7:6 20:9 3;600:0 861:2 782:8 9:1% 9:0 313:0
3T7C2F5L3 3;600:0 951:1 753:3 20:8% 12:2 5:9 3;601:3 951:1 757:7 20:3% 13:0 58:0
3T7C2F5L4 3;600:0 943:1 740:4 21:5% 9:8 11:0 3;600:7 943:1 742:3 21:3% 10:0 132:7
3T7C2F5L5 3;600:0 749:9 669:6 10:7% 11:0 7:9 3;600:2 749:9 667:5 11:0% 12:0 100:6
Avg 3;600:0 858:0 723:4 15:3% 10:7 14:2 3;600:5 858:0 720:3 15:7% 11:6 198:9
4T9C1F5L1 3;600:0 1;310:0 1;020:5 22:1% 20:2 10:1 3;600:6 1;310:0 996:3 23:9% 22:0 87:8
4T9C1F5L2 3;600:0 1;042:5 851:7 18:3% 29:4 43:7 3;600:2 1;042:5 795:2 23:7% 29:0 656:5
4T9C1F5L3 3;600:0 1;198:2 914:4 23:7% 20:4 45:3 3;600:2 1;203:3 889:6 26:1% 23:0 448:6
4T9C1F5L4 3;600:0 973:4 835:1 14:2% 21:2 11:9 3;600:1 973:4 836:1 14:1% 23:0 131:5
4T9C1F5L5 3;600:0 844:9 679:9 19:5% 19:0 30:8 3;600:1 844:9 678:8 19:7% 22:0 510:7
Avg 3;600:0 1;073:8 860:3 19:6% 22:0 28:4 3;600:2 1;074:8 839:2 21:5% 23:8 367:0
4T9C2F5L1 3;600:0 1;084:8 793:2 26:9% 21:8 53:9 3;600:2 1;084:8 769:2 29:1% 26:0 854:5
4T9C2F5L2 3;600:0 1;180:0 853:6 27:7% 11:6 15:0 3;600:0 1;180:0 817:0 30:8% 15:0 212:5
4T9C2F5L3 3;600:0 1;094:7 691:5 36:8% 15:2 12:2 3;601:3 1;106:8 707:8 36:0% 16:0 168:2
4T9C2F5L4 3;600:0 1;332:1 794:9 40:3% 11:4 18:3 3;600:2 1;332:1 769:5 42:2% 12:0 311:9
4T9C2F5L5 3;600:0 1;056:8 847:5 19:8% 15:4 27:4 3;600:4 1;059:7 852:2 19:6% 22:0 487:6
Avg 3;600:0 1;149:7 796:2 30:3% 15:1 25:4 3;600:4 1;152:7 783:2 31:5% 18:2 406:9
Table 7: Comparison of different branch-and-cut algorithms
Table 7 shows the effectiveness of our adaptive cut generation method by comparing with the solutions
from the branch-and-cut algorithm with standard cut generation. Columns in Adp B&C are the results for
the adaptive branch-and-cut algorithm, which include total CPU time, upper and lower bounds, gap, total
number of cuts added and CPU time (in seconds) used in searching the cut pool. Columns in Std B&C
are the results for the standard branch-and-cut algorithm, which are similar to the columns in Adp B&C.
Compared with the standard branch-and-cut algorithm, the adaptive branch-and-cut algorithm does not
check all of the cuts in the cut pool at every branch node. Instead, the algorithm checks a group of cuts with
higher probability if this group is more likely to generate invalid cuts. As expected, the adaptive branch-
and-cut algorithm saves time in checking the cut pool. The risk is it may miss important cuts. Results
in Table 7 show that compared with the standard branch-and-cut algorithm, the adaptive branch-and-cut
algorithm spends less than 10.0% of the time in checking the cut pool, while it can generate about 90.0%
of the number of cuts generated in the standard branch-and-cut algorithm. For instance, in the setting
3T7C2F5L, on average the adaptive branch-and-cut algorithm generates 10.7 cuts for each instance within
14.2 CPU seconds. The standard branch-and-cut algorithm needs 198.9 CPU seconds to generate 11.6 cuts.
Eventually the adaptive branch-and-cut algorithm can save CPU time to solve problems. For instance,
instance 3T7C1F5L1 needs 1,140.4 CPU seconds for the adaptive branch-and-cut algorithm, and 1,411.2
58
CPU seconds for the standard one. These results show that compared with the standard branch-and-cut
algorithm, our adaptive branch-and-cut algorithm can solve problems faster.
4.4.2 Results for Large Instances
In this section, we continue to increase the size of the problems and use the ALNS method to solve the
problems. In order to reduce the effect of the initial and final gas level, we extend the schedule to be 5 days
by repeating the customers’ demands in each day. We require that all trucks need to be back to the depot at
the end of each day. In the first experiment of this section, we show the performance of our ALNS method
by checking the convergence speed of the method.
Figure 11: Convergence of the ALNS method
Figure 12: Locations of fuel stations
59
The instance we test here has 2 fleets, each fleet includes 10 trucks, 28 customers and 7 fuel stations
need to be built within 15 candidates. Since the schedule is for 5 days, in total the ALNS method needs
to solve a routing problem for 20 trucks with about 300 customers. Figure 11 shows the convergence of the
ALNS method, which is finished in about 1 hour of CPU time. The line with squares is the total distance
with the solution from the p-median problem. We use the locations from the p-median problem in the first
step and solve the routing problem by the ALNS method without changing the locations. Therefore, in the
method only the route changes. The line with circles is the total distance solved by our heuristic method.
In this method, both the location and route change in the iterative process. Our heuristic method’s location
solution is better than the p-median’s solution in terms of total distance after the same number of iterations.
The locations of the chosen fuel stations are presented in Figure 12. The squares and circles are the locations
from the p-median and our heuristic methods, respectively. In this figure, only two out of the seven locations
are the same. This case shows the benefit of solving the location and routing problems simultaneously.
(a) Cases with different tank capacities (b) Cases with different number of fuel stations
Figure 13: Sensitivity analysis of the CTLRP
Figure 13 presents a sensitivity analysis of the CTLRP. Since the heuristic method includes randomness,
we run the solution method 5 times for each instance and record the average value. In Figure 13a, we test
an instance of the same size as the previous experiment, and only change the tank capacity from 200 to
275. Figure 13b uses the same instance, but only varies the number of fuel stations to be built from 6 to 9.
The total distance decreases as the tank capacity and number of fuel stations increase. However, the gaps
are smaller than the gap in Figure 11 between different location solutions, which means under the current
settings the locations of fuel stations have larger influence on the total distance for the fleets.
60
5 Hybrid Transit System
In this chapter, we introduce a hybrid transit system, which combines a ride sharing system with a fixed
route transit system. We present a MIP model to solve the routing problem for the shard vehicles optimally.
Since the problem is computational hard to solve, we propose a heuristic method to solve problems of large
size. Computational results are given in the last part of this chapter.
5.1 Problem Statement
In our system, participants are divided into two groups, vehicle offers and customers. Vehicle offers are
drivers who have their own cars and offer ridesharing service. Customers are passengers who do not have
cars and need to be delivered by vehicle offers. Each vehicle offer has an origin and a destination with earliest
departure time and latest return time. Each vehicle has a capacity constraint. Each customer request has
both a pickup and delivery location with time windows. If the offered vehicles are not enough to serve the
customers, other vehicles, such as taxis, will be used for serving the requests. Like offered vehicles, each
taxi also has an origin and destination, time limitation and capacity constraint. In addition, we have a
fixed-route public transit system, like buses, in the model. Buses or subways have fixed routes and travel
through stations on routes based on a regular timetable. Each customer can be delivered by vehicles or taxis
only, or vehicles combined with buses. Currently we only deal with a static scheduling problem, namely all
of the information, including vehicle offers, customer requests and bus timetables, are available when the
schedule is created.
5.2 Mixed Integer Programming Model
Assume there aren
c
customers in total and each customer has a request containing a pickup node, a delivery
node, and a time window for each node. We create a set N
c
with 2n
c
vertices to represent these customers’
pickup and delivery nodes, and number these vertices from 1 to 2n
c
. Let vertex i, 1in
c
, represent the
pickup node for customeri, then vertexn
c
+i is the corresponding delivery node. Each vertex in setN
c
has
two numbers, tc
b
i
and tc
e
i
, which are the earliest and latest access times for the vertex.
In this system, we have n
s
vehicles to be routed and scheduled. We create a set N
s
with n
s
+ 1 vertices
to represent the departure and return nodes for these vehicles. We number these vertices from 2n
c
+ 1 to
2n
c
+n
s
+ 1 and let vertex 2n
c
+j and 2n
c
+j + 1 be the departure and return nodes for vehicle j. Note
that under such a modeling method, a vertex 2n
c
+j + 1, denotes both the departure node for vehicle j + 1
61
and return node for vehicle j. In reality, these two nodes may have different locations and can be differed
by the traveling distances of arcs adjacent to vertex 2n
c
+j in our model. Letts
d
2nc+j
andts
a
2nc+j+1
denote
the earliest departure time and latest return time for vehicle j.
For the bus nodes, they are more complicated than the previous two kinds of nodes. In the real world,
everyday a bus stop may be visited by multiple buses many times. Also a bus stop may be visited by
multiple vehicles and can be pickup and drop-off nodes at the same time. In our model, we first construct
all feasible bus routes. For instance, a bus going through bus stop A, B, and C at 10:00 a.m., 10:30 a.m.,
and 11:00 a.m. can be represented by three feasible routes. Route from A to B, B to C and A to C. If
there is another bus that goes through A, B and C at different times, it can be represented by another three
routes. Suppose there are n
b
feasible bus routes in total. Since each route may be used by every customer,
we make n
c
replicates for these n
b
routes. Each route has a start node and an end node, so we create a
set N
b
with 2n
c
n
b
vertices to represent these bus routes and number these vertices from 2n
c
+n
s
+ 2 to
2n
c
+n
s
+ 1 + 2n
c
n
b
. In this model, the kth bus route’s start and end nodes for customer i are the vertices
numbered as 2n
c
+n
s
+ 1 + 2(i 1)n
b
+k and 2n
c
+n
s
+ 1 + 2(i 1)n
b
+n
b
+k. Any vertex numbered as
l in setN
b
has a time tb
l
, which is the time when the bus reaches the stop.
There are three kinds of arcs in this map, which are vehicle arcs and bus arcs and artificial arcs. Vehicle
arcs can be divided into six categories based on the properties of start and end vertices. The first is a vehicle
arc between customer vertices in set N
c
. The second is a vehicle arc from a vehicle departure vertex to
its corresponding return vertex. The third is a vehicle arc between bus vertices, which represents vehicle
traveling from one bus node to another. The following three kinds of vehicle arcs are between customer,
vehicle and bus vertices. A vehicle arc between customer and vehicle vertices represents a vehicle traveling
from its departure vertex to a customer pickup vertex, or from a customer delivery vertex to a vehicle return
vertex. A vehicle arc between customer and bus vertices represents a customer traveling between his or her
pickup/delivery node and a bus node. A vehicle arc between vehicle and bus vertices represents a vehicle
traveling from its departure vertex to an end bus vertex, or from a bus start vertex to a vehicle return
vertex. All of the vehicle arcs are defined in set A
s
. In addition to the vehicles arcs, there are bus arcs in
this map as well, whose start and end vertices are all in set N
b
, representing a feasible bus route. Let A
b
donate these bus arcs. The last kind of arcs is called artificial arcs. We useA
a
to donate the set of artificial
arcs. An artificial arc is added to the graph from a customer delivery node to end bus stop if both nodes
are in the same location and belong to the same customer. These artificial arcs guarantee no vehicles visit
the customer pickup nodes and start bus stops, or end bus stops and their corresponding customer delivery
62
nodes if they are in the same locations and the customer uses a bus route.
Table 8 lists the notations used in the model.
Parameters
n
c
=n
s
=n
b
number of customers/vehicles/bus vertices
tc
b
i
;tc
e
i
time window for customer vertex i
ts
d
i
;ts
a
i
earliest departure and latest return time for vehicle vertex i
tb
i
bus arrival time at bus vertex i
t
i;j
vehicle traveling time from vertex i to vertex j, 8(i;j)2A
s
q vehicle capacity
c
1
i;j
cost for bus arc (i;j);8(i;j)2A
b
c
2
i;j
cost for vehicle arc (i;j);8(i;j)2A
s
M a very big number
Set notation
N
c
=N
s
=N
b
set of vertices for customers/vehicles/bus
N
+
c
=N

c
set of pickup/delivery vertices for the customers
N
+
b;i
=N

b;i
set of start/end bus vertices for customer vertex i, namely for a specific customer i
N
+
b
nc
[
i=1
N
+
b;i
N

b
nc
[
i=1
N

b;i
N N
c
[N
s
[N
b
A
s
=A
b
=A
a
set of vehicle/bus/artificial arcs
A
b;i
set of bus arcs for customer vertex i, namely for a specific customer i
A A
s
[A
b
[A
a
Decision variables
y
i;j
binary variable, is 1 if bus arc (i;j) is used, 8(i;j)2A
b
x
i;j
binary variable, is 1 if vehicle arc or artificial arc (i;j) is used, 8(i;j)2A
s
[A
a
p
i
continuous variable denotes number of passengers in vehicle at vertex i;8i2N
r
i
continuous variable denotes service time at vertex i;8i2N
h
i
flow variable of vehicle after it visits vertex i;8i2N
Table 8: Notations for the exact formulation
63
5.2.1 Basic Formulation
The objective of the model is to find a minimal cost path starting from the first vehicle vertex 2n
c
+ 1 to
the last vehicle vertex 2n
c
+n
s
+ 1, while visiting all the vertices in N
c
and N
s
. There are two parts in the
objective function (95). The first part is the cost for the customers riding the bus and the second part is for
the shared vehicles.
min
X
(i;j)2A
b
c
1
i;j
y
i;j
+
X
(i;j)2As[Aa
c
2
i;j
x
i;j
(95)
Constraints (96) and (97) state that the customer and vehicle vertices must have one in and out degree
except for the first and last vehicle nodes. They ensure the solutions of the MIP are basically feasible for
the problem regardless of time and capacity and visiting sequence constraints. Constraint (98) denotes that
a customer can choose at most one bus route. Constraints (99) and (100) ensure a vehicle goes through bus
vertex i if and only if its corresponding bus arc (i;j) is used. Constraints (101) and (102) enforce similar
requirement for bus vertex j.
X
(i;j)2As[Aa
x
i;j
= 1 8j2N
c
[N
s
=f2n
c
+ 1g (96)
X
(i;j)2As[Aa
x
i;j
= 1 8i2N
c
[N
s
=f2n
c
+n
s
+ 1g (97)
X
(j;k)2A
b;i
y
j;k
1 8i2N
+
c
(98)
y
i;j

X
(k;i)2As[Aa
x
k;i
= 0 8(i;j)2A
b
(99)
y
i;j

X
(i;k)2As[Aa
x
i;k
= 0 8(i;j)2A
b
(100)
y
i;j

X
(k;j)2As[Aa
x
k;j
= 0 8(i;j)2A
b
(101)
y
i;j

X
(j;k)2As[Aa
x
j;k
= 0 8(i;j)2A
b
(102)
Constraint (103) means if a vehicle goes through a customer pickup vertex, or a bus end vertex, the
number of passengers on the vehicle will increase by 1. Constraint (104) ensures the number of customers
on the vehicle is reduced by 1 if it goes through a customer delivery vertex, or a bus start vertex. If a
vehicle goes back to its return vertex, only 1 passenger, namely the driver, should be left in the vehicle.
64
Constraint (105) enforces this requirement. Constraint (106) sets the initial number of passenger on the
vehicle as 1. Constraint (107) places an upper bound on the number of passengers in a vehicle.
p
i
+M(x
i;j
1) + 1p
j
8j2N
+
c
[N

b
; (i;j)2A
s
(103)
p
i
+M(x
i;j
1) 1p
j
8j2N

c
[N
+
b
; (i;j)2A
s
(104)
p
i
+M(x
i;j
1)p
j
8j2N
s
; (i;j)2A
s
(105)
p
i
= 1 8i2N
s
(106)
p
i
q 8i2N
+
c
[N

b
(107)
Constraints (108) means if a vehicle travels from vertex i to j via arc (i;j), the time at vertex j should
bet
i;j
later than the time at vertex i. This is a general constraint for all vertices except the vehicle vertices
who use a single vertex to represent a depart node and a return node simultaneously. Hence we use an extra
constraint (109) to demonstrate such a situation. Constraints (110) to (113) ensure the visiting time at a
vertex should be within its time window.
r
i
+M(x
i;j
1) +t
i;j
r
j
8j2N=N
s
; (i;j)2A
s
(108)
r
i
+M(x
i;j
1) +t
i;j
ts
a
j
8j2N
s
; (i;j)2A
s
(109)
tc
b
i
r
i
tc
e
i
8i2N
c
(110)
ts
d
i
r
i
8i2N
s
=f2n
c
+n
s
+ 1g (111)
tb
i
r
i
8i2N
+
b
(112)
tb
i
r
i
8i2N

b
(113)
So far we have added time and capacity constraints into the model, we still need one more group of
constraints, which are visiting sequences constraints. These constraints ensure a customer will go through
a pickup vertex, start bus vertex, end bus vertex and delivery vertex, or directly from pickup to delivery
vertex. To guarantee the customers’ visiting sequences are ordered properly, we introduce a set of variableh
i
.
Constraints (114) ensures the path is connected and the value of variable h
i
decreases along with the path.
Constraints (115) to (117) ensure the visiting sequences for all customers are correct, namely a customer’s
pickup vertex should be visited before the drop-off vertex and the delivery to a start bus vertex should be
finished before the pickup from an end bus vertex. Constraints (118) and (119) initialize the value of h
j
at
65
the beginning and end vertices of the path.
h
i
+M(x
i;j
1) 1h
j
8j2N;8(i;j)2A
s
(114)
h
i
h
n+i
8i2N
+
c
(115)
h
i
h
j
8i2N
c
;j2N
b;i
(116)
h
2nc+ns+1+i
h
2nc+ns+1+n
b
+i
8i2N
+
b
(117)
h
2nc+1
= 2n
c
+n
s
+ 2
X
(i;j)2A
b
y
i;j

X
(i;j)2Aa
x
i;j
(118)
h
2nc+ns+1
= 1 (119)
The last constraints are the domain constraints.
x
i;j
2f0; 1g 8(i;j)2A
s
[A
a
(120)
y
i;j
2f0; 1g 8(i;j)2A
b
(121)
p
i
;r
i
;h
i
0 8i2N (122)
The MIP formulation for the hybrid transit system consists of the objective function (95) and constraints
from (96) to (121).
5.2.2 Preprocessing and Valid Inequalities
Unfortunately the MIP formulation is hard to solve. In order to improve the speed, we use some strategies
to strength to formulation so that the relaxation of the MIP formulation is tightened.
The first strategy is to eliminate some variables due to feasibility. The elimination rules are mainly from
time window constraints. An arc is infeasible and the corresponding variable is eliminated if any of the
inequalities (123) to (127) are satisfied, which means the vehicle can not travel from one vertex to another
66
while meeting the time window constraint.
tc
b
i
+t
i;j
>tc
e
j
; 8i;j2N
c
(123)
ts
d
i
+t
i;j
>tc
e
j
; 8i2N
s
;j2N
c
(124)
tc
b
j
+t
i;j
>ts
d
i
; 8i2N
s
;j2N
c
(125)
tc
b
i
+t
i;j
>tb
j
; 8i2N
+
c
;j2N
+
b
(126)
tb
j
+t
i;j
>tc
e
i
; 8i2N

c
;j2N

b
(127)
Constraint (123) means a vehicle can not select a customer after another, if the first picked customer’s
earliest access time plus the direct traveling time is still later than the latest access time of the second
selected customer. Constraints (124) and (125) are related to vehicles and customers. A vehicle can not pick
up a customer if their time windows are not overlapped. Constraints (126) and (127) check the feasibility of
using a bus. A customer will not use a bus if the bus is scheduled out of the customer’s time windows.
The second strategy is to tighten the big M. In theory M can be an arbitrary large number; however
too big number will make it computationally inefficient. We select M as small as possible to keep the
constraints still valid. In constraints (103) and (104), we set M = q. In constraint (108), the minimal M
is maxr
i
minr
j
, where maxr
i
and minr
j
are from the time windows of vertex i and j. Similarly, the
minimal M of constraint (109) is ts
a
j
minr
j
. In constraint (114), we set M = 4n
c
+n
s
, which is the
maximal possible value of h
i
.
The last strategy is to add some extra constraints in to the formulation. The extra constraints are
redundant and do not reduce the solution space. However they can tight the relaxation of the problem.
t
nc+i
t
i
t
i;nc+i
; 8i2N
+
c
(128)
X
(j
0
;j)2A
b;i
x
j
0
;j

X
k2N

b;i
X
(k;k
0
)2A
b;i
x
k;k
0 = 0; 8i2N
c
;j2N
+
b;i
(129)
Constraint (128) is intuitive. A customer’s drop off time should be at least t
i;nc+i
, the direct drive time,
later than the pick up time. Constraint (129) ensures that if a customer is dropped at a bus stop, another
vehicle needs to visit one of the other stops along the bus line.
67
5.2.3 Branch-and-cut Algorithm
This hybrid transit problem is a special case of the pickup and delivery problem. The MIP problem given
above is hard to solve. Inspired by the successful implementation of the branch-and-cut algorithm in similar
problems (Cordeau, 2006; Cortés et al., 2010), we develop a branch-and-cut algorithm to solve our hybrid
transit problem.
Note that the ridesharing part of the hybrid transit system is a classical pickup and delivery problem.
Most of the cuts used for the pickup and delivery problem are also valid for the hybrid transit system. We
modified the cuts from Cordeau (2006) and keep the efficient ones in our branch-and-cut algorithm after
primary experiments. The cuts in inequalities (131) to (144) are the ones used in our model and they belong
to the subtour elimination, precedence constraints and generalized order constraints.
Subtour elimination constraints (130)
x
i;j
+x
j;i
+x
nc+i;j
+x
nc+j;i
1; 8i;j2N
+
c
(131)
x
i;nc+j
+x
nc+j;i
+x
nc+i;n+j
1; 8i;j2N
+
c
(132)
x
i;j
+x
j;i
+x
i;nc+i
+x
j;nc+i
+x
nc+i;j
+x
nc+j;i
+x
nc+j;nc+i
2; 8i;j2N
+
c
(133)
x
nc+i;nc+j
+x
nc+j;nc+i
+x
nc+i;j
+x
nc+j;i
1; 8i;j2N
+
c
(134)
x
i;nc+j
+x
nc+j;i
+x
i;j
1; 8i;j2N
+
c
(135)
x
i;nc+j
+x
nc+j;i
+x
i;nc+i
+x
nc+j;nc+i
+x
nc+i;nc+j
+x
i;j
+x
nc+i;j
2; 8i;j2N
+
c
(136)
x
nc+i;j
+ 2x
j;nc+i
+x
j;i
+x
i;nc+i
+x
nc+j;nc+i
2; 8i;j2N
+
c
(137)
x
i;nc+i
+x
nc+i;nc+j
+x
nc+j;i
+ 2x
i;nc+j
+x
i;j
2; 8i;j2N
+
c
(138)
Precedence constraints (139)
x
i;nc+j
+x
nc+j;i
+x
k;i
1; 8k2N
s
;i;j2N
+
c
(140)
x
i;nc+j
+x
nc+j;i
+x
nc+j;k
1; 8k2N
s
;i;j2N
+
c
(141)
x
i;nc+i
+x
i;nc+j
+x
nc+j;i
+x
nc+i;nc+j
+x
nc+j;nc+i
+x
k;i
2; 8k2N
s
;i;j2N
+
c
(142)
Generalized order constraints (143)
x
i;nc+j
+x
nc+j;i
+x
nc+i;j
+x
j;nc+i
1; 8i;j2N
+
c
(144)
68
5.3 Heuristic Solution Method
Since the problem is NP-hard, it is computational difficult to find an optimal solution, even for small size
problems. Therefore, we propose a heuristic method to find a near-optimal solution.
The objective function is a weighted sum of two components: total shared vehicle traveling distance
(VTD), total passenger miles on the buses (BTD), which is identical with the objective function (95) in
Section 5.2.1. The objective function is:
TotalCost =
1
VTD +
2
BTD:
Theparameters
1
;
2
aretheweightsforthefactors. SimilartoAldahaniandDessouky(2003)’sheuristic
method, we separated the hybrid transportation problem into three phases. They are:
1. Insert customer requests into shared vehicles to obtain an initial solution.
2. Improve the initial solution using Tabu search.
3. Check each customer’s feasible bus paths. Take the bus path if the total cost is reduced.
In Aldahani and Dessouky (2003)’s paper, Phase 3 is conducted before Phase 1 and Phase 2. However, after
inserting bus routes, the shared vehicles may face tighter time windows and less adjustment to the schedule
may be possible. Hence, we conduct Phase 3 after Phase 1 and 2.
In Phase 1, we first insert each request into the shared vehicles, compare the increased cost for each
vehicle, and assign the request to the vehicle with minimal increased cost. With the initial solution from
Phase 1, we use Tabu search to make improvements. The search is both within vehicles, namely switching
the visiting sequence of the nodes, and between vehicles, namely changing the vehicle which the customer
takes. A detailed description of Phases 1 and 2 can be found in Wang et al. (2016).
In Phase 3, we attempt to transfer some of the requests to use the fixed-route buses for a portion of their
trip.
1. Select feasible bus paths.
2. Randomize the sequence of requests.
3. For each customer in the sequence, delete the customer’s request in the solution from Phase 2, try to
insert a customer’s feasible bus path in the sequence instead of using a shared vehicle.
69
4. If the system’s total cost is cheaper with using a bus path, take the bus path with the lowest total
cost. Otherwise add the original customer’s request back using a shared vehicle.
5. Repeat steps 2 to 4 for n times and keep the solution with the lowest cost.
Here each bus path consists of three parts, from pickup node to beginning bus stop using a shared vehicle,
from beginning bus stop to end bus stop using a bus, and from end bus stop to delivery node using a shared
vehicle. In Step 1, the feasible bus paths must satisfy the following criteria.
Customer arrives before the bus departure time from the stop. This means there is enough time for
the shared vehicle to take the customer from the pickup node to the beginning bus node, and from the
end bus node to the delivery node.
The customer waiting time at the starting bus node can be no more than W minutes. This criterion
avoids long waiting times for the customers, which is a reasonable assumption in practice.
The ratio of the shortest traveling distance over the bus path total distance must be greater than or
equal to F 1. This criterion eliminates unexpected detour.
The ratio of the shared vehicle traveling distance over bus traveling distance must be less than or equal
to F 2. This criterion makes sure the bus is a "significant" part in the bus path.
The shortest traveling distance must be greater than F 3. If the direct distance between the pickup
node and delivery node is too short, there is no need to take a bus.
The last three criteria follow the ones in the paper by Aldahani and Dessouky (2003). In Step 2, we adopted
the same method used in Phase 1 to check every feasible bus paths.
5.4 Computational Results
In this section, we present the computational results of the exact and heuristic methods mentioned above.
First we run the proposed MIP model on small instances. We implement the branch-and-cut algorithm
based on the standard commercial solver CPLEX. The small instances are used to show the efficiency of
the preprocessing and valid inequalities presented in Section 5.2.2 and the branch-and-cut algorithm in
Section 5.2.3 in searching for an optimal solution. We use the heuristic method to solve large size problems
for a specific simulation environment.
70
All experiments are performed on a PC with Intel Core i7-4790K CPU and 32 GB RAM, running Linux
distribution Ubuntu 14.04 and CPLEX 12.5.
5.4.1 Results on Small Instances
We first create a map with 10*16 nodes. Any two adjacent nodes on the map are connected by an edge with
10 miles distance and it takes 20 minutes for a vehicle and 24 minutes for a bus to travel along an edge.
There is a horizontal bus line that goes through the map with 3 bus stops, as shown in Figure 14 where
the circles represent the bus stops. The cost parameter c
1
i;j
is $1 per mile and c
2
i;j
is $1 per mile. Vehicle
capacityq is set as 4. A bus goes through the bus line once ever day. Small instances for the MIP model are
generated by randomly picking locations for customers and vehicles on the map. We assume all vehicles are
available all day, namely the earliest departure and latest return time cover 24 hours. The time windows for
a customer’s pickup and drop-off nodes are the same and the lengths of customers’ time windows are decided
based on the shortest traveling time between these two nodes. For instance, if the shortest traveling time
between a customer’s pickup and drop-off node is 20 minutes, the time window for the customer will be 40
minutes in a case of twice shortest traveling time. Though the lengths of the time windows are controlled in
the instances, the earliest pick up time for a customer is randomly generated.
Figure 14: Bus map for small instances
Three kinds of formulations, basic, improved and branch-and-cut, are tested in this section. The basic
formulation includes the base MIP formulation and the big M selection. The improved formulation includes
everything presented in Section 5.2.2. The branch-and-cut formulation includes the improved formulation
and the branch-and-cut algorithm in Section 5.2.3. Note we stop each run after 1 CPU hour if the exact
solution procedure has not yet found the optimal solution.
71
Name
Basic Imp B&C
Time Upper Lower Gap Time Upper Lower Gap Time Upper Lower Gap
5V10C4T1 1;745 6;800 6;800 0:00% 889 6;800 6;800 0:00% 207 6;800 6;800 0:00%
5V10C4T2 605 5;000 5;000 0:00% 192 5;000 5;000 0:00% 453 5;000 5;000 0:00%
5V10C4T3 507 6;100 6;100 0:00% 1;335 6;100 6;100 0:00% 917 6;100 6;100 0:00%
5V10C4T4 176 5;750 5;750 0:00% 31 5;750 5;750 0:00% 252 5;750 5;750 0:00%
5V10C4T5 1;086 5;700 5;700 0:00% 1;290 5;700 5;700 0:00% 821 5;700 5;700 0:00%
Avg 824 5;870 5;870 0:00% 747 5;870 5;870 0:00% 530 5;870 5;870 0:00%
7V10C4T1 128 5;650 5;650 0:00% 22 5;650 5;650 0:00% 52 5;650 5;650 0:00%
7V10C4T2 3;600 6;000 5;467 8:89% 3;600 5;900 5;609 4:93% 2;714 5;800 5;800 0:00%
7V10C4T3 3;329 6;000 6;000 0:00% 1;326 6;000 6;000 0:00% 2;611 6;000 6;000 0:00%
7V10C4T4 1;889 6;050 6;050 0:00% 1;531 6;050 6;050 0:00% 1;820 6;050 6;050 0:00%
7V10C4T5 137 6;550 6;550 0:00% 147 6;550 6;550 0:00% 28 6;550 6;550 0:00%
Avg 1;817 6;050 5;943 1:78% 1;325 6;030 5;972 0:99% 1;445 6;010 6;010 0:00%
5V8C6T1 391 3;900 3;900 0:00% 288 3;900 3;900 0:00% 207 3;900 3;900 0:00%
5V8C6T2 2;088 5;650 5;650 0:00% 757 5;850 5;850 0:00% 346 5;850 5;850 0:00%
5V8C6T3 2;541 5;450 5;450 0:00% 1;825 5;450 5;450 0:00% 755 5;450 5;450 0:00%
5V8C6T4 401 4;400 4;400 0:00% 285 4;400 4;400 0:00% 234 4;400 4;400 0:00%
5V8C6T5 1;727 5;250 5;250 0:00% 1;240 5;250 5;250 0:00% 450 5;250 5;250 0:00%
Avg 1;430 4;930 4;930 0:00% 879 4;970 4;970 0:00% 398 4;970 4;970 0:00%
7V10C3T1 4 6;550 6;550 0:00% 6 6;550 6;550 0:00% 4 6;550 6;550 0:00%
7V10C3T2 1;636 6;750 6;750 0:00% 1;025 6;750 6;750 0:00% 665 6;750 6;750 0:00%
7V10C3T3 28 6;100 6;100 0:00% 14 6;100 6;100 0:00% 19 6;100 6;100 0:00%
7V10C3T4 3;600 8;000 5;324 33:44% 3;600 7;700 5;581 27:52% 3;600 7;600 5;595 26:38%
7V10C3T5 143 6;950 6;950 0:00% 63 6;950 6;950 0:00% 51 6;950 6;950 0:00%
Avg 1;082 6;870 6;335 6:69% 942 6;810 6;386 5:50% 868 6;790 6;389 5:28%
5V8C5T1 1;488 5;300 5;300 0:00% 1;273 5;300 5;300 0:00% 681 5;300 5;300 0:00%
5V8C5T2 1;232 5;400 5;400 0:00% 475 5;400 5;400 0:00% 1;144 5;400 5;400 0:00%
5V8C5T3 3;600 4;700 4;600 2:13% 1;301 4;700 4;700 0:00% 346 4;700 4;700 0:00%
5V8C5T4 1;095 5;800 5;800 0:00% 856 5;800 5;800 0:00% 1;245 5;800 5;800 0:00%
5V8C5T5 179 4;950 4;950 0:00% 18 4;950 4;950 0:00% 116 4;950 4;950 0:00%
Avg 1;519 5;230 5;210 0:43% 785 5;230 5;230 0:00% 706 5;230 5;230 0:00%
5V8C3T1 104 5;550 5;550 0:00% 63 5;550 5;550 0:00% 27 5;550 5;550 0:00%
5V8C3T2 264 6;900 6;900 0:00% 53 6;900 6;900 0:00% 56 6;900 6;900 0:00%
5V8C3T3 242 5;700 5;700 0:00% 211 5;700 5;700 0:00% 244 5;700 5;700 0:00%
5V8C3T4 60 5;950 5;950 0:00% 32 5;950 5;950 0:00% 128 5;950 5;950 0:00%
5V8C3T5 3;600 5;900 5;469 7:31% 2;519 5;900 5;900 0:00% 1;450 5;900 5;900 0:00%
Avg 854 6;000 5;914 1:46% 576 6;000 6;000 0:00% 381 6;000 6;000 0:00%
5V10C3T1 3;032 5;550 5;550 0:00% 2;105 5;550 5;550 0:00% 816 5;550 5;550 0:00%
5V10C3T2 2;687 7;500 7;500 0:00% 1;753 7;500 7;500 0:00% 638 7;500 7;500 0:00%
5V10C3T3 3;600 7;600 5;915 22:17% 3;600 6;900 6;063 12:13% 3;600 6;700 5;980 10:75%
5V10C3T4 22 6;150 6;150 0:00% 9 6;150 6;150 0:00% 13 6;150 6;150 0:00%
5V10C3T5 2;213 6;400 6;400 0:00% 1;151 6;400 6;400 0:00% 1;340 6;400 6;400 0:00%
Avg 2;311 6;640 6;303 4:43% 1;724 6;500 6;333 2:43% 1;281 6;460 6;316 2:15%
Table 9: Results comparison between different solution methods
Table 9 has three groups of columns for basic, improved and branch-and-cut formulations. Each group
includes the CPU time (in seconds), upper bound, lower bound and gap for the corresponding formulation.
In total we test 7 different settings and each setting has 5 random instances. Name 5V10C4T1 means it is
the first instance with 5 vehicles, 10 customers and time windows are set to be 4 times of the shortest time.
When comparing the average CPU time among different settings, it is clear that the CPU time increases
if there is more vehicles (5V10C4T vs 7V10C4T), or more customers (5V8C3T vs 5V10C3T). However the
influence of time windows is unclear. The average CPU time of 7V10C4T is larger than 7V10C3T, but
72
on average instances in the setting 5V8C6T uses less CPU time than those in the setting 5V8C5T. When
comparing CPU time among different formulations, the results in Table 9 show that on average, the improved
formulation is faster than the basic formulation. The branch-and-cut formulation is the fastest one among
the tested three formulations. These results confirm that the preprocessing and valid inequalities are effective
in saving CPU time, and the branch-and-cut algorithm can further improve the solution time.
5.4.2 Simulation Settings and Results
The simulation setting is similar to the setting in Section 5.4.1. The only difference is more bus lines are
added. In total, there are five bus lines in the map, which are shown in Figure 15. Each vertical bus line has
a bus stop at every node along the line in the map and each horizontal bus line has a bus stop at every other
node along the line. Buses go along these lines in both directions all day. We assume the initial headway for
all bus lines is 20 minutes.
Figure 15: Bus map for large instances
In the first group of experiments, we vary the ratios of
1
=
2
. A high value of
1
=
2
means the vehicles
are more expensive than the buses. We start the ratio from 0, which means only vehicles are available, and
until1 where buses are free for all the customers.
Ratio
20V100C3T1 20V100C3T2 20V100C3T3 20V100C3T4 20V100C3T5
VTD BTD #B VTD BTD #B VTD BTD #B VTD BTD #B VTD BTD #B
0 9;274:0 0:0 0:0 9;578:0 0:0 0:0 8;950:0 0:0 0:0 8;742:0 0:0 0:0 7;172:0 0:0 0:0
1 9;242:0 18:0 1:2 9;524:0 62:0 0:6 8;910:0 38:0 1:0 8;624:0 92:0 1:4 7;138:0 42:0 0:6
3 8;858:0 228:0 3:0 9;172:0 258:0 4:8 8;564:0 192:0 2:4 8;502:0 462:0 6:4 6;980:0 332:0 5:2
5 8;874:0 284:0 4:2 8;848:0 688:0 10:4 8;450:0 538:0 6:4 8;038:0 712:0 9:8 6;906:0 342:0 4:8
10 8;652:0 652:0 9:4 8;624:0 1;020:0 17:4 8;158:0 564:0 9:8 7;912:0 720:0 9:6 6;432:0 702:0 10:0
1 8;482:0 824:0 13:8 8;408:0 1;562:0 19:8 7;960:0 868:0 12:0 7;734:0 732:0 11:2 6;122:0 1;521:0 17:2
Table 10: Results comparison with different ratios of
1
=
2
73
Figure 16: Average results with different ratios of
1
=
2
Table 10 and Figure 16 present the detailed results with different ratios of
1
=
2
. In total there are
5 instances. Each instance has 20 vehicles and 100 customers. Considering that the heuristic algorithm
includes randomness, we run the heuristic method 5 times for each instance and record the average values of
VTD, BTD and the number of customers using buses (#B) in Table 10. Figure 16 shows the average values
of the 5 instances. When the ratio is 1, namely the costs of using vehicles and buses are equal, only about
0:6% to 1:4% of the customers use buses. More customers are likely to use buses if the cost of using vehicles
gets higher. On average, 15:2% of the customers use buses if the bus system is free. The average VTD is
over 8; 700 when the ratio is 0, and decreases to about 7; 700 when the ratio is1, saving over 10% of VTD.
The results show the ratio of
1
=
2
do have an influence on VTD and BTD.
In the second group of experiments, different numbers of available buses are tested. We control the
number of available buses by tuning the headway of the bus system. The ratio of
1
=
2
is fixed at 5 and
the headway is increased from 20 minutes to 80 minutes. When the headway is 20 minutes, there are 50
buses traveling through each bus line in a single day. When the headway is 80 minutes, the number of buses
decreases to 6 only.
Ratio
20V100C3T1 20V100C3T2 20V100C3T3 20V100C3T4 20V100C3T5
VTD BTD #B VTD BTD #B VTD BTD #B VTD BTD #B VTD BTD #B
20 8;874:0 284:0 4:2 8;848:0 688:0 10:4 8;450:0 538:0 6:4 8;038:0 712:0 9:8 6;906:0 342:0 4:8
40 8;980:0 402:0 5:0 9;050:0 648:0 8:2 8;712:0 472:0 5:8 8;420:0 692:0 8:6 7;102:0 284:0 4:0
60 9;024:0 408:0 4:8 9;502:0 82:0 0:4 8;718:0 482:0 6:0 8;418:0 722:0 7:6 7;140:0 290:0 3:8
80 9;264:0 30:0 0:8 9;492:0 78:0 0:6 8;932:0 42:0 1:4 8;622:0 100:0 2:2 7;208:0 38:0 0:8
Table 11: Results comparison with different headway
74
Figure 17: Average results with different headway
Similar to the previous group of experiments, we use the same instances and each instance is solved 5
times. Table 11 and Figure 17 present the influence of the headway on the customers’ choices. About 4:2%
to 10:4% of the customers use buses when the headway is 20 minutes. The percentage decreases to 0:6%
to 2:2% when the headway is 80 minutes. Figure 17 shows that compared with the huge drop when the
headway increases from 60 minutes to 80 minutes, the influence of headway is smaller when the headway
increases from 20 minutes to 40 minutes. The headway of the transit system will change the customers’
behavior in terms of choosing bus service.
75
6 Conclusion and Future Directions
This dissertation focuses on routing problems for fuel efficient vehicles and proposes three models belonging
to two major research directions. The first two models, the CTRPFS and the CTLRP are designed to solve
the routing problem for CNG trucks with the consideration of fuel stations, which can be categorized into the
direction of reducing emissions per traveled mile. The third model, the hybrid transit system, is developed
to combine the ridesharing service with public transit system and give a routing plan for the shared vehicles.
In this dissertation, we make contributions both in methodology and model development.
In the methodology perspective, we design a MIP model and a specific heuristic method for all three
problems. WeintroduceeffectivevariableeliminationmethodsandstrengthconstraintstotheseMIPmodels.
To further improve the solving speed, we propose an adaptive branch-and-cut algorithm for the CTRPFS
and CTLRP. Numerical experiments show that the adaptive branch-and-cut algorithm is more efficient in
the cut generation procedure compared with the standard branch-and-cut algorithm. Since all of the three
problems are NP-hard and computational difficult to solve optimally, we propose a heuristic method for each
problem. The two heuristic methods for the CTRPFS and CTLRP combine an ALNS with a local search
and MIP model. We test both methods on small and large instances to show their high solution quality
and fast searching speed. The hybrid transit system is solved by a three-phase heuristic method, which is
modified from a Tabu search method.
In the model development perspective, we get some interesting insights from our experimental results.
In the CTRPFS, the results of moderate instances show that increasing the number of CNG fueling stations
can reduce the workload of a single station, not only due to the allocation of workload, but also due to the
reduction in the total number of times trucks use a CNG fueling station. The tank capacity has an impact
on the total traveling distance, especially when the mile range is low. The Ports of Los Angeles and Long
Beach case discusses the relationship between total working time and number of trucks used in the fleet.
In the CTLRP, we first compared the solution quality from the p-median problem and our heuristic
method. The comparison confirms the benefit of simultaneously considering the location and routing prob-
lems.
In the hybrid transit system. We discuss the influence of costs and availability of buses on the customers’
behavior. Experiments show that customers tend to use buses more if the cost of using buses is lower than
shared vehicles. The availability of buses has a similar effect. More customers are likely to use buses if the
headway of buses is shorter.
More work can be done on the three problems in this dissertation. In the methodology perspective, the
76
adaptive branch-and-cut algorithm has a good performance in solving the MIP models for the CTRPFS and
CTLRP. The generalization of this algorithm is a good topic for future research. In the model development
perspective, interesting future work would consider the randomness of the parameters in the models and
change it to be stochastic vehicle routing problems. The problems we solve so far are all static. In the
CTRPFS and CTLRP, we assume the speed of trucks, fuel consumption rate and waiting time in the fuel
station are all constants. In reality, these parameters are influenced by other factors and may have large
effect on the routing and fueling decisions. In the hybrid transit system, the change of the vehicle speed
will influence the feasibility of the customers’ travel plans, so it also needs to be discussed. Another factor
calls for extra attention is about the cost splitting. In the CTLRP, the current objective is to minimize the
total traveling distance of all fleets. However, the global optimal solution may not be the best plan for an
individual fleet. How to fairly split the extra cost of building fuel stations should be resolved.
77
References
Ahmed, A., Chaudhry, A. G., Farooq, H., and Riaz, A. (2013). Short fall of compressed natural gas (CNG) and
issues of affected community; an anthropological perspective. Sci. Int (Lahore), 25(3):623–626.
Aldaihani, M. and Dessouky, M. M. (2003). Hybrid scheduling methods for paratransit operations. Computers &
Industrial Engineering, 45(1):75–96.
Altfuels Communications Group (2017). Worldwide NGV statistics, NGV journal. http://www.ngvjournal.com/
worldwide-ngv-statistics/. Accessed: 2017-06-03.
Augerat, P., Belenguer, J., Benavent, E., Corberán, A., Naddef, D., and Rinaldi, G. (1995). Computational results
with a branch and cut code for the capacitated vehicle routing problem. Technical report, RR949-M, ARTEMIS-
IMAG, Grenoble, France.
Baldacci, R., Bartolini, E., and Mingozzi, A. (2011). An exact algorithm for the pickup and delivery problem with
time windows. Operations Research, 59(2):414–426.
Baldacci, R., Mingozzi, A., Roberti, R., and Calvo, R. W. (2013). An exact algorithm for the two-echelon capacitated
vehicle routing problem. Operations Research, 61(2):298–314.
Barreto, S., Ferreira, C., Paixao, J., and Santos, B. S. (2007). Using clustering analysis in a capacitated location-
routing problem. European Journal of Operational Research, 179(3):968–977.
Bektaş, T., Demir, E., and Laporte, G. (2016). Green vehicle routing. In Green Transportation Logistics, pages
243–265. Springer.
Bektaş, T. and Laporte, G. (2011). The pollution-routing problem. Transportation Research Part B: Methodological,
45(8):1232–1250.
Belenguer, J.-M., Benavent, E., Prins, C., Prodhon, C., and Calvo, R. W. (2011). A branch-and-cut method for the
capacitated location-routing problem. Computers & Operations Research, 38(6):931–941.
Bent, R. and Van Hentenryck, P. (2006). A two-stage hybrid algorithm for pickup and delivery vehicle routing
problems with time windows. Computers & Operations Research, 33(4):875–893.
Capar, I., Kuby, M., Leon, V. J., and Tsai, Y.-J. (2013). An arc cover–path-cover formulation and strategic analysis
of alternative-fuel station locations. European Journal of Operational Research, 227(1):142–151.
Caulfield, B. (2009). Estimating the environmental benefits of ride-sharing: A case study of dublin. Transportation
Research Part D: Transport and Environment, 14(7):527–531.
78
Chen, J.-F. and Wu, T.-H. (2006). Vehicle routing problem with simultaneous deliveries and pickups. Journal of the
Operational Research Society, 57(5):579–587.
Contardo, C., Cordeau, J.-F., and Gendron, B. (2013). An exact algorithm based on cut-and-column generation for
the capacitated location-routing problem. INFORMS Journal on Computing, 26(1):88–102.
Cordeau, J.-F. (2006). A branch-and-cut algorithm for the dial-a-ride problem. Operations Research, 54(3):573–586.
Cordeau, J.-F. and Laporte, G. (2003). A tabu search heuristic for the static multi-vehicle dial-a-ride problem.
Transportation Research Part B: Methodological, 37(6):579–594.
Cortés, C. E., Matamala, M., and Contardo, C. (2010). The pickup and delivery problem with transfers: Formulation
and a branch-and-cut solution method. European Journal of Operational Research, 200(3):711–724.
Dantzig, G. B. and Ramser, J. H. (1959). The truck dispatching problem. Management Science, 6(1):80–91.
Daskin, M. S. (2011). Network and discrete location: models, algorithms, and applications. John Wiley & Sons.
Demir, E., Bektaş, T., and Laporte, G. (2014). A review of recent research on green road freight transportation.
European Journal of Operational Research, 237(3):775–793.
Dessouky, M. M. and Wang, X. (2009). Evaluation of compressed natural gas CNG sweeper operation. Metrans
Technical Report, University of Southern California, Los Angeles.
Doulabi, S. H. H. and Seifi, A. (2013). Lower and upper bounds for location-arc routing problems with vehicle
capacity constraints. European Journal of Operational Research, 224(1):189–208.
Drexl, M. and Schneider, M. (2015). A survey of variants and extensions of the location-routing problem. European
Journal of Operational Research, 241(2):283–308.
Dumas, Y., Desrosiers, J., and Soumis, F. (1991). The pickup and delivery problem with time windows. European
Journal of Operational Research, 54(1):7–22.
Dumitrescu, I., Ropke, S., Cordeau, J.-F., and Laporte, G. (2010). The traveling salesman problem with pickup and
delivery: polyhedral results and a branch-and-cut algorithm. Mathematical Programming, 121(2):269–305.
Erdoğan, S. and Miller-Hooks, E. (2012). A green vehicle routing problem. Transportation Research Part E: Logistics
and Transportation Review, 48(1):100–114.
Furuhata, M., Dessouky, M., Ordóñez, F., Brunet, M.-E., Wang, X., and Koenig, S. (2013). Ridesharing: The
state-of-the-art and future directions. Transportation Research Part B: Methodological, 57:28–46.
79
Gendron, B. and Semet, F. (2009). Formulations and relaxations for a multi-echelon capacitated location–distribution
problem. Computers & Operations Research, 36(5):1335–1355.
Ghilas, V., Demir, E., and Van Woensel, T. (2016a). The pickup and delivery problem with time windows and
scheduled lines. INFOR: Information Systems and Operational Research, 54(2):147–167.
Ghilas, V., Demir, E., and Van Woensel, T. (2016b). A scenario-based planning for the pickup and delivery problem
with time windows, scheduled lines and stochastic demands. Transportation Research Part B: Methodological,
91:34–51.
Glover, F. (1989). Tabu search-part I. ORSA Journal on computing, 1(3):190–206.
Glover, F. (1990). Tabu search-part II. ORSA Journal on computing, 2(1):4–32.
Goodchild, M. F. and Noronha, V. T. (1987). Location-allocation and impulsive shopping: the case of gasoline
retailing. Spatial analysis and location-allocation models, pages 121–136.
Hickman, J., Hassel, D., Joumard, R., Samaras, Z., and Sorenson, S. (1999). Methodology for calculating trans-
portemissionsandenergyconsumption. http://www.transport-research.info/sites/default/files/project/
documents/meet.pdf.
Hiermann, G., Puchinger, J., Ropke, S., and Hartl, R. F. (2016). The electric fleet size and mix vehicle routing
problem with time windows and recharging stations. European Journal of Operational Research, 252(3):995–1018.
Hodgson, M. J. (1990). A flow-capturing location-allocation model. Geographical Analysis, 22(3):270–279.
Intergovernmental Panel on Climate Change (2015). Climate Change 2014: Mitigation of Climate Change, volume 3.
Cambridge University Press.
Johnson, C. and Singer, M. R. (2016). Clean cities 2015 annual metrics report. Technical report, National Renewable
Energy Laboratory.
Kelley, S. and Kuby, M. (2015). Compressed natural gas (CNG) fleets in Southern California: Variations in vehicles
and route types 2. In Transportation Research Board 94th Annual Meeting, number 15-3584.
Khan, M. I., Yasmin, T., and Shakoor, A. (2015). Technical overview of compressed natural gas (CNG) as a
transportation fuel. Renewable and Sustainable Energy Reviews, 51:785–797.
Kim, J.-G. and Kuby, M. (2012). The deviation-flow refueling location model for optimizing a network of refueling
stations. international journal of hydrogen energy, 37(6):5406–5420.
Kuby, M. and Lim, S. (2005). The flow-refueling location problem for alternative-fuel vehicles. Socio-Economic
Planning Sciences, 39(2):125–145.
80
Kuby, M. J., Kelley, S. B., and Schoenemann, J. (2013). Spatial refueling patterns of alternative-fuel and gasoline
vehicle drivers in Los Angeles. Transportation Research Part D: Transport and Environment, 25:84–92.
Landrieu, A., Mati, Y., and Binder, Z. (2001). A tabu search heuristic for the single vehicle pickup and delivery
problem with time windows. Journal of Intelligent Manufacturing, 12(5-6):497–508.
Laporte, G. and Nobert, Y. (1981). An exact algorithm for minimizing routing and operating costs in depot location.
European Journal of Operational Research, 6(2):224–226.
Laughlin, M. and Burnham, A. (2016). Case study - natural gas regional transportation truck. Technical report, U.S.
Department of Energy.
Li, H. and Lim, A. (2003). A metaheuristic for the pickup and delivery problem with time windows. International
Journal on Artificial Intelligence Tools, 12(02):173–186.
LSA Associates, I. (2012). SCAG regional travel demand model and 2008 model validation. http://www.scag.ca.
gov/dataandtools/documents/validationsummaryreport_scag2008val_2012_06_05.pdf. Accessed: 2016-10-21.
Lu, Q. and Dessouky, M. (2004). An exact algorithm for the multiple vehicle pickup and delivery problem.
Transportation Science, 38(4):503–514.
Ma, L., Geng, J., Li, W., Liu, P., and Li, Z. (2013). The development of natural gas as an automotive fuel in China.
Energy policy, 62:531–539.
Mahmoudi, M. and Zhou, X. (2015). Finding optimal solutions for vehicle routing problem with pickup and delivery
services with time windows: A dynamic programming approach based on state-space-time network representations.
arXiv preprint arXiv:1507.02731.
Masson, R., Lehuédé, F., and Péton, O. (2013). An adaptive large neighborhood search for the pickup and delivery
problem with transfers. Transportation Science, 47(3):344–355.
Masson, R., Ropke, S., Lehuédé, F., and Péton, O. (2014). A branch-and-cut-and-price approach for the pickup and
delivery problem with shuttle routes. European Journal of Operational Research, 236(3):849–862.
Min, H., Jayaraman, V., and Srivastava, R. (1998). Combined location-routing problems: A synthesis and future
research directions. European Journal of Operational Research, 108(1):1–15.
Mitchell, G. (2015). Building a business case for compressed natural gas in fleet applications. National Renewable
Energy Laboratory, United States Department of Energy.
Mitrovic-Minic, S. and Laporte, G. (2006). The pickup and delivery problem with time windows and transshipment.
Infor, 44(3):217.
81
Montané, F. A. T. and Galvao, R. D. (2006). A tabu search algorithm for the vehicle routing problem with simulta-
neous pick-up and delivery service. Computers & Operations Research, 33(3):595–619.
Nagy, G. and Salhi, S. (2007). Location-routing: Issues, models and methods. European Journal of Operational
Research, 177(2):649–672.
Nanry, W. P. and Barnes, J. W. (2000). Solving the pickup and delivery problem with time windows using reactive
tabu search. Transportation Research Part B: Methodological, 34(2):107–121.
Nicholas, M., Handy, S., and Sperling, D. (2004). Using geographic information systems to evaluate siting and
networks of hydrogen stations. Transportation Research Record: Journal of the Transportation Research Board,
(1880):126–134.
Nicholas, M. and Ogden, J. (2006). Detailed analysis of urban station siting for california hydrogen highway network.
Transportation Research Record: Journal of the Transportation Research Board, (1983):121–128.
Parragh, S. N., Doerner, K. F., and Hartl, R. F. (2007). A survey on pickup and delivery problems.
Prins, C., Prodhon, C., Ruiz, A., Soriano, P., and Wolfler Calvo, R. (2007). Solving the capacitated location-routing
problem by a cooperative lagrangean relaxation-granular tabu search heuristic. Transportation Science, 41(4):470–
483.
Prodhon, C. and Prins, C. (2014). A survey of recent research on location-routing problems. European Journal of
Operational Research, 238(1):1–17.
Qian, J. and Eglese, R. (2016). Fuel emissions optimization in vehicle routing problems with time-varying speeds.
European Journal of Operational Research, 248(3):840–848.
Quadrifoglio, L., Dessouky, M. M., and Ordóñez, F. (2008). Mobility allowance shuttle transit (MAST) services: Mip
formulation and strengthening with logic constraints. European Journal of Operational Research, 185(2):481–494.
Quadrifoglio, L., Dessouky, M. M., and Palmer, K. (2007). An insertion heuristic for scheduling mobility allowance
shuttle transit (mast) services. Journal of Scheduling, 10(1):25–40.
Ropke, S. and Cordeau, J.-F. (2009). Branch and cut and price for the pickup and delivery problem with time
windows. Transportation Science, 43(3):267–286.
Ropke, S., Cordeau, J.-F., and Laporte, G. (2007). Models and branch-and-cut algorithms for pickup and delivery
problems with time windows. Networks, 49(4):258–272.
Ropke, S. and Pisinger, D. (2006). An adaptive large neighborhood search heuristic for the pickup and delivery
problem with time windows. Transportation Science, 40(4):455–472.
82
Salimifard, K. and Raeesi, R. (2014). A green routing problem: optimising CO2 emissions and costs from a bi-fuel
vehicle fleet. International Journal of Advanced Operations Management, 6(1):27–57.
Schiffer, M. and Walther, G. (2017). The electric location routing problem with time windows and partial recharging.
European Journal of Operational Research.
Schneider, M., Stenger, A., and Goeke, D. (2014). The electric vehicle-routing problem with time windows and
recharging stations. Transportation Science, 48(4):500–520.
Shahraeeni, M., Ahmed, S., Malek, K., Van Drimmelen, B., and Kjeang, E. (2015). Life cycle emissions and cost of
transportation systems: Case study on diesel and natural gas for light duty trucks in municipal fleet operations.
Journal of Natural Gas Science and Engineering, 24:26–34.
Shaw, P. (1998). Using constraint programming and local search methods to solve vehicle routing problems. In
Principles and Practice of Constraint Programming CP98, pages 417–431. Springer.
Stålhane, M., Andersson, H., Christiansen, M., Cordeau, J.-F., and Desaulniers, G. (2012). A branch-price-and-
cut method for a ship routing and scheduling problem with split loads. Computers & Operations Research,
39(12):3361–3375.
Stiglic, M., Agatz, N., Savelsbergh, M., and Gradisar, M. (2015). The benefits of meeting points in ride-sharing
systems. Transportation Research Part B: Methodological, 82:36–53.
Tong, F., Jaramillo, P., and Azevedo, I. M. (2015). Comparison of life cycle greenhouse gases from natural gas
pathways for light-duty vehicles. Energy & Fuels, 29(9):6008–6018.
Tuzun, D. and Burke, L. I. (1999). A two-phase tabu search approach to the location routing problem. European
journal of operational research, 116(1):87–99.
Upchurch, C. and Kuby, M. (2010). Comparing the p-median and flow-refueling models for locating alternative-fuel
stations. Journal of Transport Geography, 18(6):750–758.
U.S. Department of Energy (2014). Cng vehicle fueling animation. http://www.afdc.energy.gov/vehicles/cng_
tank_animation_text.html. Accessed: 2017-05-05.
U.S. Department of Energy (2017a). Alternative fueling station locator. http://www.afdc.energy.gov/locator/
stations/. Accessed: 2017-06-03.
U.S. Department of Energy (2017b). Compressed natural gas fueling stations. http://www.afdc.energy.gov/fuels/
natural_gas_cng_stations.html. Accessed: 2017-05-25.
83
U.S. Department of Energy (2017c). Los angeles clean cities coalition. https://cleancities.energy.gov/
coalitions/los-angeles. Accessed: 2017-03-22.
U.S. Department of Energy (2017d). National clean fleets partnership. https://cleancities.energy.gov/fleets/.
Accessed: 2017-03-22.
U.S. Department of Energy (2017e). Natural gas vehicle technology forum. https://cleancities.energy.gov/
natural-gas-forum/#activities. Accessed: 2017-03-22.
U.S. Environmental Protection Agency (2016). Inventory of us greenhouse gas emissions and sinks: 1990-2014.
Technical report.
U.S. Environmental Protection Agency (2017). Global greenhouse gas emissions data. https://www.epa.gov/
ghgemissions/global-greenhouse-gas-emissions-data. Accessed: 2017-03-12.
Wang, X., Dessouky, M., and Ordonez, F. (2016). A pickup and delivery problem for ridesharing considering conges-
tion. Transportation Letters, 8(5):259–269.
Wang, Y.-W. and Wang, C.-R. (2010). Locating passenger vehicle refueling stations. Transportation Research Part
E: Logistics and Transportation Review, 46(5):791–801.
Xiao, Y. and Konak, A. (2016). The heterogeneous green vehicle routing and scheduling problem with time-varying
traffic congestion. Transportation Research Part E: Logistics and Transportation Review, 88:146–166.
Xiao, Y., Zhao, Q., Kaku, I., and Xu, Y. (2012). Development of a fuel consumption optimization model for the
capacitated vehicle routing problem. Computers & Operations Research, 39(7):1419–1431.
Xu, H., Chen, Z.-L., Rajagopal, S., and Arunapuram, S. (2003). Solving a practical pickup and delivery problem.
Transportation science, 37(3):347–364.
Yang, J. and Sun, H. (2015). Battery swap station location-routing problem with capacitated electric vehicles.
Computers & Operations Research, 55:217–232.
84

Abstract (if available)

Abstract In this dissertation, we introduce three routing problems for fuel efficient vehicles. The first two models focus on exploring alternative fuels, which are CNG Truck Routing Problem with Fueling Stations (CTRPFS) and CNG Truck Location and Routing Problem (CTLRP). The last model is about changing people's behavior. We call it the hybrid transit system in which the routes of shared vehicles are scheduled together with the public transit system. Both exact and heuristic algorithms are proposed for the three models. We conduct experiments to show the efficiency of our solving methods as well as insights from real world cases.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Dynamic programming-based algorithms and heuristics for routing problems

PDF

Routing for ridesharing

PDF

Train scheduling and routing under dynamic headway control

PDF

Vehicle routing and resource allocation for health care under uncertainty

PDF

The robust vehicle routing problem

PDF

New approaches for routing courier delivery services

PDF

Package delivery with trucks and UAVs

PDF

Models and algorithms for energy efficient wireless sensor networks

PDF

Models and algorithms for pricing and routing in ride-sharing

PDF

Cost-sharing mechanism design for freight consolidation

PDF

An online cost allocation model for horizontal supply chains

PDF

Routing and inventory model for emergency response to minimize unmet demand

PDF

Train routing and timetabling algorithms for general networks

PDF

Evaluating city logistics using two-level location routing modeling and SCPM simulation

PDF

Continuous approximation formulas for location and hybrid routing/location problems

PDF

Continuous approximation for selection routing problems

PDF

Continuous approximation formulas for cumulative routing optimization problems

PDF

Novel queueing frameworks for performance analysis of urban traffic systems

PDF

Computational geometric partitioning for vehicle routing

PDF

Integration of truck scheduling and routing with parking availability

Asset Metadata

Creator Shao, Yihuan (Ethan) (author)

Core Title Routing problems for fuel efficient vehicles

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Industrial and Systems Engineering

Publication Date 09/28/2017

Defense Date 06/29/2017

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag ALNS,branch-and-cut algorithm,fuel efficient vehicles,green VRP,OAI-PMH Harvest,routing and location problem

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Dessouky, Maged (committee chair), Ordonez, Fernando (committee member), Savla, Ketan (committee member)

Creator Email shaoyihuan@gmail.com,yihuansh@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-438707

Unique identifier UC11265611

Identifier etd-ShaoYihuan-5791.pdf (filename),usctheses-c40-438707 (legacy record id)

Legacy Identifier etd-ShaoYihuan-5791.pdf

Dmrecord 438707

Document Type Dissertation

Rights Shao, Yihuan (Ethan)

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA