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

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.

This work is licensed under http://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

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. Table 1 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 residents $\left(C/P\times 100,000\right)$
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.98

Due 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 revised $k$-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. Section 2 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 Subsection 2.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

The $k$-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 of $k$-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 of $k$-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 Table 2.

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 find $k$-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 Section 4.

3.1. Problem Statement

We define the problem on a simple directed network $G\left(V,E\right)$, where $V=\left\{v_i|i=1,\cdots ,n_V\right\}$ is the set of vertices for ODs, $E=\left\{\left(v_i,v_j\right)|\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{o}\mathrm{m}\mathrm{e}{v}_i,v_j\in V\right\}$ is the set of arcs having nonnegative lengths, $\left|V\right|=n_V$, and $\left|E\right|=n_E$. 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=\left\{q\left(v_i,v_j\right)\right|v_i,v_j\in V,i<j\}$ be the set of OD pairs. For any $q\left(v_i,v_j\right)\in Q$, vehicles perform the same round trip, which can be divided into the original path $P\left(v_i,v_j\right)$ from origin $v_i$ to destination $v_j$, and the return path $P\left(v_j,v_i\right)$ from $v_j$ to $v_i$. For convenience, $L_o=\left\{P\left(v_i,v_j\right)|\mathrm{f}\mathrm{o}\mathrm{r}q\left(v_i,v_j\right)\in Q\right\}$ and $L_r=\left\{P\left(v_j,v_i\right)|\mathrm{f}\mathrm{o}\mathrm{r}q\left(v_i,v_j\right)\in Q\right\}$ are defined as the sets of original and return paths, respectively. Let $V\left(P\left(v_i,v_j\right)\right)$ be the set of constituent vertices in path $P\left(v_i,v_j\right)$; then, each path can be represented as a sequence of vertices such that $⟨v_i=v_{\left(0\right)},v_{\left(1\right)},\cdots ,v_j=v_{\left(w\right)}⟩$, where $v_{\left(r\right)}$, $r=0,1,\cdots ,w$, is the $r$-th vertex in $V\left(P\left(v_i,v_j\right)\right)$ and $\left|V\left(P\left(v_i,v_j\right)\right)\right|=w+1$ is the number of vertices in $P\left(v_i,v_j\right)$. Also, the length of $P\left(v_i,v_j\right)$ is calculated as $d\left(P\left(v_i,v_j\right)\right)={\sum }_{r=0}^{w-1}d\left(v_{\left(r\right)},v_{\left(r+1\right)}\right)$, where $d\left(v_{\left(r\right)},v_{\left(r+1\right)}\right)$ is the arc length from the $r$-th vertex to the next vertex in $P\left(v_i,v_j\right)$. Similarly, $d\left(P\left(v_j,v_i\right)\right)$ is defined as the length of the return path. If $d\left(P\left(v_i,v_j\right)\right)$ and $d\left(P\left(v_j,v_i\right)\right)$ are the shortest distances from $v_i$ to $v_j$ and from $v_j$ to $v_i$, respectively, then the corresponding original and return paths are denoted as $P_S\left(v_i,v_j\right)$ and $P_S\left(v_j,v_i\right)$. Also, we define $F_k\left(v_i,v_j\right)$ and $F_k\left(v_j,v_i\right)$ as the $k$-th shortest feasible original and return paths for OD pair $q\left(v_i,v_j\right)$. 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 as $M=\left\{m_t|t=1,\cdots ,n_M\right\}$, 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 $CS\left(F_k\left(v_i,v_j\right)\right)=⟨m_{\left(0\right)},\cdots ,m_{\left(f\right)}⟩$ as the sequence of candidate sites in $F_k\left(v_i,v_j\right)$, where $m_{\left(0\right)},\cdots ,m_{\left(f\right)}\in M$, such that $0\le f\le n_M-1$. In order to represent the distances between $v_i=v_{\left(0\right)}$, $v_j=v_{\left(w\right)}$, and site $m_{\left(l\right)}\in CS\left(F_k\left(v_i,v_j\right)\right)$, such that $0\le l\le f$, $d_{ijk}\left(v_{\left(0\right)},m_{\left(l\right)}\right)$ denotes the distance from $v_{\left(0\right)}$ to $m_{\left(l\right)}$ in $F_k\left(v_i,v_j\right)$. If $\left(v_{\left(s\right)},v_{\left(s+1\right)}\right)$ is the arc that passes through $m_{\left(l\right)}$ in $CS\left(F_k\left(v_i,v_j\right)\right)$, such that $0\le s<w$, then $d_{ijk}\left(v_{\left(0\right)},m_{\left(l\right)}\right)$ is calculated as ${\sum }_{r=0}^{s-1}d\left(v_{\left(r\right)},v_{\left(r+1\right)}\right)+d\left(v_{\left(s\right)},m_{\left(l\right)}\right)$, where $v_{\left(s\right)}$ is the vertex adjacent to $m_{\left(l\right)}$ and $d\left(v_{\left(s\right)},m_{\left(l\right)}\right)$ is the distance from $v_{\left(s\right)}$ to $m_{\left(l\right)}$. Similarly, $d_{ijk}\left(m_{\left(l\right)},v_{\left(w\right)}\right)$, which refers to the distance from $m_{\left(l\right)}$ to $v_{\left(w\right)}$, is calculated as $d\left(m_{\left(l\right)},v_{\left(s+1\right)}\right)+{\sum }_{r=s+1}^{w-1}d\left(v_{\left(r\right)},v_{\left(r+1\right)}\right)$, where $d\left(m_{\left(l\right)},v_{\left(s+1\right)}\right)$ indicates the distance from $m_{\left(l\right)}$ to $v_{\left(s+1\right)}$. Furthermore, the distance between two candidate sites $m_{\left(l\right)}$ and $m_{\left(l^{{}^{\prime }}\right)}$, where $l\ne l^{{}^{\prime }}$, is denoted as $d_{ijk}\left(m_{\left(l\right)},m_{\left(l^{{}^{\prime }}\right)}\right)$. When we select the two arcs, $\left(v_{\left(s\right)},v_{\left(s+1\right)}\right)$ and $\left(v_{\left(s^{\prime }\right)},v_{\left(s^{\prime }+1\right)}\right)$, which include $m_{\left(l\right)}$ and $m_{\left(l^{{}^{\prime }}\right)}$, respectively, such that $0\le s\le s^{\prime }<w$, $d_{ijk}\left(m_{\left(l\right)},m_{\left(l^{\prime }\right)}\right)={\sum }_{r=s}^{s^{\prime }}d\left(v_{\left(r\right)},v_{\left(r+1\right)}\right)-d\left(v_{\left(s\right)},m_{\left(l\right)}\right)-d\left(m_{\left(l^{\prime }\right)},v_{\left(s^{\prime }+1\right)}\right)$. Figure 1 shows a representation of vertex and candidate site sequences, as well as distances between vertices and candidate sites.

[figure omitted; refer to PDF]

We consider that vehicles make a complete round trip between their ODs. They have a limited driving range under free flow conditions, denoted as $R$, 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 $SS\left(F_k\left(v_i,v_j\right)\right)$, $ES\left(F_k\left(v_i,v_j\right)\right)$, and $IS\left(F_k\left(v_i,v_j\right)\right)$ as the following three sets of candidate sites in $CS\left(F_k\left(v_i,v_j\right)\right)$, which are categorized depending on their distances from $v_i$ and $v_j$:$$\begin{array}{}\text{(1)}& SS\left(F_k\left(v_i,v_j\right)\right)=\left\{m_{\left(l\right)}\in CS\left(F_k\left(v_i,v_j\right)\right)\left|d_{ijk}\left(v_{\left(0\right)},m_{\left(l\right)}\right)\right.\le \frac{R}{2}\right\},\end{array}$$$$\begin{array}{}\text{(2)}& ES\left(F_k\left(v_i,v_j\right)\right)=\left\{m_{\left(l\right)}\in CS\left(F_k\left(v_i,v_j\right)\right)\left|d_{ijk}\left({m_{\left(l\right)},v}_{\left(w\right)}\right)\right.\le \frac{R}{2}\right\},\end{array}$$$$\begin{array}{}\text{(3)}& IS\left(F_k\left(v_i,v_j\right)\right)=\left\{m_{\left(l\right)}\in CS\left(F_k\left(v_i,v_j\right)\right)\left|d_{ijk}\left(v_{\left(0\right)},m_{\left(l\right)}\right)\right.>\frac{R}{2}\mathrm{a}\mathrm{n}\mathrm{d}{d}_{ijk}\left({m_{\left(l\right)},v}_{\left(w\right)}\right)>\frac{R}{2}\right\}.\end{array}$$

Note that we consider $IS\left(F_k\left(v_i,v_j\right)\right)$ only for $F_k\left(v_i,v_j\right)$ such that $d\left(F_k\left(v_i,v_j\right)\right)>R$. For return paths, we similarly define $CS\left(F_k\left(v_j,v_i\right)\right)$, $SS\left(F_k\left(v_j,v_i\right)\right)$, $ES\left(F_k\left(v_j,v_i\right)\right)$, and $IS\left(F_k\left(v_j,v_i\right)\right)$. When a feasible set of AF refueling stations is located in $SS\left(F_k\left(v_i,v_j\right)\right)$, $ES\left(F_k\left(v_i,v_j\right)\right)$, $IS\left(F_k\left(v_i,v_j\right)\right)$, $SS\left(F_k\left(v_j,v_i\right)\right)$, $ES\left(F_k\left(v_j,v_i\right)\right)$, and $IS\left(F_k\left(v_j,v_i\right)\right)$, 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 $\beta$. 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 $\beta$ 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 $\beta$ 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\left(v_i,v_j\right),P\left(v_j,v_i\right)$	Original and return paths between interchanges $v_i$ and $v_j$
$L_o,L_r$	Sets of original and return paths
$V\left(P\left(v_i,v_j\right)\right)$	Set of constituent vertices in path $P\left(v_i,v_j\right)$
$d\left(P\left(v_i,v_j\right)\right)$	Length of $P\left(v_i,v_j\right)$
$P_S\left(v_i,v_j\right)$	Shortest original path from interchange $v_i$ to interchange $v_j$
$F_k\left(v_i,v_j\right)$	$k$ -th shortest feasible original path from interchange $v_i$ to interchange $v_j$
$M$	Set of all candidate sites for AF refueling stations
$CS\left(F_k\left(v_i,v_j\right)\right)$	Sequence of candidate sites in $F_k\left(v_i,v_j\right)$
$d_{ijk}\left(v_{\left(r\right)},m_{\left(l\right)}\right)$	Distance from the $r$-th interchange $v_{\left(r\right)}$ to the $l$-th candidate site $m_{\left(l\right)}$ in $F_k\left(v_i,v_j\right)$
$d_{ijk}\left(m_{\left(l\right)},m_{\left(l^{{}^{\prime }}\right)}\right)$	Distance between two candidate sites $m_{\left(l\right)}$ and $m_{\left(l^{{}^{\prime }}\right)}$ in $F_k\left(v_i,v_j\right)$
$SS\left(F_k\left(v_i,v_j\right)\right)$	Set of candidate sites that are located within a distance of $R/2$ from interchange $v_i$ in $F_k\left(v_i,v_j\right)$
$ES\left(F_k\left(v_i,v_j\right)\right)$	Set of candidate sites that are located within a distance of $R/2$ from interchange $v_j$ in $F_k\left(v_i,v_j\right)$
$IS\left(F_k\left(v_i,v_j\right)\right)$	Set of candidate sites that are located beyond a distance of $R/2$ from both interchanges $v_i$ and $v_j$ in $F_k\left(v_i,v_j\right)$
$R$	Vehicle’s maximum travel distance with a single refueling
$\beta$	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 vertices $\left\{v_1,\cdots ,v_{17}\right\}$, two dual-access candidate sites $\left\{m_1,m_3\right\}$, and five single-access candidate sites $\left\{m_2,m_4,m_5,m_6,m_7\right\}$ 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 $q\left(v_2,v_4\right)$, shortest paths $P_S\left(v_2,v_4\right)$ and $P_S\left(v_4,v_2\right)$ are feasible because the corresponding trips can be covered by placing a single refueling station at site $m_1$. However, for $q\left(v_1,v_3\right)$, since candidate sites are available in neither $P_S\left(v_1,v_3\right)$ nor $P_S\left(v_3,v_1\right)$, drivers need to deviate from these paths to receive refueling service in their round trips. If a refueling station is located at site $m_1$, then vertex sequence $⟨v_1=v_{\left(0\right)},v_2=v_{\left(1\right)},v_4=v_{\left(2\right)},v_2=v_{\left(3\right)},v_3=v_{\left(4\right)}⟩$ for the original path and opposite sequence $⟨v_3=v_{\left(0\right)},v_2=v_{\left(1\right)},v_4=v_{\left(2\right)},v_2=v_{\left(3\right)},v_1=v_{\left(4\right)}⟩$ for the return path are feasible paths for $q\left(v_1,v_3\right)$. Second, suppose that there are two AF refueling stations available at sites $m_3$ and $m_4$. Then, the shortest original path of $q\left(v_6,v_{11}\right)$ does not offer AF refueling availability because vehicles cannot reach $v_{11}$ 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 $v_{11}$ cannot reach the AF refueling station at site $m_3$. Thus, vehicles in $q\left(v_6,v_{11}\right)$ would need to make detours using vertex sequence $⟨v_6=v_{\left(0\right)},v_7=v_{\left(1\right)},v_8=v_{\left(2\right)},v_9=v_{\left(3\right)},v_{10}=v_{\left(4\right)},v_7=v_{\left(5\right)},v_8=v_{\left(6\right)},v_{11}=v_{\left(7\right)}⟩$ for the original path and opposite sequence $⟨v_{11}=v_{\left(0\right)},v_8=v_{\left(1\right)},v_9=v_{\left(2\right)},v_{10}=v_{\left(3\right)},v_7=v_{\left(4\right)},v_6=v_{\left(5\right)}⟩$ for the return path. These two deviation paths are feasible paths for $q\left(v_6,v_{11}\right)$. Lastly, suppose that sites $m_5$ and $m_6$ are selected as AF refueling stations. Vehicles in $q\left(v_{12},v_{17}\right)$ 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 $v_{13}$, then vertex sequence $⟨v_{12}=v_{\left(0\right)},v_{13}=v_{\left(1\right)},v_{14}=v_{\left(2\right)},v_{13}=v_{\left(3\right)},v_{15}=v_{\left(4\right)},v_{13}=v_{\left(5\right)},v_{17}=v_{\left(6\right)}⟩$ becomes a feasible original path and vertex sequence $⟨v_{17}=v_{\left(0\right)},v_{15}=v_{\left(1\right)},v_{13}=v_{\left(2\right)},v_{12}=v_{\left(3\right)}⟩$ turn out to be a feasible return path for $q\left(v_{12},v_{17}\right)$. Note that another vertex sequence $⟨v_{17}=v_{\left(0\right)},v_{13}=v_{\left(1\right)},v_{14}=v_{\left(2\right)},v_{13}=v_{\left(3\right)},v_{15}=v_{\left(4\right)},v_{13}=v_{\left(5\right)},v_{12}=v_{\left(6\right)}⟩$ 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 $⟨v_{17}=v_{\left(0\right)},v_{15}=v_{\left(1\right)},v_{13}=v_{\left(2\right)},v_{12}=v_{\left(3\right)}⟩$. From these observations, we can state four feasible path properties in the next subsection.

[figure omitted; refer to PDF]

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 existing $k$-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.

Let $F_k\left(v_i,v_j\right)$ be the $k$-th shortest feasible original path from $v_i$ to $v_j$ that contains feasible candidate sites for AF refueling stations to be able to cover trips in $P\left(v_i,v_j\right)$, and $CS\left(F_k\left(v_i,v_j\right)\right)=⟨m_{\left(0\right)},\cdots ,m_{\left(f\right)}⟩$ denote the sequence of candidate sites in $F_k\left(v_i,v_j\right)$, where $m_{\left(0\right)},\cdots ,m_{\left(f\right)}\in M$. According to the cardinality of $CS\left(F_k\left(v_i,v_j\right)\right)$, denoted as $\left|CS\left(F_k\left(v_i,v_j\right)\right)\right|$, one of the following two cases must be satisfied:

 

(a) If $\left|CS\left(F_k\left(v_i,v_j\right)\right)\right|=1$, then $SS\left(F_k\left(v_i,v_j\right)\right)\cap ES\left(F_k\left(v_i,v_j\right)\right)\ne \varnothing$.

 

(b) If $\left|CS\left(F_k\left(v_i,v_j\right)\right)\right|>1$, then $SS\left(F_k\left(v_i,v_j\right)\right)\ne \varnothing$, $ES\left(F_k\left(v_i,v_j\right)\right)\ne \varnothing$, and $d_{ijk}\left(m_{\left(l\right)},m_{\left(l+1\right)}\right)\le R$ for any two adjacent candidate sites $m_{\left(l\right)},m_{\left(l+1\right)}\in CS\left(F_k\left(v_i,v_j\right)\right)$.

Proof.

Since $F_k\left(v_i,v_j\right)$ is a feasible path, the half-full tank assumption implies that $\left|CS\left(F_k\left(v_i,v_j\right)\right)\right|\ge 1$. In case (a), there is a single candidate site available for a refueling station on $F_k\left(v_i,v_j\right)$, and by the half-full tank assumption, this site is located within a distance of $R/2$ from $v_i$ and $v_j$ to ensure that vehicles on $F_k\left(v_i,v_j\right)$ make a successful trip from $v_i$ to $v_j$. Then, by definition of $SS\left(F_k\left(v_i,v_j\right)\right)$ and $ES\left(F_k\left(v_i,v_j\right)\right)$, this site should be located in the common segment in $SS\left(F_k\left(v_i,v_j\right)\right)\cap ES\left(F_k\left(v_i,v_j\right)\right)$. This implies that $SS\left(F_k\left(v_i,v_j\right)\right)\cap ES\left(F_k\left(v_i,v_j\right)\right)\ne \varnothing$.

In case (b), there are at least two candidate sites available on $F_k\left(v_i,v_j\right)$. If either $SS\left(F_k\left(v_i,v_j\right)\right)=\varnothing$ or $ES\left(F_k\left(v_i,v_j\right)\right)=\varnothing$, then $F_k\left(v_i,v_j\right)$ cannot offer any AF refueling service. Thus, $SS\left(F_k\left(v_i,v_j\right)\right)\ne \varnothing$ and $ES\left(F_k\left(v_i,v_j\right)\right)\ne \varnothing$. For any two adjacent candidate sites $m_{\left(l\right)},m_{\left(l+1\right)}\in CS\left(F_k\left(v_i,v_j\right)\right)$, if $d_{ijk}\left(m_{\left(l\right)},m_{\left(l+1\right)}\right)>R$, then vehicles at site $m_{\left(l\right)}$ cannot reach site $m_{\left(l+1\right)}$, i.e., $F_k\left(v_i,v_j\right)$ does not have AF refueling availability. It also contradicts the definition of $F_k\left(v_i,v_j\right)$. Thus, $d_{ijk}\left(m_{\left(l\right)},m_{\left(l+1\right)}\right)\le R$ for any two adjacent candidate sites $m_{\left(l\right)},m_{\left(l+1\right)}\in CS\left(F_k\left(v_i,v_j\right)\right)$.

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 of $q\left(v_{14},v_{16}\right)$ in Figure 2. While a vertex sequence $⟨v_{16}=v_{\left(0\right)},v_{15}=v_{\left(1\right)},v_{13}=v_{\left(2\right)},v_{17}=v_{\left(3\right)},{v_{15}=v_{\left(4\right)},v_{13}=v_{\left(5\right)},v}_{14}=v_{\left(6\right)}⟩$ goes over site $m_6$ twice, vehicles can also use a vertex sequence $⟨v_{16}=v_{\left(0\right)},v_{15}=v_{\left(1\right)},v_{13}=v_{\left(2\right)},v_{14}=v_{\left(3\right)}⟩$ to travel from $v_{16}$ to $v_{14}$ with a single visit to site $m_6$. 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 $q\left(v_1,v_3\right)$ in Figure 2, a feasible vertex sequence $⟨v_1=v_{\left(0\right)},v_2=v_{\left(1\right)},v_4=v_{\left(2\right)},v_2=v_{\left(3\right)},v_3=v_{\left(4\right)}⟩$ passes through a refueling station at dual-access site $m_1$ 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.

For $q\left(v_i,v_j\right)\in Q$, suppose that $F_{k_1}\left(v_i,v_j\right)$ is an original path that satisfies Property 1 and passes through a single-access candidate site $m_{t_1}\in M$ twice:

$$\begin{array}{}\text{(4)}& CS\left(F_{k_1}\left(v_i,v_j\right)\right)=⟨m_{\left(0\right)},\cdots m_{\left(a\right)},m_{\left(a+1\right)}\cdots ,m_{\left(b\right)},m_{\left(b+1\right)},\cdots ,m_{\left(f\right)}⟩,\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{m}_{t_1}=m_{\left(a\right)}=m_{\left(b\right)}.\end{array}$$

Then, there must exist an original path, denoted as $F_{k_2}\left(v_i,v_j\right)$, that satisfies Property 1 and goes over site $m_{t_1}$ only one time, such that $CS\left(F_{k_2}\left(v_i,v_j\right)\right)<CS\left(F_{k_1}\left(v_i,v_j\right)\right)$:$$\begin{array}{}\text{(5)}& CS\left(F_{k_2}\left(v_i,v_j\right)\right)=⟨m_{\left(0\right)},\cdots m_{\left(a\right)},m_{\left(b+1\right)},\cdots ,m_{\left(f\right)}⟩.\end{array}$$

Thus, $d\left(F_{k_2}\left(v_i,v_j\right)\right)<d\left(F_{k_1}\left(v_i,v_j\right)\right)$.

Proof.

Let $m_{\left(a\right)}$ and $m_{\left(b\right)}$ be the $a$-th and $b$-th candidate sites in $CS\left(F_{k_1}\left(v_i,v_j\right)\right)=⟨m_{\left(0\right)},\cdots m_{\left(a\right)},m_{\left(a+1\right)}\cdots ,m_{\left(b\right)},m_{\left(b+1\right)},\cdots ,m_{\left(f\right)}⟩$ which are located at the same site $m_{t_1}$ such that $0\le a<b\le f$, where $f=\left|CS\left(F_{k_1}\left(v_i,v_j\right)\right)\right|-1$. By eliminating all consecutive candidate sites from $m_{\left(a+1\right)}$ to $m_{\left(b\right)}$ and the corresponding arcs in $F_{k_1}\left(v_i,v_j\right)$, we can construct path $F_{k_2}\left(v_i,v_j\right)$ with $CS\left(F_{k_2}\left(v_i,v_j\right)\right)=⟨m_{\left(0\right)},\cdots m_{\left(a\right)},m_{\left(b+1\right)},\cdots ,m_{\left(f\right)}⟩$. Then, from Property 1, since $d_{ij{k}_1}\left(m_{\left(b\right)},m_{\left(b+1\right)}\right)\le R$, vehicles at $m_{\left(b\right)}$ can reach $m_{\left(b+1\right)}$ in $F_{k_1}\left(v_i,v_j\right)$. This implies that vehicles at $m_{\left(a\right)}$ can reach $m_{\left(b+1\right)}$ in $F_{k_2}\left(v_i,v_j\right)$ because $m_{\left(a\right)}$ and $m_{\left(b\right)}$ are placed at the same site $m_{t_1}$. Thus, $CS\left(F_{k_2}\left(v_i,v_j\right)\right)$ is a feasible sequence of candidate sites for $F_{k_2}\left(v_i,v_j\right)$ and $d\left(F_{k_2}\left(v_i,v_j\right)\right)<d\left(F_{k_1}\left(v_i,v_j\right)\right)$.

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 of $q\left(v_1,v_3\right)$ from the previous network in Figure 2. Recall that a vertex sequence $⟨v_1=v_{\left(0\right)},v_2=v_{\left(1\right)},v_4=v_{\left(2\right)},v_2=v_{\left(3\right)},v_3=v_{\left(4\right)}⟩$ is one of the feasible original paths for $q\left(v_1,v_3\right)$. Now, let us consider a different vertex sequence defined by $⟨v_1=v_{\left(0\right)},v_2=v_{\left(1\right)},v_4=v_{\left(2\right)},v_5=v_{\left(3\right)},v_4=v_{\left(4\right)},v_2=v_{\left(5\right)},v_3=v_{\left(6\right)}⟩$. The candidate site sequence for this second path is $⟨m_1=m_{\left(0\right)},m_2=m_{\left(1\right)},m_1=m_{\left(2\right)}⟩$, which enables successful trips from $v_1$ to $v_3$ without running out of fuel. However, drivers would not travel further along the subpath of the second path, $⟨v_4=v_{\left(2\right)},v_5=v_{\left(3\right)},v_4=v_{\left(4\right)}⟩$ because they can reach $v_3$ without refueling at site $m_2$. They would return at $v_4$ after visiting the refueling station at site $m_1$, i.e., they would use the first feasible path, $⟨v_1=v_{\left(0\right)},v_2=v_{\left(1\right)},v_4=v_{\left(2\right)},v_2=v_{\left(3\right)},v_3=v_{\left(4\right)}⟩$, instead of the second one for the original path of $q\left(v_1,v_3\right)$. Similarly, drivers can consider two possible paths for the return path of $q\left(v_{12},v_{17}\right)$. While vehicles using a vertex sequence $⟨v_{17}=v_{\left(0\right)},v_{13}=v_{\left(1\right)},v_{14}=v_{\left(2\right)},v_{13}=v_{\left(3\right)},v_{15}=v_{\left(4\right)},v_{13}=v_{\left(5\right)},v_{12}=v_{\left(6\right)}⟩$ require each refueling stop at both sites $m_5$ and $m_6$, vehicles in a vertex sequence $⟨v_{17}=v_{\left(0\right)},v_{15}=v_{\left(1\right)},v_{13}=v_{\left(2\right)},v_{12}=v_{\left(3\right)}⟩$ need only one refueling stop at site $m_6$. Since the length of the latter path is shorter, vehicles in the return path of $q\left(v_{12},v_{17}\right)$ 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, $F_{k_1}\left(v_i,v_j\right)$ and $F_{k_2}\left(v_i,v_j\right)$, such that $d\left(F_{k_1}\left(v_i,v_j\right)\right)<d\left(F_{k_2}\left(v_i,v_j\right)\right)$, which satisfy Properties 1 and 2, $F_{k_2}\left(v_i,v_j\right)$ includes irrelevant candidate sites in set $CS\left(F_{k_2}\left(v_i,v_j\right)\right)\backslash CS\left(F_{k_1}\left(v_i,v_j\right)\right)$ when one of the following conditions is satisfied:

(a)

If $\left|CS\left(F_{k_1}\left(v_i,v_j\right)\right)\right|=1$, then (1) $CS\left(F_{k_1}\left(v_i,v_j\right)\right)\subset CS\left(F_{k_2}\left(v_i,v_j\right)\right)$, and (2) either $S\left(F_{k_1}\left(v_i,v_j\right)\right)=SS\left(F_{k_2}\left(v_i,v_j\right)\right)$ or $ES\left(F_{k_1}\left(v_i,v_j\right)\right)=ES\left(F_{k_2}\left(v_i,v_j\right)\right)$.

(b)

If $\left|CS\left(F_{k_1}\left(v_i,v_j\right)\right)\right|>1$, then (1) $CS\left(F_{k_1}\left(v_i,v_j\right)\right)\subset CS\left(F_{k_2}\left(v_i,v_j\right)\right)$, (2) $SS\left(F_{k_1}\left(v_i,v_j\right)\right)=SS\left(F_{k_2}\left(v_i,v_j\right)\right)$, (3) $ES\left(F_{k_1}\left(v_i,v_j\right)\right)=ES\left(F_{k_2}\left(v_i,v_j\right)\right)$, and (4) $0<d_{ij{k}_1}\left(m_{\left(0\right)},m_{\left(f\right)}\right)\le R$, where $CS\left(F_{k_1}\left(v_i,v_j\right)\right)=⟨m_{\left(0\right)},\cdots ,m_{\left(f\right)}⟩$.

Proof.

In case (a), suppose that $m_{\left(0\right)}\in SS\left(F_{k_1}\left(v_i,v_j\right)\right)=SS\left(F_{k_2}\left(v_i,v_j\right)\right)$. By Property 1, $m_{\left(0\right)}\in ES\left(F_{k_1}\left(v_i,v_j\right)\right)$, which implies that $d_{ij{k}_1}\left(m_{\left(0\right)},v_j\right)\le R/2$. Thus, any candidate site in $CS\left(F_{k_2}\left(v_i,v_j\right)\right)\backslash CS\left(F_{k_1}\left(v_i,v_j\right)\right)$ is irrelevant because vehicles covered by a single refueling station at site $m_{\left(0\right)}$ can complete trips from $v_i$ to $v_j$. Similarly, it can be shown when $m_{\left(0\right)}\in ES\left(F_{k_1}\left(v_i,v_j\right)\right)=ES\left(F_{k_2}\left(v_i,v_j\right)\right)$.

In case (b), by Condition (4), $d_{ij{k}_1}\left(m_{\left(a\right)},m_{\left(b\right)}\right)\le R$ such that $m_{\left(a\right)}\in SS\left(F_{k_1}\left(v_i,v_j\right)\right)$, $m_{\left(b\right)}\in ES\left(F_{k_1}\left(v_i,v_j\right)\right)$, and $0\le a<b\le f$. This implies that any additional refueling stop is unnecessary between $m_{\left(a\right)}$ and $m_{\left(b\right)}$. Then, by Conditions (1), (2), and (3), $IS\left(F_{k_2}\left(v_i,v_j\right)\right)\ne \varnothing$, which implies that $F_{k_2}\left(v_i,v_j\right)$ includes irrelevant candidate sites.

Given two feasible paths, $F_{k_1}\left(v_i,v_j\right)$ and $F_{k_2}\left(v_i,v_j\right)$ such that $d\left(F_{k_1}\left(v_i,v_j\right)\right)\le d\left(F_{k_2}\left(v_i,v_j\right)\right)$, for the same OD pair $q\left(v_i,v_j\right)\in Q$, if any subset of candidate sites that covers the trips in $F_{k_2}\left(v_i,v_j\right)$ is also able to cover the trips in $F_{k_1}\left(v_i,v_j\right)$, then we say that $F_{k_2}\left(v_i,v_j\right)$ is dominated by $F_{k_1}\left(v_i,v_j\right)$. This implies that $F_{k_2}\left(v_i,v_j\right)$ is unnecessary to be considered in our proposed model. In practice, drivers would use only $F_{k_1}\left(v_i,v_j\right)$ because $d\left(F_{k_1}\left(v_i,v_j\right)\right)\le d\left(F_{k_2}\left(v_i,v_j\right)\right)$. For example, let us take the previous example of $q\left(v_6,v_{11}\right)$ in Figure 2. Recall that $F_1\left(v_6,v_{11}\right)=⟨v_6=v_{\left(0\right)},v_7=v_{\left(1\right)},v_8=v_{\left(2\right)},v_9=v_{\left(3\right)},v_{10}=v_{\left(4\right)},v_7=v_{\left(5\right)},v_8=v_{\left(6\right)},v_{11}=v_{\left(7\right)}⟩$. Then, we can find another feasible original path, which is $F_2\left(v_6,v_{11}\right)=⟨v_6=v_{\left(0\right)},v_7=v_{\left(1\right)},v_{10}=v_{\left(2\right)},v_9=v_{\left(3\right)},v_{10}=v_{\left(4\right)},v_7=v_{\left(5\right)},v_8=v_{\left(6\right)},v_{11}=v_{\left(7\right)}⟩$ for vehicles with $R=100$. They have the same sequence of candidate sites, but vehicles may prefer $F_1\left(v_6,v_{11}\right)$ to $F_2\left(v_6,v_{11}\right)$ because $d\left(F_1\left(v_6,v_{11}\right)\right)=125<d\left(F_2\left(v_6,v_{11}\right)\right)=145$. Furthermore, in order to clarify the concept of dominated feasible paths, let us see another example of $q\left(v_{14},v_{16}\right)$ in Figure 2. For the return path, we have two feasible return paths, $F_1\left(v_{16},v_{14}\right)=⟨v_{16}=v_{\left(0\right)},v_{15}=v_{\left(1\right)},v_{13}=v_{\left(2\right)},v_{14}=v_{\left(3\right)}⟩$ and $F_2\left(v_{16},v_{14}\right)=⟨v_{16}=v_{\left(0\right)},v_{15}=v_{\left(1\right)},v_{17}=v_{\left(2\right)},v_{13}=v_{\left(3\right)},v_{14}=v_{\left(4\right)}⟩$, such that $CS\left(F_1\left(v_{16},v_{14}\right)\right)=⟨m_7=m_{\left(0\right)},m_6=m_{\left(1\right)},m_5=m_{\left(2\right)}⟩$ and $CS\left(F_2\left(v_{16},v_{14}\right)\right)=⟨m_7=m_{\left(0\right)},m_5=m_{\left(1\right)}⟩$. For $F_1\left(v_{16},v_{14}\right)$, we have five possible sets of candidate sites to cover the trip, i.e., $\left\{m_6\right\}$, $\left\{m_6,m_5\right\}$, $\left\{m_7,m_6\right\}$, $\left\{m_7,m_5\right\}$, and $\left\{m_7,m_6,m_5\right\}$. In case of $F_2\left(v_{16},v_{14}\right)$, $\left\{m_7,m_5\right\}$ is the only feasible set. That is, $CS\left(F_2\left(v_{16},v_{14}\right)\right)\subseteq CS\left(F_1\left(v_{16},v_{14}\right)\right)$. It means that $F_2\left(v_{16},v_{14}\right)$ is dominated by $F_1\left(v_{16},v_{14}\right)$ because $d\left(F_1\left(v_{16},v_{14}\right)\right)=65<d\left(F_2\left(v_{16},v_{14}\right)\right)=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.

Given $q\left(v_i,v_j\right)\in Q$, let $F_{k_1}\left(v_i,v_j\right)$ be the shortest feasible original path with $CS\left(F_{k_1}\left(v_i,v_j\right)\right)$. Suppose that there exists another feasible original path, denoted as $F_{k_2}\left(v_i,v_j\right)$, such that $CS\left(F_{k_2}\left(v_i,v_j\right)\right)\subseteq CS\left(F_{k_1}\left(v_i,v_j\right)\right)$ and $d\left(F_{k_1}\left(v_i,v_j\right)\right)\le d\left(F_{k_2}\left(v_i,v_j\right)\right)$. Then, $F_{k_2}\left(v_i,v_j\right)$ is dominated by $F_{k_1}\left(v_i,v_j\right)$.

Proof.

By definition of dominated paths, we need to show that any subset of candidate sites to cover the trips in $F_{k_2}\left(v_i,v_j\right)$ is also able to cover the trips in $F_{k_1}\left(v_i,v_j\right)$. Suppose that $FS\left(F_{k_2}\left(v_i,v_j\right)\right)$ is a set of candidate sites that covers trips in $F_{k_2}\left(v_i,v_j\right)$. This implies that $FS\left(F_{k_2}\left(v_i,v_j\right)\right)\subseteq CS\left(F_{k_2}\left(v_i,v_j\right)\right)$, and therefore, $FS\left(F_{k_2}\left(v_i,v_j\right)\right)\subseteq CS\left(F_{k_1}\left(v_i,v_j\right)\right)$ because $CS\left(F_{k_2}\left(v_i,v_j\right)\right)\subseteq CS\left(F_{k_1}\left(v_i,v_j\right)\right)$. If $\left|FS\left(F_{k_2}\left(v_i,v_j\right)\right)\right|=1$, then it is trivial to show that $FS\left(F_{k_2}\left(v_i,v_j\right)\right)$ can cover trips in $F_{k_1}\left(v_i,v_j\right)$ because $d\left(F_{k_1}\left(v_i,v_j\right)\right)\le d\left(F_{k_2}\left(v_i,v_j\right)\right)$. If $\left|FS\left(F_{k_2}\left(v_i,v_j\right)\right)\right|>1$, then we select the first and last sequences of candidate sites in $FS\left(F_{k_2}\left(v_i,v_j\right)\right)$, denoted as $m_{\left(a\right)}$ and $m_{\left(b\right)}$, respectively. When we decompose $F_{k_1}\left(v_i,v_j\right)$ into $⟨v_i=v_{\left(0\right)},\cdots ,m_{\left(a\right)}⟩$, $⟨m_{\left(a\right)},\cdots ,m_{\left(b\right)}⟩$, and $⟨m_{\left(b\right)},\cdots ,v_j=v_{\left(w\right)}⟩$, we have that $d\left(F_{k_1}\left(v_i,v_j\right)\right)=d_{ij{k}_1}\left(v_{\left(0\right)},m_{\left(a\right)}\right)+d_{ij{k}_1}\left(m_{\left(a\right)},m_{\left(b\right)}\right)+d_{ij{k}_1}\left(m_{\left(b\right)},v_{\left(w\right)}\right)$. If $d_{ij{k}_2}\left(v_{\left(0\right)},m_{\left(a\right)}\right)\le R/2<d_{ij{k}_1}\left(v_{\left(0\right)},m_{\left(a\right)}\right)$, then we can construct a shorter feasible path whose length is $d_{ij{k}_2}\left(v_{\left(0\right)},m_{\left(a\right)}\right)+d_{ij{k}_1}\left(m_{\left(a\right)},m_{\left(b\right)}\right)+d_{ij{k}_1}\left(m_{\left(b\right)},v_{\left(w\right)}\right)$. It contradicts that $F_{k_1}\left(v_i,v_j\right)$ is the shortest feasible path. Thus, $d_{ij{k}_1}\left(v_{\left(0\right)},m_{\left(a\right)}\right)\le R/2$. Similarly, it can be shown that $d_{ij{k}_1}\left(m_{\left(b\right)},v_{\left(w\right)}\right)\le R/2$. Furthermore, if $d_{ij{k}_2}\left(m_{\left(l\right)},m_{\left(l+1\right)}\right)\le R<d_{ij{k}_1}\left(m_{\left(l\right)},m_{\left(l+1\right)}\right)$ such that $a\le l<b$, then it also contradicts that $F_{k_1}\left(v_i,v_j\right)$ is the shortest feasible path. Thus, $FS\left(F_{k_2}\left(v_i,v_j\right)\right)$ can cover trips in $F_{k_1}\left(v_i,v_j\right)$ by Property 1.

3.4. $\mathit{k}$-th Shortest Nondominated Feasible Path Algorithm

In this subsection, we describe the $k$-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 $\beta$ and a maximum number $K$ of SNDF paths, the $k$-th SNDF path algorithm determines up to $K$ SNDF paths from $v_i$ to $v_j$ and another up to $K$ SNDF paths from $v_j$ to $v_i$.

The first step of the $k$-th SNDF path algorithm is to transform a simple directed network $G\left(V,E\right)$ with the set of $n_V$ vertices and the set of $n_E$ arcs into an expanded directed network, denoted as $G^{\prime }\left(V^{\prime },E^{\prime }\right)$, 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^{\prime }$. Then, the maximum cardinalities of $V^{\prime }$ and $E^{\prime }$ can be written as $\left|V^{\prime }\right|=n_V+n_{SM}+{2n}_{DM}$ and $\left|E^{\prime }\right|=n_E+2n_{SM}+{4n}_{DM}$, respectively, where $n_{SM}$ and $n_{DM}$ refer to the number of single- and dual-access candidate sites. Next, for interchange $v_j\in V^{\prime }$, $T_{v_j}$ denotes a shortest path arborescence rooted at vertex $v_j$ containing the shortest path to each $v_i$ in $V^{\prime }\backslash \left\{v_j\right\}$. The shortest path arborescence $T_{v_j}$ for interchange $v_j$ in $G^{\prime }\left(V^{\prime },E^{\prime }\right)$ can be constructed using Dijkstra’s algorithm in time $O\left(\left|E^{\prime }\right|+\left|V^{\prime }\right|\mathrm{l}\mathrm{o}\mathrm{g}\left|V^{\prime }\right|\right)$. Let $E\left(T_{v_j}\right)$ denote the set of arcs composing $T_{v_j}$. If we select arc $\left(v_a,v_b\right)\in E^{\prime }\backslash E\left(T_{v_j}\right)$ and follow the shortest path from $v_b$ to $v_j$ instead of taking the shortest path in $T_{v_j}$, then the additional distance from $T_{v_j}$, denoted as ${\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)$, is computed as: ${\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)=d\left(v_a,v_b\right)+d\left(P_S\left(v_b,v_j\right)\right)-d\left(P_S\left(v_a,v_j\right)\right)\ge 0$. Note that ${\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)>0$ for $\left(v_a,v_b\right)\in E^{\prime }\backslash E\left(T_{v_j}\right)$ and ${\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)=0$ for $\left(v_a,v_b\right)$ in $T_{v_j}$. This implies that we can calculate the length of path $P\left(v_i,v_j\right)$ by adding the values of ${\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)$ for the arcs in $P\left(v_i,v_j\right)$ to the length of $P_S\left(v_i,v_j\right)$ since ${\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)=0$ for $\left(v_a,v_b\right)$ in $P_S\left(v_i,v_j\right)$. These arcs aside from $P_S\left(v_i,v_j\right)$ do not belong to $T_{v_j}$, and Eppstein’s algorithm calls these arcs the sidetracks of $P_S\left(v_i,v_j\right)$. The sidetracks are used to generate the k-shortest path of $P\left(v_i,v_j\right)$. For $P\left(v_i,v_j\right)$, 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 $P\left(v_i,v_j\right)$ using the depth first search algorithm in time $O\left(K\mathrm{l}\mathrm{o}\mathrm{g}K\right)$. The k-th shortest path of $P\left(v_i,v_j\right)$, for $k=1,\cdots ,K$, generated through the min-heap data structure, is called a candidate path from $v_i$ to $v_j$ and is denoted as $F^{{}^{\prime }}\left(v_i,v_j\right)$.

When generating $F^{{}^{\prime }}\left(v_i,v_j\right)$, we determine whether or not $F^{{}^{\prime }}\left(v_i,v_j\right)$ is included in the set of SNDF paths for $q\left(v_i,v_j\right)$ by Properties 1, 2, 3, and 4. First, $F^{\prime }\left(v_i,v_j\right)$ needs to have AF refueling availability by Property 1. If either $SS\left(F^{\prime }\left(v_i,v_j\right)\right)\cap ES\left(F^{{}^{\prime }}\left(v_i,v_j\right)\right)\ne \varnothing$ when $\left|CS\left(F^{\prime }\left(v_i,v_j\right)\right)\right|=1$ or $SS\left(F^{\prime }\left(v_i,v_j\right)\right)\ne \varnothing$, $ES\left(F^{\prime }\left(v_i,v_j\right)\right)\ne \varnothing$, and $d_{ijk}\left(m_{\left(l\right)},m_{\left(l+1\right)}\right)\le R$ for $m_{\left(l\right)},m_{\left(l+1\right)}\in CS\left(F^{\prime }\left(v_i,v_j\right)\right)$ when $\left|CS\left(F^{\prime }\left(v_i,v_j\right)\right)\right|>1$, then $F^{\prime }\left(v_i,v_j\right)$ has AF refueling availability. Second, if $F^{\prime }\left(v_i,v_j\right)$ 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 $d\left(F^{\prime }\left(v_i,v_j\right)\right)$ by Property 2. Next, regarding Properties 3 and 4, we compare $CS\left(F^{\prime }\left(v_i,v_j\right)\right)$ with $CS\left(F_k\left(v_i,v_j\right)\right)$’s when $k>1$. If $CS\left(F_k\left(v_i,v_j\right)\right)\subset CS\left(F^{\prime }\left(v_i,v_j\right)\right)$ for $k=1,\cdots ,K-1$, we diagnose whether or not $F^{\prime }\left(v_i,v_j\right)$ has irrelevant candidate sites by Property 3. If $CS\left(F^{\prime }\left(v_i,v_j\right)\right)\subseteq CS\left(F_k\left(v_i,v_j\right)\right)$ for $k=1,\cdots ,K-1$, then $F^{\prime }\left(v_i,v_j\right)$ is removed from the set of SNDF paths for $q\left(v_i,v_j\right)$ because the associated OD pair already has shortest feasible paths which dominate $F^{\prime }\left(v_i,v_j\right)$ 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’ ${\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)$ on the root of the binary heap becomes greater than the value of $\beta d\left(P_S\left(v_i,v_j\right)\right)$. 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 $v_i$’s to a certain interchange $v_j$, we similarly apply the above procedures to other interchanges $v_j$’s. Regarding $2K$ SNDF paths in $q\left(v_i,v_j\right)$, 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: $G\left(V,E\right),M,R,\beta ,$ and $K$

 

Output: $G^{\prime }\left(V^{\prime },E^{\prime }\right)$, sets of all $K$ SNDF original paths, and $K$ SNDF return paths for all OD pairs

 

Construct an expanded directed network $G^{\prime }$ by transforming dual-candidate sites into two separate single-candidate sites on each direction and adding all sites to $G$ as vertices.

 

Set $j=1$.

 

Repeat

 

   Build a shortest path arborescence $T_{v_j}$ rooted at vertex $v_j$.

 

  Set $i=1,i\ne j$.

 

  Repeat

 

   Set $k=0$.

 

   Construct a min binary heap data structure, denoted as $H_{\left(v_i,v_j\right)}$, with sequences of sidetracked arcs using the corresponding values of ${\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)$, $\forall \left(v_a,v_b\right)\in E\left(G^{\prime }\right)$, 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 $P_S\left(v_i,v_j\right)$.

 

   Repeat

 

    Construct a candidate path $F^{\prime }\left(v_i,v_j\right)$ by inserting a sequence of sidetracked arcs on the root of $H_{\left(v_i,v_j\right)}$ into $T_{v_j}$.

 

    Construct $CS\left(F^{\prime }\left(v_i,v_j\right)\right),SS\left(F^{\prime }\left(v_i,v_j\right)\right),ES\left(F^{\prime }\left(v_i,v_j\right)\right),\mathrm{a}\mathrm{n}\mathrm{d}IS\left(F^{\prime }\left(v_i,v_j\right)\right).$

 

    If $F^{\prime }\left(v_i,v_j\right)$ is satisfied with Properties 1, 2, 3, and 4, then define $F^{\prime }\left(v_i,v_j\right)$, $SS\left(F^{\prime }\left(v_i,v_j\right)\right)$, $ES\left(F^{\prime }\left(v_i,v_j\right)\right)$ and $IS\left(F^{\prime }\left(v_i,v_j\right)\right)$ as $F_k\left(v_i,v_j\right)$, $SS\left(F_k\left(v_i,v_j\right)\right)$, $ES\left(F_k\left(v_i,v_j\right)\right)$, and $IS\left(F_k\left(v_i,v_j\right)\right)$, respectively, and set $k=k+1$.

 

    Remove the current sequence of sidetracked arcs on the root of $H_{\left(v_i,v_j\right)}$, and insert new sequences of sidetracked arcs from the min-heap data structure to the binary heap $H_{\left(v_i,v_j\right)}$ by maintaining the min-heap property.

 

   Until$k>K$, ${\sum }_{\left(v_a,v_b\right)\in AE\left(P_S\left(v_i,v_j\right)\right)}{\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)>\beta \bullet d\left(P_S\left(v_i,v_j\right)\right)$, or the min-heap data structure does not have any sidetracked arcs.

 

   Set $i=i+1$.

 

  Until$i$ reaches $n_V+1$.

 

  Set $j=j+1$.

 

Until $j$ reaches $n_V+1$.

We prove below that at most $2K$ 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 transforming $G$ into $G^{\prime }$, given $q\left(v_i,v_j\right)\in Q$, $F^{\prime }\left(v_i,v_j\right)$ is generated by Eppstein’s algorithm. The $k$-th SNDF path algorithm determines $CS\left(F^{\prime }\left(v_i,v_j\right)\right)$, and then classifies the sequence of candidate sites into $SS\left(F^{\prime }\left(v_i,v_j\right)\right)$, $ES\left(F^{\prime }\left(v_i,v_j\right)\right)$, and $IS\left(F^{\prime }\left(v_i,v_j\right)\right)$ by measuring distances between the sites, $v_i$ and $v_j$. The algorithm checks first whether or not $F^{\prime }\left(v_i,v_j\right)$ meets Properties 1 and 2. Next, if the corresponding OD pair already has SNDF paths, then $CS\left(F^{\prime }\left(v_i,v_j\right)\right)$ is compared with $CS\left(F_k\left(v_i,v_j\right)\right)$ for $k=1,\cdots ,K-1$ by Properties 3 and 4. If $CS\left(F_k\left(v_i,v_j\right)\right)\subset CS\left(F^{\prime }\left(v_i,v_j\right)\right)$ for $k=1,\cdots ,K-1$, then the algorithm checks whether or not $F^{\prime }\left(v_i,v_j\right)$ has irrelevant candidate sites. If $CS\left(F^{\prime }\left(v_i,v_j\right)\right)\subseteq CS\left(F_k\left(v_i,v_j\right)\right)$ for $k=1,\cdots ,K-1$, then $F^{\prime }\left(v_i,v_j\right)$ is not included in the set of SNDF paths for $q\left(v_i,v_j\right)$. Otherwise, the algorithm assigns $F^{\prime }\left(v_i,v_j\right)$ to the $k$-th SNDF path $F_k\left(v_i,v_j\right)$ and increases the value of $k$ by one. If $k\le K$, a next temporary path can be constructed by adding new sequences of sidetracked arcs into $T_{v_j}$. The additional distance of the next path can be calculated as ${\sum }_{\left(v_a,v_b\right)\in AE\left(P_S\left(v_i,v_j\right)\right)}{\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)$, where $AE\left(P_S\left(v_i,v_j\right)\right)$ is the set of newly added arcs and ${\mathrm{\delta }}_{v_j}\left(v_a,v_b\right)=d\left(v_a,v_b\right)+d\left(P_S\left(v_b,v_j\right)\right)-d\left(P_S\left(v_a,v_j\right)\right)$. If the additional travel distance does not exceed $\beta ·d\left(P_S\left(v_i,v_j\right)\right)$, then the above processes are repeated until $k>K$. Otherwise, the algorithm moves to the next interchange $v_i$ in the shortest path arborescence $T_{v_j}$. 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 network $G$ 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, ${\widehat{L}}_o=\left\{F_k\left(v_i,v_j\right)|k=1,\cdots ,K,\forall q\left(v_i,v_j\right)\in Q\right\}$ and ${\widehat{L}}_r=\left\{F_k\left(v_j,v_i\right)|k=1,\cdots ,K,\forall q\left(v_i,v_j\right)\in Q\right\}$ denote the sets of all $K$ SNDF original paths and $K$ SNDF return paths for all OD pairs, respectively.

Theorem 2.

In order to generate $2K$ SNDF original and return paths for all OD pairs, the complexity of the $k$-th SNDF path algorithm takes time $O\left(\left|V\right|\left(\left|E^{\prime }\right|+\left|V^{\prime }\right|\mathrm{l}\mathrm{o}\mathrm{g}\left|V^{\prime }\right|+\left|V\right|\left(K{\left|M\right|}^2+\left|E^{\prime }\right|\right)K\mathrm{l}\mathrm{o}\mathrm{g}K\right)\right)$, where $\left|V^{\prime }\right|=n_V+n_{SM}+{2n}_{DM}$ and $\left|E^{\prime }\right|=n_E+2n_{SM}+{4n}_{DM}$. Recall that $\left|V\right|=n_V$, $\left|E\right|=n_E$, and $\left|M\right|=n_M.$ Also, $n_{SM}$ and $n_{DM}$ refer to the number of single- and dual-access candidate sites, respectively.

Proof.

Time $O\left(n_{DM}\right)$ is required to transform $G\left(V,E\right)$ with the set of all candidate sites to a simple expanded directed network $G^{\prime }\left(V^{\prime },E^{\prime }\right)$ by regarding dual-access candidate sites as two distinct single-access candidate sites on arcs. Given $G^{\prime }$, the basic version of Eppstein’s algorithm takes time $O\left(\left|V^{\prime }\right|+\left|V^{\prime }\right|\mathrm{l}\mathrm{o}\mathrm{g}\left|V^{\prime }\right|+\left|V\right|K\mathrm{l}\mathrm{o}\mathrm{g}K\right)$ to find the $k$-th shortest path to $v_j$ from any other interchanges. In time $O\left(\left|E^{\prime }\right|+\left|V^{\prime }\right|\mathrm{l}\mathrm{o}\mathrm{g}\left|V^{\prime }\right|\right)$, the algorithm constructs a shortest path arborescence and a min-heap data structure to output $F^{\prime }\left(v_i,v_j\right)$. Then, in time $O\left(K\mathrm{l}\mathrm{o}\mathrm{g}K\right)$, it generates $F^{\prime }\left(v_i,v_j\right)$’s from a specific interchange to $v_j$. After we construct $CS\left(F^{\prime }\left(v_i,v_j\right)\right)$ and classify $CS\left(F^{\prime }\left(v_i,v_j\right)\right)$ into $SS\left(F^{\prime }\left(v_i,v_j\right)\right)$, $ES\left(F^{\prime }\left(v_i,v_j\right)\right)$, and $IS\left(F^{\prime }\left(v_i,v_j\right)\right)$, it takes time $O\left({\left|M\right|}^2+\left|E^{\prime }\right|\right)$ for Property 1. Next, we find out which candidate sites are repeated in $F^{\prime }\left(v_i,v_j\right)$ in time $O\left({\left|M\right|}^2\right)$ for Property 2. Also, when $k>1$, time $O\left(K{\left|M\right|}^2\right)$ is required to compare $CS\left(F^{\prime }\left(v_i,v_j\right)\right)$ with $CS\left(F_k\left(v_i,v_j\right)\right)$ 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 $\left(1+\beta \right)d\left(P_S\left(v_i,v_j\right)\right)$ in time $O\left(\left(K{\left|M\right|}^2+\left|E^{\prime }\right|\right)K\mathrm{l}\mathrm{o}\mathrm{g}K\right)$. Next, we move to other interchanges $v_i$’s, $\forall v_i\in V\left(G\right)$, $i\ne j$, which requires time $O\left(\left|V\right|\left(K{\left|M\right|}^2+\left|E^{\prime }\right|\right)K\mathrm{l}\mathrm{o}\mathrm{g}K\right)$. In summary, generating $K$ SNDF paths from other interchanges to a certain interchange $v_j$ requires time $O\left(\left|E^{\prime }\right|+\left|V^{\prime }\right|\mathrm{l}\mathrm{o}\mathrm{g}\left|V^{\prime }\right|+\left|V\right|\left(K{\left|M\right|}^2+\left|E^{\prime }\right|\right)K\mathrm{l}\mathrm{o}\mathrm{g}K\right)$. Finally, since the number of interchanges for ODs is $\left|V\right|$, the total computational complexity of the algorithm takes time $O\left(\left|V\right|\left(\left|E^{\prime }\right|+\left|V^{\prime }\right|\mathrm{l}\mathrm{o}\mathrm{g}\left|V^{\prime }\right|+\left|V\right|\left(K{\left|M\right|}^2+\left|E^{\prime }\right|\right)K\mathrm{l}\mathrm{o}\mathrm{g}K\right)\right)$ 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 on $F_k\left(v_i,v_j\right)$, when $d\left(F_k\left(v_i,v_j\right)\right)>R$, at least one AF refueling station should be located at each candidate site in $SS\left(F_k\left(v_i,v_j\right)\right)$ and $ES\left(F_k\left(v_i,v_j\right)\right)$, and then locating additional stations at candidate sites in $IS\left(F_k\left(v_i,v_j\right)\right)$ are required if the distance between the two candidate sites in $SS\left(F_k\left(v_i,v_j\right)\right)$ and $ES\left(F_k\left(v_i,v_j\right)\right)$ is greater than $R$. On the other hand, when $d\left(F_k\left(v_i,v_j\right)\right)\le R$, since the distance between the two candidate sites in $SS\left(F_k\left(v_i,v_j\right)\right)$ and $ES\left(F_k\left(v_i,v_j\right)\right)$ is less than or equal to $R$, additional AF refueling stations are unnecessary to be considered for AF refueling stations between the sites in $SS\left(F_k\left(v_i,v_j\right)\right)$ and $ES\left(F_k\left(v_i,v_j\right)\right)$. Similar covering conditions are applied to the return path. From this observation, the following types of $k$-th SNDF path are defined according to $d\left(F_k\left(v_i,v_j\right)\right)$ and $d\left(F_k\left(v_j,v_i\right)\right)$ for $F_k\left(v_i,v_j\right)$ and $F_k\left(v_j,v_i\right)$, respectively, in OD pairs $q\left(v_i,v_j\right)\in Q$:

(i)

Type 1 $k$-th SNDF original path: $F_k^{\left(1\right)}\left(v_i,v_j\right)=\left\{F_k\left(v_i,v_j\right)\in {\widehat{L}}_o|0<d\left(F_k\left(v_i,v_j\right)\right)\le R\right\}$,

(ii)

Type 1 $k$-th SNDF return path: $F_k^{\left(1\right)}\left(v_j,v_i\right)=\left\{F_k\left(v_j,v_i\right)\in {\widehat{L}}_r|0<d\left(F_k\left(v_j,v_i\right)\right)\le R\right\}$,

(iii)

Type 2 $k$-th SNDF original path: $F_k^{\left(2\right)}\left(v_i,v_j\right)=\left\{F_k\left(v_i,v_j\right)\in {\widehat{L}}_o|R<d\left(F_k\left(v_i,v_j\right)\right)\right\}$,

(iv)

Type 2 $k$-th SNDF return path: $F_k^{\left(2\right)}\left(v_j,v_i\right)=\left\{F_k\left(v_j,v_i\right)\in {\widehat{L}}_r|R<d\left(F_k\left(v_j,v_i\right)\right)\right\}$.

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

In case of Type 1 $k$-th SNDF original paths, since the distance between two candidate sites in $SS\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)$ and $ES\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)$ is at most $R$, vehicles can travel from $v_i$ to $v_j$ by visiting two refueling stations at the sites. Sometimes, if there exists a single candidate site in the common area of $SS\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)$ and $ES\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)$, 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 $SS\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)$ and $ES\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)$. Figure 3 depicts trips in $F_k^{\left(1\right)}\left(v_i,v_j\right)$ covered by AF refueling stations at either $m_{\left(1\right)}\in SS\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)\cap ES\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)$ or $m_{\left(0\right)}\in SS\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)$ and $m_{\left(2\right)}\in ES\left(F_k^{\left(1\right)}\left(v_i,v_j\right)\right)$. Similar covering conditions are also applied to Type 1 $k$-th SNDF return paths.

[figure omitted; refer to PDF]

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

When two candidate sites in each $SS\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$ and $ES\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$ 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 $SS\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$ cannot reach the other stations in $ES\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$ due to the limited driving range, then additional refueling stations at candidate sites in $IS\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$ 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 $US\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)=\left\{m_{\left(l^{{}^{\prime }}\right)}\in IS\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)\cup ES\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)\right\}$, 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 ${\alpha }_{k,m_{\left(l\right)},m_{\left(l^{\prime }\right)}}$, where $m_{\left(l\right)}\in CS\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$ and $m_{\left(l^{\prime }\right)}\in US\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$, are determined as follows:$$\begin{array}{}\text{(6)}& {\alpha }_{k,m_{\left(l\right)},m_{\left(l^{\prime }\right)}}=\left\{\begin{array}{l}1,\text{if}{d}_{ijk}\left(v_{\left(0\right)},m_{\left(l\right)}\right)<d_{ijk}\left(v_{\left(0\right)},m_{\left(l^{\prime }\right)}\right)\text{and}0<d_{ijk}\left(m_{\left(l\right)},m_{\left(l^{\prime }\right)}\right)\le R,\\ 0,\text{otherwise}.\end{array}\right.\end{array}$$

Identification coefficient ${\alpha }_{k,m_{\left(l\right)},m_{\left(l^{\prime }\right)}}=1$ indicates that site $m_{\left(l^{\prime }\right)}$ follows $m_{\left(l\right)}$ along the $k$-th path and their distance is less than or equal to $R$. In an example in Figure 4, we have $CS\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)=⟨m_{\left(0\right)},m_{\left(1\right)},m_{\left(2\right)},m_{\left(3\right)},m_{\left(4\right)}⟩$ and $US\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)=\left\{m_{\left(2\right)},m_{\left(3\right)},m_{\left(4\right)}\right\}$. If $d_{ijk}\left(m_{\left(0\right)},m_{\left(3\right)}\right)>R$, $d_{ijk}\left(m_{\left(1\right)},m_{\left(3\right)}\right)\le R$, and $d_{ijk}\left(m_{\left(2\right)},m_{\left(3\right)}\right)\le R$, then the corresponding identification coefficients for site $m_{\left(3\right)}$ are set as ${\alpha }_{k,m_{\left(0\right)},m_{\left(3\right)}}=0$, ${\alpha }_{k,m_{\left(1\right)},m_{\left(3\right)}}=1$, ${\alpha }_{k,m_{\left(2\right)},m_{\left(3\right)}}=1$, ${\alpha }_{k,m_{\left(3\right)},m_{\left(3\right)}}=0$, and ${\alpha }_{k,m_{\left(4\right)},m_{\left(3\right)}}=0$. These values mean that a refueling station must be located at either site $m_{\left(1\right)}$ or $m_{\left(2\right)}$, so as to reach site $m_{\left(3\right)}$ in $F_k^{\left(2\right)}\left(v_i,v_j\right)$. 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, $US\left(F_k^{\left(2\right)}\left(v_j,v_i\right)\right)$ is similarly defined to determine identification coefficients.

[figure omitted; refer to PDF]

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 the $p$ 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\left(v_i,v_j\right)$ : average traffic flow in round trips per time unit between interchanges $v_i$ and $v_j$,

${\alpha }_{k,m_t,m_{t^{\prime }}}$ : identification coefficient for sites $m_t$ and $m_{t^{{}^{\prime }}}$ for Type 2 $k$-th SNDF path.

4.3.2. Decision Variables.

$$\begin{array}{}\text{(7)}& x_{m_t}=\left\{\begin{array}{l}1,\text{i}\text{f}\text{r}\text{e}\text{f}\text{u}\text{e}\text{l}\text{i}\text{n}\text{g}\text{s}\text{t}\text{a}\text{t}\text{i}\text{o}\text{n}\text{l}\text{o}\text{c}\text{a}\text{t}\text{i}\text{o}\text{n}{m}_t\in M\text{i}\text{s}\text{s}\text{e}\text{l}\text{e}\text{c}\text{t}\text{e}\text{d},\\ 0,\text{o}\text{t}\text{h}\text{e}\text{r}\text{w}\text{i}\text{s}\text{e},\end{array}\right.\end{array}$$ $$\begin{array}{}\text{(8)}& y_{k,v_i,v_j}=\left\{\begin{array}{l}1,\text{i}\text{f}\text{t}\text{h}\text{e}k-\text{t}\text{h}\text{S}\text{N}\text{D}\text{F}\text{p}\text{a}\text{t}\text{h}\text{f}\text{r}\text{o}\text{m}{v}_{i}to{v}_j\text{i}\text{s}\text{c}\text{o}\text{v}\text{e}\text{r}\text{e}\text{d}\text{b}\text{y}\text{t}\text{h}\text{e}\text{c}\text{u}\text{r}\text{r}\text{e}\text{n}\text{t}\text{s}\text{e}\text{t}\text{o}\text{f}\\ \text{r}\text{e}\text{f}\text{u}\text{e}\text{l}\text{i}\text{n}\text{g}\text{s}\text{t}\text{a}\text{t}\text{i}\text{o}\text{n}\text{s},\\ 0,\text{o}\text{t}\text{h}\text{e}\text{r}\text{w}\text{i}\text{s}\text{e},\end{array}\right.\end{array}$$ $$\begin{array}{}\text{(9)}& z_{v_i,v_j}=\left\{\begin{array}{l}1,\text{i}\text{f}\text{t}\text{h}\text{e}\text{r}\text{o}\text{u}\text{n}\text{d}\text{t}\text{r}\text{i}\text{p}\text{o}\text{f}\text{t}\text{r}\text{a}\text{f}\text{f}\text{i}\text{c}\text{f}\text{l}\text{o}\text{w}f\left(v_i,v_j\right)\text{i}s\text{c}\text{o}\text{v}\text{e}\text{r}\text{e}\text{d}\text{b}\text{y}\text{t}\text{h}\text{e}\text{c}\text{u}\text{r}\text{r}\text{e}\text{n}\text{t}\text{s}\text{e}\text{t}\text{o}\text{f}\\ \text{r}\text{e}\text{f}\text{u}\text{e}\text{l}\text{i}\text{n}\text{g}\text{s}\text{t}\text{a}\text{t}\text{i}\text{o}\text{n}\text{s},\\ 0,\text{o}\text{t}\text{h}\text{e}\text{r}\text{w}\text{i}\text{s}\text{e}.\end{array}\right.\end{array}$$

Then, we propose the following a mixed-integer programming model:$$\begin{array}{}\text{(10)}& \mathrm{maximize}\sum _{q\left(v_i,v_j\right)\in Q}f\left(v_i,v_j\right)z_{v_i,v_j},\end{array}$$

subject to$$\begin{array}{}\text{(11)}& \sum _{m_t\in SS\left(F_k\left(v_i,v_j\right)\right)}{x}_{m_t}\ge y_{k,v_i,v_j},\forall F_k\left(v_i,v_j\right)\in {\widehat{L}}_o,k=1,\cdots ,K,\end{array}$$$$\begin{array}{}\text{(12)}& \sum _{m_t\in ES\left(F_k\left(v_i,v_j\right)\right)}{x}_{m_t}\ge y_{k,v_i,v_j},\forall F_k\left(v_i,v_j\right)\in {\widehat{L}}_o,k=1,\cdots ,K,\end{array}$$$$\begin{array}{}\text{(13)}& \sum _{m_t\in SS\left(F_k\left(v_j,v_i\right)\right)}{x}_{m_t}\ge y_{k,v_j,v_i},\forall F_k\left(v_j,v_i\right)\in {\widehat{L}}_r,k=1,\cdots ,K,\end{array}$$$$\begin{array}{}\text{(14)}& \sum _{m_t\in ES\left(F_k\left(v_j,v_i\right)\right)}{x}_{m_t}\ge y_{k,v_j,v_i},\forall F_k\left(v_j,v_i\right)\in {\widehat{L}}_r,k=1,\cdots ,K,\end{array}$$$$\begin{array}{}\text{(15)}& \sum _{m_t\in CS\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)}{{\alpha }_{k,m_t,m_{t^{\prime }}}x}_{m_t}\ge y_{k,v_i,v_j},\forall m_{t^{\prime }}\in US\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right),\forall F_k^{\left(2\right)}\left(v_i,v_j\right)\in {\widehat{L}}_o,k=1,\cdots ,K,\end{array}$$$$\begin{array}{}\text{(16)}& \sum _{m_t\in CS\left(F_k^{\left(2\right)}\left(v_j,v_i\right)\right)}{{\alpha }_{{k,m}_t,m_{t^{\prime }}}x}_{m_t}\ge y_{{k,v}_j,v_i},\forall m_{t^{\prime }}\in US\left(F_k^{\left(2\right)}\left(v_j,v_i\right)\right),\forall F_k^{\left(2\right)}\left(v_j,v_i\right)\in {\widehat{L}}_r,k=1,\cdots ,K,\end{array}$$$$\begin{array}{}\text{(17)}& \sum _{k=1,\cdots ,K}{y}_{k,v_i,v_j}\ge z_{v_i,v_j},\forall q\left(v_i,v_j\right)\in Q,\end{array}$$$$\begin{array}{}\text{(18)}& \sum _{k=1,\cdots ,K}{y}_{k,v_j,v_i}\ge z_{v_i,v_j},\forall q\left(v_i,v_j\right)\in Q,\end{array}$$$$\begin{array}{}\text{(19)}& \sum _{m_t\in M}{x}_{m_t}=p,\end{array}$$$$\begin{array}{}\text{(20)}& x_{m_t}\in \left\{0,1\right\},\forall m_t\in M,\end{array}$$$$\begin{array}{}\text{(21)}& y_{k,v_i,v_j}\in \left\{0,1\right\},\forall F_k\left(v_i,v_j\right)\in {\widehat{L}}_o,k=1\cdots K,\end{array}$$$$\begin{array}{}\text{(22)}& y_{k,v_j,v_i}\in \left\{0,1\right\},\forall F_k\left(v_j,v_i\right)\in {\widehat{L}}_r,k=1\cdots K,\end{array}$$$$\begin{array}{}\text{(23)}& z_{v_i,v_j}\in \left\{0,1\right\},\forall q\left(v_i,v_j\right)\in Q.\end{array}$$

In this formulation, the objective function (10) maximizes the traffic flow that is covered by the set of AF refueling stations located at the $p$ 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 $SS\left(F_k\left(v_i,v_j\right)\right)$ and another one in set of candidate sites $ES\left(F_k\left(v_i,v_j\right)\right)$. 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 ${\alpha }_{k,m_t,m_{t^{\prime }}}$ for $m_t\in CS\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$ and $m_{t^{\prime }}\in US\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$, where ${\alpha }_{k,m_t,m_{t^{\prime }}}=1$ if $d_{ijk}\left(v_i,m_t\right)<d_{ijk}\left(v_i,m_{t^{\prime }}\right)$ and $0<d_{ijk}\left(m_t,m_{t^{\prime }}\right)\le R$; otherwise, ${\alpha }_{k,m_t,m_{t^{\prime }}}=0$. This constraint set guarantees that, if a refueling station is located at $m_{t^{\prime }}$ in $US\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$, then one of sites in $CS\left(F_k^{\left(2\right)}\left(v_i,v_j\right)\right)$ should be selected to locate another station which would provide refueling for vehicles to reach the station at $m_{t^{\prime }}$. 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:$$\begin{array}{}\text{(24)}& \sum _{m_t\in M}{c}_{m_t}{x}_{m_t}\le B,\end{array}$$

where $c_{m_t}$ refers to the capital cost for building a refueling station in candidate site $m_t\in 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 the $k$-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, $\beta =0\%,50\%,$ and 100%, and three driving ranges, $R=15,20,$ and 30.

[figure omitted; refer to PDF]

All computational experiments in Section 5 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 $\beta$ (i.e., $\beta =$ 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.

Table 4 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 for $p=1$–$25$ (in sec).

$K$	$R$	$\beta =0$			$\beta =50\%$			$\beta =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 Figure 6 show the coverage of traffic flows (in %) with $\beta =$ 0%, 50%, and 100% of $d\left(P_S\left(v_i,v_j\right)\right)$ as we increase the number of refueling stations from $p=1$ to $25$. Note that we fix $R=20$ in this subsection.

[figures omitted; refer to PDF]

First, when $\beta =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 $q\left(v_3,v_{13}\right)$, i.e., $F_1\left(v_{13},v_3\right)=⟨v_{13}=v_{\left(0\right)},v_8=v_{\left(1\right)},v_4=v_{\left(2\right)},v_3=v_{\left(3\right)}⟩$ is not covered. In the case of $K=3$ and $5$ with $\beta =0\%$, the mixed-integer programming model locates one refueling station at site $m_{14}$ and then $F_3\left(v_{13},v_3\right)=⟨v_{13}=v_{\left(0\right)},v_{10}=v_{\left(1\right)},v_9=v_{\left(2\right)},v_3=v_{\left(3\right)}⟩$, which is the second SNDF return path of $q\left(v_3,v_{13}\right)$, 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 $\beta =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 $\beta =50\%$, while the improvement increases up to 13.07% when $\beta =100\%$. The $k$-th SNDF path algorithm can generate more SNDF paths under $\beta =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 $\beta =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 $\beta =0\%$, 50%, and 100% for $p=1,\cdots ,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 $\beta =0\%$, $\beta =50\%$, and $\beta =100\%$ when $R=20$.

$p$	$\beta =0\%$			$\beta =50\%$			$\beta =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 path $P_S\left(v_i,v_j\right)$ and shortest return path $P_S\left(v_j,v_i\right)$ up to $d\left(P_S\left(v_i,v_j\right)\right)\left(1+\beta \right)$ and $d\left(P_S\left(v_j,v_i\right)\right)\left(1+\beta \right)$, 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_{\beta }$, and the average deviation factor values of the vehicles that deviate from their shortest paths, denoted as $\overline{\beta }$, 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 $\beta$, we observe that the value of $C_{\beta }$ 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_{\beta }$ tends to increase when refueling stations are scarce, while the value of $C_{\beta }$ starts decreasing after the traffic flows are mostly covered. These patterns imply that the effect of $\beta$ 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 $\overline{\beta }$, we have an extreme case where the value of $\overline{\beta }$ drops rapidly when we consider $\beta =100\%$ and the number of refueling stations is small. However, in most cases, the value of $\overline{\beta }$ 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 flows $\left(C_{\beta }\right)$ and the average deviation factor $\left(\overline{\beta }\right)$ for $p=1,\cdots ,25$.

$p$	$K$ = 1							$K$ = 3							$K$ = 5
	$\beta$  = 0%	$\beta$  = 50%			$\beta$  = 100%			$\beta$  = 0%	$\beta$  = 50%			$\beta$  = 100%			$\beta$  = 0%	$\beta$  = 50%			$\beta$  = 100%
	$C$ (%)	$C$ (%)	$C_{\beta }$ (%)	$\overline{\beta }$ (%)	$C$ (%)	$C_{\beta }$ (%)	$\overline{\beta }$ (%)	$C$ (%)	$C$ (%)	$C_{\beta }$ (%)	$\overline{\beta }$ (%)	$C$ (%)	$C_{\beta }$ (%)	$\overline{\beta }$ (%)	$C$ (%)	$C$ (%)	$C_{\beta }$ (%)	$\overline{\beta }$ (%)	$C$ (%)	$C_{\beta }$ (%)	$\overline{\beta }$ (%)
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, Figures 7(a) and 7(c) provide comparisons of optimal coverages of no-detour versus detour flows when $R=15$ and 30 for $p=1,\cdots ,25$. We set the coverages of traffic flows with $\beta =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 $\beta =$ $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 $d\left(P_S\left(v_i,v_j\right)\right)$ 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 ($\beta =0\%$) and detour flows ($\beta =50\%$ or $\beta =100\%$) continuously increases until the number of refueling stations reaches to some level ($p=21$ for $\beta =50\%$, $p=22$ for $\beta =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 $\beta =0\%$ and $\beta =50\%$ is insignificant. However, if vehicles are able to deviate further from their preplanned paths, i.e., $\beta =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$.

[figures omitted; refer to PDF]

Table 7

Comparisons of the optimal coverages of no-detour ($\beta =0\%$) versus detour flows ($\beta =50\%,100\%$) and their coverage differences when $R=15,30$.

$p$	$R=15$					$R=30$
	$\beta$			Differences between 0% and 50%	Differences between 0% and 100%	$\beta$			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 Figure 8. 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 omitted; refer to PDF]

[figures omitted; refer to PDF]

Table 8 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 ($\beta =0\%$, $\beta =50\%$ with $K=1$, and $\beta =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 ($\beta =0\%$) versus detour flows ($\beta =50\%$) with $K=1,3$ when $R=300,600$.

Budget ($M)	$R=300$									$R=600$
	No-detour with $\beta =0\%$			Detour with $K=1$ and $\beta =50\%$			Detour with $K=3$ and $\beta =50\%$			No-detour with $\beta =0\%$			Detour with $K=1$ and $\beta =50\%$			Detour with $K=3$ and $\beta =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.33

First, we easily identify the effect of detour flows (both $\beta$ 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 with $R=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 Figure 9, 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 ($\beta =0\%$) and detour flows ($\beta =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 $d\left(P_S\left(v_i,v_j\right)\right)$, 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 Table 8 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 Subsection 3.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 the $k$-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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgements

The first author was supported by the Research Initiation Funds provided to new faculty by Hongik University. The second author was supported by the Settlement Research Funds (Project Number 1.190090.01) of Ulsan National Institute of Science & Technology (UNIST).

References

[1] U S Environmental Protection Agency (EPA), "Inventory of U.S. greenhouse gas emissions and sinks: 1990–2017," 2017. https://www.epa.gov/ghgemissions/inventory-us-greenhouse-gas-emissions-and-sinks , EPA Report on Greenhouse Gas Emissions, accessed date: February 22, 2019

[2] U S Department of Energy (DOE), "Alternative fuels data center: alternative fuel and advanced vehicle search," 2019. http://www.afdc.energy.gov/vehicles/search/ , DOE Report on Energy Efficiency and Renewable Energy, accessed date: February 03, 2019

[3] U S Department of Energy (DOE), "Alternative fuels data center: alternative fueling station locator," 2017. http://www.afdc.energy.gov/locator/stations/results?fuel=ELEC , DOE Report on Energy Efficiency and Renewable Energy, accessed date: February 25, 2017

[4] U S Census Bureau, "Reference – urban and rural: county classification lookup table," 2017. https://www.census.gov/geo/reference/urban-rural.html , Geography Report, accessed date: April 02, 2017

[5] J. A. Ventura, S. J. Kweon, S. W. Hwang, M. Tormay, C. Li, "Energy policy considerations in the design of an alternative-fuel refueling infrastructure to reduce GHG emissions on a transportation network," Energy Policy, vol. 111, pp. 427-439, DOI: 10.1016/j.enpol.2017.09.035, 2017.

[6] S. W. Hwang, S. J. Kweon, J. A. Ventura, "Infrastructure development for alternative fuel vehicles on a highway road system," Transportation Research Part E: Logistics and Transportation Review, vol. 77, pp. 170-183, DOI: 10.1016/j.tre.2015.02.011, 2015.

[7] S. W. Hwang, S. J. Kweon, J. A. Ventura, "Locating alternative-fuel refueling stations on a multi-class vehicle transportation network," European Journal of Operational Research, vol. 261 no. 3, pp. 941-957, DOI: 10.1016/j.ejor.2017.02.036, 2017.

[8] M. Holder, "Bikes, EVs and LNG trucks: how delivery giant UPS is steering “beyond diesel”," 2019. https://www.greenbiz.com/article/bikes-evs-and-lng-trucks-how-delivery-giant-ups-steering-beyond-diesel

[9] UPS Pressroom, "UPS invests more than $90 million in natural gas vehicles and infrastructure," 2017. https://pressroom.ups.com/pressroom/ContentDetailsViewer.page?ConceptType=PressReleases&id=1489579573572-162 , United Parcel Service, accessed date: October 09, 2019

[10] W. Hoffman, R. Pavley, "A method for the solution of the n -th best path problem," Journal of the ACM (JACM), vol. 6 no. 4, pp. 506-514, DOI: 10.1145/320998.321004, 1959.

[11] B. L. Fox, "K -th shortest paths and applications to probabilistic networks," Operations Research, vol. 23, pp. B263-B263, 1975.

[12] D. Eppstein, "Finding the k shortest paths," SIAM Journal on Computing, vol. 28 no. 2, pp. 652-673, DOI: 10.1137/S0097539795290477, 1998.

[13] V. M. Jiménez, A. Marzal, "A lazy version of eppstein’s k shortest paths algorithm," The Second International Workshop on Experimental and Efficient Algorithms, vol. 2647, pp. 179-191, 2003.

[14] H. Aljazzar, S. Leue, "K ∗ : A heuristic search algorithm for finding the k shortest paths," Artificial Intelligence, vol. 175 no. 18, pp. 2129-2154, DOI: 10.1016/j.artint.2011.07.003, 2011.

[15] G. Liu, Z. Qiu, H. Qu, L. Ji, "Computing k shortest paths using modified pulse-coupled neural network," Neurocomputing, vol. 149, pp. 1162-1176, DOI: 10.1016/j.neucom.2014.09.012, 2015.

[16] J. D. Adler, P. B. Mirchandani, G. Xue, M. Xia, "The electric vehicle shortest-walk problem with battery exchanges," Networks and Spatial Economics, vol. 16 no. 1, pp. 155-173, DOI: 10.1007/s11067-013-9221-7, 2016.

[17] E. Q. Martins, M. M. Pascoal, "A new implementation of Yen’s ranking loopless paths algorithm," Quarterly Journal of the Belgian, French and Italian Operations Research Societies, vol. 1 no. 2, pp. 121-133, DOI: 10.1007/s10288-002-0010-2, 2003.

[18] M. Pollack, "Solutions of the k th best route through a network—a review," Journal of Mathematical Analysis and Applications, vol. 3 no. 3, pp. 547-559, DOI: 10.1016/0022-247X(61)90076-2, 1961.

[19] S. Clarke, A. Krikorian, J. Rausen, "Computing the N best loopless paths in a network," Journal of the Society for Industrial and Applied Mathematics, vol. 11 no. 4, pp. 1096-1102, 1963.

[20] M. Sakarovitch, "The K shortest chains in a graph," Transportation Research, vol. 2 no. 1,DOI: 10.1016/0041-1647(68)90003-8, 1968.

[21] J. Y. Yen, "Finding the k shortest loopless paths in a network," Management Science, vol. 17 no. 11, pp. 712-716, 1971.

[22] E. L. Lawler, "A procedure for computing the k best solutions to discrete optimization problems and its application to the shortest path problem," Management Science, vol. 18 no. 7, pp. 401-405, 1972.

[23] N. Katoh, T. Ibaraki, H. Mine, "An efficient algorithm for k shortest simple paths," Networks, vol. 12 no. 4, pp. 411-427, DOI: 10.1002/net.3230120406, 1982.

[24] E. Hadjiconstantinou, N. Christofides, "An efficient implementation of an algorithm for finding K shortest simple paths," Networks, vol. 34 no. 2, pp. 88-101, DOI: 10.1002/(SICI)1097-0037(199909)34:2<88:AID-NET2>3.0.CO;2-1, 1999.

[25] E. D. Q. V. Martins, E. Queir, J. L. E. dos Santos, V. Martins, M. Margarida, M. M. B. Pascoal, "A new algorithm for ranking loopless paths," 1997. Research Report

[26] A. Perko, "Implementation of algorithms for K shortest loopless paths," Networks, vol. 16 no. 2, pp. 149-160, DOI: 10.1002/net.3230160204, 1986.

[27] J. Hershberger, S. Suri, "Vickrey prices and shortest paths: what is an edge worth?," pp. 252-259, DOI: 10.1109/SFCS.2001.959899, .

[28] J. Hershberger, M. Maxel, S. Suri, "Finding the k shortest simple paths: a new algorithm and its implementation," ACM Transactions on Algorithms (TALG), vol. 3 no. 4, pp. 851-864, 2007.

[29] W. Zeng, T. Miwa, T. Morikawa, "Application of the support vector machine and heuristic k -shortest path algorithm to determine the most eco-friendly path with a travel time constraint," Transportation Research Part D: Transport and Environment, vol. 57, pp. 458-473, DOI: 10.1016/j.trd.2017.10.001, 2017.

[30] M. Kuby, S. Lim, "The flow-refueling location problem for alternative-fuel vehicles," Socio-Economic Planning Sciences, vol. 39 no. 2, pp. 125-145, DOI: 10.1016/j.seps.2004.03.001, 2005.

[31] C. Upchurch, M. Kuby, "Comparing the p -median and flow-refueling models for locating alternative-fuel stations," Journal of Transport Geography, vol. 18 no. 6, pp. 750-758, DOI: 10.1016/j.jtrangeo.2010.06.015, 2010.

[32] S. Lim, M. Kuby, "Heuristic algorithms for siting alternative-fuel stations using the flow-refueling location model," European Journal of Operational Research, vol. 204 no. 1, pp. 51-61, DOI: 10.1016/j.ejor.2009.09.032, 2010.

[33] M. Kuby, L. Lines, R. Schultz, Z. Xie, J. G. Kim, S. Lim, "Optimization of hydrogen stations in Florida using the flow-refueling location model," International Journal of Hydrogen Energy, vol. 34 no. 15, pp. 6045-6064, DOI: 10.1016/j.ijhydene.2009.05.050, 2009.

[34] I. Capar, M. Kuby, "An efficient formulation of the flow refueling location model for alternative-fuel stations," IIE Transactions, vol. 44 no. 8, pp. 622-636, DOI: 10.1080/0740817X.2011.635175, 2012.

[35] I. Capar, M. Kuby, V. J. Leon, Y. J. Tsai, "An arc cover-path-cover formulation and strategic analysis of alternative-fuel station locations," European Journal of Operational Research, vol. 227 no. 1, pp. 142-151, DOI: 10.1016/j.ejor.2012.11.033, 2013.

[36] P. Jochem, C. Brendel, M. Reuter-Oppermann, W. Fichtner, S. Nickel, "Optimizing the allocation of fast charging infrastructure along the German autobahn," Journal of Business Economics, vol. 86 no. 5, pp. 513-535, DOI: 10.1007/s11573-015-0781-5, 2015.

[37] M. Kuby, S. Lim, "Location of alternative-fuel stations using the flow-refueling location model and dispersion of candidate sites on arcs," Networks and Spatial Economics, vol. 7 no. 2, pp. 129-152, DOI: 10.1007/s11067-006-9003-6, 2007.

[38] J. A. Ventura, S. W. Hwang, S. J. Kweon, "A continuous network location problem for a single refueling station on a tree," Computers & Operations Research, vol. 62, pp. 257-265, DOI: 10.1016/j.cor.2014.07.017, 2015.

[39] S. J. Kweon, S. W. Hwang, J. A. Ventura, "A continuous deviation-flow location problem for an alternative-fuel refueling station on a tree-like transportation network," Journal of Advanced Transportation, vol. 2017,DOI: 10.1155/2017/1705821, 2017.

[40] J. He, H. Yang, T. Q. Tang, H. J. Huang, "An optimal charging station location model with the consideration of electric vehicle’s driving range," Transportation Research Part C: Emerging Technologies, vol. 86, pp. 641-654, DOI: 10.1016/j.trc.2017.11.026, 2018.

[41] Y. W. Wang, C. C. Lin, "Locating road-vehicle refueling stations," Transportation Research Part E: Logistics and Transportation Review, vol. 45 no. 5, pp. 821-829, DOI: 10.1016/j.tre.2009.03.002, 2009.

[42] Y. W. Wang, C. R. Wang, "Locating passenger vehicle refueling stations," Transportation Research Part E: Logistics and Transportation Review, vol. 46 no. 5, pp. 791-801, DOI: 10.1016/j.tre.2009.12.001, 2010.

[43] S. A. MirHassani, R. Ebrazi, "A flexible reformulation of the refueling station location problem," Transportation Science, vol. 47 no. 4, pp. 617-628, DOI: 10.1287/trsc.1120.0430, 2013.

[44] S. H. Chung, C. Kwon, "Multi-period planning for electric car charging station locations: a case of Korean expressways," European Journal of Operational Research, vol. 242 no. 2, pp. 677-687, DOI: 10.1016/j.ejor.2014.10.029, 2015.

[45] M. Hosseini, S. A. MirHassani, "Refueling-station location problem under uncertainty," Transportation Research Part E: Logistics and Transportation Review, vol. 84, pp. 101-116, DOI: 10.1016/j.tre.2015.10.009, 2015.

[46] J. E. Kang, W. Recker, "Strategic hydrogen refueling station locations with scheduling and routing considerations of individual vehicles," Transportation Science, vol. 49 no. 4, pp. 767-783, DOI: 10.1287/trsc.2014.0519, 2014.

[47] M. Miralinaghi, B. B. Keskin, Y. Lou, A. M. Roshandeh, "Capacitated refueling station location problem with traffic deviations over multiple time periods," Networks and Spatial Economics, vol. 17 no. 1, pp. 129-151, DOI: 10.1007/s11067-016-9320-3, 2017.

[48] M. Miralinaghi, Y. Lou, B. B. Keskin, A. Zarrinmehr, R. Shabanpour, "Refueling station location problem with traffic deviation considering route choice and demand uncertainty," International Journal of Hydrogen Energy, vol. 42 no. 5, pp. 3335-3351, DOI: 10.1016/j.ijhydene.2016.12.137, 2017.

[49] M. Schiffer, G. Walther, "An adaptive large neighborhood search for the location-routing problem with intra-route facilities," Transportation Science, vol. 52 no. 2, pp. 331-352, DOI: 10.1287/trsc.2017.0746, 2018.

[50] J. Hof, M. Schneider, D. Goeke, "Solving the battery swap station location-routing problem with capacitated electric vehicles using an AVNS algorithm for vehicle-routing problems with intermediate stops," Transportation Research Part B: Methodological, vol. 97, pp. 102-112, DOI: 10.1016/j.trb.2016.11.009, 2017.

[51] A. Gagarin, P. Corcoran, "Multiple domination models for placement of electric vehicle charging stations in road networks," Computers & Operations Research, vol. 96, pp. 69-79, DOI: 10.1016/j.cor.2018.03.014, 2018.

[52] T. H. Tran, G. Nagy, T. B. T. Nguyen, N. A. Wassan, "An efficient heuristic algorithm for the alternative-fuel station location problem," European Journal of Operational Research, vol. 269 no. 1, pp. 159-170, DOI: 10.1016/j.ejor.2017.10.012, 2018.

[53] A. Langer, S. McRae, "Fueling alternatives: evidence from naturalistic driving data," 2013. Working Paper

[54] J. G. Kim, M. Kuby, "The deviation-flow refueling location model for optimizing a network of refueling stations," International Journal of Hydrogen Energy, vol. 37 no. 6, pp. 5406-5420, DOI: 10.1016/j.ijhydene.2011.08.108, 2012.

[55] J. G. Kim, M. Kuby, "A network transformation heuristic approach for the deviation flow refueling location model," Computers & Operations Research, vol. 40 no. 4, pp. 1122-1131, DOI: 10.1016/j.cor.2012.10.021, 2013.

[56] Y. Huang, S. Li, Z. S. Qian, "Optimal deployment of alternative fueling stations on transportation networks considering deviation paths," Networks and Spatial Economics, vol. 15 no. 1, pp. 183-204, DOI: 10.1007/s11067-014-9275-1, 2015.

[57] B. Yıldız, O. Arslan, O. E. Karaşan, "A branch and price approach for routing and refueling station location model," European Journal of Operational Research, vol. 248 no. 3, pp. 815-826, DOI: 10.1016/j.ejor.2015.05.021, 2016.

[58] P. Göpfert, S. Bock, "A branch&cut approach to recharging and refueling infrastructure planning," European Journal of Operational Research, vol. 279 no. 3,DOI: 10.1016/j.ejor.2019.06.031, 2019.

[59] O. Arslan, O. E. Karaşan, A. R. Mahjoub, H. Yaman, "A branch-and-cut algorithm for the alternative fuel refueling station location problem with routing," Transportation Science, vol. 53 no. 4, pp. 1107-1125, DOI: 10.1287/trsc.2018.0869, 2019.

[60] S. J. Kweon, "Location of alternative-fuel refueling stations on transportation networks considering vehicle deviations and greenhouse gas emissions," 2017. Doctoral dissertation

[61] D. Simchi-Levi, O. Berman, "A heuristic algorithm for the traveling salesman location problem on networks," Operations Research, vol. 36 no. 3, pp. 478-484, 1988.

[62] A. S. Fotheringham, M. E. O’Kelly, Spatial Interaction Models: Formulations and Applications, vol. 1, 1989.

[63] Oak Ridge National Laboratory, "Freight analysis framework version 4 freight traffic assignment," 2016. https://faf.ornl.gov/faf4/Documentation.aspx , accessed date: October 11, 2019

[64] A. Burnham, Liquefied Natural Gas. Case Study, 2013.

[65] J. Park, "Natural gas technology improvements changing attitudes," 2014. http://www.truckinginfo.com/article/story/2014/09/natural-gas.aspx , accessed date: October 10, 2019

[66] C. Upchurch, M. Kuby, S. Lim, "A model for location of capacitated alternative-fuel stations," Geographical Analysis, vol. 41 no. 1, pp. 85-106, DOI: 10.1111/j.1538-4632.2009.00744.x, 2009.