# An Alternative Fuel Refueling Station Location Model considering Detour Traffic Flows on a Highway Road System

Authors: Seong Wook Hwang; Sang Jin Kweon; Jose A. Ventura
Journal: Journal of Advanced Transportation
 (2020)

---

## Abstract

With the development of alternative fuel (AF) vehicle technologies, studies on finding the potential location of AF refueling stations in transportation networks have received considerable attention. Due to the strong limited driving range, AF vehicles for long-distance intercity trips may require multiple refueling stops at different locations on the way to their destination, which makes the AF refueling station location problem more challenging. In this paper, we consider that AF vehicles requiring multiple refueling stops at different locations during their long-distance intercity trips are capable of making detours from their preplanned paths and selecting return paths that may be different from original paths for their round trips whenever AF refueling stations are not available along the preplanned paths. These options mostly need to be considered when an AF refueling infrastructure is not fully developed on a highway system. To this end, we first propose an algorithm to generate alternative paths that may provide the multiple AF refueling stops between all origin/destination (OD) vertices. Then, a new mixed-integer programming model is proposed to locate AF refueling stations within a preselected set of candidate sites on a directed transportation network by maximizing the coverage of traffic flows along multiple paths. We first test our mathematical model with the proposed algorithm on a classical 25-vertex network with 25 candidate sites through various scenarios that consider a different number of paths for each OD pair, deviation factors, and limited driving ranges of vehicles. Then, we apply our proposed model to locate liquefied natural gas refueling stations in the state of Pennsylvania considering the construction budget. Our results show that the number of alternative paths and deviation distance available significantly affect the coverage of traffic flows at the stations as well as computational time.


---

## Body

## 1. Introduction

Reducing greenhouse gas (GHG) emissions in the transportation sector is one of the most vital steps in fighting against global warming in the United States (U.S.). According to the U.S. Environmental Protection Agency [1], the transportation sector generates the largest share of GHG emissions. In order to cut down tail-pipe emissions in the transportation sector, vehicles using alternative fuel (AF), such as biodiesel, hydrogen, electrical energy, and natural gas, have received significant attention in recent years because they emit less well-to-wheel GHG than that of vehicles using traditional fossil fuels, such as diesel and gasoline. Currently, 367 light-duty and 216 medium- and heavy-duty vehicles for the 2018 and 2019 model years are available in the U.S. AF vehicle market [2].While the public interest in using AF vehicles instead of conventional vehicles has increased, the number of public refueling stations for AF vehicles is still insufficient, especially for intercity trips between urban and rural counties. Table1 shows the 2010 census total population, number of electric charging stations, and the number of electric charging stations per 100,000 residents in the U.S. based on the Census Bureau’s urban-rural classification [3, 4]. In this table, counties are categorized into three groups according to their population density: mostly urban, mostly rural, and completely rural. It indicates that urban areas have more than twice as many electric charging stations than rural areas with a similar population size. One of the barriers to the invigoration and development of AF infrastructures in highway systems, which play a major role in intercity trips among urban and rural counties, is high construction cost. For example, the construction cost of one natural gas refueling station is at least $1.8 million in the Pennsylvania Turnpike [5].Table 1
Number of electric charging stations per 100,000 residents at urban and rural counties in the U.S.
2010 Census total population (P)
Number of electric charging stations (C)
Number of electric charging stations per 100,000 residentsC/P×100,000
Mostly urban counties
266,567,329
17,330
6.50
Mostly rural counties
36,811,523
989
2.69
Completely rural counties
5,360,991
156
2.91
Total
308,739,843
18,475
5.98Due to the sparse distribution of AF refueling stations between urban areas in highway systems, AF vehicles with a short driving range that travel long-distance intercity trips may need to use longer paths with refueling availability, including multiple refueling stops at different locations on the way to their destination, rather than their shortest paths without available refueling stations, so they can safely complete their trips. Thus, AF vehicle drivers may need to make detours to be able to refuel their vehicles.Generally, road structures in highway systems are different than those in other transportation networks. First, highway roads are divided into two pathways separated by either a raised barrier or an unpaved median. In order to offer uninterrupted traffic flow, highway roads do not have any traffic intersections, and vehicles are only able to enter/leave highway systems through entrance and exit ramps. Next, highway systems have built-in service facilities where drivers are able to take a rest and refuel their vehicles. Since highway roads are physically partitioned by a barrier or an unpaved median, some built-in service facilities, called single-access stations, can only be accessed from one side of the road, while the rest, called dual-access stations, can be accessed from both sides of the road. Hwang et al. [6, 7] and Ventura et al. [5] model this type of transportation systems as directed transportation networks.In general, concerning the design of an AF refueling infrastructure along transportation networks, a number of studies allowing repeated trips between origin/destination (OD) vertices assume that vehicles make symmetric round trips traveling along preplanned (shortest) paths between the corresponding OD vertices. It also assumes that the set of candidate station locations is a subset of vertices in the network, and therefore, all the candidate locations are dual-access sites. These assumptions imply that, if a vehicle travels from an origin to a destination on a single path by filling up at some stations, it also stops by the same stations on the return path back to the origin. These assumptions, however, have made the refueling station location problem for AF vehicles less practical since they do not reflect the characteristics of AF vehicles well at its early stage of market.By relaxing the assumptions listed above, this paper aims at determining more reliable locations of AF refueling stations in real-world applications based on the distinct features of AF vehicles traveling intercity trips on directed transportation networks. We first consider that some candidate sites for AF refueling stations are single-access and drivers may choose different return paths from original paths to be able to refuel their vehicles in both directions. Also, since several paths are available between ODs, AF refueling availability is considered as one of the AF vehicle drivers’ top priorities when they select paths to travel in highway systems. Thus, we allow AF vehicles to make nonsymmetric round trips between their ODs. We first generate multiple paths between all OD pairs through a revisedk-th shortest path algorithm considering a maximum deviation distance. Then we formulate a new mixed-integer programming model that considers a set of predetermined paths for each OD pair and a preliminary set of candidate station location sites, including single-access and dual-access sites, on a directed transportation network in which AF vehicles are able to use any of the corresponding OD paths depending on the availability of refueling service. The proposed model determines the optimal set of station locations for a given number of stations and the selected round trips for all ODs that maximize the total traffic flow covered (in round trips per time unit) by the stations. For computational experiments, we first apply the proposed model to a classical 25-vertex network with 25 candidate sites through various scenarios. We also validate our proposed model with a budget constraint to construct AF refueling stations in the state of Pennsylvania.Our model proposed in this study is applicable to various types of AF vehicles, especially to liquefied natural gas (LNG) vehicles, with their refueling station location problems. LNG vehicles are similar to the existing long-haul vehicles powered by diesel in terms of powertrain and refueling, but LNG vehicles are known to provide economic and environmental benefits; thus, LNG vehicles are well-suited for replacing the current long-haul vehicles powered by diesel on a highway road system. For example, UPS has been working to shift their high carbon-fueled vehicles to new generation of LNG vehicles since the late 20th century and has been increasing their number [8]. In the U.S., UPS has deployed their LNG vehicles mainly in Indianapolis, Chicago, Earth City, and Nashville, and plans to use these LNG vehicles in larger areas [9]. In addition to LNG vehicles, battery-electric vehicles and hydrogen fuel cell vehicles will also be applied to our proposed model once they are fully available for long haul logistics on a highway system.The remaining of this paper is organized as follows. Section2 reviews the related literature and shortly discusses the main distinctions of our research work over the existing studies as well. In Section 3, we first introduce the AF refueling station location problem with detour traffic flows on a highway road system and provide properties of feasible paths. Next, an algorithm that generates a set of multiple feasible paths for which the properties are satisfied on a given network is proposed. In Section 4, we provide covering conditions with candidate sites to cover round trips for each OD pair. Then, we propose a mixed-integer programming model to locate a given number of AF refueling stations on a directed transportation network with the objective of maximizing the total traffic flow covered, considering multiple paths for each OD pair. In Section 5, the proposed model is tested on a well-known 25-vertex network to evaluate the effects of the number of multiple paths, maximum deviation distance, and vehicle driving range on the coverage of traffic flows. The proposed model is then validated in Section 6 with an application to the state of Pennsylvania to demonstrate its performance on a large-size problem. Lastly, we provide conclusions and discuss the future work of AF refueling station location problems in Section 7.


## 2. Literature Review

In this section, we take a look at the literature relevant to this paper. First, in Subsection2.1, we review the literature addressing k-shortest path problems, and in Subsection 2.2, we examine the literature related to AF refueling station location problems. Then, we organize the main distinctions of our research work by comparing with the relevant literature in Subsection 2.3.


### 2.1.k-shortest Path Problems

Thek-shortest path problem is to find k shortest paths between two vertices in a given network in a nondecreasing order of length, where k refers to the number of shortest paths to find. The k-shortest path problem can be classified into two types according to whether paths are allowed to have cycles or required to be simple with no cycles.The first type ofk-shortest path problem, proposed by Hoffman and Pavley [10], allows repeated vertices along any path. Fox [11] develops an algorithm to apply this type of problem to probabilistic networks. Eppstein [12] uses the concept of binary heap data structure in a nondecreasing order of additional path length due to deviation to find the k shortest paths in polynomial time. Since the experimental results of Eppstein’s algorithm still take considerable time to find the k shortest paths, Jiménez and Marzal [13] present a modified version of Eppstein’s algorithm to improve its practical performance. Aljazzar and Leue [14] propose a directed search algorithm to search for the k shortest paths between two given vertices of a network. Their algorithm provides the same asymptotic worst-case complexity but uses less memory than Eppstein’s algorithm. Liu et al. [15] propose a novel k-shortest path algorithm for neural networks. Given a set of battery exchange stations, Adler et al. [16] suggest a polynomial time algorithm to solve the electric vehicle shortest-walk problem with battery exchanges considering vehicle’s limited driving range. This algorithm shows the chance of extension to the k-shortest path problem by transforming the original traffic network into the so called refueling shortest path network.The second type ofk-shortest path problem does not allow any repeated vertex along a path, and is thus called the k-shortest simple (or loopless) path problem. Since this type of problem needs an additional constraint to allow only loopless paths, it turns out to be more challenging than the first type of problem [17]. Pollack [18] introduces the concept of the k-shortest simple path problem and solves it by modifying Hoffman and Pavley’s method [10], so as to avoid paths containing repeated vertices. Clarke et al. [19] present a branch-and-bound procedure to find the k-shortest simple paths, but this method requires a significant computational effort and storage requirements in the main memory. Sakarovitch [20] first identifies several shortest paths that may contain repeated vertices by using the efficient version of Hoffman and Pavley’s method [10], and then picks up k-shortest simple paths among them. By applying a procedure that partitions a path into two subpaths, Yen [21] develops an efficient algorithm to find the k-shortest simple path. Lawler [22] generalizes a procedure that can reduce the amount of storage required in solving the k-shortest simple path problem. Katoh et al. [23] present an improved version of Yen’s algorithm that solves the problem efficiently in an undirected network. The practical performance of Yen’s algorithm is comparatively analyzed with other k-shortest path algorithms [24, 25]. An implementation of Yen’s algorithm is also studied to improve its practical performance [17, 26]. Hershberger and Suri [27, 28] suggest the efficient replacement path algorithm for finding k shortest simple paths in a directed network. Zeng et al. [29] present a heuristic k-shortest path algorithm that is based on Yen’s algorithm when determining an eco-friendly path that results in minimum carbon dioxide emissions from light-duty vehicles.


### 2.2. Alternative Fuel Refueling Station Location Problems

The maximal covering location and the set-covering location approaches are two main streams of research for addressing AF refueling station location problems. Kuby and Lim [30] are one of the first applying a maximal covering location model to solve an AF refueling station location problem. They introduce the flow-refueling location model (FRLM) to find the optimal location of refueling stations for AF vehicles by considering their limited driving range per refueling with the objective of maximizing the total traffic flow covered. Upchurch and Kuby [31] show that the FRLM identifies more stable locations for AF refueling stations than the p-median model, especially in a statewide case study. In general, the FRLM requires a significant computational effort to pregenerate all feasible location combinations of refueling stations that allow vehicles to make round trips between ODs. Lim and Kuby [32] propose three heuristic versions for the FRLM. Kuby et al. [33] apply two of them to locate hydrogen refueling stations in Florida. Capar and Kuby [34] present a new formulation for the FRLM that skips the pregeneration of all feasible combinations on every path. Capar et al. [35] suggest an arc-cover-path-cover model that focuses on the arcs comprising each path, so as to solve the problem efficiently without the pregeneration of all feasible location combinations of stations on each path for the FRLM. Jochem et al. [36] apply the Capar et al. [35] model to allocate charging stations in the German autobahn. While most AF refueling station location problems assume that the AF refueling stations can only be located at the vertices, Kuby and Lim [37] and Ventura et al. [38] consider additional candidate sites along arcs. Kweon et al. [39] extend the approach suggested in Ventura et al. [38] to locate a refueling station anywhere along a tree network to the case where a portion of drivers are willing to deviate to receive refueling service. He et al. [40] propose a bilevel model to solve the optimal locations of electric charging stations, through taking the driving range limitation of an electric vehicle, the battery charging time required, and the situation in which some electric vehicle drivers possibly charge at home into account.While Kuby and Lim [30] and the following studies solve the AF refueling station location problem using a maximal covering location problem, Wang and Lin [41] solve this problem using a set-covering model with the objective of minimizing the total cost of locating stations to cover all the traffic flow on a given transportation network. Wang and Wang [42] integrate Wang and Lin’s model [41] into the classic set-covering model considering vertex-based and flow-based demands for the AF refueling station location problem. Since Wang and Lin’s model [41] requires a significant computational effort to evaluate the effect of the limited vehicle driving range on the number of charging stations needed for achieving multiple origin-destination intercity travel with electric vehicles on Taiwan, MirHassani and Ebrazi [43] propose a novel approach by using the conservation of flow law, which is able to solve large-scale problems. Chung and Kwon [44] extend MirHassani and Ebrazi’s model [43] to a multiperiod planning problem for allocating charging stations in the Korea Expressway. Hosseini and MirHassani [45] use the idea of MirHassani and Ebrazi’s model [43] to propose a two-stage stochastic mixed integer programming model for the refueling station location problem, where the traffic flow of AF vehicles is uncertain and portable AF refueling stations are considered. Kang and Recker [46] use the idea of the set-covering problem to locate hydrogen refueling stations with the assumption that at most one refueling stop per trip is required in a city. Assuming that vehicles only require one refueling stop per trip, two refueling station location models, including the versions that consider limited capacity of refueling stations [47], and refueling demand uncertainty and driver route choice behavior [48], are developed to minimize the total cost imposed on a planner and drivers over multiple time periods. Using the Adaptive Large Neighborhood Search algorithm [49] and the Adaptive Variable Neighborhood Search algorithm [50], locations for battery swap stations and electric vehicle routes are determined to provide services with the objective of minimizing the sum of the station construction cost and routing cost. To minimize the total cost to locate electric vehicle charging stations in road networks, Gagarin and Corcoran [51] suggest a novel approach that searches for the dominating set of locations among the candidate locations whose distance is below a certain threshold from a given driver. Using a parallel computing strategy, Tran et al. [52] propose an efficient heuristic algorithm for location of AF refueling stations based on the solution of a sequence of subproblems.In general, drivers often deviate from their original paths of shortest travel time or distance to be able to refuel their vehicles [53]. Both the maximal covering location and the set-covering location approaches have been extended to consider driver’s deviation options under a variety of situations. Kim and Kuby [54] address the deviation version of the FRLM (DFRLM) where drivers are able to deviate from their paths of shortest length between ODs, and Kim and Kuby [55] propose a heuristic algorithm for the DFRLM to solve large-scale problems. Huang et al. [56] develop a new model, called the multipath refueling location model (MPRLM), by considering multiple deviation paths between ODs. For intracity trips that require at most one refueling stop, Miralinaghi et al. [47, 48] suggest deviation versions of the set-covering location problem to find potential locations of AF refueling stations. For intercity trips of large-scale problems, Yıldız et al. [57] use a branch and price algorithm, which does not need pregeneration of path generation, for the AF refueling station location problem adding the routing aspect of drivers. Most recently, Göpfert and Bock [58] and Arslan et al. [59] suggest novel projection and branch and cut methods in dealing with the deviation version of the refueling infrastructure planning, so as to extend the computational efficiency even further to solve large-size problem instances with less computational effort.


### 2.3. Main Distinctions of Our Research Work

Comparing with the literature listed above, we shortly provide the main distinctions of our research work over the existing studies as follows:(i)
Our problem is an extension of Hwang et al. [6] problem to consider potential deviation paths on directed transportation networks, such as highway network systems. This leads to consider (1) the mixed set of single-access and dual-access candidate sites to locate AF refueling stations, and (2) nonsymmetric round trips between ODs, where return paths are allowed to be different from original paths for refueling services in both directions.(ii)
Our study is well-suited for the AF refueling station location problem, specially with LNG vehicles traveling long-distance intercity round trips. Some number of recent studies, including Miralinaghi et al. [47, 48], limits this type of problem only for AF vehicles with intracity trips, which needs at most one refueling stop per trip. On the other hand, we apply covering condition procedures depending on LNG vehicle driving range, so as for LNG vehicles traveling intercity trips to allow multiple refueling stops at different locations.(iii)
In our research work, paths are not fixed for every OD pair. Instead, paths are flexible to consider detour traffic flows. While Kim and Kuby [54, 55] and Huang et al. [56] apply Hoffman and Pavley’s algorithm [10] and Yen’s algorithm [21] to take deviation paths into account, we develop a new algorithm based on Eppstein’s algorithm [12] to rigorously and efficiently find k shortest paths allowing repeated vertices along paths within the tolerance (i.e., maximum deviation distance) of the driver. This helps reduce the computational effort in solving our proposed mixed-integer programming model.(iv)
We conduct computational experiments on a classical 25-vertex network (with 300 OD pairs) and a large network at Pennsylvania state (with 2,211 OD pairs) (1) to evaluate the performance of our research work; and (2) to analyze the effects of number of alternative paths including deviation paths, maximum deviation distance, and vehicle driving range on the optimal location as well as the corresponding total traffic flow covered.The major differences between the proposed model and the existing studies that are directly relevant to ours are summarized in Table2.Table 2
Comparative summary between the proposed model and the existing studies directly relevant to our research work.
Research work
Proposed model
Kim and Kuby [54]
Kim and Kuby [55]
Huang et al. [56]
Miralinaghi et al. [47]
Miralinaghi et al. [48]
Hwang et al. [6]
Hwang et al. [7]
Ventura et al. [5]
Deviation
Allowed
Allowed
Allowed
Allowed
Allowed
Allowed
Not allowed
Not allowed
Not allowed
Method to find deviation paths
Eppstein’s algorithm [12] to find k-shortest paths allowing repeated vertices along paths
Hoffman and Pavley’s algorithm [10] to find k-shortest paths allowing repeated vertices along paths
Greedy-adding and Greedy-adding with substitution heuristics to findk-shortest paths allowing repeated vertices along paths
Yen’s algorithm [21] to find k-shortest simple paths which do not allow repeated vertices along paths
The link-based user equilibrium condition for network traffic to find a deviation path containing one refueling stop
The link-based user equilibrium framework for network traffic to find a deviation path containing one refueling stop
N/A
N/A
N/A
Directed network?
Directed
Undirected
Undirected
Undirected
Undirected
Undirected
Directed
Directed
Directed
Candidate sites
Dual- and single-access
Dual-access
Dual-access
Dual-access
Dual-access
Dual-access
Dual- and single-access
Dual- and single-access
Dual- and single-access
Multiple refueling stops or one refueling stop
Multiple refueling stops
Multiple refueling stops
Multiple refueling stops
Multiple refueling stops
One refueling stop
One refueling stop
Multiple refueling stops
Multiple refueling stops
Multiple refueling stops
Main constraints
(i) Number of stations
(i) Number of stations
(i) Number of stations
(i) Satisfy all travel demand
(i) Satisfy all multiperiod travel demand
(i) Satisfy all travel demand when it is uncertain
(i) Number of stations
(i) Number of stations
(i) Number of stations
(ii) Number of deviation paths
(ii) Vehicle driving range
(ii) Vehicle driving range
(ii) Vehicle driving range
(ii) Capacity of stations
(ii) Vehicle driving range
(ii) Vehicle driving range
(ii) Vehicle driving range
(iii) Deviation factor
(iii) Flow volume decay with the increase of deviation
(iii) Flow volume decay with the increase of deviation
(iii) Flow
(iii) Multiple vehicle classes
(iii) Different capital cost to build stations by region
(iv) Vehicle driving range
(iv) Different fuel tank levels at OD pairs
(v) Different capital cost to build stations by type and/or region
Objective(s)
Maximize the total traffic flow covered, considering multiple paths including shortest and deviation paths for each OD pair
Maximize the total flows refueled on deviation paths
Maximize the total flow refueled when deviation is possible
Minimize the total cost of locating AF stations when multiple paths including shortest paths and deviation paths are considered between an OD pair
Minimize the sum of the total station construction and station operational costs for the multiperiod travel demand
Minimize the sum of the total annualized construction costs with the cost of the total system travel time when travel demand is uncertain
Maximize the total traffic flow covered by the stations
Maximize the total traffic flow covered by the stations on a multiclass vehicle transportation network with different fuel tank levels
(i) Maximize the total vehicle-miles traveled per year covered by the stations
(ii) Minimize the capital cost for constructing the refueling infrastructure
Applied to
(i) 25-vertex network
(i) 25-vertex network
(i) 25-vertex network
(i) 25-vertex network
(i) Sioux-Falls network
(i) Sioux-Falls network
(i) Pennsylvania Turnpike
(i) Pennsylvania Turnpike
(i) Pennsylvania Turnpike
(ii) The state of Pennsylvania
(ii) The state of Florida
(ii) Sioux-Falls network
(ii) Mashhad road network
Nonsymmetric round trips?
Allowed
Not allowed
Not allowed
Not allowed
Not allowed
Not allowed
Not allowed
Not allowed
Not allowed


## 2.1.k-shortest Path Problems

Thek-shortest path problem is to find k shortest paths between two vertices in a given network in a nondecreasing order of length, where k refers to the number of shortest paths to find. The k-shortest path problem can be classified into two types according to whether paths are allowed to have cycles or required to be simple with no cycles.The first type ofk-shortest path problem, proposed by Hoffman and Pavley [10], allows repeated vertices along any path. Fox [11] develops an algorithm to apply this type of problem to probabilistic networks. Eppstein [12] uses the concept of binary heap data structure in a nondecreasing order of additional path length due to deviation to find the k shortest paths in polynomial time. Since the experimental results of Eppstein’s algorithm still take considerable time to find the k shortest paths, Jiménez and Marzal [13] present a modified version of Eppstein’s algorithm to improve its practical performance. Aljazzar and Leue [14] propose a directed search algorithm to search for the k shortest paths between two given vertices of a network. Their algorithm provides the same asymptotic worst-case complexity but uses less memory than Eppstein’s algorithm. Liu et al. [15] propose a novel k-shortest path algorithm for neural networks. Given a set of battery exchange stations, Adler et al. [16] suggest a polynomial time algorithm to solve the electric vehicle shortest-walk problem with battery exchanges considering vehicle’s limited driving range. This algorithm shows the chance of extension to the k-shortest path problem by transforming the original traffic network into the so called refueling shortest path network.The second type ofk-shortest path problem does not allow any repeated vertex along a path, and is thus called the k-shortest simple (or loopless) path problem. Since this type of problem needs an additional constraint to allow only loopless paths, it turns out to be more challenging than the first type of problem [17]. Pollack [18] introduces the concept of the k-shortest simple path problem and solves it by modifying Hoffman and Pavley’s method [10], so as to avoid paths containing repeated vertices. Clarke et al. [19] present a branch-and-bound procedure to find the k-shortest simple paths, but this method requires a significant computational effort and storage requirements in the main memory. Sakarovitch [20] first identifies several shortest paths that may contain repeated vertices by using the efficient version of Hoffman and Pavley’s method [10], and then picks up k-shortest simple paths among them. By applying a procedure that partitions a path into two subpaths, Yen [21] develops an efficient algorithm to find the k-shortest simple path. Lawler [22] generalizes a procedure that can reduce the amount of storage required in solving the k-shortest simple path problem. Katoh et al. [23] present an improved version of Yen’s algorithm that solves the problem efficiently in an undirected network. The practical performance of Yen’s algorithm is comparatively analyzed with other k-shortest path algorithms [24, 25]. An implementation of Yen’s algorithm is also studied to improve its practical performance [17, 26]. Hershberger and Suri [27, 28] suggest the efficient replacement path algorithm for finding k shortest simple paths in a directed network. Zeng et al. [29] present a heuristic k-shortest path algorithm that is based on Yen’s algorithm when determining an eco-friendly path that results in minimum carbon dioxide emissions from light-duty vehicles.


## 2.2. Alternative Fuel Refueling Station Location Problems

The maximal covering location and the set-covering location approaches are two main streams of research for addressing AF refueling station location problems. Kuby and Lim [30] are one of the first applying a maximal covering location model to solve an AF refueling station location problem. They introduce the flow-refueling location model (FRLM) to find the optimal location of refueling stations for AF vehicles by considering their limited driving range per refueling with the objective of maximizing the total traffic flow covered. Upchurch and Kuby [31] show that the FRLM identifies more stable locations for AF refueling stations than the p-median model, especially in a statewide case study. In general, the FRLM requires a significant computational effort to pregenerate all feasible location combinations of refueling stations that allow vehicles to make round trips between ODs. Lim and Kuby [32] propose three heuristic versions for the FRLM. Kuby et al. [33] apply two of them to locate hydrogen refueling stations in Florida. Capar and Kuby [34] present a new formulation for the FRLM that skips the pregeneration of all feasible combinations on every path. Capar et al. [35] suggest an arc-cover-path-cover model that focuses on the arcs comprising each path, so as to solve the problem efficiently without the pregeneration of all feasible location combinations of stations on each path for the FRLM. Jochem et al. [36] apply the Capar et al. [35] model to allocate charging stations in the German autobahn. While most AF refueling station location problems assume that the AF refueling stations can only be located at the vertices, Kuby and Lim [37] and Ventura et al. [38] consider additional candidate sites along arcs. Kweon et al. [39] extend the approach suggested in Ventura et al. [38] to locate a refueling station anywhere along a tree network to the case where a portion of drivers are willing to deviate to receive refueling service. He et al. [40] propose a bilevel model to solve the optimal locations of electric charging stations, through taking the driving range limitation of an electric vehicle, the battery charging time required, and the situation in which some electric vehicle drivers possibly charge at home into account.While Kuby and Lim [30] and the following studies solve the AF refueling station location problem using a maximal covering location problem, Wang and Lin [41] solve this problem using a set-covering model with the objective of minimizing the total cost of locating stations to cover all the traffic flow on a given transportation network. Wang and Wang [42] integrate Wang and Lin’s model [41] into the classic set-covering model considering vertex-based and flow-based demands for the AF refueling station location problem. Since Wang and Lin’s model [41] requires a significant computational effort to evaluate the effect of the limited vehicle driving range on the number of charging stations needed for achieving multiple origin-destination intercity travel with electric vehicles on Taiwan, MirHassani and Ebrazi [43] propose a novel approach by using the conservation of flow law, which is able to solve large-scale problems. Chung and Kwon [44] extend MirHassani and Ebrazi’s model [43] to a multiperiod planning problem for allocating charging stations in the Korea Expressway. Hosseini and MirHassani [45] use the idea of MirHassani and Ebrazi’s model [43] to propose a two-stage stochastic mixed integer programming model for the refueling station location problem, where the traffic flow of AF vehicles is uncertain and portable AF refueling stations are considered. Kang and Recker [46] use the idea of the set-covering problem to locate hydrogen refueling stations with the assumption that at most one refueling stop per trip is required in a city. Assuming that vehicles only require one refueling stop per trip, two refueling station location models, including the versions that consider limited capacity of refueling stations [47], and refueling demand uncertainty and driver route choice behavior [48], are developed to minimize the total cost imposed on a planner and drivers over multiple time periods. Using the Adaptive Large Neighborhood Search algorithm [49] and the Adaptive Variable Neighborhood Search algorithm [50], locations for battery swap stations and electric vehicle routes are determined to provide services with the objective of minimizing the sum of the station construction cost and routing cost. To minimize the total cost to locate electric vehicle charging stations in road networks, Gagarin and Corcoran [51] suggest a novel approach that searches for the dominating set of locations among the candidate locations whose distance is below a certain threshold from a given driver. Using a parallel computing strategy, Tran et al. [52] propose an efficient heuristic algorithm for location of AF refueling stations based on the solution of a sequence of subproblems.In general, drivers often deviate from their original paths of shortest travel time or distance to be able to refuel their vehicles [53]. Both the maximal covering location and the set-covering location approaches have been extended to consider driver’s deviation options under a variety of situations. Kim and Kuby [54] address the deviation version of the FRLM (DFRLM) where drivers are able to deviate from their paths of shortest length between ODs, and Kim and Kuby [55] propose a heuristic algorithm for the DFRLM to solve large-scale problems. Huang et al. [56] develop a new model, called the multipath refueling location model (MPRLM), by considering multiple deviation paths between ODs. For intracity trips that require at most one refueling stop, Miralinaghi et al. [47, 48] suggest deviation versions of the set-covering location problem to find potential locations of AF refueling stations. For intercity trips of large-scale problems, Yıldız et al. [57] use a branch and price algorithm, which does not need pregeneration of path generation, for the AF refueling station location problem adding the routing aspect of drivers. Most recently, Göpfert and Bock [58] and Arslan et al. [59] suggest novel projection and branch and cut methods in dealing with the deviation version of the refueling infrastructure planning, so as to extend the computational efficiency even further to solve large-size problem instances with less computational effort.


## 2.3. Main Distinctions of Our Research Work

Comparing with the literature listed above, we shortly provide the main distinctions of our research work over the existing studies as follows:(i)
Our problem is an extension of Hwang et al. [6] problem to consider potential deviation paths on directed transportation networks, such as highway network systems. This leads to consider (1) the mixed set of single-access and dual-access candidate sites to locate AF refueling stations, and (2) nonsymmetric round trips between ODs, where return paths are allowed to be different from original paths for refueling services in both directions.(ii)
Our study is well-suited for the AF refueling station location problem, specially with LNG vehicles traveling long-distance intercity round trips. Some number of recent studies, including Miralinaghi et al. [47, 48], limits this type of problem only for AF vehicles with intracity trips, which needs at most one refueling stop per trip. On the other hand, we apply covering condition procedures depending on LNG vehicle driving range, so as for LNG vehicles traveling intercity trips to allow multiple refueling stops at different locations.(iii)
In our research work, paths are not fixed for every OD pair. Instead, paths are flexible to consider detour traffic flows. While Kim and Kuby [54, 55] and Huang et al. [56] apply Hoffman and Pavley’s algorithm [10] and Yen’s algorithm [21] to take deviation paths into account, we develop a new algorithm based on Eppstein’s algorithm [12] to rigorously and efficiently find k shortest paths allowing repeated vertices along paths within the tolerance (i.e., maximum deviation distance) of the driver. This helps reduce the computational effort in solving our proposed mixed-integer programming model.(iv)
We conduct computational experiments on a classical 25-vertex network (with 300 OD pairs) and a large network at Pennsylvania state (with 2,211 OD pairs) (1) to evaluate the performance of our research work; and (2) to analyze the effects of number of alternative paths including deviation paths, maximum deviation distance, and vehicle driving range on the optimal location as well as the corresponding total traffic flow covered.The major differences between the proposed model and the existing studies that are directly relevant to ours are summarized in Table2.Table 2
Comparative summary between the proposed model and the existing studies directly relevant to our research work.
Research work
Proposed model
Kim and Kuby [54]
Kim and Kuby [55]
Huang et al. [56]
Miralinaghi et al. [47]
Miralinaghi et al. [48]
Hwang et al. [6]
Hwang et al. [7]
Ventura et al. [5]
Deviation
Allowed
Allowed
Allowed
Allowed
Allowed
Allowed
Not allowed
Not allowed
Not allowed
Method to find deviation paths
Eppstein’s algorithm [12] to find k-shortest paths allowing repeated vertices along paths
Hoffman and Pavley’s algorithm [10] to find k-shortest paths allowing repeated vertices along paths
Greedy-adding and Greedy-adding with substitution heuristics to findk-shortest paths allowing repeated vertices along paths
Yen’s algorithm [21] to find k-shortest simple paths which do not allow repeated vertices along paths
The link-based user equilibrium condition for network traffic to find a deviation path containing one refueling stop
The link-based user equilibrium framework for network traffic to find a deviation path containing one refueling stop
N/A
N/A
N/A
Directed network?
Directed
Undirected
Undirected
Undirected
Undirected
Undirected
Directed
Directed
Directed
Candidate sites
Dual- and single-access
Dual-access
Dual-access
Dual-access
Dual-access
Dual-access
Dual- and single-access
Dual- and single-access
Dual- and single-access
Multiple refueling stops or one refueling stop
Multiple refueling stops
Multiple refueling stops
Multiple refueling stops
Multiple refueling stops
One refueling stop
One refueling stop
Multiple refueling stops
Multiple refueling stops
Multiple refueling stops
Main constraints
(i) Number of stations
(i) Number of stations
(i) Number of stations
(i) Satisfy all travel demand
(i) Satisfy all multiperiod travel demand
(i) Satisfy all travel demand when it is uncertain
(i) Number of stations
(i) Number of stations
(i) Number of stations
(ii) Number of deviation paths
(ii) Vehicle driving range
(ii) Vehicle driving range
(ii) Vehicle driving range
(ii) Capacity of stations
(ii) Vehicle driving range
(ii) Vehicle driving range
(ii) Vehicle driving range
(iii) Deviation factor
(iii) Flow volume decay with the increase of deviation
(iii) Flow volume decay with the increase of deviation
(iii) Flow
(iii) Multiple vehicle classes
(iii) Different capital cost to build stations by region
(iv) Vehicle driving range
(iv) Different fuel tank levels at OD pairs
(v) Different capital cost to build stations by type and/or region
Objective(s)
Maximize the total traffic flow covered, considering multiple paths including shortest and deviation paths for each OD pair
Maximize the total flows refueled on deviation paths
Maximize the total flow refueled when deviation is possible
Minimize the total cost of locating AF stations when multiple paths including shortest paths and deviation paths are considered between an OD pair
Minimize the sum of the total station construction and station operational costs for the multiperiod travel demand
Minimize the sum of the total annualized construction costs with the cost of the total system travel time when travel demand is uncertain
Maximize the total traffic flow covered by the stations
Maximize the total traffic flow covered by the stations on a multiclass vehicle transportation network with different fuel tank levels
(i) Maximize the total vehicle-miles traveled per year covered by the stations
(ii) Minimize the capital cost for constructing the refueling infrastructure
Applied to
(i) 25-vertex network
(i) 25-vertex network
(i) 25-vertex network
(i) 25-vertex network
(i) Sioux-Falls network
(i) Sioux-Falls network
(i) Pennsylvania Turnpike
(i) Pennsylvania Turnpike
(i) Pennsylvania Turnpike
(ii) The state of Pennsylvania
(ii) The state of Florida
(ii) Sioux-Falls network
(ii) Mashhad road network
Nonsymmetric round trips?
Allowed
Not allowed
Not allowed
Not allowed
Not allowed
Not allowed
Not allowed
Not allowed
Not allowed


## 3. The AF Refueling Station Location Problem with Detour Traffic on a Highway Road System

In this section, we first introduce the AF refueling station location problem with detour traffic flows on a highway road system, where AF vehicles are able to make detours for refueling and to select different paths between original and return trips. Next, we provide three small instances to describe concepts of feasible paths on a directed simple network in this problem. Also, based on the examples, we derive four properties of feasible paths. Then, an algorithm is presented to determine multiple paths for which the four properties are satisfied for all ODs. This paper aims to locate a given number of AF refueling stations on a directed transportation network so as to maximize the coverage of detour traffic flows by considering multiple paths between ODs. To this end, a mixed-integer programming model will be presented in Section4.


### 3.1. Problem Statement

We define the problem on a simple directed networkGV,E, where V=vi|i=1,⋯,nV  is the set of vertices for ODs, E=vi,vj|forsomevi,vj∈V is the set of arcs having nonnegative lengths, V=nV, and E=nE. The road network in this problem has neither loops nor multiple parallel road segments in the same direction between any pair of vertices. Let Q=qvi,vjvi,vj∈V,i<j} be the set of OD pairs. For any qvi,vj∈Q, vehicles perform the same round trip, which can be divided into the original path Pvi,vj from origin vi to destination vj, and the return path Pvj,vi from vj to vi. For convenience, Lo=Pvi,vj|forqvi,vj∈Q and Lr=Pvj,vi|forqvi,vj∈Q are defined as the sets of original and return paths, respectively. Let VPvi,vj be the set of constituent vertices in path Pvi,vj; then, each path can be represented as a sequence of vertices such that vi=v0,v1,⋯,vj=vw, where vr, r=0,1,⋯,w, is the r-th vertex in VPvi,vj and VPvi,vj=w+1 is the number of vertices in Pvi,vj. Also, the length of Pvi,vj is calculated as dPvi,vj=∑r=0w−1dvr,vr+1, where dvr,vr+1 is the arc length from the r-th vertex to the next vertex in Pvi,vj. Similarly, dPvj,vi is defined as the length of the return path. If dPvi,vj and dPvj,vi are the shortest distances from vi to vj and from vj to vi, respectively, then the corresponding original and return paths are denoted as PSvi,vj and PSvj,vi. Also, we define Fkvi,vj and Fkvj,vi as the k-th shortest feasible original and return paths for OD pair qvi,vj. In Subsections 3.2 and 3.3, we will discuss in detail the properties of feasible paths.Next, the locations of candidate sites for AF refueling stations, denoted asM=mt|t=1,⋯,nM, are assumed to be predetermined, where some candidate sites can only be accessed from one side of the road (i.e., single-access) and the rest can be accessed from both sides of the road (i.e., dual-access). We define CSFkvi,vj=m0,⋯,mf as the sequence of candidate sites in Fkvi,vj, where m0,⋯,mf∈M, such that 0≤f≤nM−1. In order to represent the distances between vi=v0, vj=vw, and site ml∈CSFkvi,vj, such that 0≤l≤f, dijkv0,ml denotes the distance from v0 to ml in Fkvi,vj. If vs,vs+1 is the arc that passes through ml in CSFkvi,vj, such that 0≤s<w, then dijkv0,ml is calculated as ∑r=0s−1dvr,vr+1+dvs,ml, where vs is the vertex adjacent to ml and dvs,ml is the distance from vs to ml. Similarly, dijkml,vw, which refers to the distance from ml to vw, is calculated as dml,vs+1+∑r=s+1w−1dvr,vr+1, where dml,vs+1 indicates the distance from ml to vs+1. Furthermore, the distance between two candidate sites ml and ml′, where l≠l′, is denoted as dijkml,ml′. When we select the two arcs, vs,vs+1 and vs′,vs′+1, which include ml and ml′, respectively, such that 0≤s≤s′<w, dijkml,ml′=∑r=ss′dvr,vr+1−dvs,ml−dml′,vs′+1. Figure 1 shows a representation of vertex and candidate site sequences, as well as distances between vertices and candidate sites.Figure 1
A sequence of vertices inFkvi,vj with a sequence of candidate sites and distances between vertices and candidate sites.We consider that vehicles make a complete round trip between their ODs. They have a limited driving range under free flow conditions, denoted asR, which refers to the maximum travel distance with a single refueling. Since vehicles’ home locations or final destinations are generally far away from highway interchanges, we assume they have at least a half-full tank when they enter and exit the road network. From now on, we call this assumption the half-full tank assumption. The half-full tank assumption was first introduced by Kuby and Lim [30], considering that data about the actual fuel tank level of vehicles when entering and exiting a highway road system are hard to obtain or likely to be inaccurately estimated. This assumption has been then followed by the existing literature. The half-full tank assumption originally aims at ensuring that vehicles can repeat round trips several times without running out of fuel during their round trip. That is, a vehicle accessing the last refueling station and reaching the origin/destination with at least a half-full tank is able to make the original/return trip from the origin/destination with a half-full tank and access the same station without running out of fuel. The half-full tank assumption in our study similarly makes vehicles refueled at stations positioned within a distance of R/2 from their origin interchanges and again at stations placed within the same distance from their destination interchanges in both the original and return trips. In this respect, we define SSFkvi,vj, ESFkvi,vj, and ISFkvi,vj as the following three sets of candidate sites in CSFkvi,vj, which are categorized depending on their distances from vi and vj:(1)SSFkvi,vj=ml∈CSFkvi,vjdijkv0,ml≤R2,(2)ESFkvi,vj=ml∈CSFkvi,vjdijkml,vw≤R2,(3)ISFkvi,vj=ml∈CSFkvi,vj  dijkv0,ml>R2anddijkml,vw>R2.Note that we considerISFkvi,vj only for Fkvi,vj such that dFkvi,vj>R. For return paths, we similarly define CSFkvj,vi, SSFkvj,vi, ESFkvj,vi, and ISFkvj,vi. When a feasible set of AF refueling stations is located in SSFkvi,vj, ESFkvi,vj, ISFkvi,vj, SSFkvj,vi, ESFkvj,vi, and ISFkvj,vi, the corresponding paths can be covered by the refueling stations. In Section 4, we will discuss in detail the covering conditions that depend on path lengths. The half-full tank assumption can be relaxed using Hwang et al.’s [7] model to consider different fuel tank levels of vehicles at their origins and destinations when detours are available on directed transportation networks, but it is left for future research to mainly focus on addressing the refueling station location problem with deviation options on a highway road system in this study.In general, drivers deviate from their preferred paths, e.g., the least time or shortest distance paths, in as short a distance as possible if the preferred paths do not offer refueling availability. Also, a GPS navigation system provides a certain number of routes to travelers who detour from their familiar paths. In this respect, we first consider that vehicles can deviate from their shortest path up to a maximum deviation distance, which is calculated by multiplying the length of the shortest path by a positive deviation factor, denoted asβ. Then, a predetermined maximum number of paths, denoted as K, whose length does not exceed the maximum deviation distance is nominated for each driving direction. This implies that, if a shortest distance path does not provide refueling service, then drivers can select up to K alternative paths for the same OD pair depending on their refueling availability and the limited additional deviation travel distance. Table 3 summarizes the relevant notation and parameters. Note that data about the actual values of deviation factor β and maximum number of paths K for each OD pair are difficult to obtain or likely to be inaccurately predicted. Besides, different drivers would have different standards for the proper values of β and K, as well as different vehicle driving ranges R. Thus, this study does not fix their values. Instead, in Section 5 we change their values for a given number of refueling stations to demonstrate the coupled effects of these parameters on the coverage of OD traffic flows.Table 3
Definition of notations and parameters.
G
Simple directed network
V
Set of vertices for OD interchange pairs
E
Set of arcs
Q
Set of OD interchange pairs
P
v
i
,
v
j
,
P
v
j
,
v
i
Original and return paths between interchangesvi and vj
L
o
,
L
r
Sets of original and return paths
V
P
v
i
,
v
j
Set of constituent vertices in pathPvi,vj
d
P
v
i
,
v
j
Length ofPvi,vj
P
S
v
i
,
v
j
Shortest original path from interchangevi to interchange vj
F
k
v
i
,
v
j
k-th shortest feasible original path from interchange vi to interchange vj
M
Set of all candidate sites for AF refueling stations
C
S
F
k
v
i
,
v
j
Sequence of candidate sites inFkvi,vj
d
i
j
k
v
r
,
m
l
Distance from ther-th interchange vr to the l-th candidate site ml in Fkvi,vj
d
i
j
k
m
l
,
m
l
′
Distance between two candidate sitesml and ml′ in Fkvi,vj
S
S
F
k
v
i
,
v
j
Set of candidate sites that are located within a distance ofR/2 from interchange vi in Fkvi,vj
E
S
F
k
v
i
,
v
j
Set of candidate sites that are located within a distance ofR/2 from interchange vj in Fkvi,vj
I
S
F
k
v
i
,
v
j
Set of candidate sites that are located beyond a distance ofR/2 from both interchanges vi and vj in Fkvi,vj
R
Vehicle’s maximum travel distance with a single refueling
β
Positive deviation factor
K
Predetermined maximum number of paths


### 3.2. Three Small Instances

The concept of feasible path on a directed transportation network is illustrated with three small disconnected networks with seventeen verticesv1,⋯,v17, two dual-access candidate sites m1,m3, and five single-access candidate sites m2,m4,m5,m6,m7 in Figure 2. Suppose we are trying to determine whether vehicles with R=100 can perform round trips without running out of fuel. First, since vehicles are considered to have at least a half-full tank at theirs ODs, for qv2,v4, shortest paths PSv2,v4 and PSv4,v2 are feasible because the corresponding trips can be covered by placing a single refueling station at site m1. However, for qv1,v3, since candidate sites are available in neither PSv1,v3 nor PSv3,v1, drivers need to deviate from these paths to receive refueling service in their round trips. If a refueling station is located at site m1, then vertex sequence v1=v0,v2=v1,v4=v2,v2=v3,v3=v4 for the original path and opposite sequence v3=v0,v2=v1,v4=v2,v2=v3,v1=v4 for the return path are feasible paths for qv1,v3. Second, suppose that there are two AF refueling stations available at sites m3 and m4. Then, the shortest original path of qv6,v11 does not offer AF refueling availability because vehicles cannot reach v11 with at least a half-full tank when they exit the network. For the shortest return path, vehicles that reenter the network with a half-full tank at v11 cannot reach the AF refueling station at site m3. Thus, vehicles in qv6,v11 would need to make detours using vertex sequence v6=v0,v7=v1,v8=v2,v9=v3,v10=v4,v7=v5,v8=v6,v11=v7 for the original path and opposite sequence v11=v0,v8=v1,v9=v2,v10=v3,v7=v4,v6=v5 for the return path. These two deviation paths are feasible paths for qv6,v11. Lastly, suppose that sites m5 and m6 are selected as AF refueling stations. Vehicles in qv12,v17 cannot be refueled by these two AF refueling stations if they use the shortest original and return paths. However, if vehicles make multiple cycles at vertex v13, then vertex sequence v12=v0,v13=v1,v14=v2,v13=v3,v15=v4,v13=v5,v17=v6 becomes a feasible original path and vertex sequence v17=v0,v15=v1,v13=v2,v12=v3 turn out to be a feasible return path for qv12,v17. Note that another vertex sequence v17=v0,v13=v1,v14=v2,v13=v3,v15=v4,v13=v5,v12=v6 also provides AF refueling availability for the return path, but drivers would not choose this vertex sequence because this path is longer than the path defined by the previous vertex sequence v17=v0,v15=v1,v13=v2,v12=v3. From these observations, we can state four feasible path properties in the next subsection.Figure 2
Three detour examples of a directed transportation network.


### 3.3. Four Properties of Feasible Paths

We derive four properties of feasible paths on a directed transportation network. Based on these four properties, an existingk-th shortest path algorithm will be modified in the next subsection to generate multiple feasible paths for all OD pairs. Since we consider round trips, the four properties are applicable to both feasible original and return paths. The first property represents that feasible paths have a proper sequence of candidate sites to be selected for AF refueling stations which cover trips in feasible paths.Property 1.
LetFkvi,vj be the k-th shortest feasible original path from vi to vj that contains feasible candidate sites for AF refueling stations to be able to cover trips in Pvi,vj, and CSFkvi,vj=m0,⋯,mf denote the sequence of candidate sites in Fkvi,vj, where m0,⋯,mf∈M. According to the cardinality of CSFkvi,vj, denoted as CSFkvi,vj, one of the following two cases must be satisfied:(a) IfCSFkvi,vj=1, then SSFkvi,vj∩ESFkvi,vj≠∅.
(b) IfCSFkvi,vj>1, then SSFkvi,vj≠∅, ESFkvi,vj≠∅, and dijkml,ml+1≤R for any two adjacent candidate sites ml,ml+1∈CSFkvi,vj.Proof.
SinceFkvi,vj is a feasible path, the half-full tank assumption implies that CSFkvi,vj≥1. In case (a), there is a single candidate site available for a refueling station on Fkvi,vj, and by the half-full tank assumption, this site is located within a distance of R/2 from vi and vj to ensure that vehicles on Fkvi,vj make a successful trip from vi to vj. Then, by definition of SSFkvi,vj and ESFkvi,vj, this site should be located in the common segment in SSFkvi,vj∩ESFkvi,vj. This implies that SSFkvi,vj∩ESFkvi,vj≠∅.
In case (b), there are at least two candidate sites available onFkvi,vj. If either SSFkvi,vj=∅ or ESFkvi,vj=∅, then Fkvi,vj cannot offer any AF refueling service. Thus, SSFkvi,vj≠∅ and ESFkvi,vj≠∅. For any two adjacent candidate sites ml,ml+1∈CSFkvi,vj, if dijkml,ml+1>R, then vehicles at site ml cannot reach site ml+1, i.e., Fkvi,vj does not have AF refueling availability. It also contradicts the definition of Fkvi,vj. Thus, dijkml,ml+1≤R for any two adjacent candidate sites ml,ml+1∈CSFkvi,vj.In practice, vehicles do not need to visit the same refueling station located at a single-access site several times. For example, let us consider two return path alternatives ofqv14,v16 in Figure 2. While a vertex sequence v16=v0,v15=v1,v13=v2,v17=v3,v15=v4,v13=v5,v14=v6 goes over site m6 twice, vehicles can also use a vertex sequence v16=v0,v15=v1,v13=v2,v14=v3 to travel from v16 to v14 with a single visit to site m6. That is, a shorter feasible path is always preferred to a longer feasible path that goes through a single-access candidate multiple times. On the other hand, it is unavoidable for vehicles to revisit the same refueling station when the station is placed at a dual-access candidate site. From one of the original paths of qv1,v3 in Figure 2, a feasible vertex sequence v1=v0,v2=v1,v4=v2,v2=v3,v3=v4 passes through a refueling station at dual-access site m1 twice, but in opposite directions. In this respect, prior to checking whether paths go over candidate sites multiple times, we treat a dual-access candidate site as two distinct single-access sites. From this observation, we derive the second property for candidate sites of feasible paths.Property 2.
Forqvi,vj∈Q, suppose that Fk1vi,vj is an original path that satisfies Property 1 and passes through a single-access candidate site mt1∈M twice:(4)
C
S
F
k
1
v
i
,
v
j
=
m
0
,
⋯
m
a
,
m
a
+
1
⋯
,
m
b
,
m
b
+
1
,
⋯
,
m
f
,
w
h
e
r
e
m
t
1
=
m
a
=
m
b
.Then, there must exist an original path, denoted asFk2vi,vj, that satisfies Property 1 and goes over site mt1 only one time, such that CSFk2vi,vj<CSFk1vi,vj:(5)CSFk2vi,vj=m0,⋯ma,mb+1,⋯,mf.Thus,dFk2vi,vj<dFk1vi,vj.Proof.
Letma and mb be the a-th and b-th candidate sites in CSFk1vi,vj=m0,⋯ma,ma+1⋯,mb,mb+1,⋯,mf which are located at the same site mt1 such that 0≤a<b≤f, where f=CSFk1vi,vj−1. By eliminating all consecutive candidate sites from ma+1 to mb and the corresponding arcs in Fk1vi,vj, we can construct path Fk2vi,vj with CSFk2vi,vj=m0,⋯ma,mb+1,⋯,mf. Then, from Property 1, since dijk1mb,mb+1≤R, vehicles at mb can reach mb+1 in Fk1vi,vj. This implies that vehicles at ma can reach mb+1 in Fk2vi,vj because ma and mb are placed at the same site mt1. Thus, CSFk2vi,vj is a feasible sequence of candidate sites for Fk2vi,vj and dFk2vi,vj<dFk1vi,vj.When drivers decide to make detours from their preplanned paths for refueling, they would be reluctant to travel along unnecessary paths to reach available refueling stations. For instance, let us take an original path ofqv1,v3 from the previous network in Figure 2. Recall that a vertex sequence v1=v0,v2=v1,v4=v2,v2=v3,v3=v4 is one of the feasible original paths for qv1,v3. Now, let us consider a different vertex sequence defined by v1=v0,v2=v1,v4=v2,v5=v3,v4=v4,v2=v5,v3=v6. The candidate site sequence for this second path is m1=m0,m2=m1,m1=m2, which enables successful trips from v1 to v3 without running out of fuel. However, drivers would not travel further along the subpath of the second path, v4=v2,v5=v3,v4=v4 because they can reach v3 without refueling at site m2. They would return at v4 after visiting the refueling station at site m1, i.e., they would use the first feasible path, v1=v0,v2=v1,v4=v2,v2=v3,v3=v4, instead of the second one for the original path of qv1,v3. Similarly, drivers can consider two possible paths for the return path of qv12,v17. While vehicles using a vertex sequence v17=v0,v13=v1,v14=v2,v13=v3,v15=v4,v13=v5,v12=v6 require each refueling stop at both sites m5 and m6, vehicles in a vertex sequence v17=v0,v15=v1,v13=v2,v12=v3 need only one refueling stop at site m6. Since the length of the latter path is shorter, vehicles in the return path of qv12,v17 would use the second vertex sequence instead of the first one. Based on these observations, we conclude that drivers would not make detours along a vertex sequence which requires unnecessary travel to reach available refueling stations. From now on, we call a site in unnecessary subpaths an irrelevant candidate site, so that feasible paths can be represented by vertex sequences without irrelevant candidate sites. Property 3 provides conditions under which there exist irrelevant candidate sites in paths.Property 3.
Given two paths,Fk1vi,vj and Fk2vi,vj, such that dFk1vi,vj<dFk2vi,vj, which satisfy Properties 1 and 2, Fk2vi,vj includes irrelevant candidate sites in set CSFk2vi,vj\CSFk1vi,vj when one of the following conditions is satisfied:(a)
IfCSFk1vi,vj=1, then (1) CSFk1vi,vj⊂CSFk2vi,vj, and (2) either SFk1vi,vj=SSFk2vi,vj or ESFk1vi,vj=ESFk2vi,vj.
(b)
IfCSFk1vi,vj>1, then (1) CSFk1vi,vj⊂CSFk2vi,vj, (2) SSFk1vi,vj=SSFk2vi,vj, (3) ESFk1vi,vj=ESFk2vi,vj, and (4) 0<dijk1m0,mf≤R, where CSFk1vi,vj=m0,⋯,mf.Proof.
In case (a), suppose thatm0∈SSFk1vi,vj=SSFk2vi,vj. By Property 1, m0∈ESFk1vi,vj, which implies that dijk1m0,vj≤R/2. Thus, any candidate site in CSFk2vi,vj\CSFk1vi,vj is irrelevant because vehicles covered by a single refueling station at site m0 can complete trips from vi to vj. Similarly, it can be shown when m0∈ESFk1vi,vj=ESFk2vi,vj.
In case (b), by Condition (4),dijk1ma,mb≤R such that ma∈SSFk1vi,vj, mb∈ESFk1vi,vj, and 0≤a<b≤f. This implies that any additional refueling stop is unnecessary between ma and mb. Then, by Conditions (1), (2), and (3), ISFk2vi,vj≠∅, which implies that Fk2vi,vj includes irrelevant candidate sites.Given two feasible paths,Fk1vi,vj and Fk2vi,vj such that dFk1vi,vj≤dFk2vi,vj, for the same OD pair qvi,vj∈Q, if any subset of candidate sites that covers the trips in Fk2vi,vj is also able to cover the trips in Fk1vi,vj, then we say that Fk2vi,vj is dominated by Fk1vi,vj. This implies that Fk2vi,vj is unnecessary to be considered in our proposed model. In practice, drivers would use only Fk1vi,vj because dFk1vi,vj≤dFk2vi,vj. For example, let us take the previous example of qv6,v11 in Figure 2. Recall that F1v6,v11=v6=v0,v7=v1,v8=v2,v9=v3,v10=v4,v7=v5,v8=v6,v11=v7. Then, we can find another feasible original path, which is F2v6,v11=v6=v0,v7=v1,v10=v2,v9=v3,v10=v4,v7=v5,v8=v6,v11=v7 for vehicles with R=100. They have the same sequence of candidate sites, but vehicles may prefer F1v6,v11 to F2v6,v11 because dF1v6,v11=125<dF2v6,v11=145. Furthermore, in order to clarify the concept of dominated feasible paths, let us see another example of qv14,v16 in Figure 2. For the return path, we have two feasible return paths, F1v16,v14=v16=v0,v15=v1,v13=v2,v14=v3 and F2v16,v14=v16=v0,v15=v1,v17=v2,v13=v3,v14=v4, such that CSF1v16,v14=m7=m0,m6=m1,m5=m2 and CSF2v16,v14=m7=m0,m5=m1. For F1v16,v14, we have five possible sets of candidate sites to cover the trip, i.e., m6, m6,m5, m7,m6, m7,m5, and m7,m6,m5. In case of F2v16,v14, m7,m5 is the only feasible set. That is, CSF2v16,v14⊆CSF1v16,v14. It means that F2v16,v14 is dominated by F1v16,v14 because dF1v16,v14=65<dF2v16,v14=100. In this respect, Property 4 shows that the shortest feasible paths can dominate other feasible paths when the candidate sites of the shortest feasible paths include the candidate sites of the other paths.Property 4.
Givenqvi,vj∈Q, let Fk1vi,vj be the shortest feasible original path with CSFk1vi,vj. Suppose that there exists another feasible original path, denoted as Fk2vi,vj, such that CSFk2vi,vj⊆CSFk1vi,vj and dFk1vi,vj≤dFk2vi,vj. Then, Fk2vi,vj is dominated by Fk1vi,vj.Proof.
By definition of dominated paths, we need to show that any subset of candidate sites to cover the trips inFk2vi,vj is also able to cover the trips in Fk1vi,vj. Suppose that FSFk2vi,vj is a set of candidate sites that covers trips in Fk2vi,vj. This implies that FSFk2vi,vj⊆CSFk2vi,vj, and therefore, FSFk2vi,vj⊆CSFk1vi,vj because CSFk2vi,vj⊆CSFk1vi,vj. If FSFk2vi,vj=1, then it is trivial to show that FSFk2vi,vj can cover trips in Fk1vi,vj because dFk1vi,vj≤dFk2vi,vj. If FSFk2vi,vj>1, then we select the first and last sequences of candidate sites in FSFk2vi,vj, denoted as ma and mb, respectively. When we decompose Fk1vi,vj into vi=v0,⋯,ma, ma,⋯,mb, and mb,⋯,vj=vw, we have that dFk1vi,vj=dijk1v0,ma+dijk1ma,mb+dijk1mb,vw. If dijk2v0,ma≤R/2<dijk1v0,ma, then we can construct a shorter feasible path whose length is dijk2v0,ma+dijk1ma,mb+dijk1mb,vw. It contradicts that Fk1vi,vj is the shortest feasible path. Thus, dijk1v0,ma≤R/2. Similarly, it can be shown that dijk1mb,vw≤R/2. Furthermore, if dijk2ml,ml+1≤R<dijk1ml,ml+1 such that a≤l<b, then it also contradicts that Fk1vi,vj is the shortest feasible path. Thus, FSFk2vi,vj can cover trips in Fk1vi,vj by Property 1.


### 3.4.k-th Shortest Nondominated Feasible Path Algorithm

In this subsection, we describe thek-th shortest nondominated feasible (SNDF) path algorithm to generate multiple paths that (1) offer AF refueling availability, (2) contain unique single-access candidate sites, (3) do not have irrelevant candidate sites, and (4) are not dominated by other paths for all OD pairs. We extend the k-th shortest path algorithm developed by Eppstein [12] to first construct candidate paths between ODs for SNDF paths, and then eliminate some of the paths that violate the four properties discussed on the previous subsection. Finally, given a deviation factor β and a maximum number K of SNDF paths, the k-th SNDF path algorithm determines up to K SNDF paths from vi to vj and another up to K SNDF paths from vj to vi.The first step of thek-th SNDF path algorithm is to transform a simple directed network GV,E with the set of nV vertices and the set of nE arcs into an expanded directed network, denoted as G′V′,E′, by transforming each dual-access candidate site in G into two distinct single-access candidate sites, one in each side of the road, and setting all candidate sites in G as vertices in G′. Then, the maximum cardinalities of V′ and E′ can be written as V′=nV+nSM+2nDM and E′=nE+2nSM+4nDM, respectively, where nSM and nDM refer to the number of single- and dual-access candidate sites. Next, for interchange vj∈V′, Tvj denotes a shortest path arborescence rooted at vertex vj containing the shortest path to each vi in V′\vj. The shortest path arborescence Tvj for interchange vj in G′V′,E′ can be constructed using Dijkstra’s algorithm in time OE′+V′logV′. Let ETvj denote the set of arcs composing Tvj. If we select arc va,vb∈E′\ETvj and follow the shortest path from vb to vj instead of taking the shortest path in Tvj, then the additional distance from Tvj, denoted as δvjva,vb, is computed as: δvjva,vb=dva,vb+dPSvb,vj−dPSva,vj≥0. Note that δvjva,vb>0 for va,vb∈E′\ETvj and δvjva,vb=0 for va,vb in Tvj. This implies that we can calculate the length of path Pvi,vj by adding the values of δvjva,vb for the arcs in Pvi,vj to the length of PSvi,vj since δvjva,vb=0 for va,vb in PSvi,vj. These arcs aside from PSvi,vj do not belong to Tvj, and Eppstein’s algorithm calls these arcs the sidetracks of PSvi,vj. The sidetracks are used to generate the k-shortest path of Pvi,vj. For Pvi,vj, a special tree-based data structure, called the min-heap data structure, in which it is a completely binary tree with data structure in a nondecreasing order of values of sidetracks, is applied to build K shortest paths of Pvi,vj using the depth first search algorithm in time OKlogK. The k-th shortest path of Pvi,vj, for k=1,⋯,K, generated through the min-heap data structure, is called a candidate path from vi to vj and is denoted as F′vi,vj.When generatingF′vi,vj, we determine whether or not F′vi,vj is included in the set of SNDF paths for qvi,vj by Properties 1, 2, 3, and 4. First, F′vi,vj needs to have AF refueling availability by Property 1. If either SSF′vi,vj∩ESF′vi,vj≠∅ when CSF′vi,vj=1 or SSF′vi,vj≠∅, ESF′vi,vj≠∅, and dijkml,ml+1≤R for ml,ml+1∈CSF′vi,vj when CSF′vi,vj>1, then F′vi,vj has AF refueling availability. Second, if F′vi,vj repeats some single-access candidate sites, then the path is excluded from the set of SNDF paths because we also construct a similar feasible path without repeating these candidate sites whose length is shorter than dF′vi,vj by Property 2. Next, regarding Properties 3 and 4, we compare CSF′vi,vj with CSFkvi,vj’s when k>1. If CSFkvi,vj⊂CSF′vi,vj for k=1,⋯,K−1, we diagnose whether or not F′vi,vj has irrelevant candidate sites by Property 3. If CSF′vi,vj⊆CSFkvi,vj for k=1,⋯,K−1, then F′vi,vj is removed from the set of SNDF paths for qvi,vj because the associated OD pair already has shortest feasible paths which dominate F′vi,vj by Property 4. Until reaching the maximum number, K, of SNDF paths, Eppstein’s algorithm updates a binary heap that includes sets of sidetracked arcs to find the next paths. However, the k-th SNDF path algorithm can even stop before reaching K if the sum of sidetracked arcs’ δvjva,vb on the root of the binary heap becomes greater than the value of βdPSvi,vj. This represents that the length of a newly constructed path cannot exceed a maximum deviation distance of drivers. After determining at most K SNDF paths from interchanges vi’s to a certain interchange vj, we similarly apply the above procedures to other interchanges vj’s. Regarding 2K SNDF paths in qvi,vj, Algorithm 1 describes the pseudocode of the k-th SNDF path algorithm for all OD pairs.Algorithm 1: Overview of the k-th SNDF path algorithm for constructing at most 2K SNDF paths.
Algorithm 1:k-th SNDF Path Algorithm
Input:GV,E,M,R,β, and K
Output:G′V′,E′, sets of all K SNDF original paths, and K SNDF return paths for all OD pairs
Construct an expanded directed networkG′ by transforming dual-candidate sites into two separate single-candidate sites on each direction and adding all sites to G as vertices.
Setj=1.
Repeat
Build a shortest path arborescenceTvj rooted at vertex vj.
Seti=1,i≠j.
Repeat
Setk=0.
Construct a min binary heap data structure, denoted asHvi,vj, with sequences of sidetracked arcs using the corresponding values of δvjva,vb, ∀va,vb∈EG′, through Eppstein’s algorithm. Note that the binary heap initially has the empty sequence of sidetracked arcs because the empty sequence of sidetrack arcs represents PSvi,vj.
Repeat
Construct a candidate pathF′vi,vj by inserting a sequence of sidetracked arcs on the root of Hvi,vj into Tvj.
ConstructCSF′vi,vj,SSF′vi,vj,ESF′vi,vj,andISF′vi,vj.
IfF′vi,vj is satisfied with Properties 1, 2, 3, and 4, then define F′vi,vj, SSF′vi,vj, ESF′vi,vj and ISF′vi,vj as Fkvi,vj, SSFkvi,vj, ESFkvi,vj, and ISFkvi,vj, respectively, and set k=k+1.
Remove the current sequence of sidetracked arcs on the root ofHvi,vj, and insert new sequences of sidetracked arcs from the min-heap data structure to the binary heap Hvi,vj by maintaining the min-heap property.
Untilk>K, ∑va,vb∈AEPSvi,vjδvjva,vb>β∙dPSvi,vj, or the min-heap data structure does not have any sidetracked arcs.
Seti=i+1.
Untili reaches nV+1.
Setj=j+1.
Until
j reaches nV+1.We prove below that at most2K SNDF paths (i.e., up to K original paths and another up to K SNDF return paths) for each OD pair can be determined by the k-th SNDF path algorithm.Theorem 1.
Given a value for the maximum number of paths,K, the k-th SNDF path algorithm generates at most 2K SNDF paths for each OD pair in G.Proof.
After transformingG into G′, given qvi,vj∈Q, F′vi,vj is generated by Eppstein’s algorithm. The k-th SNDF path algorithm determines CSF′vi,vj, and then classifies the sequence of candidate sites into SSF′vi,vj, ESF′vi,vj, and ISF′vi,vj by measuring distances between the sites, vi and vj. The algorithm checks first whether or not F′vi,vj meets Properties 1 and 2. Next, if the corresponding OD pair already has SNDF paths, then CSF′vi,vj is compared with CSFkvi,vj for k=1,⋯,K−1 by Properties 3 and 4. If CSFkvi,vj⊂CSF′vi,vj for k=1,⋯,K−1, then the algorithm checks whether or not F′vi,vj has irrelevant candidate sites. If CSF′vi,vj⊆CSFkvi,vj for k=1,⋯,K−1, then F′vi,vj is not included in the set of SNDF paths for qvi,vj. Otherwise, the algorithm assigns F′vi,vj to the k-th SNDF path Fkvi,vj and increases the value of k by one. If k≤K, a next temporary path can be constructed by adding new sequences of sidetracked arcs into Tvj. The additional distance of the next path can be calculated as ∑va,vb∈AEPSvi,vjδvjva,vb, where AEPSvi,vj is the set of newly added arcs and δvjva,vb=dva,vb+dPSvb,vj−dPSva,vj. If the additional travel distance does not exceed β·dPSvi,vj, then the above processes are repeated until k>K. Otherwise, the algorithm moves to the next interchange vi in the shortest path arborescence Tvj. After determining at most K SNDF paths, we repeat the algorithm by building a new shortest path arborescence for other interchanges. Therefore, for each OD pair, this algorithm determines up to K SNDF paths from origin to destination and up to K other SNDF paths from destination to origin, which leads to at most 2K SNDF paths.Since the original networkG has single- and dual-access candidate sites for AF refueling stations, vehicles can select different return paths from their SNDF original paths when they come back to their origin interchanges. This means that nonsymmetric round trips can be considered for each OD pair in our proposed model. In practice, drivers would make different detours from their SNDF original paths if the original paths do not provide AF refueling service in the opposite direction. For convenience, L^o=Fkvi,vj|k=1,⋯,K,∀qvi,vj∈Q and L^r=Fkvj,vi|k=1,⋯,K,∀qvi,vj∈Q denote the sets of all K SNDF original paths and K SNDF return paths for all OD pairs, respectively.Theorem 2.
In order to generate2K SNDF original and return paths for all OD pairs, the complexity of the k-th SNDF path algorithm takes time OVE′+V′logV′+VKM2+E′KlogK, where V′=nV+nSM+2nDM and E′=nE+2nSM+4nDM. Recall that V=nV, E=nE, and M=nM. Also, nSM and nDM refer to the number of single- and dual-access candidate sites, respectively.Proof.
TimeOnDM is required to transform GV,E with the set of all candidate sites to a simple expanded directed network G′V′,E′ by regarding dual-access candidate sites as two distinct single-access candidate sites on arcs. Given G′, the basic version of Eppstein’s algorithm takes time OV′+V′logV′+VKlogK to find the k-th shortest path to vj from any other interchanges. In time OE′+V′logV′, the algorithm constructs a shortest path arborescence and a min-heap data structure to output F′vi,vj. Then, in time OKlogK, it generates F′vi,vj’s from a specific interchange to vj. After we construct CSF′vi,vj and classify CSF′vi,vj into SSF′vi,vj, ESF′vi,vj, and ISF′vi,vj, it takes time OM2+E′ for Property 1. Next, we find out which candidate sites are repeated in F′vi,vj in time OM2 for Property 2. Also, when k>1, time OKM2 is required to compare CSF′vi,vj with CSFkvi,vj for Properties 3 and 1. These processes are repeated until either K feasible paths are generated or the length of a newly constructed path does not exceed 1+βdPSvi,vj in time OKM2+E′KlogK. Next, we move to other interchanges vi’s, ∀vi∈VG, i≠j, which requires time OVKM2+E′KlogK. In summary, generating K SNDF paths from other interchanges to a certain interchange vj requires time OE′+V′logV′+VKM2+E′KlogK. Finally, since the number of interchanges for ODs is V, the total computational complexity of the algorithm takes time OVE′+V′logV′+VKM2+E′KlogK in order to generate 2K SNDF original and return paths in all OD pairs.


## 3.1. Problem Statement

We define the problem on a simple directed networkGV,E, where V=vi|i=1,⋯,nV  is the set of vertices for ODs, E=vi,vj|forsomevi,vj∈V is the set of arcs having nonnegative lengths, V=nV, and E=nE. The road network in this problem has neither loops nor multiple parallel road segments in the same direction between any pair of vertices. Let Q=qvi,vjvi,vj∈V,i<j} be the set of OD pairs. For any qvi,vj∈Q, vehicles perform the same round trip, which can be divided into the original path Pvi,vj from origin vi to destination vj, and the return path Pvj,vi from vj to vi. For convenience, Lo=Pvi,vj|forqvi,vj∈Q and Lr=Pvj,vi|forqvi,vj∈Q are defined as the sets of original and return paths, respectively. Let VPvi,vj be the set of constituent vertices in path Pvi,vj; then, each path can be represented as a sequence of vertices such that vi=v0,v1,⋯,vj=vw, where vr, r=0,1,⋯,w, is the r-th vertex in VPvi,vj and VPvi,vj=w+1 is the number of vertices in Pvi,vj. Also, the length of Pvi,vj is calculated as dPvi,vj=∑r=0w−1dvr,vr+1, where dvr,vr+1 is the arc length from the r-th vertex to the next vertex in Pvi,vj. Similarly, dPvj,vi is defined as the length of the return path. If dPvi,vj and dPvj,vi are the shortest distances from vi to vj and from vj to vi, respectively, then the corresponding original and return paths are denoted as PSvi,vj and PSvj,vi. Also, we define Fkvi,vj and Fkvj,vi as the k-th shortest feasible original and return paths for OD pair qvi,vj. In Subsections 3.2 and 3.3, we will discuss in detail the properties of feasible paths.Next, the locations of candidate sites for AF refueling stations, denoted asM=mt|t=1,⋯,nM, are assumed to be predetermined, where some candidate sites can only be accessed from one side of the road (i.e., single-access) and the rest can be accessed from both sides of the road (i.e., dual-access). We define CSFkvi,vj=m0,⋯,mf as the sequence of candidate sites in Fkvi,vj, where m0,⋯,mf∈M, such that 0≤f≤nM−1. In order to represent the distances between vi=v0, vj=vw, and site ml∈CSFkvi,vj, such that 0≤l≤f, dijkv0,ml denotes the distance from v0 to ml in Fkvi,vj. If vs,vs+1 is the arc that passes through ml in CSFkvi,vj, such that 0≤s<w, then dijkv0,ml is calculated as ∑r=0s−1dvr,vr+1+dvs,ml, where vs is the vertex adjacent to ml and dvs,ml is the distance from vs to ml. Similarly, dijkml,vw, which refers to the distance from ml to vw, is calculated as dml,vs+1+∑r=s+1w−1dvr,vr+1, where dml,vs+1 indicates the distance from ml to vs+1. Furthermore, the distance between two candidate sites ml and ml′, where l≠l′, is denoted as dijkml,ml′. When we select the two arcs, vs,vs+1 and vs′,vs′+1, which include ml and ml′, respectively, such that 0≤s≤s′<w, dijkml,ml′=∑r=ss′dvr,vr+1−dvs,ml−dml′,vs′+1. Figure 1 shows a representation of vertex and candidate site sequences, as well as distances between vertices and candidate sites.Figure 1
A sequence of vertices inFkvi,vj with a sequence of candidate sites and distances between vertices and candidate sites.We consider that vehicles make a complete round trip between their ODs. They have a limited driving range under free flow conditions, denoted asR, which refers to the maximum travel distance with a single refueling. Since vehicles’ home locations or final destinations are generally far away from highway interchanges, we assume they have at least a half-full tank when they enter and exit the road network. From now on, we call this assumption the half-full tank assumption. The half-full tank assumption was first introduced by Kuby and Lim [30], considering that data about the actual fuel tank level of vehicles when entering and exiting a highway road system are hard to obtain or likely to be inaccurately estimated. This assumption has been then followed by the existing literature. The half-full tank assumption originally aims at ensuring that vehicles can repeat round trips several times without running out of fuel during their round trip. That is, a vehicle accessing the last refueling station and reaching the origin/destination with at least a half-full tank is able to make the original/return trip from the origin/destination with a half-full tank and access the same station without running out of fuel. The half-full tank assumption in our study similarly makes vehicles refueled at stations positioned within a distance of R/2 from their origin interchanges and again at stations placed within the same distance from their destination interchanges in both the original and return trips. In this respect, we define SSFkvi,vj, ESFkvi,vj, and ISFkvi,vj as the following three sets of candidate sites in CSFkvi,vj, which are categorized depending on their distances from vi and vj:(1)SSFkvi,vj=ml∈CSFkvi,vjdijkv0,ml≤R2,(2)ESFkvi,vj=ml∈CSFkvi,vjdijkml,vw≤R2,(3)ISFkvi,vj=ml∈CSFkvi,vj  dijkv0,ml>R2anddijkml,vw>R2.Note that we considerISFkvi,vj only for Fkvi,vj such that dFkvi,vj>R. For return paths, we similarly define CSFkvj,vi, SSFkvj,vi, ESFkvj,vi, and ISFkvj,vi. When a feasible set of AF refueling stations is located in SSFkvi,vj, ESFkvi,vj, ISFkvi,vj, SSFkvj,vi, ESFkvj,vi, and ISFkvj,vi, the corresponding paths can be covered by the refueling stations. In Section 4, we will discuss in detail the covering conditions that depend on path lengths. The half-full tank assumption can be relaxed using Hwang et al.’s [7] model to consider different fuel tank levels of vehicles at their origins and destinations when detours are available on directed transportation networks, but it is left for future research to mainly focus on addressing the refueling station location problem with deviation options on a highway road system in this study.In general, drivers deviate from their preferred paths, e.g., the least time or shortest distance paths, in as short a distance as possible if the preferred paths do not offer refueling availability. Also, a GPS navigation system provides a certain number of routes to travelers who detour from their familiar paths. In this respect, we first consider that vehicles can deviate from their shortest path up to a maximum deviation distance, which is calculated by multiplying the length of the shortest path by a positive deviation factor, denoted asβ. Then, a predetermined maximum number of paths, denoted as K, whose length does not exceed the maximum deviation distance is nominated for each driving direction. This implies that, if a shortest distance path does not provide refueling service, then drivers can select up to K alternative paths for the same OD pair depending on their refueling availability and the limited additional deviation travel distance. Table 3 summarizes the relevant notation and parameters. Note that data about the actual values of deviation factor β and maximum number of paths K for each OD pair are difficult to obtain or likely to be inaccurately predicted. Besides, different drivers would have different standards for the proper values of β and K, as well as different vehicle driving ranges R. Thus, this study does not fix their values. Instead, in Section 5 we change their values for a given number of refueling stations to demonstrate the coupled effects of these parameters on the coverage of OD traffic flows.Table 3
Definition of notations and parameters.
G
Simple directed network
V
Set of vertices for OD interchange pairs
E
Set of arcs
Q
Set of OD interchange pairs
P
v
i
,
v
j
,
P
v
j
,
v
i
Original and return paths between interchangesvi and vj
L
o
,
L
r
Sets of original and return paths
V
P
v
i
,
v
j
Set of constituent vertices in pathPvi,vj
d
P
v
i
,
v
j
Length ofPvi,vj
P
S
v
i
,
v
j
Shortest original path from interchangevi to interchange vj
F
k
v
i
,
v
j
k-th shortest feasible original path from interchange vi to interchange vj
M
Set of all candidate sites for AF refueling stations
C
S
F
k
v
i
,
v
j
Sequence of candidate sites inFkvi,vj
d
i
j
k
v
r
,
m
l
Distance from ther-th interchange vr to the l-th candidate site ml in Fkvi,vj
d
i
j
k
m
l
,
m
l
′
Distance between two candidate sitesml and ml′ in Fkvi,vj
S
S
F
k
v
i
,
v
j
Set of candidate sites that are located within a distance ofR/2 from interchange vi in Fkvi,vj
E
S
F
k
v
i
,
v
j
Set of candidate sites that are located within a distance ofR/2 from interchange vj in Fkvi,vj
I
S
F
k
v
i
,
v
j
Set of candidate sites that are located beyond a distance ofR/2 from both interchanges vi and vj in Fkvi,vj
R
Vehicle’s maximum travel distance with a single refueling
β
Positive deviation factor
K
Predetermined maximum number of paths


## 3.2. Three Small Instances

The concept of feasible path on a directed transportation network is illustrated with three small disconnected networks with seventeen verticesv1,⋯,v17, two dual-access candidate sites m1,m3, and five single-access candidate sites m2,m4,m5,m6,m7 in Figure 2. Suppose we are trying to determine whether vehicles with R=100 can perform round trips without running out of fuel. First, since vehicles are considered to have at least a half-full tank at theirs ODs, for qv2,v4, shortest paths PSv2,v4 and PSv4,v2 are feasible because the corresponding trips can be covered by placing a single refueling station at site m1. However, for qv1,v3, since candidate sites are available in neither PSv1,v3 nor PSv3,v1, drivers need to deviate from these paths to receive refueling service in their round trips. If a refueling station is located at site m1, then vertex sequence v1=v0,v2=v1,v4=v2,v2=v3,v3=v4 for the original path and opposite sequence v3=v0,v2=v1,v4=v2,v2=v3,v1=v4 for the return path are feasible paths for qv1,v3. Second, suppose that there are two AF refueling stations available at sites m3 and m4. Then, the shortest original path of qv6,v11 does not offer AF refueling availability because vehicles cannot reach v11 with at least a half-full tank when they exit the network. For the shortest return path, vehicles that reenter the network with a half-full tank at v11 cannot reach the AF refueling station at site m3. Thus, vehicles in qv6,v11 would need to make detours using vertex sequence v6=v0,v7=v1,v8=v2,v9=v3,v10=v4,v7=v5,v8=v6,v11=v7 for the original path and opposite sequence v11=v0,v8=v1,v9=v2,v10=v3,v7=v4,v6=v5 for the return path. These two deviation paths are feasible paths for qv6,v11. Lastly, suppose that sites m5 and m6 are selected as AF refueling stations. Vehicles in qv12,v17 cannot be refueled by these two AF refueling stations if they use the shortest original and return paths. However, if vehicles make multiple cycles at vertex v13, then vertex sequence v12=v0,v13=v1,v14=v2,v13=v3,v15=v4,v13=v5,v17=v6 becomes a feasible original path and vertex sequence v17=v0,v15=v1,v13=v2,v12=v3 turn out to be a feasible return path for qv12,v17. Note that another vertex sequence v17=v0,v13=v1,v14=v2,v13=v3,v15=v4,v13=v5,v12=v6 also provides AF refueling availability for the return path, but drivers would not choose this vertex sequence because this path is longer than the path defined by the previous vertex sequence v17=v0,v15=v1,v13=v2,v12=v3. From these observations, we can state four feasible path properties in the next subsection.Figure 2
Three detour examples of a directed transportation network.


## 3.3. Four Properties of Feasible Paths

We derive four properties of feasible paths on a directed transportation network. Based on these four properties, an existingk-th shortest path algorithm will be modified in the next subsection to generate multiple feasible paths for all OD pairs. Since we consider round trips, the four properties are applicable to both feasible original and return paths. The first property represents that feasible paths have a proper sequence of candidate sites to be selected for AF refueling stations which cover trips in feasible paths.Property 1.
LetFkvi,vj be the k-th shortest feasible original path from vi to vj that contains feasible candidate sites for AF refueling stations to be able to cover trips in Pvi,vj, and CSFkvi,vj=m0,⋯,mf denote the sequence of candidate sites in Fkvi,vj, where m0,⋯,mf∈M. According to the cardinality of CSFkvi,vj, denoted as CSFkvi,vj, one of the following two cases must be satisfied:(a) IfCSFkvi,vj=1, then SSFkvi,vj∩ESFkvi,vj≠∅.
(b) IfCSFkvi,vj>1, then SSFkvi,vj≠∅, ESFkvi,vj≠∅, and dijkml,ml+1≤R for any two adjacent candidate sites ml,ml+1∈CSFkvi,vj.Proof.
SinceFkvi,vj is a feasible path, the half-full tank assumption implies that CSFkvi,vj≥1. In case (a), there is a single candidate site available for a refueling station on Fkvi,vj, and by the half-full tank assumption, this site is located within a distance of R/2 from vi and vj to ensure that vehicles on Fkvi,vj make a successful trip from vi to vj. Then, by definition of SSFkvi,vj and ESFkvi,vj, this site should be located in the common segment in SSFkvi,vj∩ESFkvi,vj. This implies that SSFkvi,vj∩ESFkvi,vj≠∅.
In case (b), there are at least two candidate sites available onFkvi,vj. If either SSFkvi,vj=∅ or ESFkvi,vj=∅, then Fkvi,vj cannot offer any AF refueling service. Thus, SSFkvi,vj≠∅ and ESFkvi,vj≠∅. For any two adjacent candidate sites ml,ml+1∈CSFkvi,vj, if dijkml,ml+1>R, then vehicles at site ml cannot reach site ml+1, i.e., Fkvi,vj does not have AF refueling availability. It also contradicts the definition of Fkvi,vj. Thus, dijkml,ml+1≤R for any two adjacent candidate sites ml,ml+1∈CSFkvi,vj.In practice, vehicles do not need to visit the same refueling station located at a single-access site several times. For example, let us consider two return path alternatives ofqv14,v16 in Figure 2. While a vertex sequence v16=v0,v15=v1,v13=v2,v17=v3,v15=v4,v13=v5,v14=v6 goes over site m6 twice, vehicles can also use a vertex sequence v16=v0,v15=v1,v13=v2,v14=v3 to travel from v16 to v14 with a single visit to site m6. That is, a shorter feasible path is always preferred to a longer feasible path that goes through a single-access candidate multiple times. On the other hand, it is unavoidable for vehicles to revisit the same refueling station when the station is placed at a dual-access candidate site. From one of the original paths of qv1,v3 in Figure 2, a feasible vertex sequence v1=v0,v2=v1,v4=v2,v2=v3,v3=v4 passes through a refueling station at dual-access site m1 twice, but in opposite directions. In this respect, prior to checking whether paths go over candidate sites multiple times, we treat a dual-access candidate site as two distinct single-access sites. From this observation, we derive the second property for candidate sites of feasible paths.Property 2.
Forqvi,vj∈Q, suppose that Fk1vi,vj is an original path that satisfies Property 1 and passes through a single-access candidate site mt1∈M twice:(4)
C
S
F
k
1
v
i
,
v
j
=
m
0
,
⋯
m
a
,
m
a
+
1
⋯
,
m
b
,
m
b
+
1
,
⋯
,
m
f
,
w
h
e
r
e
m
t
1
=
m
a
=
m
b
.Then, there must exist an original path, denoted asFk2vi,vj, that satisfies Property 1 and goes over site mt1 only one time, such that CSFk2vi,vj<CSFk1vi,vj:(5)CSFk2vi,vj=m0,⋯ma,mb+1,⋯,mf.Thus,dFk2vi,vj<dFk1vi,vj.Proof.
Letma and mb be the a-th and b-th candidate sites in CSFk1vi,vj=m0,⋯ma,ma+1⋯,mb,mb+1,⋯,mf which are located at the same site mt1 such that 0≤a<b≤f, where f=CSFk1vi,vj−1. By eliminating all consecutive candidate sites from ma+1 to mb and the corresponding arcs in Fk1vi,vj, we can construct path Fk2vi,vj with CSFk2vi,vj=m0,⋯ma,mb+1,⋯,mf. Then, from Property 1, since dijk1mb,mb+1≤R, vehicles at mb can reach mb+1 in Fk1vi,vj. This implies that vehicles at ma can reach mb+1 in Fk2vi,vj because ma and mb are placed at the same site mt1. Thus, CSFk2vi,vj is a feasible sequence of candidate sites for Fk2vi,vj and dFk2vi,vj<dFk1vi,vj.When drivers decide to make detours from their preplanned paths for refueling, they would be reluctant to travel along unnecessary paths to reach available refueling stations. For instance, let us take an original path ofqv1,v3 from the previous network in Figure 2. Recall that a vertex sequence v1=v0,v2=v1,v4=v2,v2=v3,v3=v4 is one of the feasible original paths for qv1,v3. Now, let us consider a different vertex sequence defined by v1=v0,v2=v1,v4=v2,v5=v3,v4=v4,v2=v5,v3=v6. The candidate site sequence for this second path is m1=m0,m2=m1,m1=m2, which enables successful trips from v1 to v3 without running out of fuel. However, drivers would not travel further along the subpath of the second path, v4=v2,v5=v3,v4=v4 because they can reach v3 without refueling at site m2. They would return at v4 after visiting the refueling station at site m1, i.e., they would use the first feasible path, v1=v0,v2=v1,v4=v2,v2=v3,v3=v4, instead of the second one for the original path of qv1,v3. Similarly, drivers can consider two possible paths for the return path of qv12,v17. While vehicles using a vertex sequence v17=v0,v13=v1,v14=v2,v13=v3,v15=v4,v13=v5,v12=v6 require each refueling stop at both sites m5 and m6, vehicles in a vertex sequence v17=v0,v15=v1,v13=v2,v12=v3 need only one refueling stop at site m6. Since the length of the latter path is shorter, vehicles in the return path of qv12,v17 would use the second vertex sequence instead of the first one. Based on these observations, we conclude that drivers would not make detours along a vertex sequence which requires unnecessary travel to reach available refueling stations. From now on, we call a site in unnecessary subpaths an irrelevant candidate site, so that feasible paths can be represented by vertex sequences without irrelevant candidate sites. Property 3 provides conditions under which there exist irrelevant candidate sites in paths.Property 3.
Given two paths,Fk1vi,vj and Fk2vi,vj, such that dFk1vi,vj<dFk2vi,vj, which satisfy Properties 1 and 2, Fk2vi,vj includes irrelevant candidate sites in set CSFk2vi,vj\CSFk1vi,vj when one of the following conditions is satisfied:(a)
IfCSFk1vi,vj=1, then (1) CSFk1vi,vj⊂CSFk2vi,vj, and (2) either SFk1vi,vj=SSFk2vi,vj or ESFk1vi,vj=ESFk2vi,vj.
(b)
IfCSFk1vi,vj>1, then (1) CSFk1vi,vj⊂CSFk2vi,vj, (2) SSFk1vi,vj=SSFk2vi,vj, (3) ESFk1vi,vj=ESFk2vi,vj, and (4) 0<dijk1m0,mf≤R, where CSFk1vi,vj=m0,⋯,mf.Proof.
In case (a), suppose thatm0∈SSFk1vi,vj=SSFk2vi,vj. By Property 1, m0∈ESFk1vi,vj, which implies that dijk1m0,vj≤R/2. Thus, any candidate site in CSFk2vi,vj\CSFk1vi,vj is irrelevant because vehicles covered by a single refueling station at site m0 can complete trips from vi to vj. Similarly, it can be shown when m0∈ESFk1vi,vj=ESFk2vi,vj.
In case (b), by Condition (4),dijk1ma,mb≤R such that ma∈SSFk1vi,vj, mb∈ESFk1vi,vj, and 0≤a<b≤f. This implies that any additional refueling stop is unnecessary between ma and mb. Then, by Conditions (1), (2), and (3), ISFk2vi,vj≠∅, which implies that Fk2vi,vj includes irrelevant candidate sites.Given two feasible paths,Fk1vi,vj and Fk2vi,vj such that dFk1vi,vj≤dFk2vi,vj, for the same OD pair qvi,vj∈Q, if any subset of candidate sites that covers the trips in Fk2vi,vj is also able to cover the trips in Fk1vi,vj, then we say that Fk2vi,vj is dominated by Fk1vi,vj. This implies that Fk2vi,vj is unnecessary to be considered in our proposed model. In practice, drivers would use only Fk1vi,vj because dFk1vi,vj≤dFk2vi,vj. For example, let us take the previous example of qv6,v11 in Figure 2. Recall that F1v6,v11=v6=v0,v7=v1,v8=v2,v9=v3,v10=v4,v7=v5,v8=v6,v11=v7. Then, we can find another feasible original path, which is F2v6,v11=v6=v0,v7=v1,v10=v2,v9=v3,v10=v4,v7=v5,v8=v6,v11=v7 for vehicles with R=100. They have the same sequence of candidate sites, but vehicles may prefer F1v6,v11 to F2v6,v11 because dF1v6,v11=125<dF2v6,v11=145. Furthermore, in order to clarify the concept of dominated feasible paths, let us see another example of qv14,v16 in Figure 2. For the return path, we have two feasible return paths, F1v16,v14=v16=v0,v15=v1,v13=v2,v14=v3 and F2v16,v14=v16=v0,v15=v1,v17=v2,v13=v3,v14=v4, such that CSF1v16,v14=m7=m0,m6=m1,m5=m2 and CSF2v16,v14=m7=m0,m5=m1. For F1v16,v14, we have five possible sets of candidate sites to cover the trip, i.e., m6, m6,m5, m7,m6, m7,m5, and m7,m6,m5. In case of F2v16,v14, m7,m5 is the only feasible set. That is, CSF2v16,v14⊆CSF1v16,v14. It means that F2v16,v14 is dominated by F1v16,v14 because dF1v16,v14=65<dF2v16,v14=100. In this respect, Property 4 shows that the shortest feasible paths can dominate other feasible paths when the candidate sites of the shortest feasible paths include the candidate sites of the other paths.Property 4.
Givenqvi,vj∈Q, let Fk1vi,vj be the shortest feasible original path with CSFk1vi,vj. Suppose that there exists another feasible original path, denoted as Fk2vi,vj, such that CSFk2vi,vj⊆CSFk1vi,vj and dFk1vi,vj≤dFk2vi,vj. Then, Fk2vi,vj is dominated by Fk1vi,vj.Proof.
By definition of dominated paths, we need to show that any subset of candidate sites to cover the trips inFk2vi,vj is also able to cover the trips in Fk1vi,vj. Suppose that FSFk2vi,vj is a set of candidate sites that covers trips in Fk2vi,vj. This implies that FSFk2vi,vj⊆CSFk2vi,vj, and therefore, FSFk2vi,vj⊆CSFk1vi,vj because CSFk2vi,vj⊆CSFk1vi,vj. If FSFk2vi,vj=1, then it is trivial to show that FSFk2vi,vj can cover trips in Fk1vi,vj because dFk1vi,vj≤dFk2vi,vj. If FSFk2vi,vj>1, then we select the first and last sequences of candidate sites in FSFk2vi,vj, denoted as ma and mb, respectively. When we decompose Fk1vi,vj into vi=v0,⋯,ma, ma,⋯,mb, and mb,⋯,vj=vw, we have that dFk1vi,vj=dijk1v0,ma+dijk1ma,mb+dijk1mb,vw. If dijk2v0,ma≤R/2<dijk1v0,ma, then we can construct a shorter feasible path whose length is dijk2v0,ma+dijk1ma,mb+dijk1mb,vw. It contradicts that Fk1vi,vj is the shortest feasible path. Thus, dijk1v0,ma≤R/2. Similarly, it can be shown that dijk1mb,vw≤R/2. Furthermore, if dijk2ml,ml+1≤R<dijk1ml,ml+1 such that a≤l<b, then it also contradicts that Fk1vi,vj is the shortest feasible path. Thus, FSFk2vi,vj can cover trips in Fk1vi,vj by Property 1.


## 3.4.k-th Shortest Nondominated Feasible Path Algorithm

In this subsection, we describe thek-th shortest nondominated feasible (SNDF) path algorithm to generate multiple paths that (1) offer AF refueling availability, (2) contain unique single-access candidate sites, (3) do not have irrelevant candidate sites, and (4) are not dominated by other paths for all OD pairs. We extend the k-th shortest path algorithm developed by Eppstein [12] to first construct candidate paths between ODs for SNDF paths, and then eliminate some of the paths that violate the four properties discussed on the previous subsection. Finally, given a deviation factor β and a maximum number K of SNDF paths, the k-th SNDF path algorithm determines up to K SNDF paths from vi to vj and another up to K SNDF paths from vj to vi.The first step of thek-th SNDF path algorithm is to transform a simple directed network GV,E with the set of nV vertices and the set of nE arcs into an expanded directed network, denoted as G′V′,E′, by transforming each dual-access candidate site in G into two distinct single-access candidate sites, one in each side of the road, and setting all candidate sites in G as vertices in G′. Then, the maximum cardinalities of V′ and E′ can be written as V′=nV+nSM+2nDM and E′=nE+2nSM+4nDM, respectively, where nSM and nDM refer to the number of single- and dual-access candidate sites. Next, for interchange vj∈V′, Tvj denotes a shortest path arborescence rooted at vertex vj containing the shortest path to each vi in V′\vj. The shortest path arborescence Tvj for interchange vj in G′V′,E′ can be constructed using Dijkstra’s algorithm in time OE′+V′logV′. Let ETvj denote the set of arcs composing Tvj. If we select arc va,vb∈E′\ETvj and follow the shortest path from vb to vj instead of taking the shortest path in Tvj, then the additional distance from Tvj, denoted as δvjva,vb, is computed as: δvjva,vb=dva,vb+dPSvb,vj−dPSva,vj≥0. Note that δvjva,vb>0 for va,vb∈E′\ETvj and δvjva,vb=0 for va,vb in Tvj. This implies that we can calculate the length of path Pvi,vj by adding the values of δvjva,vb for the arcs in Pvi,vj to the length of PSvi,vj since δvjva,vb=0 for va,vb in PSvi,vj. These arcs aside from PSvi,vj do not belong to Tvj, and Eppstein’s algorithm calls these arcs the sidetracks of PSvi,vj. The sidetracks are used to generate the k-shortest path of Pvi,vj. For Pvi,vj, a special tree-based data structure, called the min-heap data structure, in which it is a completely binary tree with data structure in a nondecreasing order of values of sidetracks, is applied to build K shortest paths of Pvi,vj using the depth first search algorithm in time OKlogK. The k-th shortest path of Pvi,vj, for k=1,⋯,K, generated through the min-heap data structure, is called a candidate path from vi to vj and is denoted as F′vi,vj.When generatingF′vi,vj, we determine whether or not F′vi,vj is included in the set of SNDF paths for qvi,vj by Properties 1, 2, 3, and 4. First, F′vi,vj needs to have AF refueling availability by Property 1. If either SSF′vi,vj∩ESF′vi,vj≠∅ when CSF′vi,vj=1 or SSF′vi,vj≠∅, ESF′vi,vj≠∅, and dijkml,ml+1≤R for ml,ml+1∈CSF′vi,vj when CSF′vi,vj>1, then F′vi,vj has AF refueling availability. Second, if F′vi,vj repeats some single-access candidate sites, then the path is excluded from the set of SNDF paths because we also construct a similar feasible path without repeating these candidate sites whose length is shorter than dF′vi,vj by Property 2. Next, regarding Properties 3 and 4, we compare CSF′vi,vj with CSFkvi,vj’s when k>1. If CSFkvi,vj⊂CSF′vi,vj for k=1,⋯,K−1, we diagnose whether or not F′vi,vj has irrelevant candidate sites by Property 3. If CSF′vi,vj⊆CSFkvi,vj for k=1,⋯,K−1, then F′vi,vj is removed from the set of SNDF paths for qvi,vj because the associated OD pair already has shortest feasible paths which dominate F′vi,vj by Property 4. Until reaching the maximum number, K, of SNDF paths, Eppstein’s algorithm updates a binary heap that includes sets of sidetracked arcs to find the next paths. However, the k-th SNDF path algorithm can even stop before reaching K if the sum of sidetracked arcs’ δvjva,vb on the root of the binary heap becomes greater than the value of βdPSvi,vj. This represents that the length of a newly constructed path cannot exceed a maximum deviation distance of drivers. After determining at most K SNDF paths from interchanges vi’s to a certain interchange vj, we similarly apply the above procedures to other interchanges vj’s. Regarding 2K SNDF paths in qvi,vj, Algorithm 1 describes the pseudocode of the k-th SNDF path algorithm for all OD pairs.Algorithm 1: Overview of the k-th SNDF path algorithm for constructing at most 2K SNDF paths.
Algorithm 1:k-th SNDF Path Algorithm
Input:GV,E,M,R,β, and K
Output:G′V′,E′, sets of all K SNDF original paths, and K SNDF return paths for all OD pairs
Construct an expanded directed networkG′ by transforming dual-candidate sites into two separate single-candidate sites on each direction and adding all sites to G as vertices.
Setj=1.
Repeat
Build a shortest path arborescenceTvj rooted at vertex vj.
Seti=1,i≠j.
Repeat
Setk=0.
Construct a min binary heap data structure, denoted asHvi,vj, with sequences of sidetracked arcs using the corresponding values of δvjva,vb, ∀va,vb∈EG′, through Eppstein’s algorithm. Note that the binary heap initially has the empty sequence of sidetracked arcs because the empty sequence of sidetrack arcs represents PSvi,vj.
Repeat
Construct a candidate pathF′vi,vj by inserting a sequence of sidetracked arcs on the root of Hvi,vj into Tvj.
ConstructCSF′vi,vj,SSF′vi,vj,ESF′vi,vj,andISF′vi,vj.
IfF′vi,vj is satisfied with Properties 1, 2, 3, and 4, then define F′vi,vj, SSF′vi,vj, ESF′vi,vj and ISF′vi,vj as Fkvi,vj, SSFkvi,vj, ESFkvi,vj, and ISFkvi,vj, respectively, and set k=k+1.
Remove the current sequence of sidetracked arcs on the root ofHvi,vj, and insert new sequences of sidetracked arcs from the min-heap data structure to the binary heap Hvi,vj by maintaining the min-heap property.
Untilk>K, ∑va,vb∈AEPSvi,vjδvjva,vb>β∙dPSvi,vj, or the min-heap data structure does not have any sidetracked arcs.
Seti=i+1.
Untili reaches nV+1.
Setj=j+1.
Until
j reaches nV+1.We prove below that at most2K SNDF paths (i.e., up to K original paths and another up to K SNDF return paths) for each OD pair can be determined by the k-th SNDF path algorithm.Theorem 1.
Given a value for the maximum number of paths,K, the k-th SNDF path algorithm generates at most 2K SNDF paths for each OD pair in G.Proof.
After transformingG into G′, given qvi,vj∈Q, F′vi,vj is generated by Eppstein’s algorithm. The k-th SNDF path algorithm determines CSF′vi,vj, and then classifies the sequence of candidate sites into SSF′vi,vj, ESF′vi,vj, and ISF′vi,vj by measuring distances between the sites, vi and vj. The algorithm checks first whether or not F′vi,vj meets Properties 1 and 2. Next, if the corresponding OD pair already has SNDF paths, then CSF′vi,vj is compared with CSFkvi,vj for k=1,⋯,K−1 by Properties 3 and 4. If CSFkvi,vj⊂CSF′vi,vj for k=1,⋯,K−1, then the algorithm checks whether or not F′vi,vj has irrelevant candidate sites. If CSF′vi,vj⊆CSFkvi,vj for k=1,⋯,K−1, then F′vi,vj is not included in the set of SNDF paths for qvi,vj. Otherwise, the algorithm assigns F′vi,vj to the k-th SNDF path Fkvi,vj and increases the value of k by one. If k≤K, a next temporary path can be constructed by adding new sequences of sidetracked arcs into Tvj. The additional distance of the next path can be calculated as ∑va,vb∈AEPSvi,vjδvjva,vb, where AEPSvi,vj is the set of newly added arcs and δvjva,vb=dva,vb+dPSvb,vj−dPSva,vj. If the additional travel distance does not exceed β·dPSvi,vj, then the above processes are repeated until k>K. Otherwise, the algorithm moves to the next interchange vi in the shortest path arborescence Tvj. After determining at most K SNDF paths, we repeat the algorithm by building a new shortest path arborescence for other interchanges. Therefore, for each OD pair, this algorithm determines up to K SNDF paths from origin to destination and up to K other SNDF paths from destination to origin, which leads to at most 2K SNDF paths.Since the original networkG has single- and dual-access candidate sites for AF refueling stations, vehicles can select different return paths from their SNDF original paths when they come back to their origin interchanges. This means that nonsymmetric round trips can be considered for each OD pair in our proposed model. In practice, drivers would make different detours from their SNDF original paths if the original paths do not provide AF refueling service in the opposite direction. For convenience, L^o=Fkvi,vj|k=1,⋯,K,∀qvi,vj∈Q and L^r=Fkvj,vi|k=1,⋯,K,∀qvi,vj∈Q denote the sets of all K SNDF original paths and K SNDF return paths for all OD pairs, respectively.Theorem 2.
In order to generate2K SNDF original and return paths for all OD pairs, the complexity of the k-th SNDF path algorithm takes time OVE′+V′logV′+VKM2+E′KlogK, where V′=nV+nSM+2nDM and E′=nE+2nSM+4nDM. Recall that V=nV, E=nE, and M=nM. Also, nSM and nDM refer to the number of single- and dual-access candidate sites, respectively.Proof.
TimeOnDM is required to transform GV,E with the set of all candidate sites to a simple expanded directed network G′V′,E′ by regarding dual-access candidate sites as two distinct single-access candidate sites on arcs. Given G′, the basic version of Eppstein’s algorithm takes time OV′+V′logV′+VKlogK to find the k-th shortest path to vj from any other interchanges. In time OE′+V′logV′, the algorithm constructs a shortest path arborescence and a min-heap data structure to output F′vi,vj. Then, in time OKlogK, it generates F′vi,vj’s from a specific interchange to vj. After we construct CSF′vi,vj and classify CSF′vi,vj into SSF′vi,vj, ESF′vi,vj, and ISF′vi,vj, it takes time OM2+E′ for Property 1. Next, we find out which candidate sites are repeated in F′vi,vj in time OM2 for Property 2. Also, when k>1, time OKM2 is required to compare CSF′vi,vj with CSFkvi,vj for Properties 3 and 1. These processes are repeated until either K feasible paths are generated or the length of a newly constructed path does not exceed 1+βdPSvi,vj in time OKM2+E′KlogK. Next, we move to other interchanges vi’s, ∀vi∈VG, i≠j, which requires time OVKM2+E′KlogK. In summary, generating K SNDF paths from other interchanges to a certain interchange vj requires time OE′+V′logV′+VKM2+E′KlogK. Finally, since the number of interchanges for ODs is V, the total computational complexity of the algorithm takes time OVE′+V′logV′+VKM2+E′KlogK in order to generate 2K SNDF original and return paths in all OD pairs.


## 4. Model Formulation

Due to the limited driving range and the amount of fuel remaining at ODs, AF refueling stations should be properly placed at candidate sites along paths to ensure vehicles safely travel from one point to another without running out of fuel. In other words, in order to cover trips onFkvi,vj, when dFkvi,vj>R, at least one AF refueling station should be located at each candidate site in SSFkvi,vj and ESFkvi,vj, and then locating additional stations at candidate sites in ISFkvi,vj are required if the distance between the two candidate sites in SSFkvi,vj and ESFkvi,vj is greater than R. On the other hand, when dFkvi,vj≤R, since the distance between the two candidate sites in SSFkvi,vj and ESFkvi,vj is less than or equal to R, additional AF refueling stations are unnecessary to be considered for AF refueling stations between the sites in SSFkvi,vj and ESFkvi,vj. Similar covering conditions are applied to the return path. From this observation, the following types of k-th SNDF path are defined according to dFkvi,vj and dFkvj,vi for Fkvi,vj and Fkvj,vi, respectively, in OD pairs qvi,vj∈Q:(i)
Type 1k-th SNDF original path: Fk1vi,vj=Fkvi,vj∈L^o|0<dFkvi,vj≤R,(ii)
Type 1k-th SNDF return path: Fk1vj,vi=Fkvj,vi∈L^r|0<dFkvj,vi≤R,(iii)
Type 2k-th SNDF original path: Fk2vi,vj=Fkvi,vj∈L^o|R<dFkvi,vj,(iv)
Type 2k-th SNDF return path: Fk2vj,vi=Fkvj,vi∈L^r|R<dFkvj,vi.


### 4.1. Coverage for Type 1k-th SNDF Path

In case of Type 1k-th SNDF original paths, since the distance between two candidate sites in SSFk1vi,vj and ESFk1vi,vj is at most R, vehicles can travel from vi to vj by visiting two refueling stations at the sites. Sometimes, if there exists a single candidate site in the common area of SSFk1vi,vj and ESFk1vi,vj, then locating a refueling station at this common area is sufficient to cover trips on Type 1 k-th SNDF original paths. Thus, for proper AF refueling availability on Type 1 k-th SNDF original paths, AF refueling stations should be located at either one common or two distinct candidate sites in SSFk1vi,vj and ESFk1vi,vj. Figure 3 depicts trips in Fk1vi,vj covered by AF refueling stations at either m1∈SSFk1vi,vj∩ESFk1vi,vj or m0∈SSFk1vi,vj and m2∈ESFk1vi,vj. Similar covering conditions are also applied to Type 1 k-th SNDF return paths.Figure 3
An example of Type 1k-th SNDF original path.


### 4.2. Coverage for Type 2k-th SNDF Path

When two candidate sites in eachSSFk2vi,vj and ESFk2vi,vj are available to locate AF refueling stations on Type 2 SNDF original paths, the distance between the two sites can be longer than R because the length of Type 2 k-th SNDF original paths is greater than R. That is, if vehicles at a station in SSFk2vi,vj cannot reach the other stations in ESFk2vi,vj due to the limited driving range, then additional refueling stations at candidate sites in ISFk2vi,vj should be located to ensure that the distance between consecutive stations is less than or equal to R. In order to consider this covering condition for Type 2 k-th SNDF original paths, we define identification coefficients that indicate whether a vehicle that is refueled at a station located in a candidate site can reach another station located at a subsequent candidate site to receive refueling services. For candidate sites in USFk2vi,vj=ml′∈ISFk2vi,vj∪ESFk2vi,vj, we determine identification coefficients and use them as constraint coefficients for covering conditions of Type 2 k-th SNDF original paths in our proposed model.The values of identification coefficientsαk,ml,ml′, where ml∈CSFk2vi,vj and ml′∈USFk2vi,vj, are determined as follows:(6)αk,ml,ml′=1,ifdijkv0,ml<dijkv0,ml′and0<dijkml,ml′≤R,0,otherwise.Identification coefficientαk,ml,ml′=1 indicates that site ml′ follows ml along the k-th path and their distance is less than or equal to R. In an example in Figure 4, we have CSFk2vi,vj=m0,m1,m2,m3,m4 and USFk2vi,vj=m2,m3,m4. If dijkm0,m3>R, dijkm1,m3≤R, and dijkm2,m3≤R, then the corresponding identification coefficients for site m3 are set as αk,m0,m3=0, αk,m1,m3=1, αk,m2,m3=1, αk,m3,m3=0, and αk,m4,m3=0. These values mean that a refueling station must be located at either site m1 or m2, so as to reach site m3 in Fk2vi,vj. Note that Hwang et al. [6] apply a similar procedure that uses identification coefficients to setup coverage restrictions to a unique path for each OD pair. For Type 2 k-th SNDF return paths, USFk2vj,vi is similarly defined to determine identification coefficients.Figure 4
An example of Type 2k-th SNDF original path.


### 4.3. Mixed-Integer Programming Model

Given a predetermined maximum number of feasible paths and predetermined locations of candidate sites for AF refueling stations on a directed transportation network, a mixed-integer programming model is formulated to select thep optimal station locations that maximize the traffic flow covered (in round trips per time unit) along multiple feasible original and return paths between ODs. We first introduce two types of parameters and three types of decision variables as follows:


#### 4.3.1. Parameters

f
v
i
,
v
j: average traffic flow in round trips per time unit between interchanges vi and vj,α
k
,
m
t
,
m
t
′: identification coefficient for sites mt and mt′ for Type 2 k-th SNDF path.


#### 4.3.2. Decision Variables.

(7)
x
m
t
=
1
,
i
f
r
e
f
u
e
l
i
n
g
s
t
a
t
i
o
n
l
o
c
a
t
i
o
n
m
t
∈
M
i
s
s
e
l
e
c
t
e
d
,
0
,
o
t
h
e
r
w
i
s
e
,(8)
y
k
,
v
i
,
v
j
=
1
,
i
f
t
h
e
k
−
t
h
S
N
D
F
p
a
t
h
f
r
o
m
v
i
t
o
v
j
i
s
c
o
v
e
r
e
d
b
y
t
h
e
c
u
r
r
e
n
t
s
e
t
o
f
r
e
f
u
e
l
i
n
g
s
t
a
t
i
o
n
s
,
0
,
o
t
h
e
r
w
i
s
e
,(9)
z
v
i
,
v
j
=
1
,
i
f
t
h
e
r
o
u
n
d
t
r
i
p
o
f
t
r
a
f
f
i
c
f
l
o
w
f
v
i
,
v
j
i
s
c
o
v
e
r
e
d
b
y
t
h
e
c
u
r
r
e
n
t
s
e
t
o
f
r
e
f
u
e
l
i
n
g
s
t
a
t
i
o
n
s
,
0
,
o
t
h
e
r
w
i
s
e
.Then, we propose the following a mixed-integer programming model:(10)maximize∑qvi,vj∈Qfvi,vjzvi,vj,subject to(11)∑mt∈SSFkvi,vjxmt≥yk,vi,vj,∀Fkvi,vj∈L^o,k=1,⋯,K,(12)∑mt∈ESFkvi,vjxmt≥yk,vi,vj,∀Fkvi,vj∈L^o,k=1,⋯,K,(13)∑mt∈SSFkvj,vixmt≥yk,vj,vi,∀Fkvj,vi∈L^r,k=1,⋯,K,(14)∑mt∈ESFkvj,vixmt≥yk,vj,vi,∀Fkvj,vi∈L^r,k=1,⋯,K,(15)∑mt∈CSFk2vi,vjαk,mt,mt′xmt≥yk,vi,vj,∀mt′∈USFk2vi,vj,∀Fk2vi,vj∈L^o,k=1,⋯,K,(16)∑mt∈CSFk2vj,viαk,mt,mt′xmt≥yk,vj,vi,∀mt′∈USFk2vj,vi,∀Fk2vj,vi∈L^r,k=1,⋯,K,(17)∑k=1,⋯,Kyk,vi,vj≥zvi,vj,∀qvi,vj∈Q,(18)∑k=1,⋯,Kyk,vj,vi≥zvi,vj,∀qvi,vj∈Q,(19)∑mt∈Mxmt=p,(20)xmt∈0,1,∀mt∈M,(21)yk,vi,vj∈0,1,∀Fkvi,vj∈L^o,k=1⋯K,(22)yk,vj,vi∈0,1,∀Fkvj,vi∈L^r,k=1⋯K,(23)zvi,vj∈0,1,∀qvi,vj∈Q.In this formulation, the objective function (10) maximizes the traffic flow that is covered by the set of AF refueling stations located at thep selected candidate sites on the transportation network. Constraint sets (11) and (12) guarantee that trips along K SNDF original paths of Types 1 and 2 receive refueling service from one station in set of candidate sites SSFkvi,vj and another one in set of candidate sites ESFkvi,vj. The k-th SNDF return paths have similar refueling conditions in Constraint sets (13) and (14). Next, coefficients in Constraint set (15) are determined by identification coefficients αk,mt,mt′ for mt∈CSFk2vi,vj and mt′∈USFk2vi,vj, where αk,mt,mt′=1 if dijkvi,mt<dijkvi,mt′ and 0<dijkmt,mt′≤R; otherwise, αk,mt,mt′=0. This constraint set guarantees that, if a refueling station is located at mt′ in USFk2vi,vj, then one of sites in CSFk2vi,vj should be selected to locate another station which would provide refueling for vehicles to reach the station at mt′. Similarly, Constraint set (16) is applied to candidate sites on Type 2 SNDF return paths. In addition, since round trips in OD pairs are considered in this problem, all possible combinations of original and return paths are determined by Constraint sets (17) and (18), respectively. These constraint sets ensure vehicles on a directed transportation network to make nonsymmetric round trips if necessary for refueling service. Constraint (19) allows exactly p AF refueling stations to be located at the predetermined candidate site locations. Lastly, the three kinds of binary decision variables are defined in Constraint sets (20) to (23).The main stakeholders of the proposed model that maximizes the covered traffic flows can be shareholders, investors, employers and employees, and any other private entities whose stake is directly or indirectly tied to this objective function, assuming that the revenue is proportional to the traffic flow covered by a given number of stations; that is, the more traffic flow covered by the stations, the higher revenue the stakeholders expect. If we assume that the construction cost does not significantly differ by station type or region, then Constraint (19) can play a role as the budget constraint. This assumption, however, would not be always true in practice [60]. So, we can relax this assumption by replacing Constraint (19) by the following constraint:(24)∑mt∈Mcmtxmt≤B,wherecmt refers to the capital cost for building a refueling station in candidate site mt∈M, and B the refueling infrastructure budget. Constraint (24) allows different capital costs by station type (e.g., single-access vs. dual-access) or by region and forces the model to determine candidate sites for refueling stations within the budget. This constraint is used to build AF refueling stations in a statewide network to showcase that the proposed model is well-suited for solving practical problems in Section 6.


## 4.1. Coverage for Type 1k-th SNDF Path

In case of Type 1k-th SNDF original paths, since the distance between two candidate sites in SSFk1vi,vj and ESFk1vi,vj is at most R, vehicles can travel from vi to vj by visiting two refueling stations at the sites. Sometimes, if there exists a single candidate site in the common area of SSFk1vi,vj and ESFk1vi,vj, then locating a refueling station at this common area is sufficient to cover trips on Type 1 k-th SNDF original paths. Thus, for proper AF refueling availability on Type 1 k-th SNDF original paths, AF refueling stations should be located at either one common or two distinct candidate sites in SSFk1vi,vj and ESFk1vi,vj. Figure 3 depicts trips in Fk1vi,vj covered by AF refueling stations at either m1∈SSFk1vi,vj∩ESFk1vi,vj or m0∈SSFk1vi,vj and m2∈ESFk1vi,vj. Similar covering conditions are also applied to Type 1 k-th SNDF return paths.Figure 3
An example of Type 1k-th SNDF original path.


## 4.2. Coverage for Type 2k-th SNDF Path

When two candidate sites in eachSSFk2vi,vj and ESFk2vi,vj are available to locate AF refueling stations on Type 2 SNDF original paths, the distance between the two sites can be longer than R because the length of Type 2 k-th SNDF original paths is greater than R. That is, if vehicles at a station in SSFk2vi,vj cannot reach the other stations in ESFk2vi,vj due to the limited driving range, then additional refueling stations at candidate sites in ISFk2vi,vj should be located to ensure that the distance between consecutive stations is less than or equal to R. In order to consider this covering condition for Type 2 k-th SNDF original paths, we define identification coefficients that indicate whether a vehicle that is refueled at a station located in a candidate site can reach another station located at a subsequent candidate site to receive refueling services. For candidate sites in USFk2vi,vj=ml′∈ISFk2vi,vj∪ESFk2vi,vj, we determine identification coefficients and use them as constraint coefficients for covering conditions of Type 2 k-th SNDF original paths in our proposed model.The values of identification coefficientsαk,ml,ml′, where ml∈CSFk2vi,vj and ml′∈USFk2vi,vj, are determined as follows:(6)αk,ml,ml′=1,ifdijkv0,ml<dijkv0,ml′and0<dijkml,ml′≤R,0,otherwise.Identification coefficientαk,ml,ml′=1 indicates that site ml′ follows ml along the k-th path and their distance is less than or equal to R. In an example in Figure 4, we have CSFk2vi,vj=m0,m1,m2,m3,m4 and USFk2vi,vj=m2,m3,m4. If dijkm0,m3>R, dijkm1,m3≤R, and dijkm2,m3≤R, then the corresponding identification coefficients for site m3 are set as αk,m0,m3=0, αk,m1,m3=1, αk,m2,m3=1, αk,m3,m3=0, and αk,m4,m3=0. These values mean that a refueling station must be located at either site m1 or m2, so as to reach site m3 in Fk2vi,vj. Note that Hwang et al. [6] apply a similar procedure that uses identification coefficients to setup coverage restrictions to a unique path for each OD pair. For Type 2 k-th SNDF return paths, USFk2vj,vi is similarly defined to determine identification coefficients.Figure 4
An example of Type 2k-th SNDF original path.


## 4.3. Mixed-Integer Programming Model

Given a predetermined maximum number of feasible paths and predetermined locations of candidate sites for AF refueling stations on a directed transportation network, a mixed-integer programming model is formulated to select thep optimal station locations that maximize the traffic flow covered (in round trips per time unit) along multiple feasible original and return paths between ODs. We first introduce two types of parameters and three types of decision variables as follows:


### 4.3.1. Parameters

f
v
i
,
v
j: average traffic flow in round trips per time unit between interchanges vi and vj,α
k
,
m
t
,
m
t
′: identification coefficient for sites mt and mt′ for Type 2 k-th SNDF path.


### 4.3.2. Decision Variables.

(7)
x
m
t
=
1
,
i
f
r
e
f
u
e
l
i
n
g
s
t
a
t
i
o
n
l
o
c
a
t
i
o
n
m
t
∈
M
i
s
s
e
l
e
c
t
e
d
,
0
,
o
t
h
e
r
w
i
s
e
,(8)
y
k
,
v
i
,
v
j
=
1
,
i
f
t
h
e
k
−
t
h
S
N
D
F
p
a
t
h
f
r
o
m
v
i
t
o
v
j
i
s
c
o
v
e
r
e
d
b
y
t
h
e
c
u
r
r
e
n
t
s
e
t
o
f
r
e
f
u
e
l
i
n
g
s
t
a
t
i
o
n
s
,
0
,
o
t
h
e
r
w
i
s
e
,(9)
z
v
i
,
v
j
=
1
,
i
f
t
h
e
r
o
u
n
d
t
r
i
p
o
f
t
r
a
f
f
i
c
f
l
o
w
f
v
i
,
v
j
i
s
c
o
v
e
r
e
d
b
y
t
h
e
c
u
r
r
e
n
t
s
e
t
o
f
r
e
f
u
e
l
i
n
g
s
t
a
t
i
o
n
s
,
0
,
o
t
h
e
r
w
i
s
e
.Then, we propose the following a mixed-integer programming model:(10)maximize∑qvi,vj∈Qfvi,vjzvi,vj,subject to(11)∑mt∈SSFkvi,vjxmt≥yk,vi,vj,∀Fkvi,vj∈L^o,k=1,⋯,K,(12)∑mt∈ESFkvi,vjxmt≥yk,vi,vj,∀Fkvi,vj∈L^o,k=1,⋯,K,(13)∑mt∈SSFkvj,vixmt≥yk,vj,vi,∀Fkvj,vi∈L^r,k=1,⋯,K,(14)∑mt∈ESFkvj,vixmt≥yk,vj,vi,∀Fkvj,vi∈L^r,k=1,⋯,K,(15)∑mt∈CSFk2vi,vjαk,mt,mt′xmt≥yk,vi,vj,∀mt′∈USFk2vi,vj,∀Fk2vi,vj∈L^o,k=1,⋯,K,(16)∑mt∈CSFk2vj,viαk,mt,mt′xmt≥yk,vj,vi,∀mt′∈USFk2vj,vi,∀Fk2vj,vi∈L^r,k=1,⋯,K,(17)∑k=1,⋯,Kyk,vi,vj≥zvi,vj,∀qvi,vj∈Q,(18)∑k=1,⋯,Kyk,vj,vi≥zvi,vj,∀qvi,vj∈Q,(19)∑mt∈Mxmt=p,(20)xmt∈0,1,∀mt∈M,(21)yk,vi,vj∈0,1,∀Fkvi,vj∈L^o,k=1⋯K,(22)yk,vj,vi∈0,1,∀Fkvj,vi∈L^r,k=1⋯K,(23)zvi,vj∈0,1,∀qvi,vj∈Q.In this formulation, the objective function (10) maximizes the traffic flow that is covered by the set of AF refueling stations located at thep selected candidate sites on the transportation network. Constraint sets (11) and (12) guarantee that trips along K SNDF original paths of Types 1 and 2 receive refueling service from one station in set of candidate sites SSFkvi,vj and another one in set of candidate sites ESFkvi,vj. The k-th SNDF return paths have similar refueling conditions in Constraint sets (13) and (14). Next, coefficients in Constraint set (15) are determined by identification coefficients αk,mt,mt′ for mt∈CSFk2vi,vj and mt′∈USFk2vi,vj, where αk,mt,mt′=1 if dijkvi,mt<dijkvi,mt′ and 0<dijkmt,mt′≤R; otherwise, αk,mt,mt′=0. This constraint set guarantees that, if a refueling station is located at mt′ in USFk2vi,vj, then one of sites in CSFk2vi,vj should be selected to locate another station which would provide refueling for vehicles to reach the station at mt′. Similarly, Constraint set (16) is applied to candidate sites on Type 2 SNDF return paths. In addition, since round trips in OD pairs are considered in this problem, all possible combinations of original and return paths are determined by Constraint sets (17) and (18), respectively. These constraint sets ensure vehicles on a directed transportation network to make nonsymmetric round trips if necessary for refueling service. Constraint (19) allows exactly p AF refueling stations to be located at the predetermined candidate site locations. Lastly, the three kinds of binary decision variables are defined in Constraint sets (20) to (23).The main stakeholders of the proposed model that maximizes the covered traffic flows can be shareholders, investors, employers and employees, and any other private entities whose stake is directly or indirectly tied to this objective function, assuming that the revenue is proportional to the traffic flow covered by a given number of stations; that is, the more traffic flow covered by the stations, the higher revenue the stakeholders expect. If we assume that the construction cost does not significantly differ by station type or region, then Constraint (19) can play a role as the budget constraint. This assumption, however, would not be always true in practice [60]. So, we can relax this assumption by replacing Constraint (19) by the following constraint:(24)∑mt∈Mcmtxmt≤B,wherecmt refers to the capital cost for building a refueling station in candidate site mt∈M, and B the refueling infrastructure budget. Constraint (24) allows different capital costs by station type (e.g., single-access vs. dual-access) or by region and forces the model to determine candidate sites for refueling stations within the budget. This constraint is used to build AF refueling stations in a statewide network to showcase that the proposed model is well-suited for solving practical problems in Section 6.


## 4.3.1. Parameters

f
v
i
,
v
j: average traffic flow in round trips per time unit between interchanges vi and vj,α
k
,
m
t
,
m
t
′: identification coefficient for sites mt and mt′ for Type 2 k-th SNDF path.


## 4.3.2. Decision Variables.

(7)
x
m
t
=
1
,
i
f
r
e
f
u
e
l
i
n
g
s
t
a
t
i
o
n
l
o
c
a
t
i
o
n
m
t
∈
M
i
s
s
e
l
e
c
t
e
d
,
0
,
o
t
h
e
r
w
i
s
e
,(8)
y
k
,
v
i
,
v
j
=
1
,
i
f
t
h
e
k
−
t
h
S
N
D
F
p
a
t
h
f
r
o
m
v
i
t
o
v
j
i
s
c
o
v
e
r
e
d
b
y
t
h
e
c
u
r
r
e
n
t
s
e
t
o
f
r
e
f
u
e
l
i
n
g
s
t
a
t
i
o
n
s
,
0
,
o
t
h
e
r
w
i
s
e
,(9)
z
v
i
,
v
j
=
1
,
i
f
t
h
e
r
o
u
n
d
t
r
i
p
o
f
t
r
a
f
f
i
c
f
l
o
w
f
v
i
,
v
j
i
s
c
o
v
e
r
e
d
b
y
t
h
e
c
u
r
r
e
n
t
s
e
t
o
f
r
e
f
u
e
l
i
n
g
s
t
a
t
i
o
n
s
,
0
,
o
t
h
e
r
w
i
s
e
.Then, we propose the following a mixed-integer programming model:(10)maximize∑qvi,vj∈Qfvi,vjzvi,vj,subject to(11)∑mt∈SSFkvi,vjxmt≥yk,vi,vj,∀Fkvi,vj∈L^o,k=1,⋯,K,(12)∑mt∈ESFkvi,vjxmt≥yk,vi,vj,∀Fkvi,vj∈L^o,k=1,⋯,K,(13)∑mt∈SSFkvj,vixmt≥yk,vj,vi,∀Fkvj,vi∈L^r,k=1,⋯,K,(14)∑mt∈ESFkvj,vixmt≥yk,vj,vi,∀Fkvj,vi∈L^r,k=1,⋯,K,(15)∑mt∈CSFk2vi,vjαk,mt,mt′xmt≥yk,vi,vj,∀mt′∈USFk2vi,vj,∀Fk2vi,vj∈L^o,k=1,⋯,K,(16)∑mt∈CSFk2vj,viαk,mt,mt′xmt≥yk,vj,vi,∀mt′∈USFk2vj,vi,∀Fk2vj,vi∈L^r,k=1,⋯,K,(17)∑k=1,⋯,Kyk,vi,vj≥zvi,vj,∀qvi,vj∈Q,(18)∑k=1,⋯,Kyk,vj,vi≥zvi,vj,∀qvi,vj∈Q,(19)∑mt∈Mxmt=p,(20)xmt∈0,1,∀mt∈M,(21)yk,vi,vj∈0,1,∀Fkvi,vj∈L^o,k=1⋯K,(22)yk,vj,vi∈0,1,∀Fkvj,vi∈L^r,k=1⋯K,(23)zvi,vj∈0,1,∀qvi,vj∈Q.In this formulation, the objective function (10) maximizes the traffic flow that is covered by the set of AF refueling stations located at thep selected candidate sites on the transportation network. Constraint sets (11) and (12) guarantee that trips along K SNDF original paths of Types 1 and 2 receive refueling service from one station in set of candidate sites SSFkvi,vj and another one in set of candidate sites ESFkvi,vj. The k-th SNDF return paths have similar refueling conditions in Constraint sets (13) and (14). Next, coefficients in Constraint set (15) are determined by identification coefficients αk,mt,mt′ for mt∈CSFk2vi,vj and mt′∈USFk2vi,vj, where αk,mt,mt′=1 if dijkvi,mt<dijkvi,mt′ and 0<dijkmt,mt′≤R; otherwise, αk,mt,mt′=0. This constraint set guarantees that, if a refueling station is located at mt′ in USFk2vi,vj, then one of sites in CSFk2vi,vj should be selected to locate another station which would provide refueling for vehicles to reach the station at mt′. Similarly, Constraint set (16) is applied to candidate sites on Type 2 SNDF return paths. In addition, since round trips in OD pairs are considered in this problem, all possible combinations of original and return paths are determined by Constraint sets (17) and (18), respectively. These constraint sets ensure vehicles on a directed transportation network to make nonsymmetric round trips if necessary for refueling service. Constraint (19) allows exactly p AF refueling stations to be located at the predetermined candidate site locations. Lastly, the three kinds of binary decision variables are defined in Constraint sets (20) to (23).The main stakeholders of the proposed model that maximizes the covered traffic flows can be shareholders, investors, employers and employees, and any other private entities whose stake is directly or indirectly tied to this objective function, assuming that the revenue is proportional to the traffic flow covered by a given number of stations; that is, the more traffic flow covered by the stations, the higher revenue the stakeholders expect. If we assume that the construction cost does not significantly differ by station type or region, then Constraint (19) can play a role as the budget constraint. This assumption, however, would not be always true in practice [60]. So, we can relax this assumption by replacing Constraint (19) by the following constraint:(24)∑mt∈Mcmtxmt≤B,wherecmt refers to the capital cost for building a refueling station in candidate site mt∈M, and B the refueling infrastructure budget. Constraint (24) allows different capital costs by station type (e.g., single-access vs. dual-access) or by region and forces the model to determine candidate sites for refueling stations within the budget. This constraint is used to build AF refueling stations in a statewide network to showcase that the proposed model is well-suited for solving practical problems in Section 6.


## 5. Computational Experiments

The proposed mixed-integer programming model with thek-th SNDF path algorithm is applied to a well-known directed transportation network with 25 vertices, 43 arcs, and 300 OD pairs provided by Simchi-Levi and Berman [61]. Figure 5 shows the arc lengths next to each arc. It also displays 25 randomly selected candidate station locations, including 10 single- and 15 dual-access sites, to allocate the AF refueling stations in the test network. A single-access candidate site is symbolized by a triangular shape on one side of the arc, and a dual-access candidate site is depicted by a diamond shape in the middle of the arc. To analyze the effects of number of SNDF paths, maximum deviation distance, and vehicle driving range on the OD flow coverage, we consider three route options, K=1,3, and 5, three deviation factors, β=0%,50%, and 100%, and three driving ranges, R=15,20, and 30.Figure 5
Test network with 300 OD pairs and 25 candidate sites [61].All computational experiments in Section5 were conducted on an Intel i5 2.2 GHZ Dual-Core laptop with 12 GB RAM. Our computational experiments have three procedures. In the first procedure that is denoted as P1, to determine 2K SNDF original and return paths for all OD pairs, we ran the k-th SNDF path algorithm with the three limited driving ranges R (i.e., R= 15, 20, and 30) and the three deviation factors β (i.e., β= 0%, 50%, and 100%) using MATLAB R2013a. Note that the same deviation factors were applied to both original and return paths. Then, in the second procedure, denoted as P2, the mixed-integer programming model with the SNDF paths and set of candidate sites was programmed in MATLAB to generate CPLEX format files that follow the syntax rules of CPLEX. In the last procedure, denoted as P3, all the problems were solved by version 12.4 of CPLEX.Table4 provides the computational times for the three solution procedures above, P1, P2, and P3. In general, the overall computational time increases as the number of SNDF paths increases and the deviation factor is larger. Specifically, the large number of SNDF paths affects the computational time because it is directly related to the problem size. In the k-th SNDF path algorithm, 2K SNDF original and return paths are generated for all OD pairs, so that the number of decision variables and constraints in the mathematical model for this optimization application largely increase as K increases. However, the results show that, the longer the limited driving range, the shorter the computational time. This can be explained by the relationship between the limited driving range and Type 2 SNDF paths. That is, if the limited driving range is sufficiently large, the paths that satisfy the four properties of feasible paths are quickly determined and the number of constraints for Type 2 SNDF paths decreases, so that the computational time also decreases.Table 4
Computational time forp=1–25 (in sec).
K
R
β
=
0
β
=
50
%
β
=
100
%
P1
P2
P3
P1
P2
P3
P1
P2
P3
1
15
24.22
0.33
0.55
30.28
0.41
0.88
33.16
0.45
0.94
20
23.88
0.30
0.67
24.07
0.34
0.75
24.42
0.32
0.80
30
23.82
0.24
0.65
23.83
0.24
0.75
23.88
0.26
0.69
3
15
24.31
0.40
0.68
36.42
0.99
2.97
45.05
0.98
3.68
20
24.21
0.38
0.86
28.03
0.68
2.72
31.74
0.74
3.74
30
24.11
0.35
0.85
26.26
0.44
2.03
27.42
0.48
2.82
5
15
24.16
0.48
0.75
43.09
1.49
4.79
56.30
1.61
8.40
20
24.04
0.44
0.95
33.46
1.07
4.42
39.98
1.10
8.40
30
24.25
0.40
1.22
30.61
0.68
2.76
33.53
0.69
3.86


### 5.1. Effect of the Number of SNDF Paths

We first analyze the effect of the different number of SNDF paths on the coverage of OD traffic flows for each deviation factor. The three graphs in Figure6 show the coverage of traffic flows (in %) with β= 0%, 50%, and 100% of dPSvi,vj as we increase the number of refueling stations from p=1 to 25. Note that we fix R=20 in this subsection.Figure 6
Tradeoff between number of stations and coverage for different number of SNDF paths. (a)β=0%. (b) β=50%. (c) β=100%.
(a)
(b)
(c)First, whenβ=0%, the coverage is not significantly affected by the values of K, but a slight improvement is obtained between K=1 and K=3,5 when the number of refueling stations increases from p=3 to 15. For example, when p=8 and K=1, the first SNDF return path of qv3,v13, i.e., F1v13,v3=v13=v0,v8=v1,v4=v2,v3=v3 is not covered. In the case of K=3 and 5 with β=0%, the mixed-integer programming model locates one refueling station at site m14 and then F3v13,v3=v13=v0,v10=v1,v9=v2,v3=v3, which is the second SNDF return path of qv3,v13, is covered by the one refueling station. This indicates that some OD pairs have multiple shortest paths which satisfy the properties of feasible paths, so that all these paths can be selected in our model. Note that K=3 is enough for multiple SNDF paths to be considered in the network when β=0% because there is no difference in coverage between K=3 and K=5 from p=1 to 25. Second, as we allow drivers to take paths with long detour distances, the effect of the value of K on the coverage of traffic flows by a given number of refueling stations is significantly noticeable. In particular, the coverage improvement between K=1 and K=5 ranges from 0% to 8.31% when β=50%, while the improvement increases up to 13.07% when β=100%. The k-th SNDF path algorithm can generate more SNDF paths under β=100%, so that more paths can be considered to locate the optimal sets of refueling stations in the mixed-integer programming model. Then, more traffic flows using 2K SNDF original and return paths can be covered by a given number of stations. Note that the coverage improvement between different number of SNDF paths fades away slowly as the refueling network becomes mature. Lastly, in all cases, OD pairs are not fully covered when we locate refueling stations at all candidate sites. For example, considering K=5 with β=100% for both original and return paths, which is the maximum deviation factor in the experiments, 70.58% of all traffic flows can be refueled with p=25. This is because some refueling stations at single-access candidate sites can provide refueling service only to traffic flows in the same driving direction. Table 5 summarizes the coverage of traffic flows with β=0%, 50%, and 100% for p=1,⋯,25 when we allow a different number of SNDF paths for each OD pair.Table 5
Summary of coverages for different number of SNDF paths withβ=0%, β=50%, and β=100% when R=20.
p
β
=
0
%
β
=
50
%
β
=
100
%
K
=
1
K
=
3
K
=
5
K
=
1
K
=
3
K
=
5
K
=
1
K
=
3
K
=
5
1
17.33
17.33
17.33
17.33
17.33
17.33
17.33
18.03
18.87
2
30.14
30.50
30.50
30.14
30.75
32.33
30.15
34.80
39.74
3
35.09
35.45
35.45
35.09
37.41
39.96
36.82
43.79
49.89
4
39.26
39.91
39.91
39.72
42.81
45.35
42.58
50.65
55.37
5
43.57
44.23
44.23
44.18
46.78
49.80
47.27
56.12
59.82
6
45.27
46.73
46.73
46.50
49.91
52.87
51.96
59.95
64.05
7
46.32
48.43
48.43
48.23
53.35
55.44
53.84
63.19
66.28
8
47.36
49.65
49.65
49.76
55.26
58.00
55.58
64.76
67.61
9
48.27
50.41
50.41
51.23
57.74
59.54
57.11
66.86
68.74
10
49.32
50.98
50.98
52.92
58.98
60.45
58.53
68.48
69.47
11
50.36
51.52
51.52
54.45
59.92
61.08
60.27
69.13
69.90
12
50.98
51.86
51.86
55.92
60.52
61.47
61.80
69.81
70.17
13
51.50
51.93
51.93
56.96
61.11
61.85
62.89
70.20
70.47
14
51.86
52.05
52.05
58.17
61.70
61.95
64.19
70.39
70.52
15
52.05
52.12
52.12
59.11
61.94
62.13
65.31
70.49
70.56
16
52.12
52.15
52.15
59.74
62.13
62.19
66.46
70.53
70.58
17
52.15
52.15
52.15
60.70
62.19
62.24
67.57
70.55
70.58
18
52.15
52.16
52.16
61.29
62.24
62.28
68.53
70.57
70.58
19
52.16
52.16
52.16
61.74
62.28
62.31
69.27
70.58
70.58
20
52.16
52.16
52.16
62.07
62.31
62.34
69.92
70.58
70.58
21
52.16
52.16
52.16
62.26
62.34
62.34
70.27
70.58
70.58
22
52.16
52.16
52.16
62.30
62.34
62.34
70.43
70.58
70.58
23
52.16
52.16
52.16
62.34
62.34
62.34
70.51
70.58
70.58
24
52.16
52.16
52.16
62.34
62.34
62.34
70.55
70.58
70.58
25
52.16
52.16
52.16
62.34
62.34
62.34
70.58
70.58
70.58


### 5.2. Effect of the Maximum Deviation Distance

This subsection analyzes how much traffic flow makes detours from their preplanned paths for refueling and how long is the average deviation distance of covered detour traffic flows. For each OD pair, vehicles are able to deviate from shortest original pathPSvi,vj and shortest return path PSvj,vi up to dPSvi,vj1+β and dPSvj,vi1+β, respectively, in order to reach AF refueling stations. Table 6 shows the proportions of the covered detour traffic flows with respect to the total flow covered, denoted as Cβ, and the average deviation factor values of the vehicles that deviate from their shortest paths, denoted as β¯, when the number of refueling stations ranges from p= 1 to 25. Note that we fix R=20 and consider K= 1, 3, and 5 in the table. For example, given p=10, when vehicles are allowed to deviate to alternative paths whose lengths are at most 50% longer than those of the corresponding shortest paths and K=1, the covered detour flows account for only 12.47% of the total covered flows and their average deviation factor is 19.99% of their shortest path lengths. First, regardless of the value of β, we observe that the value of Cβ generally increases from p= 1 to 25 when we consider only one SNDF path between ODs. This represents that the deviation factor largely affects the coverage if vehicles have few options to select their route regardless of the number of refueling stations. However, with K= 3 and 5, the value of Cβ tends to increase when refueling stations are scarce, while the value of Cβ starts decreasing after the traffic flows are mostly covered. These patterns imply that the effect of β on the traffic flow covered gradually diminishes as the number of refueling stations increases, especially when vehicles have multiple path options between OD pairs. Next, regarding β¯, we have an extreme case where the value of β¯ drops rapidly when we consider β=100% and the number of refueling stations is small. However, in most cases, the value of β¯ does not change significantly as the coverage of traffic flows increases. This result can be explained by the predetermined locations of candidate sites. In other words, even though we allow trips along many feasible paths with long detour distances between OD pairs, the coverage proportion of detour flows is usually small with respect to all traffic flows because the candidate sites are predetermined before establishing AF refueling stations.Table 6
Coverage proportion of detour flows among the total covered flowsCβ and the average deviation factor β¯ for p=1,⋯,25.
p
K= 1
K= 3
K= 5
β = 0%
β = 50%
β = 100%
β = 0%
β = 50%
β = 100%
β = 0%
β = 50%
β = 100%
C (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C
β (%)
β
¯ (%)
1
17.33
17.33
0.00
0.00
17.33
0.00
0.00
17.33
17.33
0.00
0.00
18.03
100.00
57.61
21.23
21.23
0.00
0.00
32.28
21.24
63.64
2
30.14
30.14
0.00
0.00
30.15
0.05
85.71
30.50
30.75
1.75
20.23
34.80
12.36
65.61
33.98
34.49
6.29
24.86
42.70
22.99
52.05
3
35.09
35.09
0.00
0.00
36.82
19.04
73.97
35.45
37.41
8.96
21.22
43.79
21.75
45.41
39.76
45.54
12.61
23.47
54.89
29.39
45.08
4
39.26
39.72
1.14
28.57
42.58
18.86
69.07
39.91
42.81
8.89
21.64
50.65
32.34
46.22
44.61
51.61
13.62
23.42
62.44
29.91
45.61
5
43.57
44.18
1.39
20.81
47.27
18.13
63.62
44.23
46.78
13.39
20.69
56.12
30.14
46.84
48.90
54.87
17.69
24.75
65.64
30.31
37.56
6
45.27
46.50
6.30
24.15
51.96
18.04
53.43
46.73
49.91
16.38
22.14
59.95
33.04
38.36
51.08
57.19
18.03
24.77
68.18
35.41
36.81
7
46.32
48.23
6.15
27.32
53.84
18.24
54.83
48.43
53.35
14.90
21.70
63.19
34.09
34.61
52.95
58.53
13.25
25.41
69.17
39.12
36.65
8
47.36
49.76
5.96
27.32
55.58
17.78
53.74
49.65
55.26
26.08
24.42
64.76
36.89
36.83
53.72
59.52
29.08
25.86
69.75
38.73
35.38
9
48.27
51.23
5.79
27.32
57.11
17.27
53.74
50.41
57.74
24.84
24.80
66.86
38.71
34.47
54.33
60.10
27.06
25.16
70.24
39.43
37.71
10
49.32
52.92
12.47
19.99
58.53
23.50
37.31
50.98
58.98
25.49
23.38
68.48
39.69
36.11
54.90
61.09
27.30
26.90
70.45
40.35
38.14
11
50.36
54.45
12.12
19.99
60.27
22.92
38.49
51.52
59.92
25.75
23.74
69.13
40.22
35.37
55.29
61.67
28.05
26.30
70.51
33.06
39.27
12
50.98
55.92
11.80
19.99
61.80
22.32
38.49
51.86
60.52
27.28
23.81
69.81
41.92
35.15
55.85
61.91
27.33
26.53
70.54
45.81
35.71
13
51.50
56.96
11.59
19.99
62.89
22.78
41.57
51.93
61.11
18.86
24.09
70.20
40.88
33.65
56.24
62.11
20.87
25.77
70.57
43.36
35.80
14
51.86
58.17
13.25
21.01
64.19
23.92
41.63
52.05
61.70
20.74
24.24
70.39
41.90
34.89
56.51
62.21
21.54
24.22
70.58
38.06
33.65
15
52.05
59.11
13.87
19.02
65.31
25.10
42.99
52.12
61.94
20.46
23.84
70.49
36.65
31.89
56.58
62.26
20.76
24.75
70.58
38.41
33.36
16
52.12
59.74
14.04
21.57
66.46
25.28
38.95
52.15
62.13
23.44
24.39
70.53
38.38
31.42
56.70
62.31
22.37
23.28
70.58
40.33
31.71
17
52.15
60.70
14.55
19.72
67.57
26.41
40.26
52.15
62.19
23.45
23.31
70.55
39.53
32.32
56.77
62.35
23.33
23.33
70.58
45.25
35.78
18
52.15
61.29
15.36
18.39
68.53
26.87
36.59
52.16
62.24
22.96
23.30
70.57
38.36
30.40
56.80
62.35
21.15
23.26
70.58
36.90
31.69
19
52.16
61.74
15.88
17.50
69.27
27.87
34.25
52.16
62.28
22.70
20.51
70.58
36.36
29.92
56.80
62.35
23.85
22.99
70.58
34.54
34.33
20
52.16
62.07
16.33
18.53
69.92
28.66
32.56
52.16
62.31
22.78
20.30
70.58
35.14
32.48
56.80
62.35
22.90
21.92
70.58
32.92
33.56
21
52.16
62.26
16.28
18.53
70.27
29.14
33.38
52.16
62.34
22.05
19.65
70.58
34.94
35.54
56.80
62.35
21.38
20.93
70.58
32.04
33.56
22
52.16
62.30
16.33
18.31
70.43
29.38
34.06
52.16
62.34
21.72
20.35
70.58
30.93
31.49
56.80
62.35
21.17
20.66
70.58
30.93
31.46
23
52.16
62.34
16.32
19.15
70.51
29.49
35.22
52.16
62.34
21.46
20.20
70.58
28.45
35.63
56.80
62.35
21.46
20.43
70.58
28.45
35.63
24
52.16
62.34
16.33
18.35
70.55
29.55
34.43
52.16
62.34
24.73
19.52
70.58
27.95
34.50
56.80
62.35
23.61
21.22
70.58
29.67
33.39
25
52.16
62.34
16.33
18.35
70.58
29.55
34.16
52.16
62.34
16.33
18.35
70.58
26.10
34.16
56.80
62.35
16.33
18.35
70.58
26.10
34.16


### 5.3. Effect of the Limited Driving Range

In order to observe the effect of the limited driving range on the coverage for different deviation factors, Figures7(a) and 7(c) provide comparisons of optimal coverages of no-detour versus detour flows when R=15 and 30 for p=1,⋯,25. We set the coverages of traffic flows with β=0%, i.e., the coverage of no-detour flows, as the reference value of comparisons. Then, Figures 7(b) and 7(d) show their differences between the coverages of no-detour and detour flows for β= 50% and 100%. Note that we fix K=3 in all the cases. In general, compared to the coverage of no-detour flows, the improvement in the coverage of detour flows is easily observed for both R=15 and R=30. For example, at most an additional 15% of traffic flows can be covered if vehicles with R=15 are allowed to travel an additional distance of up to 50% of dPSvi,vj for all OD pairs. The coverage improvement rises as the deviation factor increases for both R=15 and R=30. On the other hand, the trends of improvement for R=15 and R=30 have different patterns. When R=15, the coverage difference between no-detour (β=0%) and detour flows (β=50% or β=100%) continuously increases until the number of refueling stations reaches to some level (p=21 for β=50%, p=22 for β=100%), and then the coverage difference is steady. This implies that, if an additional refueling station is added to a network, the coverge of vehicles with short driving ranges can increase significantly when a long deviation distance is allowable. Next, when R=30, the coverage difference between β=0% and β=50% is insignificant. However, if vehicles are able to deviate further from their preplanned paths, i.e., β=100%, the trends of improvement fluctuate according to the number of refueling stations, especially for p= 1–8. This shows that a significant coverage improvement can be expected by adding one refueling station at a time on a network where the infrastructure of AF refueling stations is insufficient if vehicles with a long driving range are allowed to make adequate detours. Table 7 summarizes comparisons of optimal coverages of no-detour versus detour flows and their coverage differences for R=15, 30 when K=3.Figure 7
Comparisons of the optimal coverages of no-detour versus detour flows and their coverage differences whenR=15,30 with K=3. (a) Coverages of no-detour versus detour flows with β=50%. (b) Coverage differences between no-detour and detour flows for R=15,30 and β=50%. (c) Coverages of no-detour versus detour flows with β=100%. (d) Coverage differences between no-detour and detour flows for R=15,30 and β=100%.
(a)
(b)
(c)
(d)Table 7
Comparisons of the optimal coverages of no-detour (β=0%) versus detour flows (β=50%,100%) and their coverage differences when R=15,30.
p
R
=
15
R
=
30
β
Differences between 0% and 50%
Differences between 0% and 100%
β
Differences between 0% and 50%
Differences between 0% and 100%
0%
50%
100%
0%
50%
100%
1
10.38
10.38
10.88
0.00
0.50
21.23
21.23
28.12
0.00
6.89
2
16.80
16.88
19.93
0.08
3.13
33.98
34.18
38.56
0.20
4.58
3
20.68
21.21
27.73
0.53
7.05
39.76
44.29
49.07
4.53
9.31
4
24.47
25.02
32.14
0.55
7.67
44.61
50.35
55.26
5.74
10.65
5
27.28
29.44
37.14
2.16
9.86
48.90
53.71
60.00
4.81
11.10
6
28.64
32.59
41.74
3.95
13.10
51.08
56.44
63.33
5.36
12.25
7
29.44
34.82
44.54
5.38
15.10
52.95
57.88
65.91
4.93
12.96
8
29.94
37.13
47.15
7.19
17.21
53.72
59.09
67.70
5.37
13.98
9
30.42
39.32
49.41
8.90
18.99
54.33
59.62
69.25
5.29
14.92
10
30.80
41.70
51.22
10.90
20.42
54.90
60.73
69.77
5.83
14.87
11
31.09
43.43
52.43
12.34
21.34
55.29
61.26
70.04
5.97
14.75
12
31.31
44.22
53.23
12.91
21.92
55.85
61.54
70.34
5.69
14.49
13
31.42
45.12
53.86
13.70
22.44
56.24
61.78
70.51
5.54
14.27
14
31.63
45.54
54.45
13.91
22.82
56.51
62.01
70.54
5.50
14.03
15
31.72
45.83
54.80
14.11
23.08
56.58
62.19
70.58
5.61
14.00
16
31.79
46.22
55.10
14.43
23.31
56.70
62.26
70.58
5.56
13.88
17
31.82
46.42
55.24
14.60
23.42
56.77
62.31
70.58
5.54
13.81
18
31.82
46.59
55.36
14.77
23.54
56.80
62.35
70.58
5.55
13.78
19
31.82
46.74
55.38
14.92
23.56
56.80
62.35
70.58
5.55
13.78
20
31.82
46.78
55.45
14.96
23.63
56.80
62.35
70.58
5.55
13.78
21
31.82
46.81
55.47
14.99
23.65
56.80
62.35
70.58
5.55
13.78
22
31.82
46.82
55.48
15.00
23.66
56.80
62.35
70.58
5.55
13.78
23
31.82
46.82
55.49
15.00
23.67
56.80
62.35
70.58
5.55
13.78
24
31.82
46.82
55.49
15.00
23.67
56.80
62.35
70.58
5.55
13.78
25
31.82
46.82
55.49
15.00
23.67
56.80
62.35
70.58
5.55
13.78


## 5.1. Effect of the Number of SNDF Paths

We first analyze the effect of the different number of SNDF paths on the coverage of OD traffic flows for each deviation factor. The three graphs in Figure6 show the coverage of traffic flows (in %) with β= 0%, 50%, and 100% of dPSvi,vj as we increase the number of refueling stations from p=1 to 25. Note that we fix R=20 in this subsection.Figure 6
Tradeoff between number of stations and coverage for different number of SNDF paths. (a)β=0%. (b) β=50%. (c) β=100%.
(a)
(b)
(c)First, whenβ=0%, the coverage is not significantly affected by the values of K, but a slight improvement is obtained between K=1 and K=3,5 when the number of refueling stations increases from p=3 to 15. For example, when p=8 and K=1, the first SNDF return path of qv3,v13, i.e., F1v13,v3=v13=v0,v8=v1,v4=v2,v3=v3 is not covered. In the case of K=3 and 5 with β=0%, the mixed-integer programming model locates one refueling station at site m14 and then F3v13,v3=v13=v0,v10=v1,v9=v2,v3=v3, which is the second SNDF return path of qv3,v13, is covered by the one refueling station. This indicates that some OD pairs have multiple shortest paths which satisfy the properties of feasible paths, so that all these paths can be selected in our model. Note that K=3 is enough for multiple SNDF paths to be considered in the network when β=0% because there is no difference in coverage between K=3 and K=5 from p=1 to 25. Second, as we allow drivers to take paths with long detour distances, the effect of the value of K on the coverage of traffic flows by a given number of refueling stations is significantly noticeable. In particular, the coverage improvement between K=1 and K=5 ranges from 0% to 8.31% when β=50%, while the improvement increases up to 13.07% when β=100%. The k-th SNDF path algorithm can generate more SNDF paths under β=100%, so that more paths can be considered to locate the optimal sets of refueling stations in the mixed-integer programming model. Then, more traffic flows using 2K SNDF original and return paths can be covered by a given number of stations. Note that the coverage improvement between different number of SNDF paths fades away slowly as the refueling network becomes mature. Lastly, in all cases, OD pairs are not fully covered when we locate refueling stations at all candidate sites. For example, considering K=5 with β=100% for both original and return paths, which is the maximum deviation factor in the experiments, 70.58% of all traffic flows can be refueled with p=25. This is because some refueling stations at single-access candidate sites can provide refueling service only to traffic flows in the same driving direction. Table 5 summarizes the coverage of traffic flows with β=0%, 50%, and 100% for p=1,⋯,25 when we allow a different number of SNDF paths for each OD pair.Table 5
Summary of coverages for different number of SNDF paths withβ=0%, β=50%, and β=100% when R=20.
p
β
=
0
%
β
=
50
%
β
=
100
%
K
=
1
K
=
3
K
=
5
K
=
1
K
=
3
K
=
5
K
=
1
K
=
3
K
=
5
1
17.33
17.33
17.33
17.33
17.33
17.33
17.33
18.03
18.87
2
30.14
30.50
30.50
30.14
30.75
32.33
30.15
34.80
39.74
3
35.09
35.45
35.45
35.09
37.41
39.96
36.82
43.79
49.89
4
39.26
39.91
39.91
39.72
42.81
45.35
42.58
50.65
55.37
5
43.57
44.23
44.23
44.18
46.78
49.80
47.27
56.12
59.82
6
45.27
46.73
46.73
46.50
49.91
52.87
51.96
59.95
64.05
7
46.32
48.43
48.43
48.23
53.35
55.44
53.84
63.19
66.28
8
47.36
49.65
49.65
49.76
55.26
58.00
55.58
64.76
67.61
9
48.27
50.41
50.41
51.23
57.74
59.54
57.11
66.86
68.74
10
49.32
50.98
50.98
52.92
58.98
60.45
58.53
68.48
69.47
11
50.36
51.52
51.52
54.45
59.92
61.08
60.27
69.13
69.90
12
50.98
51.86
51.86
55.92
60.52
61.47
61.80
69.81
70.17
13
51.50
51.93
51.93
56.96
61.11
61.85
62.89
70.20
70.47
14
51.86
52.05
52.05
58.17
61.70
61.95
64.19
70.39
70.52
15
52.05
52.12
52.12
59.11
61.94
62.13
65.31
70.49
70.56
16
52.12
52.15
52.15
59.74
62.13
62.19
66.46
70.53
70.58
17
52.15
52.15
52.15
60.70
62.19
62.24
67.57
70.55
70.58
18
52.15
52.16
52.16
61.29
62.24
62.28
68.53
70.57
70.58
19
52.16
52.16
52.16
61.74
62.28
62.31
69.27
70.58
70.58
20
52.16
52.16
52.16
62.07
62.31
62.34
69.92
70.58
70.58
21
52.16
52.16
52.16
62.26
62.34
62.34
70.27
70.58
70.58
22
52.16
52.16
52.16
62.30
62.34
62.34
70.43
70.58
70.58
23
52.16
52.16
52.16
62.34
62.34
62.34
70.51
70.58
70.58
24
52.16
52.16
52.16
62.34
62.34
62.34
70.55
70.58
70.58
25
52.16
52.16
52.16
62.34
62.34
62.34
70.58
70.58
70.58


## 5.2. Effect of the Maximum Deviation Distance

This subsection analyzes how much traffic flow makes detours from their preplanned paths for refueling and how long is the average deviation distance of covered detour traffic flows. For each OD pair, vehicles are able to deviate from shortest original pathPSvi,vj and shortest return path PSvj,vi up to dPSvi,vj1+β and dPSvj,vi1+β, respectively, in order to reach AF refueling stations. Table 6 shows the proportions of the covered detour traffic flows with respect to the total flow covered, denoted as Cβ, and the average deviation factor values of the vehicles that deviate from their shortest paths, denoted as β¯, when the number of refueling stations ranges from p= 1 to 25. Note that we fix R=20 and consider K= 1, 3, and 5 in the table. For example, given p=10, when vehicles are allowed to deviate to alternative paths whose lengths are at most 50% longer than those of the corresponding shortest paths and K=1, the covered detour flows account for only 12.47% of the total covered flows and their average deviation factor is 19.99% of their shortest path lengths. First, regardless of the value of β, we observe that the value of Cβ generally increases from p= 1 to 25 when we consider only one SNDF path between ODs. This represents that the deviation factor largely affects the coverage if vehicles have few options to select their route regardless of the number of refueling stations. However, with K= 3 and 5, the value of Cβ tends to increase when refueling stations are scarce, while the value of Cβ starts decreasing after the traffic flows are mostly covered. These patterns imply that the effect of β on the traffic flow covered gradually diminishes as the number of refueling stations increases, especially when vehicles have multiple path options between OD pairs. Next, regarding β¯, we have an extreme case where the value of β¯ drops rapidly when we consider β=100% and the number of refueling stations is small. However, in most cases, the value of β¯ does not change significantly as the coverage of traffic flows increases. This result can be explained by the predetermined locations of candidate sites. In other words, even though we allow trips along many feasible paths with long detour distances between OD pairs, the coverage proportion of detour flows is usually small with respect to all traffic flows because the candidate sites are predetermined before establishing AF refueling stations.Table 6
Coverage proportion of detour flows among the total covered flowsCβ and the average deviation factor β¯ for p=1,⋯,25.
p
K= 1
K= 3
K= 5
β = 0%
β = 50%
β = 100%
β = 0%
β = 50%
β = 100%
β = 0%
β = 50%
β = 100%
C (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C (%)
C
β (%)
β
¯ (%)
C (%)
C
β (%)
β
¯ (%)
1
17.33
17.33
0.00
0.00
17.33
0.00
0.00
17.33
17.33
0.00
0.00
18.03
100.00
57.61
21.23
21.23
0.00
0.00
32.28
21.24
63.64
2
30.14
30.14
0.00
0.00
30.15
0.05
85.71
30.50
30.75
1.75
20.23
34.80
12.36
65.61
33.98
34.49
6.29
24.86
42.70
22.99
52.05
3
35.09
35.09
0.00
0.00
36.82
19.04
73.97
35.45
37.41
8.96
21.22
43.79
21.75
45.41
39.76
45.54
12.61
23.47
54.89
29.39
45.08
4
39.26
39.72
1.14
28.57
42.58
18.86
69.07
39.91
42.81
8.89
21.64
50.65
32.34
46.22
44.61
51.61
13.62
23.42
62.44
29.91
45.61
5
43.57
44.18
1.39
20.81
47.27
18.13
63.62
44.23
46.78
13.39
20.69
56.12
30.14
46.84
48.90
54.87
17.69
24.75
65.64
30.31
37.56
6
45.27
46.50
6.30
24.15
51.96
18.04
53.43
46.73
49.91
16.38
22.14
59.95
33.04
38.36
51.08
57.19
18.03
24.77
68.18
35.41
36.81
7
46.32
48.23
6.15
27.32
53.84
18.24
54.83
48.43
53.35
14.90
21.70
63.19
34.09
34.61
52.95
58.53
13.25
25.41
69.17
39.12
36.65
8
47.36
49.76
5.96
27.32
55.58
17.78
53.74
49.65
55.26
26.08
24.42
64.76
36.89
36.83
53.72
59.52
29.08
25.86
69.75
38.73
35.38
9
48.27
51.23
5.79
27.32
57.11
17.27
53.74
50.41
57.74
24.84
24.80
66.86
38.71
34.47
54.33
60.10
27.06
25.16
70.24
39.43
37.71
10
49.32
52.92
12.47
19.99
58.53
23.50
37.31
50.98
58.98
25.49
23.38
68.48
39.69
36.11
54.90
61.09
27.30
26.90
70.45
40.35
38.14
11
50.36
54.45
12.12
19.99
60.27
22.92
38.49
51.52
59.92
25.75
23.74
69.13
40.22
35.37
55.29
61.67
28.05
26.30
70.51
33.06
39.27
12
50.98
55.92
11.80
19.99
61.80
22.32
38.49
51.86
60.52
27.28
23.81
69.81
41.92
35.15
55.85
61.91
27.33
26.53
70.54
45.81
35.71
13
51.50
56.96
11.59
19.99
62.89
22.78
41.57
51.93
61.11
18.86
24.09
70.20
40.88
33.65
56.24
62.11
20.87
25.77
70.57
43.36
35.80
14
51.86
58.17
13.25
21.01
64.19
23.92
41.63
52.05
61.70
20.74
24.24
70.39
41.90
34.89
56.51
62.21
21.54
24.22
70.58
38.06
33.65
15
52.05
59.11
13.87
19.02
65.31
25.10
42.99
52.12
61.94
20.46
23.84
70.49
36.65
31.89
56.58
62.26
20.76
24.75
70.58
38.41
33.36
16
52.12
59.74
14.04
21.57
66.46
25.28
38.95
52.15
62.13
23.44
24.39
70.53
38.38
31.42
56.70
62.31
22.37
23.28
70.58
40.33
31.71
17
52.15
60.70
14.55
19.72
67.57
26.41
40.26
52.15
62.19
23.45
23.31
70.55
39.53
32.32
56.77
62.35
23.33
23.33
70.58
45.25
35.78
18
52.15
61.29
15.36
18.39
68.53
26.87
36.59
52.16
62.24
22.96
23.30
70.57
38.36
30.40
56.80
62.35
21.15
23.26
70.58
36.90
31.69
19
52.16
61.74
15.88
17.50
69.27
27.87
34.25
52.16
62.28
22.70
20.51
70.58
36.36
29.92
56.80
62.35
23.85
22.99
70.58
34.54
34.33
20
52.16
62.07
16.33
18.53
69.92
28.66
32.56
52.16
62.31
22.78
20.30
70.58
35.14
32.48
56.80
62.35
22.90
21.92
70.58
32.92
33.56
21
52.16
62.26
16.28
18.53
70.27
29.14
33.38
52.16
62.34
22.05
19.65
70.58
34.94
35.54
56.80
62.35
21.38
20.93
70.58
32.04
33.56
22
52.16
62.30
16.33
18.31
70.43
29.38
34.06
52.16
62.34
21.72
20.35
70.58
30.93
31.49
56.80
62.35
21.17
20.66
70.58
30.93
31.46
23
52.16
62.34
16.32
19.15
70.51
29.49
35.22
52.16
62.34
21.46
20.20
70.58
28.45
35.63
56.80
62.35
21.46
20.43
70.58
28.45
35.63
24
52.16
62.34
16.33
18.35
70.55
29.55
34.43
52.16
62.34
24.73
19.52
70.58
27.95
34.50
56.80
62.35
23.61
21.22
70.58
29.67
33.39
25
52.16
62.34
16.33
18.35
70.58
29.55
34.16
52.16
62.34
16.33
18.35
70.58
26.10
34.16
56.80
62.35
16.33
18.35
70.58
26.10
34.16


## 5.3. Effect of the Limited Driving Range

In order to observe the effect of the limited driving range on the coverage for different deviation factors, Figures7(a) and 7(c) provide comparisons of optimal coverages of no-detour versus detour flows when R=15 and 30 for p=1,⋯,25. We set the coverages of traffic flows with β=0%, i.e., the coverage of no-detour flows, as the reference value of comparisons. Then, Figures 7(b) and 7(d) show their differences between the coverages of no-detour and detour flows for β= 50% and 100%. Note that we fix K=3 in all the cases. In general, compared to the coverage of no-detour flows, the improvement in the coverage of detour flows is easily observed for both R=15 and R=30. For example, at most an additional 15% of traffic flows can be covered if vehicles with R=15 are allowed to travel an additional distance of up to 50% of dPSvi,vj for all OD pairs. The coverage improvement rises as the deviation factor increases for both R=15 and R=30. On the other hand, the trends of improvement for R=15 and R=30 have different patterns. When R=15, the coverage difference between no-detour (β=0%) and detour flows (β=50% or β=100%) continuously increases until the number of refueling stations reaches to some level (p=21 for β=50%, p=22 for β=100%), and then the coverage difference is steady. This implies that, if an additional refueling station is added to a network, the coverge of vehicles with short driving ranges can increase significantly when a long deviation distance is allowable. Next, when R=30, the coverage difference between β=0% and β=50% is insignificant. However, if vehicles are able to deviate further from their preplanned paths, i.e., β=100%, the trends of improvement fluctuate according to the number of refueling stations, especially for p= 1–8. This shows that a significant coverage improvement can be expected by adding one refueling station at a time on a network where the infrastructure of AF refueling stations is insufficient if vehicles with a long driving range are allowed to make adequate detours. Table 7 summarizes comparisons of optimal coverages of no-detour versus detour flows and their coverage differences for R=15, 30 when K=3.Figure 7
Comparisons of the optimal coverages of no-detour versus detour flows and their coverage differences whenR=15,30 with K=3. (a) Coverages of no-detour versus detour flows with β=50%. (b) Coverage differences between no-detour and detour flows for R=15,30 and β=50%. (c) Coverages of no-detour versus detour flows with β=100%. (d) Coverage differences between no-detour and detour flows for R=15,30 and β=100%.
(a)
(b)
(c)
(d)Table 7
Comparisons of the optimal coverages of no-detour (β=0%) versus detour flows (β=50%,100%) and their coverage differences when R=15,30.
p
R
=
15
R
=
30
β
Differences between 0% and 50%
Differences between 0% and 100%
β
Differences between 0% and 50%
Differences between 0% and 100%
0%
50%
100%
0%
50%
100%
1
10.38
10.38
10.88
0.00
0.50
21.23
21.23
28.12
0.00
6.89
2
16.80
16.88
19.93
0.08
3.13
33.98
34.18
38.56
0.20
4.58
3
20.68
21.21
27.73
0.53
7.05
39.76
44.29
49.07
4.53
9.31
4
24.47
25.02
32.14
0.55
7.67
44.61
50.35
55.26
5.74
10.65
5
27.28
29.44
37.14
2.16
9.86
48.90
53.71
60.00
4.81
11.10
6
28.64
32.59
41.74
3.95
13.10
51.08
56.44
63.33
5.36
12.25
7
29.44
34.82
44.54
5.38
15.10
52.95
57.88
65.91
4.93
12.96
8
29.94
37.13
47.15
7.19
17.21
53.72
59.09
67.70
5.37
13.98
9
30.42
39.32
49.41
8.90
18.99
54.33
59.62
69.25
5.29
14.92
10
30.80
41.70
51.22
10.90
20.42
54.90
60.73
69.77
5.83
14.87
11
31.09
43.43
52.43
12.34
21.34
55.29
61.26
70.04
5.97
14.75
12
31.31
44.22
53.23
12.91
21.92
55.85
61.54
70.34
5.69
14.49
13
31.42
45.12
53.86
13.70
22.44
56.24
61.78
70.51
5.54
14.27
14
31.63
45.54
54.45
13.91
22.82
56.51
62.01
70.54
5.50
14.03
15
31.72
45.83
54.80
14.11
23.08
56.58
62.19
70.58
5.61
14.00
16
31.79
46.22
55.10
14.43
23.31
56.70
62.26
70.58
5.56
13.88
17
31.82
46.42
55.24
14.60
23.42
56.77
62.31
70.58
5.54
13.81
18
31.82
46.59
55.36
14.77
23.54
56.80
62.35
70.58
5.55
13.78
19
31.82
46.74
55.38
14.92
23.56
56.80
62.35
70.58
5.55
13.78
20
31.82
46.78
55.45
14.96
23.63
56.80
62.35
70.58
5.55
13.78
21
31.82
46.81
55.47
14.99
23.65
56.80
62.35
70.58
5.55
13.78
22
31.82
46.82
55.48
15.00
23.66
56.80
62.35
70.58
5.55
13.78
23
31.82
46.82
55.49
15.00
23.67
56.80
62.35
70.58
5.55
13.78
24
31.82
46.82
55.49
15.00
23.67
56.80
62.35
70.58
5.55
13.78
25
31.82
46.82
55.49
15.00
23.67
56.80
62.35
70.58
5.55
13.78


## 6. Case Study: The State of Pennsylvania

In this section, we apply the proposed model with different capital costs to build LNG refueling stations for long haul logistics on a highway system by their station type (single-access vs. dual-access) into the state of Pennsylvania with 2,211 OD pairs, so as to demonstrate the performance in a large-size problem.From the Topologically Integrated Geographic Encoding and Referencing (TIGER) database, we first make a map of Pennsylvania with center points of 67 counties and centroid-to-centroid routes, as shown in Figure8. Then, we randomly add 50 of single- and 50 of dual-access candidate locations for LNG refueling stations. We use $2.75 M and $3.40 M as the construction cost of single- and dual-access LNG refueling stations, respectively, which are derived as the median values of the estimated construction costs to develop LNG refueling infrastructure in the Pennsylvania Turnpike by Ventura et al. [5]. The OD flows between counties are estimated by a simple gravity spatial interaction model [62]. We use the Freight Analysis Framework Version 4 long distance truck volume for the year of 2045 in the state of Pennsylvania, which was developed by the Federal Highway Administration in 2016 [63]. Also, to consider a conservative range for current and older models of LNG trucks with a single fuel tank [64] and improving technologies of dual fuel tank systems [65], we set R=300 and 600 for LNG trucks. Lastly, we conduct this experiment by using Intel(R) i7-6700 K CPU 4.00 GHz and 16 GB RAM desktop PC.Figure 8
A road network in the state of Pennsylvania with 100 candidate sites for LNG stations.Figure 9
Trade-off between the optimal coverage and the construction budget for no-detour (β=0%) versus detour flows (β=50%) with K=1,3 when R=300,600. (a) R=300. (b) R=600.
(a)
(b)Table8 provides comparisons of number of selected stations and optimal coverage for no-detour versus 50-percent detour flows with K=1 and 3 when R=300 and 600. Also, Figure 9 shows trade-off between budget and optimal coverage for the two LNG truck driving ranges (R=300,600) and three route options (β=0%, β=50% with K=1, and β=50% with K=3). We apply our coverage maximization model with the budget constraint to the state of Pennsylvania using the budget range from $10 M to $200 M in $10 M increments.Table 8
Comparisons of the number of stations and the optimal coverage for no-detour (β=0%) versus detour flows (β=50%) with K=1,3 when R=300,600.
Budget ($M)
R
=
300
R
=
600
No-detour withβ=0%
Detour withK=1 and β=50%
Detour withK=3 and β=50%
No-detour withβ=0%
Detour withK=1 and β=50%
Detour withK=3 and β=50%
Single (#)
Dual (#)
Coverage (%)
Single (#)
Dual (#)
Coverage (%)
Single (#)
Dual (#)
Coverage (%)
Single (#)
Dual (#)
Coverage (%)
Single (#)
Dual (#)
Coverage (%)
Single (#)
Dual (#)
Coverage (%)
10
0
2
19.99
1
2
21.95
1
2
27.91
0
2
21.14
1
2
23.10
1
2
29.40
20
1
5
31.23
1
5
33.55
1
5
42.95
1
5
33.11
1
5
35.15
1
5
45.21
30
1
8
39.71
1
8
45.31
1
8
54.56
1
8
41.56
1
8
44.18
1
8
56.57
40
0
11
46.70
0
11
49.21
0
11
62.34
0
11
48.48
0
11
51.22
0
11
64.13
50
0
14
52.65
0
14
55.74
2
13
68.29
0
14
54.65
0
14
57.68
0
14
69.94
60
2
16
57.58
2
16
60.89
2
16
72.62
0
17
59.34
0
17
62.61
2
16
73.98
70
0
20
61.11
3
18
64.81
3
18
75.21
0
20
63.06
0
20
66.56
0
20
76.63
80
0
23
64.34
3
21
68.33
3
21
77.90
0
23
66.21
0
23
69.92
0
23
79.01
90
3
24
66.98
3
24
71.62
4
23
80.04
0
26
68.41
0
26
72.63
0
26
80.86
100
4
26
68.54
4
26
73.75
4
26
81.79
0
29
70.06
4
26
74.59
0
29
82.36
110
4
29
70.10
4
29
75.49
4
29
83.85
2
30
71.42
4
29
76.45
4
29
83.60
120
4
32
72.00
9
28
77.06
4
32
84.45
2
33
72.63
4
32
77.95
4
32
84.60
130
4
35
72.53
10
30
78.54
4
35
85.16
4
35
73.74
10
30
79.25
4
35
85.29
140
5
37
73.31
10
33
79.89
5
37
85.66
5
37
74.44
10
33
80.64
5
37
85.80
150
5
40
73.97
11
35
81.04
5
40
86.16
5
40
75.05
10
36
81.79
5
40
86.27
160
6
42
74.39
16
34
82.47
6
42
86.48
6
42
75.45
12
37
82.73
7
41
86.57
170
11
41
74.80
16
37
83.07
11
41
86.97
7
44
75.73
16
37
83.58
11
41
86.81
180
11
44
75.12
17
39
83.88
12
43
86.96
11
44
76.03
16
40
84.34
12
43
87.00
190
14
44
75.32
18
41
84.57
15
43
87.05
14
44
76.21
17
42
84.95
14
44
87.06
200
18
44
75.44
22
41
85.22
17
45
87.07
18
44
76.33
22
41
85.47
17
45
87.08
Computational time (in sec)
P1
P2
P3
P1
P2
P3
P1
P2
P3
P1
P2
P3
P1
P2
P3
P1
P2
P3
367.61
17.97
5.13
332.54
18.56
6.69
352.87
49.09
71.09
345.90
18.11
5.16
325.82
18.72
6.42
338.66
49.18
50.33First, we easily identify the effect of detour flows (bothβ and K) on the optimal coverage for different budgets. For example, the smallest budget of $10 M allows only construction of two dual-access LNG refueling stations with the optimal coverages of 19.99% and 21.14% for R=300 and 600, respectively, when trucks merely use shortest paths between each OD. When we allow trucks to make detours between ODs, the corresponding coverages increase to 27.91% and 29.40% for the both limited driving ranges. This implies, when we do not have enough budget to construct dual-access LNG refueling stations, construction of additional single-access stations can be beneficial to trucks making detours between ODs.Next, we observe that we can spend a large amount of budget effectively when trucks have more route options between ODs. If we have more than $120 M and trucks withR=300 have three SNDF paths (i.e., K=3), more dual-access LNG refueling stations are selected instead of single-access LNG refueling stations comparing to the case of trucks with a single SNDF path (i.e., K=1). It means that since trucks can have many alternative paths to make round trips between ODs, the proposed model with the budget constraint effectively selects more dual-access LNG refueling stations on their routes and improves the optimal coverage for a given budget. The case of trucks with R=600 also shows a similar situation.Next, as shown in Figure9, the marginal optimal coverage resulting from additional refueling stations generally decreases with the amount of construction budget available. In other words, initially, a large percentage of the optimal coverage can be achieved with a relatively small construction budget. After that, the rate of increase of optimal coverage declines as the construction budget increases. One of the interesting things here is that the optimal coverage gap between no-detour and detour flows, however, increases as the construction budget increases. For example, for R=300, when the budget is initially $10 M, the optimal coverage gap between no-detour (β=0%) and detour flows (β=50%) with K=3 is 7.91%; when the budget later becomes $200 M, this gap increases to 11.63%. This result implies that, even if the LNG refueling infrastructure becomes mature (e.g., the budget is around $200 M or more), the effect of detour flows on the optimal coverage would be even more noticeable. The next interesting thing is that the effect of the different number of SNDF paths on the optimal coverage declines as the construction budget increases. That is, when LNG trucks are allowed to deviate from their shortest path up to 50% of dPSvi,vj, the optimal coverage gap between detour flows with K=1 and detour flows with K=3 declines as the construction budget increases. For example, for R=300, at the budget of initial $10 M, this optimal coverage gap is 5.96%; later at the budget of $200 M, the gap declines to 1.85%. It supports our result analysis discussed in Subsection 5.1 on the effect of the number of SNDF paths. Another interesting one here is that the effect of driving range of R=600 on the optimal coverage also declines as the construction budget increases, especially when considering detour options. In other words, comparing Figures 9(a) and 9(b) show that the optimal coverage gap between R=300 and 600 for detour flows with K=3 declines overall as the budget increases. For example, at the budget of initial $10 M, this optimal coverage gap is 1.49%; later at the budget greater than or equal to $120 M, this gap declines to less than or equal to 0.15%. Thus, if the LNG refueling infrastructure becomes mature (e.g., if the budget is around $120 M or more) and if drivers allow reasonable detour options, using LNG trucks with driving range R=600 would not lead to a significant improvement on the LNG refueling coverage.Lastly, the last row of Table8 presents the computational times for P1, P2, and P3 when we run the case study of the state of Pennsylvania. As the computational times of the small instances in Section 5 are shown, the large number of SNDF paths increases the overall computational time. Also, we verify the relationship between the limited driving range and the computational time; that is, the longer the limited driving range, the shorter the computational time because, when the limited driving range is long, the paths that meet the four properties of feasible paths are determined quickly and the number of constraints for Type 2 SNDF paths decreases.


## 7. Conclusions and Future Research

Most path-based demand models available in the literature assume that AF vehicles travel along a unique path between ODs to complete their round trips. In this research work, we have considered that, if preplanned paths for each OD pair cannot offer any AF refueling service, then alternative paths can be used by AF vehicles to receive refueling service. With the four feasible path properties derived in Subsection3.3, we have proposed an algorithm that generates K SNDF paths in each path direction for all OD pairs. A new mixed-integer programming model has been also proposed to optimally locate AF refueling stations at candidate sites with the objective of maximizing the coverage of traffic flows using multiple SNDF paths on a directed transportation network. In order to test our proposed model with the k-th SNDF path algorithm, we have run two experiments: (i) a classical 25-vertex network containing 300 OD pairs with 10 single- and 15 dual-access candidate sites, and (ii) a statewide network of Pennsylvania state containing 2,211 OD pairs with 50 single- and 50 dual-access candidate sites. For our experiment on the classical 25-vertex network, we have analyzed the effects of number of SNDF paths, deviation distance, and vehicle driving range on the coverage of no-detour and detour flows for a given number of AF refueling stations. Computational experiments have shown that the number of SNDF paths and maximum deviation factor largely affect the coverage of traffic flows, especially when the network does not have enough AF refueling stations to cover most of the traffic flows. Next, for our experiment on the case study of Pennsylvania state, we have also analyzed the effects of detour flows, station type (single-access vs. dual-access), and longer driving range for different construction budget levels for the LNG refueling infrastructure to present more insights into our society.For future research, our proposed model can be applied to a variety of real-world problems to locate AF refueling stations in highways. To reduce the complexity of thek-th SNDF path algorithm, we can also attempt to develop a more efficient version of the k-th SNDF path algorithm for large-scale networks. Upchurch et al. [66] consider available fueling capacity of AF refueling stations to cover traffic flows. Similarly, our proposed model can be extended to contemplate the limited capacity of AF refueling stations. Lastly, we generalize our solution approach to address the continuous version of the AF refueling station location problem where stations can be located anywhere along the network. As we have discussed the effect of predetermined candidate sites on the coverage in Subsection 5.2, if we were able to locate stations anywhere along the network, the suboptimality of the optimal solutions of the discrete version of the problem could be significantly improved.



---

*Source: 9473831-2020-02-20.xml*