# Optimal Facility Location Model Based on Genetic Simulated Annealing Algorithm for Siting Urban Refueling Stations

Authors: Dawei Chen; Shun Zhou; Yuanchang Xie; Xuhong Li
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper analyzes the impact factors and principles of siting urban refueling stations and proposes a three-stage method. The main objective of the method is to minimize refueling vehicles’ detour time. The first stage aims at identifying the most frequently traveled road segments for siting refueling stations. The second stage focuses on adding additional refueling stations to serve vehicles whose demands are not directly satisfied by the refueling stations identified in the first stage. The last stage further adjusts and optimizes the refueling station plan generated by the first two stages. A genetic simulated annealing algorithm is proposed to solve the optimization problem in the second stage and the results are compared to those from the genetic algorithm. A case study is also conducted to demonstrate the effectiveness of the proposed method and algorithm. The results indicate the proposed method can provide practical and effective solutions that help planners and government agencies make informed refueling station location decisions.


---

## Body

## 1. Introduction

With the rapid economic development in China, automobile ownership has increased drastically in the past decade. Due to this, there is an urgent need for upgrading transportation infrastructure to accommodate the newly generated traffic. Refueling stations, either gasoline stations or electric vehicle charging stations, are a very important type of service facility and are essential for the appropriate operations of urban transportation systems. The locations of refueling stations can directly affect the efficacy of urban road traffic system. Inappropriate location designs can substantially increase vehicle travel distances for finding refueling stations and cause unnecessary congestion and emissions. Therefore, investigating the optimal refueling station locations has important practical significance for developing countries such as China. This research is also important for some developed countries that are investing heavily in alternative-fueled vehicles (e.g., electric vehicles), where needs for optimally locating electric vehicle charging stations exist.Some existing refueling station location models are based on qualitative analysis. These methods can generally be categorized into two groups: maximum coverage and maximum profit methods [1]. The maximum coverage methods are mainly used for high population density and congested areas. Their objective is to maximize the coverage of refueling stations and reduce the unnecessary travel distance caused by vehicles searching for refueling stations. The maximum profit methods are designed for low population density areas. Their primary goal is to ensure that refueling stations at selected locations can have adequate demand to maintain certain level of profit [2]. Both methods cannot provide quantitative results to demonstrate the economic and traffic impacts of different refueling station location plans.Several related quantitative researches have also been conducted. Current et al. [3] investigated a network design problem. They formulated it as a Maximum Covering/Shortest Path (MCSP) model with two objectives. The first objective was to find the shortest path between any pair of start and end points, and the second objective was to maximize the number of demand nodes covered by the shortest path. In their research, if a node is on the shortest path or within certain distance of any nodes on the shortest path, then it is considered as being covered. Bapna et al. [4] extended the MCSP formulation to a Maximum Covering/Shortest Spanning Subgraph Problem (MC3SP) that considers multiple pairs of start and end points. The focus of their research was to optimally locate unleaded gasoline stations to facilitate trips among large and populated cities. Building upon the MCSP and MC3SP, Y.-W. Wang and C.-R. Wang [5] formulated a flow-based set covering location model. Similar to [3, 4], they assumed the refueling demand of a city to be concentrated on the center of the city. In other words, a city was considered as a single demand point. All these models [3–5] tried to cover as many nodes or population along the shortest path as possible while keeping the total cost low. However, a more important issue in refueling station location modeling is to cover more OD (origin and destination) flows with refueling demand. Also different from these studies, the focus of our research is to investigate refueling station location problems within urban areas rather than for intercity trips.Goodchild and Noronha [6] proposed a combined model to site a finite number of refueling stations to maximize the refueling market share of a firm. As pointed out by Kuby and Lim [7], a major problem of this combined model is that it used link flow data instead of OD flow data as the input. Therefore, OD flows may be double-counted as they pass through several road links. Hodgson [8] conducted one of the pioneer studies on siting refueling stations and developed a Flow Capturing Location Model (FCLM). The goal of the FCLM was to capture more OD flows (instead of more link flows as in [6]) on the shortest paths between OD pairs. This can avoid double-counting OD flows captured by multiple refueling stations along their shortest paths. A number of improvements have been done based on the FCLM. Hodgson and Rosing [9] extended the FCLM and proposed a hybrid flow capturing plus p-median model with two conflicting objectives, which were maximizing the coverage of OD flows and minimizing demand-weighted travel distance. Later on, Kuby and Lim [7] proposed a Flow-Refueling Location Model (FRLM) based on the FCLM. A major difference between FCLM and FRLM is that a facility combination concept was proposed. OD flows were covered by a combination of facilities instead of individual facilities in the FRLM model. Recently, Lim and Kuby proposed a number of heuristic algorithms for solving the developed FRLM [10, 11]. These studies [7–10] all considered a single shortest path connecting each OD pair. Their common objective was to optimally locate a prespecified amount of refueling stations such that the total number of shortest paths being covered could be maximized. In the real world, traffic flow between an OD pair rarely takes a single shortest path. Travelers typically follow the user equilibrium (UE) principle and try to minimize their individual travel times. This usually results in many different routes between an OD pair. Unfortunately, this diverse route choice behavior is not considered in these refueling station location studies [7–10].In the following subsection, a comprehensive analysis of factors affecting refueling station location planning is presented. Based upon this analysis, a three-stage method is proposed to optimally locate refueling stations. Instead of focusing on the regional level and considering each city as a single demand point [3–5], this study centers on urban areas and the entire city is divided into small traffic analysis zones (TAZs). In addition, the model introduced in this study explicitly considers the diverse route choice behavior between origins and destinations and can better reflect real world traffic scenarios.


## 2. Principles for Siting Urban Refueling Stations

### 2.1. Principles of Location Design

Siting urban refueling stations should consider a comprehensive set of factors. Government agencies are mostly concerned about congestion, safety, environment, and land use factors, while refueling station operators are more concerned about demand and profit. Besides, fuel price difference at various refueling stations is also an important factor to consider. However, in China gasoline prices are decided by the government and there is almost no difference in terms of price among different gas stations. Therefore, the price competition among different refueling stations/companies is not considered in this study. In general, siting refueling stations should (a) minimize their impacts on urban traffic, (b) be on main roads to ensure profitability, (c) keep an appropriate distance between refueling stations, (d) be far from residential and other sensitive areas, and (e) be based on urban land use planning.Principle (a) is from the perspective of minimizing the influence of refueling stations on urban road traffic. For example, a vehicle with refueling need should not have to detour too far for refueling. Intuitively, this requires refueling stations to be located on or near the paths of most vehicles, so that they do not need to make a detour or just need to make a very short detour.Principle (b) considers the economic benefit of refueling stations. The number of refueling vehicles that may potentially be served by a refueling station should be taken into consideration during the location design. Ideally, this number should be within a reasonable range to ensure that a refueling station’s profit is above certain level and the refueling demand does not exceed its capacity.Principle (c) takes the service ranges of refueling stations into consideration. It is to make sure that refueling stations are not too concentrated in certain areas. Keeping an appropriate distance between refueling stations sometimes is necessary due to reasons such as safety.Principle (d) is mainly about the safety and environmental concerns. Based on this principle, siting refueling stations should consider the safety of refueling station facilities and equipment, refueling vehicles, travelers, and nearby residents. Also, refueling stations should be far away from some environmentally and socially sensitive areas.Last but certainly not least, principle (e) optimizes refueling station locations from the urban planning point of view. Specifically, urban refueling station locations should take into account land use development and future land use planning. Land use directly affects trip productions and attractions. Consequently, it has major impacts on the distribution of future refueling demand and refueling station locations.From the standpoint of model development for optimally siting refueling stations, principle (a) can be considered in the objective function, and principles (b) through (e) can be either included as constraints or used to determine candidate locations for future refueling stations.


### 2.2. Three-Stage Model

Based on the analysis in the previous section, if principle (a) is adopted as the objective, there will be two cases in the refueling station location planning. The first case is that refueling stations are located on the paths of some refueling vehicles. Hence, these vehicles do not need to detour for refueling. The other case is that refueling stations are not the paths of the remaining refueling vehicles and these vehicles have to detour for refueling. In light of these two cases, siting refueling stations can be implemented in three stages. In the first stage, candidate refueling stations on road segments traveled by most refueling vehicles are selected. In the second stage, the focus is on refueling vehicles that are not covered in the first stage. The objective is to minimize the total detour distance of these vehicles by properly selecting additional refueling stations. The third stage is based on the results of the first and second stages. It considers the redistribution of the paths for refueling vehicles due to traffic congestion and refueling station capacities. The redistribution result is used to fine-tune the location plan generated in the first and second stages.Specifically, the first stage starts with assigning all vehicles (including refueling vehicles) to the road network. Based on the traffic assignment results, the paths of refueling vehicles are identified. From these paths, road segments that cover the highest refueling vehicle OD demand are selected for constructing refueling stations. The refueling ODs covered/served by these road segments are then removed from the original refueling vehicle OD. The regular vehicle OD and the unsatisfied refueling vehicle OD are assigned again to the road network to find additional candidate locations.At the end of the first stage, the demands of some refueling ODs may still be unmet. Therefore, additional refueling stations need to be constructed and some refueling vehicles have to detour. When selecting locations for these additional refueling stations, it is natural to consider the objective of minimizing the total detour distance.In the first and second stages, refueling vehicles are assigned to the road network in an incremental manner and refueling station locations are added gradually. In the third stage, all vehicles (including refueling vehicles) are assigned to the network simultaneously given the locations selected in the first and second stages. This most likely will result in a traffic flow pattern (e.g., choice of stations for refueling vehicles) that is different from those in the previous two stages. For instance, some previously overloaded refueling stations (based on the result of the second stage) may end up with not having enough demand. For those without enough demand in the second stage, they may not have adequate capacity to satisfy the refueling demand. For the lack of demand scenario, the corresponding refueling stations can be simply removed/unselected. For the lack of capacity scenario, additional refueling stations need to be added. To add new refueling stations, the first step is to find refueling vehicle ODs being assigned to the overloaded refueling stations. Next, road segments on these OD pairs’ paths and the refueling traffic volume on each of these segments are identified. Based on the principles discussed in the previous section, these segments’ refueling volumes are then sorted in a descending order and used to select new refueling stations. During this adjustment process, at least one new station that helps with diverting refueling demand should be selected. The new location plan may need to be fed into the second stage to be further improved. This three-stage model is summarized in Figure1.Figure 1
Flowchart for the three stages.


## 2.1. Principles of Location Design

Siting urban refueling stations should consider a comprehensive set of factors. Government agencies are mostly concerned about congestion, safety, environment, and land use factors, while refueling station operators are more concerned about demand and profit. Besides, fuel price difference at various refueling stations is also an important factor to consider. However, in China gasoline prices are decided by the government and there is almost no difference in terms of price among different gas stations. Therefore, the price competition among different refueling stations/companies is not considered in this study. In general, siting refueling stations should (a) minimize their impacts on urban traffic, (b) be on main roads to ensure profitability, (c) keep an appropriate distance between refueling stations, (d) be far from residential and other sensitive areas, and (e) be based on urban land use planning.Principle (a) is from the perspective of minimizing the influence of refueling stations on urban road traffic. For example, a vehicle with refueling need should not have to detour too far for refueling. Intuitively, this requires refueling stations to be located on or near the paths of most vehicles, so that they do not need to make a detour or just need to make a very short detour.Principle (b) considers the economic benefit of refueling stations. The number of refueling vehicles that may potentially be served by a refueling station should be taken into consideration during the location design. Ideally, this number should be within a reasonable range to ensure that a refueling station’s profit is above certain level and the refueling demand does not exceed its capacity.Principle (c) takes the service ranges of refueling stations into consideration. It is to make sure that refueling stations are not too concentrated in certain areas. Keeping an appropriate distance between refueling stations sometimes is necessary due to reasons such as safety.Principle (d) is mainly about the safety and environmental concerns. Based on this principle, siting refueling stations should consider the safety of refueling station facilities and equipment, refueling vehicles, travelers, and nearby residents. Also, refueling stations should be far away from some environmentally and socially sensitive areas.Last but certainly not least, principle (e) optimizes refueling station locations from the urban planning point of view. Specifically, urban refueling station locations should take into account land use development and future land use planning. Land use directly affects trip productions and attractions. Consequently, it has major impacts on the distribution of future refueling demand and refueling station locations.From the standpoint of model development for optimally siting refueling stations, principle (a) can be considered in the objective function, and principles (b) through (e) can be either included as constraints or used to determine candidate locations for future refueling stations.


## 2.2. Three-Stage Model

Based on the analysis in the previous section, if principle (a) is adopted as the objective, there will be two cases in the refueling station location planning. The first case is that refueling stations are located on the paths of some refueling vehicles. Hence, these vehicles do not need to detour for refueling. The other case is that refueling stations are not the paths of the remaining refueling vehicles and these vehicles have to detour for refueling. In light of these two cases, siting refueling stations can be implemented in three stages. In the first stage, candidate refueling stations on road segments traveled by most refueling vehicles are selected. In the second stage, the focus is on refueling vehicles that are not covered in the first stage. The objective is to minimize the total detour distance of these vehicles by properly selecting additional refueling stations. The third stage is based on the results of the first and second stages. It considers the redistribution of the paths for refueling vehicles due to traffic congestion and refueling station capacities. The redistribution result is used to fine-tune the location plan generated in the first and second stages.Specifically, the first stage starts with assigning all vehicles (including refueling vehicles) to the road network. Based on the traffic assignment results, the paths of refueling vehicles are identified. From these paths, road segments that cover the highest refueling vehicle OD demand are selected for constructing refueling stations. The refueling ODs covered/served by these road segments are then removed from the original refueling vehicle OD. The regular vehicle OD and the unsatisfied refueling vehicle OD are assigned again to the road network to find additional candidate locations.At the end of the first stage, the demands of some refueling ODs may still be unmet. Therefore, additional refueling stations need to be constructed and some refueling vehicles have to detour. When selecting locations for these additional refueling stations, it is natural to consider the objective of minimizing the total detour distance.In the first and second stages, refueling vehicles are assigned to the road network in an incremental manner and refueling station locations are added gradually. In the third stage, all vehicles (including refueling vehicles) are assigned to the network simultaneously given the locations selected in the first and second stages. This most likely will result in a traffic flow pattern (e.g., choice of stations for refueling vehicles) that is different from those in the previous two stages. For instance, some previously overloaded refueling stations (based on the result of the second stage) may end up with not having enough demand. For those without enough demand in the second stage, they may not have adequate capacity to satisfy the refueling demand. For the lack of demand scenario, the corresponding refueling stations can be simply removed/unselected. For the lack of capacity scenario, additional refueling stations need to be added. To add new refueling stations, the first step is to find refueling vehicle ODs being assigned to the overloaded refueling stations. Next, road segments on these OD pairs’ paths and the refueling traffic volume on each of these segments are identified. Based on the principles discussed in the previous section, these segments’ refueling volumes are then sorted in a descending order and used to select new refueling stations. During this adjustment process, at least one new station that helps with diverting refueling demand should be selected. The new location plan may need to be fed into the second stage to be further improved. This three-stage model is summarized in Figure1.Figure 1
Flowchart for the three stages.


## 3. Optimization Model for Siting Urban Refueling Stations

### 3.1. Variables and Refueling Demand

#### 3.1.1. Definitions of Variables and Parameters

Given a road networkG(N,A), N is the set of all nodes in G(N,A), and A is the set of all road sections in G(N,A), am(m=1,2,…,n)∈A. Suppose all major and minor arterial segments, cm(m=1,2,…,n)∈C, are considered as candidate locations for refueling stations. Before the models are introduced, the following variables and parameters used in the models are defined below.(1)
xai is the traffic flow of road segment ai;(2)
tai is the travel time of road segment ai;(3)
ζ is a binary variable. 1 stands for constructing a station and 0 for not;(4)
M is the minimum refueling demand for constructing a refueling station;(5)
S is the service area of the road network;(6)
p(p=1,2,…,n) is the set of candidate locations for constructing refueling stations;(7)
bp is the cost of constructing a refueling station at candidate location p, and Bmax⁡ stands for the maximum average cost;(8)
tpp′ is the travel time between refueling stations p and p′ and should be in the range of tdmin⁡,tdmax⁡;(9)
ρ stands for the density of refueling stations and its range should be ρmin⁡,ρmax⁡;(10)
Qij is the travel demand between origin i and destination j and qij is the corresponding refueling demand;(11)
fkij is the traffic volume on path k between origin i and destination j.


#### 3.1.2. Vehicle Refueling Demand

Except for the vehicle refueling demandqij, most of the variables and parameters defined in the previous section either are given or can be readily calculated. The vehicle refueling demand however is closely related to the vehicle travel demand (Qij) and should be determined carefully.For different OD pairs, even the same vehicle’s average refueling frequencies can be different. This is mainly due to two factors. First, different OD pairs have different average travel distances, which naturally result in different fuel consumptions. Second, vehicles from different OD pairs may experience different levels of congestion, which will also affect their fuel consumptions. On a particular day, travel demand of an OD pair is separated into demand with refueling needs and demand without refueling needs. The demand between origini and destination j with refueling needs is estimated using (1). This study assumes that travelers do not have a special preference for certain refueling stations. This is a reasonable assumption for planning purpose, since in China gasoline prices are the same at different stations and it is hard to tell which refueling station(s) will be more attractive in the future. Consider(1)qij=Dijuij=Qij·Av·uijCgas,∀i,j.In (1), Dij is the demand function between i and j. This monotonic function depends on the travel distance between i and j and has an upper limit; uij is the average travel distance between i and j in kilometers; A(v) is the average vehicle fuel consumption measured in liter/100 km. It is a function of the actual driving speed v in km/hr, which can be estimated using the method in [11]; and Cgas is the average fuel tank size measured in liters.


### 3.2. First-Stage Model Based on FCLM

The first-stage model is shown in (2) below and illustrated in Figure 2. It is built upon the Flow-Capturing Location-Allocation Model (FCLM) proposed by Hodgson [8] and the idea for the first stage discussed previously. One has(2)Maximum∑ai∈Aqaiζais.t.qai≥M(ai∈P)ζai=0ai∉P1ai∈P.Figure 2
First-stage location model based on FCLM.


### 3.3. Location Models for the Second Stage

After completing the first stage, the refueling demands of some vehicles are still unmet. When adding additional refueling stations (also called public refueling stations in this study) to address this issue, our objective is to minimize the total detour distance for those ODs with unmet refueling demands. Specifically, a bilevel location model is proposed that reflects the concerns of both government agencies and travelers. In the upper-level model, the total detour time is used as the objective function. Constraints considered include travel time between refueling stations, average construction cost, and refueling station density. The lower-level model assigns the unmet refueling ODs to the network using Stochastic User Equilibrium (SUE). This assignment is based on the resultant road network of the first stage. This bilevel second-stage model is also shown in Figure3. The third stage does not have a specific model. It is simply to fine-tune the result of the second-stage model.Figure 3
Location models of the second stage.Step 1 (upper-level location model).
Consider(3)min⁡z(x)=∑i∑jqij(tij′-t-ij)(4)s.t.tdmin⁡<tpp′<tdmax⁡p,p′∈P(5)∑p∈Pζpbp∑p∈Pζp≤Bmax⁡(6)ρmin⁡≤∑p∈PζpS≤ρmax⁡.In (3), tij′ is the average travel time of refueling vehicles for OD pair (i,j), t¯ij stands for the average travel time of vehicles that do not need refueling, and the remaining variables have all been defined previously.Step 2 (lower-level location model).
Consider(7)min⁡Zx,q=∑ai∈A∫0xaitaiωdω-∑ip∫0qipDip-1ωdω-∑pj∫0qpjDpj-1ωdω-∑ij∫0qij′Dij-1ωdωs.t.∑afkij=qij′,∑afkip=qip,∑afkpj=qpj,∀i,j,pfkij,fkip,fkpj≥0,∀k,i,j,pqij′,qip,qpj≥0,∀i,j,p,p∈P,(8)tipopt+tpoptj<tiP′+tP′s,popt∈P,P′={p∣p=P-popt}.
In constraint (8), popt represents the refueling station closest to an OD pair among all candidate refueling stations P and P′ is the set of refueling stations excluding popt.


## 3.1. Variables and Refueling Demand

### 3.1.1. Definitions of Variables and Parameters

Given a road networkG(N,A), N is the set of all nodes in G(N,A), and A is the set of all road sections in G(N,A), am(m=1,2,…,n)∈A. Suppose all major and minor arterial segments, cm(m=1,2,…,n)∈C, are considered as candidate locations for refueling stations. Before the models are introduced, the following variables and parameters used in the models are defined below.(1)
xai is the traffic flow of road segment ai;(2)
tai is the travel time of road segment ai;(3)
ζ is a binary variable. 1 stands for constructing a station and 0 for not;(4)
M is the minimum refueling demand for constructing a refueling station;(5)
S is the service area of the road network;(6)
p(p=1,2,…,n) is the set of candidate locations for constructing refueling stations;(7)
bp is the cost of constructing a refueling station at candidate location p, and Bmax⁡ stands for the maximum average cost;(8)
tpp′ is the travel time between refueling stations p and p′ and should be in the range of tdmin⁡,tdmax⁡;(9)
ρ stands for the density of refueling stations and its range should be ρmin⁡,ρmax⁡;(10)
Qij is the travel demand between origin i and destination j and qij is the corresponding refueling demand;(11)
fkij is the traffic volume on path k between origin i and destination j.


### 3.1.2. Vehicle Refueling Demand

Except for the vehicle refueling demandqij, most of the variables and parameters defined in the previous section either are given or can be readily calculated. The vehicle refueling demand however is closely related to the vehicle travel demand (Qij) and should be determined carefully.For different OD pairs, even the same vehicle’s average refueling frequencies can be different. This is mainly due to two factors. First, different OD pairs have different average travel distances, which naturally result in different fuel consumptions. Second, vehicles from different OD pairs may experience different levels of congestion, which will also affect their fuel consumptions. On a particular day, travel demand of an OD pair is separated into demand with refueling needs and demand without refueling needs. The demand between origini and destination j with refueling needs is estimated using (1). This study assumes that travelers do not have a special preference for certain refueling stations. This is a reasonable assumption for planning purpose, since in China gasoline prices are the same at different stations and it is hard to tell which refueling station(s) will be more attractive in the future. Consider(1)qij=Dijuij=Qij·Av·uijCgas,∀i,j.In (1), Dij is the demand function between i and j. This monotonic function depends on the travel distance between i and j and has an upper limit; uij is the average travel distance between i and j in kilometers; A(v) is the average vehicle fuel consumption measured in liter/100 km. It is a function of the actual driving speed v in km/hr, which can be estimated using the method in [11]; and Cgas is the average fuel tank size measured in liters.


## 3.1.1. Definitions of Variables and Parameters

Given a road networkG(N,A), N is the set of all nodes in G(N,A), and A is the set of all road sections in G(N,A), am(m=1,2,…,n)∈A. Suppose all major and minor arterial segments, cm(m=1,2,…,n)∈C, are considered as candidate locations for refueling stations. Before the models are introduced, the following variables and parameters used in the models are defined below.(1)
xai is the traffic flow of road segment ai;(2)
tai is the travel time of road segment ai;(3)
ζ is a binary variable. 1 stands for constructing a station and 0 for not;(4)
M is the minimum refueling demand for constructing a refueling station;(5)
S is the service area of the road network;(6)
p(p=1,2,…,n) is the set of candidate locations for constructing refueling stations;(7)
bp is the cost of constructing a refueling station at candidate location p, and Bmax⁡ stands for the maximum average cost;(8)
tpp′ is the travel time between refueling stations p and p′ and should be in the range of tdmin⁡,tdmax⁡;(9)
ρ stands for the density of refueling stations and its range should be ρmin⁡,ρmax⁡;(10)
Qij is the travel demand between origin i and destination j and qij is the corresponding refueling demand;(11)
fkij is the traffic volume on path k between origin i and destination j.


## 3.1.2. Vehicle Refueling Demand

Except for the vehicle refueling demandqij, most of the variables and parameters defined in the previous section either are given or can be readily calculated. The vehicle refueling demand however is closely related to the vehicle travel demand (Qij) and should be determined carefully.For different OD pairs, even the same vehicle’s average refueling frequencies can be different. This is mainly due to two factors. First, different OD pairs have different average travel distances, which naturally result in different fuel consumptions. Second, vehicles from different OD pairs may experience different levels of congestion, which will also affect their fuel consumptions. On a particular day, travel demand of an OD pair is separated into demand with refueling needs and demand without refueling needs. The demand between origini and destination j with refueling needs is estimated using (1). This study assumes that travelers do not have a special preference for certain refueling stations. This is a reasonable assumption for planning purpose, since in China gasoline prices are the same at different stations and it is hard to tell which refueling station(s) will be more attractive in the future. Consider(1)qij=Dijuij=Qij·Av·uijCgas,∀i,j.In (1), Dij is the demand function between i and j. This monotonic function depends on the travel distance between i and j and has an upper limit; uij is the average travel distance between i and j in kilometers; A(v) is the average vehicle fuel consumption measured in liter/100 km. It is a function of the actual driving speed v in km/hr, which can be estimated using the method in [11]; and Cgas is the average fuel tank size measured in liters.


## 3.2. First-Stage Model Based on FCLM

The first-stage model is shown in (2) below and illustrated in Figure 2. It is built upon the Flow-Capturing Location-Allocation Model (FCLM) proposed by Hodgson [8] and the idea for the first stage discussed previously. One has(2)Maximum∑ai∈Aqaiζais.t.qai≥M(ai∈P)ζai=0ai∉P1ai∈P.Figure 2
First-stage location model based on FCLM.


## 3.3. Location Models for the Second Stage

After completing the first stage, the refueling demands of some vehicles are still unmet. When adding additional refueling stations (also called public refueling stations in this study) to address this issue, our objective is to minimize the total detour distance for those ODs with unmet refueling demands. Specifically, a bilevel location model is proposed that reflects the concerns of both government agencies and travelers. In the upper-level model, the total detour time is used as the objective function. Constraints considered include travel time between refueling stations, average construction cost, and refueling station density. The lower-level model assigns the unmet refueling ODs to the network using Stochastic User Equilibrium (SUE). This assignment is based on the resultant road network of the first stage. This bilevel second-stage model is also shown in Figure3. The third stage does not have a specific model. It is simply to fine-tune the result of the second-stage model.Figure 3
Location models of the second stage.Step 1 (upper-level location model).
Consider(3)min⁡z(x)=∑i∑jqij(tij′-t-ij)(4)s.t.tdmin⁡<tpp′<tdmax⁡p,p′∈P(5)∑p∈Pζpbp∑p∈Pζp≤Bmax⁡(6)ρmin⁡≤∑p∈PζpS≤ρmax⁡.In (3), tij′ is the average travel time of refueling vehicles for OD pair (i,j), t¯ij stands for the average travel time of vehicles that do not need refueling, and the remaining variables have all been defined previously.Step 2 (lower-level location model).
Consider(7)min⁡Zx,q=∑ai∈A∫0xaitaiωdω-∑ip∫0qipDip-1ωdω-∑pj∫0qpjDpj-1ωdω-∑ij∫0qij′Dij-1ωdωs.t.∑afkij=qij′,∑afkip=qip,∑afkpj=qpj,∀i,j,pfkij,fkip,fkpj≥0,∀k,i,j,pqij′,qip,qpj≥0,∀i,j,p,p∈P,(8)tipopt+tpoptj<tiP′+tP′s,popt∈P,P′={p∣p=P-popt}.
In constraint (8), popt represents the refueling station closest to an OD pair among all candidate refueling stations P and P′ is the set of refueling stations excluding popt.


## 4. Algorithms

### 4.1. Overview

This study relies on two core models: the first-stage model based on the FCLM and the second-stage model for adding and fine-tuning refueling stations. The first-stage model is relatively simple. Thus, this paper will briefly discuss the algorithm for the first-stage model and will focus on developing the algorithm for solving the second-stage model.The traffic assignment in the first-stage model is based on SUE assignment. From the result of each assignment, a road segment with the highest refueling traffic volume is selected. This segment is selected for constructing a refueling station if its refueling volume is greater than a prespecified value. The corresponding OD will be excluded from further considerations. This process repeats until no qualified segments can be found for constructing refueling stations.The upper-level of the second-stage model is a nonlinear program and the lower-level is a SUE problem. Heuristic algorithms are very popular for solving nonlinear programming problems. For instance, genetic algorithm (GA) is often used for finding global optimal solutions, and simulated annealing (SA) has good local search capability. This paper proposes an integrated genetic simulated annealing (GSA) algorithm to solve the upper-level model, generating feasible solutions for the lower-level model. The lower-level model basically performs SUE traffic assignments and provides the upper-level model withxai,tai. In this study, both GA and GSA algorithms are considered and their results are compared.


### 4.2. Genetic Simulated Annealing (GSA) Algorithm

GA was motivated by Darwin’s theory of natural selection, while SA was inspired by the process of heating and cooling of materials in a controlled setting. By integrating these two algorithms, we hope to further improve GA by enhancing its local search with SA.The proposed algorithm consists of an inner loop and an outer loop, which are based on GA and SA, respectively. The inner loop generates binary vectors (e.g., chromosomes) representing various refueling station location solutions and evaluates their fitness. The algorithm then switches to the outer loop and applies SA to further improve the fitness of these chromosomes locally. Upon finishing the local search, the algorithm returns to the inner loop and applies crossover and mutation operators to avoid being stuck at local minimum solutions. The crossover and mutation operators, initial temperature of annealing, and the annealing rate are all important factors that may affect the algorithm’s performance and should be carefully considered. The following steps describe the proposed algorithm in more detail.Step 1.
Since this is essentially a facility location problem, a binary vector of lengthn is used to represent each refueling station location solution. Each vector element (or each gene) represents a candidate location. 1 means the corresponding location is selected to construct a refueling station, and 0 means it is not selected. Pseudorandom number generator (PRNG) is used to generate initial refueling station location solution set Pop(k), k=0. The solution set size is Popsize, meaning it includes Popsize solutions. Also, the initial temperature is set to T0; final temperature is Te; cooling rate is TR; and the maximum number of iterations for the outer loop is set to Oitermax.Step 2.
This step calculates the fitness values for individual solutions/chromosomes. Based on the results, the fittest chromosomei and its fitness value F(i) are determined. This chromosome is the one that has the best objective function value (i.e., fitness value) among all chromosomes in the current population. Therefore, we set Sbest=i and Fmin⁡=F(i). Also, fi(Tk) is used to denote the target value for chromosome i under temperature Tk, where i∈Pop(k) is a refueling station location solution in current solution set.Step 3.
The entire algorithm will stop upon reaching the final temperatureTe and it will output the best chromosome Sbest and the optimal fitness value Fmin⁡. Otherwise it will randomly select a state j∈N(i) from the neighborhood of each chromosome i∈Pop(k) in the current solution set. The new state j will be either accepted or rejected according to the acceptance probability pij=min⁡1,exp⁡-fj(Tk)-fi(Tk)/Tk. fi(Tk) and fj(Tk) are the target values of states i and j, respectively. This step takes Popsize iterations and generates a new solution set Newpop1 for refueling stations.Step 4.
This step calculates the fitness valuesfi(Tk)=exp⁡-fi(Tk)-fmin⁡/Tk of all solutions in Newpop1, where fmin⁡ is the minimum fitness value among all solutions. At the end, Popsize solutions/chromosomes from Newpop1 are selected with probabilities inversely proportional to their fitness values to form a new solution set Newpop2.Step 5.
Crossover and mutation are applied toNewpop2, which generate crosspop and mutpop, respectively. The crossover and mutation processes are described below.(i) Crossover Operator. For crossover, two parent chromosomes (representing two refueling station location solutions) A1 and A2 are selected to form a new chromosome. Each gene of the new chromosome is selected from the corresponding gene in either A1 or A2 based on probabilities inversely proportional to their fitness values. Suppose the fitness values of A1 and A2 are f1 and f2, respectively; the crossover process is carried out as follows.(a)
Calculate the probabilities for genes to be selected fromA1 and A2(9)ForA1:P1=f1f1+f2,ForA2:P2=f2f1+f2.(b)
Generate a sequenceαi(i=1,2,…,k) with each element being a random number within [0,1]. Compare each element in this sequence with P1 or P2 to determine whether a gene from A1 or A2 should be selected. For example, if P1<P2, for αi<P1, the ith gene from A1 will be selected for the ith gene in the new solution; otherwise, the ith gene from A2 will be selected.(ii) Mutation Operator. Assume a parent solution selected for mutation is X=x1x2x3⋯xn. xk(k=1,2,3…,n) is the mutation point and it is defined on the domain of [Ak,Bk]. A new solution is generated based on mutation as follows:(10)xk′=xk+f(t,Bk-xk)rand(0,1)=0xk+f(t,xk-Ak)rand(0,1)=1f(t,y)=y(1-rβ(1-t/T)),where r is a random real number uniformly distributed between 0 and 1, t is the current number of iterations, T sets the maximum number of iterations, and β is a prespecified parameter and can be adjusted based on t. It is set to a large value at the beginning and reduces gradually as t increases. rand(0,1) randomly selects 0 or 1 as the output. Finally, the mutated chromosomes and their parent chromosomes X are both included in the solution set.The mutation operator is applied to chromosomes with the highest and lowest fitness values. Applying it to the fittest chromosome is for accelerating the speed for local search, while applying it to the least-fit chromosome is to increase the diversity of the population and to avoid premature convergences.Step 6.
SetPop(k+1)=mutpop and find chromosome i that minimizes fi(Tk). If fi(Tk)≤fmin⁡, set Sbest=i and Fmin⁡=fi(Tk). Finally, set Tk=T0exp⁡-C·k1/N and k=k+1 and go to Step 3.


## 4.1. Overview

This study relies on two core models: the first-stage model based on the FCLM and the second-stage model for adding and fine-tuning refueling stations. The first-stage model is relatively simple. Thus, this paper will briefly discuss the algorithm for the first-stage model and will focus on developing the algorithm for solving the second-stage model.The traffic assignment in the first-stage model is based on SUE assignment. From the result of each assignment, a road segment with the highest refueling traffic volume is selected. This segment is selected for constructing a refueling station if its refueling volume is greater than a prespecified value. The corresponding OD will be excluded from further considerations. This process repeats until no qualified segments can be found for constructing refueling stations.The upper-level of the second-stage model is a nonlinear program and the lower-level is a SUE problem. Heuristic algorithms are very popular for solving nonlinear programming problems. For instance, genetic algorithm (GA) is often used for finding global optimal solutions, and simulated annealing (SA) has good local search capability. This paper proposes an integrated genetic simulated annealing (GSA) algorithm to solve the upper-level model, generating feasible solutions for the lower-level model. The lower-level model basically performs SUE traffic assignments and provides the upper-level model withxai,tai. In this study, both GA and GSA algorithms are considered and their results are compared.


## 4.2. Genetic Simulated Annealing (GSA) Algorithm

GA was motivated by Darwin’s theory of natural selection, while SA was inspired by the process of heating and cooling of materials in a controlled setting. By integrating these two algorithms, we hope to further improve GA by enhancing its local search with SA.The proposed algorithm consists of an inner loop and an outer loop, which are based on GA and SA, respectively. The inner loop generates binary vectors (e.g., chromosomes) representing various refueling station location solutions and evaluates their fitness. The algorithm then switches to the outer loop and applies SA to further improve the fitness of these chromosomes locally. Upon finishing the local search, the algorithm returns to the inner loop and applies crossover and mutation operators to avoid being stuck at local minimum solutions. The crossover and mutation operators, initial temperature of annealing, and the annealing rate are all important factors that may affect the algorithm’s performance and should be carefully considered. The following steps describe the proposed algorithm in more detail.Step 1.
Since this is essentially a facility location problem, a binary vector of lengthn is used to represent each refueling station location solution. Each vector element (or each gene) represents a candidate location. 1 means the corresponding location is selected to construct a refueling station, and 0 means it is not selected. Pseudorandom number generator (PRNG) is used to generate initial refueling station location solution set Pop(k), k=0. The solution set size is Popsize, meaning it includes Popsize solutions. Also, the initial temperature is set to T0; final temperature is Te; cooling rate is TR; and the maximum number of iterations for the outer loop is set to Oitermax.Step 2.
This step calculates the fitness values for individual solutions/chromosomes. Based on the results, the fittest chromosomei and its fitness value F(i) are determined. This chromosome is the one that has the best objective function value (i.e., fitness value) among all chromosomes in the current population. Therefore, we set Sbest=i and Fmin⁡=F(i). Also, fi(Tk) is used to denote the target value for chromosome i under temperature Tk, where i∈Pop(k) is a refueling station location solution in current solution set.Step 3.
The entire algorithm will stop upon reaching the final temperatureTe and it will output the best chromosome Sbest and the optimal fitness value Fmin⁡. Otherwise it will randomly select a state j∈N(i) from the neighborhood of each chromosome i∈Pop(k) in the current solution set. The new state j will be either accepted or rejected according to the acceptance probability pij=min⁡1,exp⁡-fj(Tk)-fi(Tk)/Tk. fi(Tk) and fj(Tk) are the target values of states i and j, respectively. This step takes Popsize iterations and generates a new solution set Newpop1 for refueling stations.Step 4.
This step calculates the fitness valuesfi(Tk)=exp⁡-fi(Tk)-fmin⁡/Tk of all solutions in Newpop1, where fmin⁡ is the minimum fitness value among all solutions. At the end, Popsize solutions/chromosomes from Newpop1 are selected with probabilities inversely proportional to their fitness values to form a new solution set Newpop2.Step 5.
Crossover and mutation are applied toNewpop2, which generate crosspop and mutpop, respectively. The crossover and mutation processes are described below.(i) Crossover Operator. For crossover, two parent chromosomes (representing two refueling station location solutions) A1 and A2 are selected to form a new chromosome. Each gene of the new chromosome is selected from the corresponding gene in either A1 or A2 based on probabilities inversely proportional to their fitness values. Suppose the fitness values of A1 and A2 are f1 and f2, respectively; the crossover process is carried out as follows.(a)
Calculate the probabilities for genes to be selected fromA1 and A2(9)ForA1:P1=f1f1+f2,ForA2:P2=f2f1+f2.(b)
Generate a sequenceαi(i=1,2,…,k) with each element being a random number within [0,1]. Compare each element in this sequence with P1 or P2 to determine whether a gene from A1 or A2 should be selected. For example, if P1<P2, for αi<P1, the ith gene from A1 will be selected for the ith gene in the new solution; otherwise, the ith gene from A2 will be selected.(ii) Mutation Operator. Assume a parent solution selected for mutation is X=x1x2x3⋯xn. xk(k=1,2,3…,n) is the mutation point and it is defined on the domain of [Ak,Bk]. A new solution is generated based on mutation as follows:(10)xk′=xk+f(t,Bk-xk)rand(0,1)=0xk+f(t,xk-Ak)rand(0,1)=1f(t,y)=y(1-rβ(1-t/T)),where r is a random real number uniformly distributed between 0 and 1, t is the current number of iterations, T sets the maximum number of iterations, and β is a prespecified parameter and can be adjusted based on t. It is set to a large value at the beginning and reduces gradually as t increases. rand(0,1) randomly selects 0 or 1 as the output. Finally, the mutated chromosomes and their parent chromosomes X are both included in the solution set.The mutation operator is applied to chromosomes with the highest and lowest fitness values. Applying it to the fittest chromosome is for accelerating the speed for local search, while applying it to the least-fit chromosome is to increase the diversity of the population and to avoid premature convergences.Step 6.
SetPop(k+1)=mutpop and find chromosome i that minimizes fi(Tk). If fi(Tk)≤fmin⁡, set Sbest=i and Fmin⁡=fi(Tk). Finally, set Tk=T0exp⁡-C·k1/N and k=k+1 and go to Step 3.


## 5. Case Study

### 5.1. Description of Case Study Problem

#### 5.1.1. Road Network and Parameters

To validate the proposed model and algorithm, a case study was conducted using the urban road network in Figures4 and 5. This network consisted of 128 nodes and 234 segments. In Figure 5, green lines represent major arterials, yellow lines are for minor arterials, and blue lines denote secondary roads. The numbers next to each segment are its capacity and free-flow travel time. The entire study area was about 90 km2. In order to mitigate traffic congestion in downtown, no refueling stations were allowed within the red circle in Figure 5. Also, it was assumed that no refueling stations can be built on secondary roads.Figure 4
Road network with node and segment numbers.Figure 5
Segment capacity and travel time.Additionally, the following assumptions were made: the maximum average construction cost for a refueling station is 3 million Chinese yuan; the density of refueling stations is between 0.1 and 0.3 stations/km2; the minimum distance between two refueling stations is 4 km; and each station should serve at least 300 vehicles in order to be profitable and no more than 500 vehicles due to capacity constraint.


#### 5.1.2. Vehicle Refueling Demand

It was assumed that each node is the centroid of a Traffic Analysis Zone (TAZ). The travel demand between two TAZs was estimated byfij=PiPj/dij1.5, where Pi and Pj stand for the trip production of nodes i and the trip attraction of node j, respectively. Both Pi and Pj were assumed to be 5000 trips/day for all nodes. dij is the distance between two nodes.Considering peak-hour factor and mode split, the peak-hour travel demand between two TAZs was estimated to be 38,537 pcu/h. By applying (1), the corresponding refueling demand was estimated to be 6,775 pcu/h. Due to page limit, detailed steps for the above calculations are not included in this paper. Also, based on the refueling station density and capacity constraints, it was assumed that between 14 and 23 refueling stations should be constructed.


### 5.2. Modeling Results

#### 5.2.1. Results of the First-Stage Model

TransCAD 5.0 was used to solve the first-stage model. As shown in Table1, the first-stage model finished after 8 iterations. During this stage, the capacity of refueling stations was not considered. In other words, we assumed that a refueling station can satisfy all passing vehicles’ refueling demand.Table 1
Result of the first-stage model.
Number of iterations
Segment number
Segment direction
Segment flow (pcu/h)
Refueling demand (pcu/h)
1
125
112 → 121
4637
788
2
162
133 → 132
2902
610
3
111
105 → 106
3966
552
4
103
101 → 102
2341
421
5
56
75 → 74
2612
418
6
154
128 → 129
1983
333
7
98
98 → 107
1853
308
8
194
154 → 202
2141
278Table1 shows that the refueling demand of Segment 194 is 278 pcu/h, which is less than 300 pcu/h and does not meet the minimum demand requirement for constructing a refueling station. Therefore, the first-stage model stopped after the 8th iteration. The result in Table 1 suggests that there is no obvious correlation between segment flow and segment refueling demand. The result of the first-stage model is also shown in Figure 6.Figure 6
Locations selected for refueling station during the first stage.The stations selected during the first stage can directly satisfy a refueling demand of 3,430 pcu/h, which is about half of the total demand, meaning that approximately 50% of the refueling vehicles do not need to detour. In the following section, the unmet refueling demand will be discussed.


#### 5.2.2. Result of the Second-Stage Model

(i) Model Initialization. As described previously, the second-stage model first generated Popsize=100 chromosomes and each chromosome consisted of 151 genes. Each gene is a binary variable, representing whether the corresponding segment is selected for constructing a refueling station.(ii) Parameter Setting. To prevent premature convergence, the following parameters were used: initial temperature in the annealing process T0=1000, terminating temperature Te=0, cooling speed TR=0.95, maximum number of iterations for each temperature was 100, crossover rate was 0.40~0.99, mutation rate was 0.0001~0.1, crossover probability was 0.9, mutation probability was 0.04, and the maximum number of iterations Oitermax=2000.The constraints in the second-stage model were considered by introducing a penalty. All solutions not meeting constraints (4)~(6) were punished by a multiplication of 10. The elitist strategy was adopted to select chromosomes to be included in the next generation, meaning chromosomes were selected with probabilities inversely proportional to their fitness values.(iii) Modeling Result. VB and GISDK in TransCAD were used to implement the second-stage algorithm. Among them, VB was used for coding the GSA algorithm, and GISDK was used for implementing the SUE assignment algorithm. After 876 cooling iterations, the temperature was 3.06E-17 and an optimal solution was obtained. The average run time was 24 minutes and 46 seconds. The obtained shortest total detour time for refueling vehicles was 6,688 minutes.Table2 lists the shortest total detour time and the detour times for different chromosomes at each iteration step. The iterative process is also illustrated in Figure 7, in which the horizontal axis stands for temperature and the vertical axis stands for fitness value. Figure 7 shows that the optimal fitness curve goes up and down at the beginning and tends to stabilize and level out after 876 iterations. This optimal fitness curve is not monotonically decreasing because the simulated annealing module sometimes accepts chromosomes with suboptimal (i.e., worse) fitness values with a small probability. The newly added stations in the second stage are shown in Figure 8. These stations provide service to the remaining 49.37% refueling vehicles.Table 2
Result of the second-stage model.
Number of iterations
Shortest total detour time (min)
Detour times of other solutions (min)
1
26752
27078
28244
27927
28470
⋯
28154
28637
2
26066
27129
28141
28923
28228
⋯
27130
26639
3
25045
28991
26479
27130
27910
⋯
28192
28087
4
25384
28136
28467
26486
27130
⋯
25631
26177
5
25218
25596
25922
26140
26551
⋯
26841
25676
6
23806
24469
25921
26685
24534
⋯
25170
25191
7
22362
22796
22344
22867
22845
⋯
25917
22948
8
22269
23082
23505
23600
22796
⋯
22796
25692
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
872
6991
8729
9455
7086
7138
⋯
7283
7136
873
6782
8007
8795
6926
6792
⋯
8869
6936
874
6688
6688
6688
8360
8024
⋯
6688
6688
875
6688
6688
6688
8125
7727
⋯
6688
6688
876
6688
6688
6688
6688
6688
⋯
6688
6688Figure 7
Iterative process of the second-stage model.Figure 8
Refueling stations selected during the second stage.(iv) Comparison of Algorithms. In this section, GA and the GSA algorithms are compared. The result is reported in Table 3. As can be seen, on average it takes the GA more iterations (1226) than the GSA algorithm to obtain the same solution. Also, the GA requires longer CPU time. For a larger and more complex road network, the advantages of the GSA algorithm could be more significant.Table 3
Algorithm performance comparison.
Algorithm
Average number of iterations
Average CPU time
Genetic simulated annealing algorithm
876
24 min 46 s
Genetic algorithm
1226
34 min 39 s(v) Sensitivity Analysis of Model Parameters. To assess the sensitivity of the GSA algorithm to model parameters, we experimented with various population sizes, crossover rates, and mutation rates. The modeling results are summarized and compared in Table 4. It seems that these parameters do have an important impact on model performance and should be decided carefully based on test runs.Table 4
Algorithm parameter sensitivity analysis results.
Parameter
Population size
Crossover rate
Mutation rate
Iterations
Calculating time
Bypass time
Optimization scheme
100
0.9
0.04
876
24′46′′
6688′
Test scheme 1
90
0.9
0.04
704
20′15′′
6847′
Test scheme 2
110
0.9
0.04
991
32′51′′
6688′
Test scheme 3
100
0.85
0.04
657
18′34′′
7351′
Test scheme 4
100
0.95
0.04
1324
40′42′′
6688′
Test scheme 5
100
0.9
0.03
754
21′23′′
7023′
Test scheme 6
100
0.9
0.05
1521
51′34′′
6688′


#### 5.2.3. Result of the Third Stage

The stations selected in the previous two stages are shown in Figure9, since these stations were selected in an incremental manner. If we assign the refueling traffic to the network considering all stations in Figure 9 simultaneously, the refueling vehicles will be redistributed among these stations as shown in Table 5, which shows that three refueling stations now have demand greater than their capacities (i.e., 500 pcu/h) and another three stations are unprofitable (demand less than 300 pcu/h).Table 5
Refueling traffic volumes of each refueling station based on the location plan generated by the second-stage model.
Station number
Refueling demand
Number 2
605
Number 5
577
Number 6
535
Number 1
476
Number 3
468
Number 4
417
Number 7
378
Number 8
347
Number 14
332
Number 9
329
Number 10
327
Number 17
321
Number 11
318
Number 16
313
Number 12
316
Number 13
277
Number 15
241
Number 18
198Figure 9
Preliminary location plan for refueling stations.The purpose of the third stage was to adjust the solutions generated in the previous two stages according to the method described in Figure1. The result is shown in Table 6. As we adjusted/relocated refueling stations, the standard deviation (STDEV) of the numbers of refueling vehicles served by different refueling stations decreased. However, the total detour time increased. After 3 adjustments, the numbers of refueling vehicles served by each refueling station all fell in the range of [300,500] pcu/h as shown in Table 7. If we continue to relocate refueling stations, the total detour time of refueling vehicles most likely will keep increasing. Therefore, the third stage stopped here, since all selected refueling stations had a reasonable demand and the total detour time was still low. The final solution to this case study is shown in Figure 10.Table 6
Result of the third stage.
Number of adjustments
Adjusting station
Adjusting method
STDEV
Total detour time/min
0
—
—
114.7
6688
1
Number 15
link157 → link185
94.8
7058
2
Number 18
link165 → link168
71.3
7367
3
Number 13
Link34 → link13
62.1
7520Table 7
Final refueling traffic volumes at each station.
Station number
Refueling volume
Number 1
481
Number 2
473
Number 3
473
Number 5
450
Number 6
443
Number 4
421
Number 7
382
Number 8
364
Number 9
352
Number 14
342
Number 15
334
Number 10
333
Number 17
326
Number 13
324
Number 11
321
Number 18
320
Number 12
319
Number 16
317Figure 10
Final solution to the case study.The refueling stations selected during the first stage were directly on the paths of refueling vehicles, and the average detour time of these refueling vehicles was 0 minutes; the average detour time of detouring refueling vehicles served by stations selected during the second stage was approximately 2 minutes; the average detour time of new detouring refueling vehicles due to the adjustments in the third stage was 2.5 minutes; the average detour time of all detouring refueling vehicles was 2.1 minutes; and the average detouring time of all refueling vehicles was 1.1 minutes.


## 5.1. Description of Case Study Problem

### 5.1.1. Road Network and Parameters

To validate the proposed model and algorithm, a case study was conducted using the urban road network in Figures4 and 5. This network consisted of 128 nodes and 234 segments. In Figure 5, green lines represent major arterials, yellow lines are for minor arterials, and blue lines denote secondary roads. The numbers next to each segment are its capacity and free-flow travel time. The entire study area was about 90 km2. In order to mitigate traffic congestion in downtown, no refueling stations were allowed within the red circle in Figure 5. Also, it was assumed that no refueling stations can be built on secondary roads.Figure 4
Road network with node and segment numbers.Figure 5
Segment capacity and travel time.Additionally, the following assumptions were made: the maximum average construction cost for a refueling station is 3 million Chinese yuan; the density of refueling stations is between 0.1 and 0.3 stations/km2; the minimum distance between two refueling stations is 4 km; and each station should serve at least 300 vehicles in order to be profitable and no more than 500 vehicles due to capacity constraint.


### 5.1.2. Vehicle Refueling Demand

It was assumed that each node is the centroid of a Traffic Analysis Zone (TAZ). The travel demand between two TAZs was estimated byfij=PiPj/dij1.5, where Pi and Pj stand for the trip production of nodes i and the trip attraction of node j, respectively. Both Pi and Pj were assumed to be 5000 trips/day for all nodes. dij is the distance between two nodes.Considering peak-hour factor and mode split, the peak-hour travel demand between two TAZs was estimated to be 38,537 pcu/h. By applying (1), the corresponding refueling demand was estimated to be 6,775 pcu/h. Due to page limit, detailed steps for the above calculations are not included in this paper. Also, based on the refueling station density and capacity constraints, it was assumed that between 14 and 23 refueling stations should be constructed.


## 5.1.1. Road Network and Parameters

To validate the proposed model and algorithm, a case study was conducted using the urban road network in Figures4 and 5. This network consisted of 128 nodes and 234 segments. In Figure 5, green lines represent major arterials, yellow lines are for minor arterials, and blue lines denote secondary roads. The numbers next to each segment are its capacity and free-flow travel time. The entire study area was about 90 km2. In order to mitigate traffic congestion in downtown, no refueling stations were allowed within the red circle in Figure 5. Also, it was assumed that no refueling stations can be built on secondary roads.Figure 4
Road network with node and segment numbers.Figure 5
Segment capacity and travel time.Additionally, the following assumptions were made: the maximum average construction cost for a refueling station is 3 million Chinese yuan; the density of refueling stations is between 0.1 and 0.3 stations/km2; the minimum distance between two refueling stations is 4 km; and each station should serve at least 300 vehicles in order to be profitable and no more than 500 vehicles due to capacity constraint.


## 5.1.2. Vehicle Refueling Demand

It was assumed that each node is the centroid of a Traffic Analysis Zone (TAZ). The travel demand between two TAZs was estimated byfij=PiPj/dij1.5, where Pi and Pj stand for the trip production of nodes i and the trip attraction of node j, respectively. Both Pi and Pj were assumed to be 5000 trips/day for all nodes. dij is the distance between two nodes.Considering peak-hour factor and mode split, the peak-hour travel demand between two TAZs was estimated to be 38,537 pcu/h. By applying (1), the corresponding refueling demand was estimated to be 6,775 pcu/h. Due to page limit, detailed steps for the above calculations are not included in this paper. Also, based on the refueling station density and capacity constraints, it was assumed that between 14 and 23 refueling stations should be constructed.


## 5.2. Modeling Results

### 5.2.1. Results of the First-Stage Model

TransCAD 5.0 was used to solve the first-stage model. As shown in Table1, the first-stage model finished after 8 iterations. During this stage, the capacity of refueling stations was not considered. In other words, we assumed that a refueling station can satisfy all passing vehicles’ refueling demand.Table 1
Result of the first-stage model.
Number of iterations
Segment number
Segment direction
Segment flow (pcu/h)
Refueling demand (pcu/h)
1
125
112 → 121
4637
788
2
162
133 → 132
2902
610
3
111
105 → 106
3966
552
4
103
101 → 102
2341
421
5
56
75 → 74
2612
418
6
154
128 → 129
1983
333
7
98
98 → 107
1853
308
8
194
154 → 202
2141
278Table1 shows that the refueling demand of Segment 194 is 278 pcu/h, which is less than 300 pcu/h and does not meet the minimum demand requirement for constructing a refueling station. Therefore, the first-stage model stopped after the 8th iteration. The result in Table 1 suggests that there is no obvious correlation between segment flow and segment refueling demand. The result of the first-stage model is also shown in Figure 6.Figure 6
Locations selected for refueling station during the first stage.The stations selected during the first stage can directly satisfy a refueling demand of 3,430 pcu/h, which is about half of the total demand, meaning that approximately 50% of the refueling vehicles do not need to detour. In the following section, the unmet refueling demand will be discussed.


### 5.2.2. Result of the Second-Stage Model

(i) Model Initialization. As described previously, the second-stage model first generated Popsize=100 chromosomes and each chromosome consisted of 151 genes. Each gene is a binary variable, representing whether the corresponding segment is selected for constructing a refueling station.(ii) Parameter Setting. To prevent premature convergence, the following parameters were used: initial temperature in the annealing process T0=1000, terminating temperature Te=0, cooling speed TR=0.95, maximum number of iterations for each temperature was 100, crossover rate was 0.40~0.99, mutation rate was 0.0001~0.1, crossover probability was 0.9, mutation probability was 0.04, and the maximum number of iterations Oitermax=2000.The constraints in the second-stage model were considered by introducing a penalty. All solutions not meeting constraints (4)~(6) were punished by a multiplication of 10. The elitist strategy was adopted to select chromosomes to be included in the next generation, meaning chromosomes were selected with probabilities inversely proportional to their fitness values.(iii) Modeling Result. VB and GISDK in TransCAD were used to implement the second-stage algorithm. Among them, VB was used for coding the GSA algorithm, and GISDK was used for implementing the SUE assignment algorithm. After 876 cooling iterations, the temperature was 3.06E-17 and an optimal solution was obtained. The average run time was 24 minutes and 46 seconds. The obtained shortest total detour time for refueling vehicles was 6,688 minutes.Table2 lists the shortest total detour time and the detour times for different chromosomes at each iteration step. The iterative process is also illustrated in Figure 7, in which the horizontal axis stands for temperature and the vertical axis stands for fitness value. Figure 7 shows that the optimal fitness curve goes up and down at the beginning and tends to stabilize and level out after 876 iterations. This optimal fitness curve is not monotonically decreasing because the simulated annealing module sometimes accepts chromosomes with suboptimal (i.e., worse) fitness values with a small probability. The newly added stations in the second stage are shown in Figure 8. These stations provide service to the remaining 49.37% refueling vehicles.Table 2
Result of the second-stage model.
Number of iterations
Shortest total detour time (min)
Detour times of other solutions (min)
1
26752
27078
28244
27927
28470
⋯
28154
28637
2
26066
27129
28141
28923
28228
⋯
27130
26639
3
25045
28991
26479
27130
27910
⋯
28192
28087
4
25384
28136
28467
26486
27130
⋯
25631
26177
5
25218
25596
25922
26140
26551
⋯
26841
25676
6
23806
24469
25921
26685
24534
⋯
25170
25191
7
22362
22796
22344
22867
22845
⋯
25917
22948
8
22269
23082
23505
23600
22796
⋯
22796
25692
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
872
6991
8729
9455
7086
7138
⋯
7283
7136
873
6782
8007
8795
6926
6792
⋯
8869
6936
874
6688
6688
6688
8360
8024
⋯
6688
6688
875
6688
6688
6688
8125
7727
⋯
6688
6688
876
6688
6688
6688
6688
6688
⋯
6688
6688Figure 7
Iterative process of the second-stage model.Figure 8
Refueling stations selected during the second stage.(iv) Comparison of Algorithms. In this section, GA and the GSA algorithms are compared. The result is reported in Table 3. As can be seen, on average it takes the GA more iterations (1226) than the GSA algorithm to obtain the same solution. Also, the GA requires longer CPU time. For a larger and more complex road network, the advantages of the GSA algorithm could be more significant.Table 3
Algorithm performance comparison.
Algorithm
Average number of iterations
Average CPU time
Genetic simulated annealing algorithm
876
24 min 46 s
Genetic algorithm
1226
34 min 39 s(v) Sensitivity Analysis of Model Parameters. To assess the sensitivity of the GSA algorithm to model parameters, we experimented with various population sizes, crossover rates, and mutation rates. The modeling results are summarized and compared in Table 4. It seems that these parameters do have an important impact on model performance and should be decided carefully based on test runs.Table 4
Algorithm parameter sensitivity analysis results.
Parameter
Population size
Crossover rate
Mutation rate
Iterations
Calculating time
Bypass time
Optimization scheme
100
0.9
0.04
876
24′46′′
6688′
Test scheme 1
90
0.9
0.04
704
20′15′′
6847′
Test scheme 2
110
0.9
0.04
991
32′51′′
6688′
Test scheme 3
100
0.85
0.04
657
18′34′′
7351′
Test scheme 4
100
0.95
0.04
1324
40′42′′
6688′
Test scheme 5
100
0.9
0.03
754
21′23′′
7023′
Test scheme 6
100
0.9
0.05
1521
51′34′′
6688′


### 5.2.3. Result of the Third Stage

The stations selected in the previous two stages are shown in Figure9, since these stations were selected in an incremental manner. If we assign the refueling traffic to the network considering all stations in Figure 9 simultaneously, the refueling vehicles will be redistributed among these stations as shown in Table 5, which shows that three refueling stations now have demand greater than their capacities (i.e., 500 pcu/h) and another three stations are unprofitable (demand less than 300 pcu/h).Table 5
Refueling traffic volumes of each refueling station based on the location plan generated by the second-stage model.
Station number
Refueling demand
Number 2
605
Number 5
577
Number 6
535
Number 1
476
Number 3
468
Number 4
417
Number 7
378
Number 8
347
Number 14
332
Number 9
329
Number 10
327
Number 17
321
Number 11
318
Number 16
313
Number 12
316
Number 13
277
Number 15
241
Number 18
198Figure 9
Preliminary location plan for refueling stations.The purpose of the third stage was to adjust the solutions generated in the previous two stages according to the method described in Figure1. The result is shown in Table 6. As we adjusted/relocated refueling stations, the standard deviation (STDEV) of the numbers of refueling vehicles served by different refueling stations decreased. However, the total detour time increased. After 3 adjustments, the numbers of refueling vehicles served by each refueling station all fell in the range of [300,500] pcu/h as shown in Table 7. If we continue to relocate refueling stations, the total detour time of refueling vehicles most likely will keep increasing. Therefore, the third stage stopped here, since all selected refueling stations had a reasonable demand and the total detour time was still low. The final solution to this case study is shown in Figure 10.Table 6
Result of the third stage.
Number of adjustments
Adjusting station
Adjusting method
STDEV
Total detour time/min
0
—
—
114.7
6688
1
Number 15
link157 → link185
94.8
7058
2
Number 18
link165 → link168
71.3
7367
3
Number 13
Link34 → link13
62.1
7520Table 7
Final refueling traffic volumes at each station.
Station number
Refueling volume
Number 1
481
Number 2
473
Number 3
473
Number 5
450
Number 6
443
Number 4
421
Number 7
382
Number 8
364
Number 9
352
Number 14
342
Number 15
334
Number 10
333
Number 17
326
Number 13
324
Number 11
321
Number 18
320
Number 12
319
Number 16
317Figure 10
Final solution to the case study.The refueling stations selected during the first stage were directly on the paths of refueling vehicles, and the average detour time of these refueling vehicles was 0 minutes; the average detour time of detouring refueling vehicles served by stations selected during the second stage was approximately 2 minutes; the average detour time of new detouring refueling vehicles due to the adjustments in the third stage was 2.5 minutes; the average detour time of all detouring refueling vehicles was 2.1 minutes; and the average detouring time of all refueling vehicles was 1.1 minutes.


## 5.2.1. Results of the First-Stage Model

TransCAD 5.0 was used to solve the first-stage model. As shown in Table1, the first-stage model finished after 8 iterations. During this stage, the capacity of refueling stations was not considered. In other words, we assumed that a refueling station can satisfy all passing vehicles’ refueling demand.Table 1
Result of the first-stage model.
Number of iterations
Segment number
Segment direction
Segment flow (pcu/h)
Refueling demand (pcu/h)
1
125
112 → 121
4637
788
2
162
133 → 132
2902
610
3
111
105 → 106
3966
552
4
103
101 → 102
2341
421
5
56
75 → 74
2612
418
6
154
128 → 129
1983
333
7
98
98 → 107
1853
308
8
194
154 → 202
2141
278Table1 shows that the refueling demand of Segment 194 is 278 pcu/h, which is less than 300 pcu/h and does not meet the minimum demand requirement for constructing a refueling station. Therefore, the first-stage model stopped after the 8th iteration. The result in Table 1 suggests that there is no obvious correlation between segment flow and segment refueling demand. The result of the first-stage model is also shown in Figure 6.Figure 6
Locations selected for refueling station during the first stage.The stations selected during the first stage can directly satisfy a refueling demand of 3,430 pcu/h, which is about half of the total demand, meaning that approximately 50% of the refueling vehicles do not need to detour. In the following section, the unmet refueling demand will be discussed.


## 5.2.2. Result of the Second-Stage Model

(i) Model Initialization. As described previously, the second-stage model first generated Popsize=100 chromosomes and each chromosome consisted of 151 genes. Each gene is a binary variable, representing whether the corresponding segment is selected for constructing a refueling station.(ii) Parameter Setting. To prevent premature convergence, the following parameters were used: initial temperature in the annealing process T0=1000, terminating temperature Te=0, cooling speed TR=0.95, maximum number of iterations for each temperature was 100, crossover rate was 0.40~0.99, mutation rate was 0.0001~0.1, crossover probability was 0.9, mutation probability was 0.04, and the maximum number of iterations Oitermax=2000.The constraints in the second-stage model were considered by introducing a penalty. All solutions not meeting constraints (4)~(6) were punished by a multiplication of 10. The elitist strategy was adopted to select chromosomes to be included in the next generation, meaning chromosomes were selected with probabilities inversely proportional to their fitness values.(iii) Modeling Result. VB and GISDK in TransCAD were used to implement the second-stage algorithm. Among them, VB was used for coding the GSA algorithm, and GISDK was used for implementing the SUE assignment algorithm. After 876 cooling iterations, the temperature was 3.06E-17 and an optimal solution was obtained. The average run time was 24 minutes and 46 seconds. The obtained shortest total detour time for refueling vehicles was 6,688 minutes.Table2 lists the shortest total detour time and the detour times for different chromosomes at each iteration step. The iterative process is also illustrated in Figure 7, in which the horizontal axis stands for temperature and the vertical axis stands for fitness value. Figure 7 shows that the optimal fitness curve goes up and down at the beginning and tends to stabilize and level out after 876 iterations. This optimal fitness curve is not monotonically decreasing because the simulated annealing module sometimes accepts chromosomes with suboptimal (i.e., worse) fitness values with a small probability. The newly added stations in the second stage are shown in Figure 8. These stations provide service to the remaining 49.37% refueling vehicles.Table 2
Result of the second-stage model.
Number of iterations
Shortest total detour time (min)
Detour times of other solutions (min)
1
26752
27078
28244
27927
28470
⋯
28154
28637
2
26066
27129
28141
28923
28228
⋯
27130
26639
3
25045
28991
26479
27130
27910
⋯
28192
28087
4
25384
28136
28467
26486
27130
⋯
25631
26177
5
25218
25596
25922
26140
26551
⋯
26841
25676
6
23806
24469
25921
26685
24534
⋯
25170
25191
7
22362
22796
22344
22867
22845
⋯
25917
22948
8
22269
23082
23505
23600
22796
⋯
22796
25692
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
872
6991
8729
9455
7086
7138
⋯
7283
7136
873
6782
8007
8795
6926
6792
⋯
8869
6936
874
6688
6688
6688
8360
8024
⋯
6688
6688
875
6688
6688
6688
8125
7727
⋯
6688
6688
876
6688
6688
6688
6688
6688
⋯
6688
6688Figure 7
Iterative process of the second-stage model.Figure 8
Refueling stations selected during the second stage.(iv) Comparison of Algorithms. In this section, GA and the GSA algorithms are compared. The result is reported in Table 3. As can be seen, on average it takes the GA more iterations (1226) than the GSA algorithm to obtain the same solution. Also, the GA requires longer CPU time. For a larger and more complex road network, the advantages of the GSA algorithm could be more significant.Table 3
Algorithm performance comparison.
Algorithm
Average number of iterations
Average CPU time
Genetic simulated annealing algorithm
876
24 min 46 s
Genetic algorithm
1226
34 min 39 s(v) Sensitivity Analysis of Model Parameters. To assess the sensitivity of the GSA algorithm to model parameters, we experimented with various population sizes, crossover rates, and mutation rates. The modeling results are summarized and compared in Table 4. It seems that these parameters do have an important impact on model performance and should be decided carefully based on test runs.Table 4
Algorithm parameter sensitivity analysis results.
Parameter
Population size
Crossover rate
Mutation rate
Iterations
Calculating time
Bypass time
Optimization scheme
100
0.9
0.04
876
24′46′′
6688′
Test scheme 1
90
0.9
0.04
704
20′15′′
6847′
Test scheme 2
110
0.9
0.04
991
32′51′′
6688′
Test scheme 3
100
0.85
0.04
657
18′34′′
7351′
Test scheme 4
100
0.95
0.04
1324
40′42′′
6688′
Test scheme 5
100
0.9
0.03
754
21′23′′
7023′
Test scheme 6
100
0.9
0.05
1521
51′34′′
6688′


## 5.2.3. Result of the Third Stage

The stations selected in the previous two stages are shown in Figure9, since these stations were selected in an incremental manner. If we assign the refueling traffic to the network considering all stations in Figure 9 simultaneously, the refueling vehicles will be redistributed among these stations as shown in Table 5, which shows that three refueling stations now have demand greater than their capacities (i.e., 500 pcu/h) and another three stations are unprofitable (demand less than 300 pcu/h).Table 5
Refueling traffic volumes of each refueling station based on the location plan generated by the second-stage model.
Station number
Refueling demand
Number 2
605
Number 5
577
Number 6
535
Number 1
476
Number 3
468
Number 4
417
Number 7
378
Number 8
347
Number 14
332
Number 9
329
Number 10
327
Number 17
321
Number 11
318
Number 16
313
Number 12
316
Number 13
277
Number 15
241
Number 18
198Figure 9
Preliminary location plan for refueling stations.The purpose of the third stage was to adjust the solutions generated in the previous two stages according to the method described in Figure1. The result is shown in Table 6. As we adjusted/relocated refueling stations, the standard deviation (STDEV) of the numbers of refueling vehicles served by different refueling stations decreased. However, the total detour time increased. After 3 adjustments, the numbers of refueling vehicles served by each refueling station all fell in the range of [300,500] pcu/h as shown in Table 7. If we continue to relocate refueling stations, the total detour time of refueling vehicles most likely will keep increasing. Therefore, the third stage stopped here, since all selected refueling stations had a reasonable demand and the total detour time was still low. The final solution to this case study is shown in Figure 10.Table 6
Result of the third stage.
Number of adjustments
Adjusting station
Adjusting method
STDEV
Total detour time/min
0
—
—
114.7
6688
1
Number 15
link157 → link185
94.8
7058
2
Number 18
link165 → link168
71.3
7367
3
Number 13
Link34 → link13
62.1
7520Table 7
Final refueling traffic volumes at each station.
Station number
Refueling volume
Number 1
481
Number 2
473
Number 3
473
Number 5
450
Number 6
443
Number 4
421
Number 7
382
Number 8
364
Number 9
352
Number 14
342
Number 15
334
Number 10
333
Number 17
326
Number 13
324
Number 11
321
Number 18
320
Number 12
319
Number 16
317Figure 10
Final solution to the case study.The refueling stations selected during the first stage were directly on the paths of refueling vehicles, and the average detour time of these refueling vehicles was 0 minutes; the average detour time of detouring refueling vehicles served by stations selected during the second stage was approximately 2 minutes; the average detour time of new detouring refueling vehicles due to the adjustments in the third stage was 2.5 minutes; the average detour time of all detouring refueling vehicles was 2.1 minutes; and the average detouring time of all refueling vehicles was 1.1 minutes.


## 6. Conclusions

Urban refueling stations are an important component of urban transportation systems. Based on the potential impacts of urban refueling stations, several principles of siting urban refueling stations have been proposed. Built upon these principles and impact factors, this paper proposed a three-stage method to optimally site urban refueling stations, developed a genetic simulated annealing algorithm for solving the three-stage model, and conducted a case study to verify the effectiveness of the model and algorithm. The case study results suggest that (i) an optimal urban refueling station plan developed using the three-stage method can significantly shorten the detour time of refueling vehicles. Consequently, this can help reduce the negative impacts of detouring refueling vehicles on urban traffic, which in many cities is already highly congested; and (ii) the three-stage model takes the user equilibrium principle into consideration and can generate reasonable and practical facility location solutions.



---

*Source: 981370-2015-09-07.xml*