# A Novel Convex Hull Coverage Algorithm for the Deployment of Electric Taxi Swap Stations Based on Urban Traffic Flow

Authors: Xudong Zhu; Lijun Yao; Yuxin Mao; Renjie Wang; Hanxiao Shi; Ji Qi
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

Electric taxis have been adopted as a new energy public transportation tool as opposed to traditional taxis in modern city. Designing an efficient swapping station deployment scheme has become an important issue to improve the endurance capability of electric taxis. In this study, based on real operation trajectory data from 3997 taxis in Suzhou city, the battery swapping demand of taxis based on urban traffic flow is obtained to construct a network coverage deployment model under rule constraints, where the main optimization goal is to minimize the number of swapping stations and load balancing. According to this model, a traffic drive planning algorithm based on computational geometry is presented. The experimental results illustrate that the deployment scheme obtained by the proposed algorithm is significantly optimized in terms of the deployment cost and service load and has a lower algorithm time complexity than the typical unified deployment scheme. Therefore, the proposed method can be applied to improve the operating efficiency of the urban electric taxi system.


---

## Body

## 1. Introduction

With the development of the global economy, the harm to the environment caused by traditional diesel locomotives has become increasingly serious. Based on this, cities around the world have begun to vigorously promote new energy vehicles, especially electric taxis, in order to improve the urban public transportation system and alleviate the energy and environmental crisis brought by the transportation sector. However, first the current technical barriers need to be addressed; the actual endurance of electric vehicles is approximately 150 to 250 kilometers, and a taxi needs to travel 400 to 600 kilometers a day. At present, the basic energy supply methods for electric vehicles mainly include the charging mode and battery swapping mode. The charging mode is too expensive for taxis, so the battery swapping mode has become the first choice for energy supply. The deployment plan of the battery swapping station affects the operating cost and efficiency of the electric taxi system. In this paper, under the conditions of meeting the demand of urban taxis for battery swapping, the main goal is to minimize the deployment cost of the power station. At the same time, it requires a service load balance of the deployment station and the utilization rate of the station, by considering the overall urban traffic flow. In order to discover relevant influencing factors, Yunet al. [1] developed models to identify and compare the factors that influence PHEV user’s charging location choices (home, workplace, and public stations). Chen et al. [2] developed binomial logit models to analyze the spatiotemporal characteristics of multimode travelers by combining the taxi FCD, the metro smartcard data, and the GPS trajectories of Mobike, one of the most popular shared bicycles in China, 2018. Sun and Ding [3] proposed a two-level growth model (GM) to identify the potential multilevel factors that may affect online ride-hailing service demand.In the field of facility network planning, Liang (2016) and Liu et al. [4, 5] established a multiobjective optimization model with time windows based on dynamic network traffic flow and adopted a two-stage algorithm to establish the location and capacity model of electric vehicle charging and swapping stations. In the field of facility deployment, early academic research [6–11] directly adopted location optimization methods based on oil and gas demand models. Hakimi (1964) [12] proposed the P-center problem and the P-median problem; Gavranovi (2014) et al. [13] designed the charging station network model by limiting the number and service volume of a certain area for the P-median problem. Toregas (1971) et al. [14], aimed at minimizing the total cost of the service facility construction and on the premise of meeting demand, first proposed the set coverage model (SCLM); Church and Revelle (1974) [15], based on Toregas, established the maximum coverage location model (MCLM) to maximize the coverage of a limited number of service facilities; Current (1985) et al. [16] aimed to maximize user coverage and minimize travel paths and proposed multiobjective programming location problem (MCSPP); Bapna (2002) et al. [17] proposed the MC3SP model based on the (MCSPP) model. In the early stage, the goal was to minimize the construction cost of service facilities. In the later stage, the goal was to minimize vehicle owners’ refilling costs and maximize coverage of user demand to optimize vehicle refilling service quality. In the field of spatial discrete optimization, Kuby and Lim (2005) [18] proposed an independent optimization model to determine the candidate locations of stations as well as the number and scope of service stations. Ma et al. [19, 20] proposed a coverage location problem based on time satisfaction and solved the model. Quan et al. [21, 22] used the particle swarm algorithm to solve the optimal deployment model of the substation. Ge et al. [23–26] established models with different goals, designed the relevant parameters of the Voronoi diagram, and used the particle swarm algorithm to solve the model. There is little research at present on the deployment of electric taxi swap stations, and most of the existing literature does not fully consider the surrounding environment of the candidate station or the load balance of the swap station.In this paper, we focus on the origin and destination of taxis, and other conventional origins and destinations are solved by public transportation and private vehicles. Based on this, we analyze the real operation trajectory data of 3997 taxis in Suzhou and obtain the taxi battery swapping demand based on urban traffic flow. Under the constraint of deployment rules, such as the service coverage and service load of a single swapping station, the main optimization objective is to minimize the number of swapping stations and load balancing, and the graph-covering algorithm for convex hull coverage is designed in combination with the abstract road network topology to obtain the approximate optimal solution for swapping station deployment, with a low algorithm time complexity. We have taken the urban traffic flow factor into consideration when constructing the network model. In addition, this paper also considers the differences between different areas of a city, divides the deployment area into available areas and empty areas, determines the minimum coverage that surrounds the empty areas, and gives an optimal deployment scheme to reduce the ineffective coverage rate and improve the utilization rate of the swapping station. The experimental results show that, compared with the unified deployment scheme, the deployment scheme obtained by the algorithm described in this paper is significantly optimized in terms of the deployment cost and service load, which provides a reference for improving the operating efficiency of the urban electric taxi system.The first section of this paper describes and models the deployment problem of electric taxi swap stations. The second section introduces the graph-covering algorithm of convex hull coverage by considering traffic flow in detail. The third section considers the improved deployment scheme for urban empty areas. Section4 analyzes the correctness and performance of the algorithm. In Section 5, taxis in Suzhou are selected for example verification, and Section 6 gives the conclusions of this paper, as well as future research directions.


## 2. Problem and Modeling

### 2.1. Problem Description and NP-C Proof

#### 2.1.1. Problem Description

The main problem facing the deployment of electric taxi swap stations is how to ensure that the service load of each swap station is balanced under the premise of minimizing the cost of station construction. This paper needs to design a deployment scheme under the constraints of the service load and service radius and use the graph-covering road network graphGV,E to complete station deployment. First, the service load capacity of any swapping station should be close to, but not exceed, the maximum load limit L; second, for the sake of simplicity, the service coverage of the swapping station in this paper is the coverage of the circle centered on the swapping station, and the service radius of any swapping station should strictly not exceed the maximum service radius R. The service load of swapping station sk comes from the total battery swapping demand SLsk of taxis traveling on all road sections within its service radius SRsk. The service load of the swapping station is largely related to the traffic load of the section within the service radius, is a function of time, and is closely related to the battery status of the taxi on the section. To simplify the calculation, this paper temporarily ignores the time factor and understands that the service load of a swapping station is proportional to the total traffic load of all sections within its service radius.In summary, the deployment problem of electric taxi swap stations is described as follows: given a road network graphGV,E, find the smallest station set BS=s1;s2;s3;…;sBS∀k∈1,2,3,…,BS,SRsk≤R&SLsk≤L&MinVarSLs1;SLs2;…;SLsBS so that the service coverage of the stations in the set BS can completely cover the convex hull of the vertex set V. This problem is a high-dimensional nonlinear programming problem, and it can be proven that this problem is a strictly NP-complete problem.


#### 2.1.2. Proof of an NP-Complete Problem

If the location of each road section in the city, the average battery swapping demand of the road section, the service radius limit, and the service load limit of the swapping station are known, Theorem1 can be obtained.Theorem 1.
The coverage problem based on the above framework definition is NP-complete.Proof.
The known minimum vertex cover (VC) problem is NP-complete. By modifying some conditions in the problem, this paper reduces this problem to a VC problem, which proves that it is also an NP-complete problem.
Generate ann×m lattice in the traffic network graph GV,E, where the distance between adjacent points is j, and the lattice just covers the convex hull of V. Assuming that the cost of constructing and operating any swapping station is fixed and the same, the service radius of the swapping station is limited to R, the service load is limited to L, and the traffic load of section v is TLv. The undirected graph G′V′,E′ can then be obtained through the following steps:Step 1: initializeV′=∅. For each vertex v in GV,E, find the nearest point v′ in the lattice. If v′∉V′, add v′ to V′.Step 2: forv′∈V′ and v∈V, if the Euclidean distance EDv′,v≤R and the total traffic load of v′∑v∈VTLv′,v≤L, then edge v′,v is added to edge set E′.Step 3: forv′∈V′, set the vertex v in maxEDv′,v|v′,v∈E″ to vmax, for ∀v″=V′/v′, if EDv′,v″<EDv′,vmax, add v′,v″ to the edge set E′.Step 4: for∀v∈V, add v to V′ in turn.
Then, the covering problem is transformed based on the undirected graphG′V′,E′, finding the minimum vertex subset V″ from V′, where any edge of E′ has at least one vertex in V″. This problem is a VC problem. From the above assumptions, we can see that the VC problem is a subproblem of the covering problem, so the original problem is NP-complete.


### 2.2. Problem Analysis

In the early study of urban traffic calculations [27], we learned that the traffic flow of the urban road network has obvious patterns and that the pattern and traffic load of road sections are determined by such factors as regional function, surrounding population, road level, transportation planning, and economic conditions. Because this paper covers all the taxis’ swapping demand in a load balancing manner, it is first necessary to determine the traffic load pattern of the road network. From a statistical point of view, the average demand for the swapping station is directly proportional to the average traffic load of the road section.


#### 2.2.1. Road Traffic Load Analysis

If the location of the road section is fixed, the area function, surrounding population, road level, and transportation planning will also be determined and will not change over a long period of time. Based on this, this paper establishes a model of the relationship between the traffic mode and the location of the road section. This paper thus defines the traffic load of sectionvi in the road network on day d as TLvid. By using the Kolmogorov-Smirnov test (K-S test), reference [28] proposes that the time statistics of vehicles arriving at road section vi are similar to the exponential distribution; that is, the traffic load statistics of road section vi in a fixed time period obey a Poisson distribution. Therefore, the probability distribution function of TLvi is(1)PTLviX=k=e−λviλvikk!,where TLvi represents the average traffic load statistics of road section vi in a day and λvi represents the incidence of taxis randomly arriving at road section vi in a day. λvi is calculated using the trajectory data of all taxis in Suzhou in 366 days in 2012:(2)λvi=∑d=1366TLvid366.Then, calculate and count the average traffic loadTLvi of road section vi by the following formula:(3)TLvi=EX=∑k=0∞ke−λviλvikk!=λvi.


#### 2.2.2. Service Load Analysis of the Swapping Station

As mentioned earlier, the service loadSLsk of swapping station sk is determined by the total traffic load of all sections within the service radius, namely,(4)TLsk=∑EDskvi≤SRskTLvi.Although the total traffic loadTLsk cannot be directly used as the service load SLsk of the swapping station, because they are still linearly proportional, this paper defines the service load SLsk of the swapping station sk as(5)SLsk=TLsk∗TSL∑k=1BSTLsk,where TSL represents the total service load of all swapping stations, that is, the total demand for battery swapping of all taxis in the entire city, and TV is the number of taxis. Assuming that each taxi needs to change electricity x times a day on average, TSL can be calculated as(6)TSL=TV∗x.Under normal circumstances, and based on the aforementioned facts and assumptions in this paper, a taxi only needs to swap the battery once a day to meet the daily travel demand.


### 2.3. Building a Network Model

The road network consists of road sections and intersections. The center of road sectioni is abstracted as node vi to obtain road section set V=v1,v2,…,vn. If vi and vj intersect, then connect vi and vj to generate edge evi,vj∈E, thereby obtaining the road network topology graph GV,E, which is located in the latitude and longitude coordinate system.Define the taxi set asTV=n0,n1,n2,…,nTV−1, where nj represents the j-th electric taxi; the trajectory of taxi nj on day d is TRnjd=vk1,vk2,vk3,…,vki, and vki1≤i≤l represents a certain road section in the road network; fvinjd represents the frequency of taxi nj passing through section vi on day d, and the traffic load of section vi on day d is TLvid=∑j=1TV−1fvinjd. The objective function and constraint conditions based on the above conditions are as follows:(7)Mins1,s2,…,sBSBS,VarSLs1,SLs2,…,SLsBSs.t.SRsk≤R,k∈1,2,…,BS,SLsk≤L,k∈1,2,…,BS,∑k=1BSSLsk=TSL,∀vi∈V,∃k∈1,2,…,BS,EDskvi≤SRsk,where BS is the total number of deployed swapping stations; Var is the variance function; SRsk≤R indicates that the service radius of swapping station sk is less than or equal to the maximum service radius of the station; SLsk≤L indicates that the service load of swapping station sk is less than or equal to the maximum load of the station; TSL is the total service load of all swapping stations; and EDsk,vi is the Euclidean distance between swapping station sk and section vi. In this model, we have taken the urban traffic flow factor into consideration.


## 2.1. Problem Description and NP-C Proof

### 2.1.1. Problem Description

The main problem facing the deployment of electric taxi swap stations is how to ensure that the service load of each swap station is balanced under the premise of minimizing the cost of station construction. This paper needs to design a deployment scheme under the constraints of the service load and service radius and use the graph-covering road network graphGV,E to complete station deployment. First, the service load capacity of any swapping station should be close to, but not exceed, the maximum load limit L; second, for the sake of simplicity, the service coverage of the swapping station in this paper is the coverage of the circle centered on the swapping station, and the service radius of any swapping station should strictly not exceed the maximum service radius R. The service load of swapping station sk comes from the total battery swapping demand SLsk of taxis traveling on all road sections within its service radius SRsk. The service load of the swapping station is largely related to the traffic load of the section within the service radius, is a function of time, and is closely related to the battery status of the taxi on the section. To simplify the calculation, this paper temporarily ignores the time factor and understands that the service load of a swapping station is proportional to the total traffic load of all sections within its service radius.In summary, the deployment problem of electric taxi swap stations is described as follows: given a road network graphGV,E, find the smallest station set BS=s1;s2;s3;…;sBS∀k∈1,2,3,…,BS,SRsk≤R&SLsk≤L&MinVarSLs1;SLs2;…;SLsBS so that the service coverage of the stations in the set BS can completely cover the convex hull of the vertex set V. This problem is a high-dimensional nonlinear programming problem, and it can be proven that this problem is a strictly NP-complete problem.


### 2.1.2. Proof of an NP-Complete Problem

If the location of each road section in the city, the average battery swapping demand of the road section, the service radius limit, and the service load limit of the swapping station are known, Theorem1 can be obtained.Theorem 1.
The coverage problem based on the above framework definition is NP-complete.Proof.
The known minimum vertex cover (VC) problem is NP-complete. By modifying some conditions in the problem, this paper reduces this problem to a VC problem, which proves that it is also an NP-complete problem.
Generate ann×m lattice in the traffic network graph GV,E, where the distance between adjacent points is j, and the lattice just covers the convex hull of V. Assuming that the cost of constructing and operating any swapping station is fixed and the same, the service radius of the swapping station is limited to R, the service load is limited to L, and the traffic load of section v is TLv. The undirected graph G′V′,E′ can then be obtained through the following steps:Step 1: initializeV′=∅. For each vertex v in GV,E, find the nearest point v′ in the lattice. If v′∉V′, add v′ to V′.Step 2: forv′∈V′ and v∈V, if the Euclidean distance EDv′,v≤R and the total traffic load of v′∑v∈VTLv′,v≤L, then edge v′,v is added to edge set E′.Step 3: forv′∈V′, set the vertex v in maxEDv′,v|v′,v∈E″ to vmax, for ∀v″=V′/v′, if EDv′,v″<EDv′,vmax, add v′,v″ to the edge set E′.Step 4: for∀v∈V, add v to V′ in turn.
Then, the covering problem is transformed based on the undirected graphG′V′,E′, finding the minimum vertex subset V″ from V′, where any edge of E′ has at least one vertex in V″. This problem is a VC problem. From the above assumptions, we can see that the VC problem is a subproblem of the covering problem, so the original problem is NP-complete.


## 2.1.1. Problem Description

The main problem facing the deployment of electric taxi swap stations is how to ensure that the service load of each swap station is balanced under the premise of minimizing the cost of station construction. This paper needs to design a deployment scheme under the constraints of the service load and service radius and use the graph-covering road network graphGV,E to complete station deployment. First, the service load capacity of any swapping station should be close to, but not exceed, the maximum load limit L; second, for the sake of simplicity, the service coverage of the swapping station in this paper is the coverage of the circle centered on the swapping station, and the service radius of any swapping station should strictly not exceed the maximum service radius R. The service load of swapping station sk comes from the total battery swapping demand SLsk of taxis traveling on all road sections within its service radius SRsk. The service load of the swapping station is largely related to the traffic load of the section within the service radius, is a function of time, and is closely related to the battery status of the taxi on the section. To simplify the calculation, this paper temporarily ignores the time factor and understands that the service load of a swapping station is proportional to the total traffic load of all sections within its service radius.In summary, the deployment problem of electric taxi swap stations is described as follows: given a road network graphGV,E, find the smallest station set BS=s1;s2;s3;…;sBS∀k∈1,2,3,…,BS,SRsk≤R&SLsk≤L&MinVarSLs1;SLs2;…;SLsBS so that the service coverage of the stations in the set BS can completely cover the convex hull of the vertex set V. This problem is a high-dimensional nonlinear programming problem, and it can be proven that this problem is a strictly NP-complete problem.


## 2.1.2. Proof of an NP-Complete Problem

If the location of each road section in the city, the average battery swapping demand of the road section, the service radius limit, and the service load limit of the swapping station are known, Theorem1 can be obtained.Theorem 1.
The coverage problem based on the above framework definition is NP-complete.Proof.
The known minimum vertex cover (VC) problem is NP-complete. By modifying some conditions in the problem, this paper reduces this problem to a VC problem, which proves that it is also an NP-complete problem.
Generate ann×m lattice in the traffic network graph GV,E, where the distance between adjacent points is j, and the lattice just covers the convex hull of V. Assuming that the cost of constructing and operating any swapping station is fixed and the same, the service radius of the swapping station is limited to R, the service load is limited to L, and the traffic load of section v is TLv. The undirected graph G′V′,E′ can then be obtained through the following steps:Step 1: initializeV′=∅. For each vertex v in GV,E, find the nearest point v′ in the lattice. If v′∉V′, add v′ to V′.Step 2: forv′∈V′ and v∈V, if the Euclidean distance EDv′,v≤R and the total traffic load of v′∑v∈VTLv′,v≤L, then edge v′,v is added to edge set E′.Step 3: forv′∈V′, set the vertex v in maxEDv′,v|v′,v∈E″ to vmax, for ∀v″=V′/v′, if EDv′,v″<EDv′,vmax, add v′,v″ to the edge set E′.Step 4: for∀v∈V, add v to V′ in turn.
Then, the covering problem is transformed based on the undirected graphG′V′,E′, finding the minimum vertex subset V″ from V′, where any edge of E′ has at least one vertex in V″. This problem is a VC problem. From the above assumptions, we can see that the VC problem is a subproblem of the covering problem, so the original problem is NP-complete.


## 2.2. Problem Analysis

In the early study of urban traffic calculations [27], we learned that the traffic flow of the urban road network has obvious patterns and that the pattern and traffic load of road sections are determined by such factors as regional function, surrounding population, road level, transportation planning, and economic conditions. Because this paper covers all the taxis’ swapping demand in a load balancing manner, it is first necessary to determine the traffic load pattern of the road network. From a statistical point of view, the average demand for the swapping station is directly proportional to the average traffic load of the road section.


### 2.2.1. Road Traffic Load Analysis

If the location of the road section is fixed, the area function, surrounding population, road level, and transportation planning will also be determined and will not change over a long period of time. Based on this, this paper establishes a model of the relationship between the traffic mode and the location of the road section. This paper thus defines the traffic load of sectionvi in the road network on day d as TLvid. By using the Kolmogorov-Smirnov test (K-S test), reference [28] proposes that the time statistics of vehicles arriving at road section vi are similar to the exponential distribution; that is, the traffic load statistics of road section vi in a fixed time period obey a Poisson distribution. Therefore, the probability distribution function of TLvi is(1)PTLviX=k=e−λviλvikk!,where TLvi represents the average traffic load statistics of road section vi in a day and λvi represents the incidence of taxis randomly arriving at road section vi in a day. λvi is calculated using the trajectory data of all taxis in Suzhou in 366 days in 2012:(2)λvi=∑d=1366TLvid366.Then, calculate and count the average traffic loadTLvi of road section vi by the following formula:(3)TLvi=EX=∑k=0∞ke−λviλvikk!=λvi.


### 2.2.2. Service Load Analysis of the Swapping Station

As mentioned earlier, the service loadSLsk of swapping station sk is determined by the total traffic load of all sections within the service radius, namely,(4)TLsk=∑EDskvi≤SRskTLvi.Although the total traffic loadTLsk cannot be directly used as the service load SLsk of the swapping station, because they are still linearly proportional, this paper defines the service load SLsk of the swapping station sk as(5)SLsk=TLsk∗TSL∑k=1BSTLsk,where TSL represents the total service load of all swapping stations, that is, the total demand for battery swapping of all taxis in the entire city, and TV is the number of taxis. Assuming that each taxi needs to change electricity x times a day on average, TSL can be calculated as(6)TSL=TV∗x.Under normal circumstances, and based on the aforementioned facts and assumptions in this paper, a taxi only needs to swap the battery once a day to meet the daily travel demand.


## 2.2.1. Road Traffic Load Analysis

If the location of the road section is fixed, the area function, surrounding population, road level, and transportation planning will also be determined and will not change over a long period of time. Based on this, this paper establishes a model of the relationship between the traffic mode and the location of the road section. This paper thus defines the traffic load of sectionvi in the road network on day d as TLvid. By using the Kolmogorov-Smirnov test (K-S test), reference [28] proposes that the time statistics of vehicles arriving at road section vi are similar to the exponential distribution; that is, the traffic load statistics of road section vi in a fixed time period obey a Poisson distribution. Therefore, the probability distribution function of TLvi is(1)PTLviX=k=e−λviλvikk!,where TLvi represents the average traffic load statistics of road section vi in a day and λvi represents the incidence of taxis randomly arriving at road section vi in a day. λvi is calculated using the trajectory data of all taxis in Suzhou in 366 days in 2012:(2)λvi=∑d=1366TLvid366.Then, calculate and count the average traffic loadTLvi of road section vi by the following formula:(3)TLvi=EX=∑k=0∞ke−λviλvikk!=λvi.


## 2.2.2. Service Load Analysis of the Swapping Station

As mentioned earlier, the service loadSLsk of swapping station sk is determined by the total traffic load of all sections within the service radius, namely,(4)TLsk=∑EDskvi≤SRskTLvi.Although the total traffic loadTLsk cannot be directly used as the service load SLsk of the swapping station, because they are still linearly proportional, this paper defines the service load SLsk of the swapping station sk as(5)SLsk=TLsk∗TSL∑k=1BSTLsk,where TSL represents the total service load of all swapping stations, that is, the total demand for battery swapping of all taxis in the entire city, and TV is the number of taxis. Assuming that each taxi needs to change electricity x times a day on average, TSL can be calculated as(6)TSL=TV∗x.Under normal circumstances, and based on the aforementioned facts and assumptions in this paper, a taxi only needs to swap the battery once a day to meet the daily travel demand.


## 2.3. Building a Network Model

The road network consists of road sections and intersections. The center of road sectioni is abstracted as node vi to obtain road section set V=v1,v2,…,vn. If vi and vj intersect, then connect vi and vj to generate edge evi,vj∈E, thereby obtaining the road network topology graph GV,E, which is located in the latitude and longitude coordinate system.Define the taxi set asTV=n0,n1,n2,…,nTV−1, where nj represents the j-th electric taxi; the trajectory of taxi nj on day d is TRnjd=vk1,vk2,vk3,…,vki, and vki1≤i≤l represents a certain road section in the road network; fvinjd represents the frequency of taxi nj passing through section vi on day d, and the traffic load of section vi on day d is TLvid=∑j=1TV−1fvinjd. The objective function and constraint conditions based on the above conditions are as follows:(7)Mins1,s2,…,sBSBS,VarSLs1,SLs2,…,SLsBSs.t.SRsk≤R,k∈1,2,…,BS,SLsk≤L,k∈1,2,…,BS,∑k=1BSSLsk=TSL,∀vi∈V,∃k∈1,2,…,BS,EDskvi≤SRsk,where BS is the total number of deployed swapping stations; Var is the variance function; SRsk≤R indicates that the service radius of swapping station sk is less than or equal to the maximum service radius of the station; SLsk≤L indicates that the service load of swapping station sk is less than or equal to the maximum load of the station; TSL is the total service load of all swapping stations; and EDsk,vi is the Euclidean distance between swapping station sk and section vi. In this model, we have taken the urban traffic flow factor into consideration.


## 3. Convex Hull Covering Algorithm

### 3.1. Overview of the Convex Hull Covering Algorithm

Define the convex hullCHG as the smallest convex set covering the road section set V in graph GV,E. Its boundary BCHG is a convex polygon whose vertices are all from the road section set V, and the line connecting the two vertices with the farthest distance on the boundary is the convex hull diameter. The main idea of the convex hull coverage algorithm is to gradually cover all of the demand points by using a circle coverage method, where the center of each circle is the station location of the swapping station, and the coverage covered by the circle is the service coverage of the station. The deployment algorithm described in this paper is shown in Algorithm 1: (1) calculate the convex hull CHG and its diameter of the road network GV,E using taxi trajectory data and road network data; (2) calculate the smallest set of covering circles covering the boundary of the convex hull with the load of the service coverage as the constraint condition; (3) generate the inner convex hull according to the results of the circle coverage; (4) iteratively execute steps (2) and (3) until all areas are covered. The set of the center of the covering circle is the deployment scheme of the swapping station, the details of which are described as follows .Figure 1
Schematic diagram of selecting the next circle center.Figure 2
Schematic diagram of the convex hull coverage.Algorithm 1: The BSSP algorithm
(1)
procedure produce convex hull (G)(2)
based on traffic flow generate (GV,E)(3)
Graham algorithm (G)(4)
end procedure(5)
procedure cover (G)(6)
Preparata algorithm (CHG)(7)
deploy the first stations0,G(8)
cover (G1)(9)
end procedure


### 3.2. Convex Hull Covering Algorithm Steps

Step 1: initializek=0, GkVk,Ek=G, call the Graham algorithm [29] to obtain the convex hull CHGk, and find the convex hull boundary BCHGk.Step 2: call the Preparata algorithm [30] according to the convex hull CHGk to obtain the convex hull diameter vikvjk¯,vik,vjk∈V.Step 3: generate the first covering circle around the convex hull. Take the pointo0 on the convex hull diameter vikvjk¯ to maximize SRs0=EDvi,o0 and SLs0≤L.The service loadSLs0 is calculated using formulas (4) and (5), and SRs0 and o0 represent the service radius and candidate station of s0, respectively. According to the properties of the convex hull, it can be proven that there must be two intersection points between the circle Co0,SRs0 and the boundary BCHGk of the convex hull. Let these two intersection points be q0 and qb, respectively. The specific process of generating the first covering circle is as follows.Step 4: if the recently drawn circle isCoi+1,SRSi+1, point qi and point qi+1 are the intersection points of the circle and the convex hull boundary BCHGk, respectively. Draw the next circle by using the following iterative algorithm.Step 4.1: initializeSRsi+1=R.Step 4.2: calculate the positionoi+1 of the center of the next circle according to the service radius SRsi+1. The detailed process is as follows.Assuming that the newly drawn circle isCoi,SRSi, which intersects the convex hull boundary at qi and intersects the circle Coi+1,SRsi+1 to be drawn at qi and qi′; then the center discriminant equation (8) can be obtained. However, knowing only the coordinates of point qi and the service radius SRsi+1, an infinite number of solutions for the center oi+1 can be obtained.(8)EDqi′,oi+1=EDqi,oi+1=EDqi+1,oi+1=SRsi+1.In this paper, in order to obtain the unique solution of the circle center coordinates, the solution with the highest coverage rate is selected from the solutions of the countless circle centeroi+1. Considering the two situations that may be encountered when drawing the next circle, mark them as PART I and PART II, respectively, and mark the arched areas covered by each situation as I, II, III, and IV with shadows, as shown in Figure 1. For PART I, the total area of the shadow area can be obtained by separately calculating the area of each arch and then summing them. If x=2R/c is defined and the arc corresponding to the arch is a minor arc, R is the radius of the circle, and c is the chord length of the arch, then the area S of an arch is calculated, as shown in the following equation:(9)S=c2∗14∗x2∗arcsin1x−x2−14,x∈1,∞.If equation (9) is directly used to calculate the areas of the four arched regions, then the final equation expression is very complicated. By analyzing the calculation formula of the arch area, it can be determined that if the length of the chord is given, then the final area of the arch is determined by the parameter x. Since the new center of the circle must be in the convex hull, the shaded area IV in PART I must be less than half of the circle oi+1. Under this premise, when the chord length c is a given constant, then the final arch area S decreases monotonically with respect to x. To simplify the calculation, this paper approximates the problem of minimizing the shadow area into a problem of minimizing the sum of chord length LA. For PART I in Figure 1, the judgment formula of LA minimization is adopted:(10)LA=EDqi−1′qi′SRsi+EDqi′qiSRsi+EDqi′qiSRsi+1+EDqiqi+1SRsi+1.Among them, pointqi−1′ and point qi are fixed, and point qi′ and point qi+1 are determined by the next circle. This paper derives formula (10) for minimizing LA, and it can be seen that the three chord lengths in the shadow area can be calculated in the following way.Assuming that the center of the circleoi+1 is a,b, the radius is r1, the center of the circle oi is c,d, the radius is r2, and the linear equation of the chord is qiqi+1 (where the convex hull boundary intersecting the circle is Ax+By+C=0), obviously the linear equation of the common chord qi′qi formed by the intersection of the two circles is 2xc−a+2yd−b+a2−c2+b2−d2−r12+r22=0. Among them, the linear equation where the common chord is located is obtained by subtracting the circle oi+1 and the circle oi. Suppose the equation of the line where the chord qi−1′qi′ is located is A1x+B1y+C1=0. First, according to the distance formula from the point to the straight line, the distance from the center of the circle a,b to the line Ax+By+C=0 is d1, the distance from the center of the circle c,d to the common chord qi′qi is d2, and the distance from the center c,d to the distance of the straight line A1x+B1y+C1=0 is d3. Then, according to the Pythagorean theorem, the lengths of the three chord lengths in the four shaded areas can be calculated separately and represented by L1, L2, L3. Then, the total chord length LA=L1+L2+L3 is obtained. Finally, the obtained three chord lengths are substituted into the LA minimization judgment formula (10) to obtain the minimum chord length.Figure1 PART II demonstrates the situation where the next circle drawn hits the corner of the convex hull. On the one hand, it is assumed that the covering circle is much smaller than the convex hull and that there are very few corners in a circle; on the other hand, most of the corners in the convex hull are obtuse angles. Therefore, the shaded area IV in PART II in Figure 1 can be approximately calculated as the area of the arcuate area formed by the intersection of circle Coi+1,SRsi+1 and line segment qiqi+1¯. Although shadow area IV in PART II is no longer an arched area at this time, it can still be considered that its area is proportional to the distance between the line segments qiqi+1¯, so equation (11) is used to simplify this special case.In Figure1, if circle Coi+1,SRsi+1 is the second circle in the convex hull, then circle Coi−1,SRsi−1 and shadow area I do not exist. At this time, this paper cannot directly calculate LA through equation (9). However, for the second circle this paper can still minimize the other three shaded areas II, III, and IV. In this case, the formula for calculating LA is as follows:(11)LA=EDqi′qiSRsi+EDqi′qiSRsi+1+EDqiqi+1SRsi+1.When calculating the position of the next circle center, we must strictly avoid the intersection of two nonadjacent circles. In this paper, the Euclidean distance between the two circle centersoi+1 and oi−1 is greater than the sum of their radii to meet the following conditions:(12)EDοi−1οi+1>SRsi−1+SRsi+1.Step 4.3: if swapping stationsi+1 is deployed in oi+1, then its current service radius and service load are SRsi,SLsi, and the service radius SRsi+1 is updated through iteration (13), where L is the maximum load of the swapping station. Obviously, the service radius SRsi can be continuously updated according to the relationship between the current service load SLsi and L, but the update process may cross an excessively large step and continue to iterate at this time. If α<SRsi+1′/SRsi+1≤1 is satisfied, where α is a fixed constant close to 1, then swapping station si+1 is determined to be deployed in oi+1. Otherwise, SRsi+1=SRsi+1′, so jump to Step 4.2.(13)SRsi′=SRsi∗LSLsi,SRsi=SRsi′.CircleCoi+1,SRsi+1 intersects the convex hull boundary BCHGk at qi and qi+1 and intersects circle Coi,SRsi at qi and qi′. As shown in Figure 2, Step 4 is repeated until EDqi+1,qb≤2R, where qi+1 is the intersection of circle Coi+1,SRsi+1 and convex hull boundary BCHGk.CircleCoi+1,SRsi+1 intersects the convex hull boundary BCHGk at qi and qi+1 and intersects circle Coi,SRsi at qi and qi′. As shown in Figure 2, Step 4 is repeated until EDqi+1,qb≤2R, where qi+1 is the intersection of circle Coi+1,SRsi+1 and convex hull boundary BCHGk.Step 5: draw a circleCob,SRSb with a radius SRsb=R and pass through points qi+1 and qb. Suppose that a swapping station sb is deployed at ob and calculate the service load SLsb of the station. If SLsb≤L, then ob is a candidate location for the deployment of swapping station sb. Otherwise, call an iterative program similar to Step 4 to adjust the position of the service radius SRsb and the center ob. The specific process of generating the first and next covering circles is shown in Algorithms 2 and 3.Algorithm 2: Deploying the first station.
(1)
procedureDEPLOY THE FIRST STATIONsi,G(2)
SRsi←R(3)
vivj¯←DiameterCHG(4)
Selectoifromvivj¯where EDvioi=SRsi(5)
Locationsi←oi(6)
ifSLsi>Lthen(7)
repeat(8)
SRsi←ReductionSRsi(9)
Selectoifromvivj¯where EDvioi=SRsi(10)
Locationsi←oi(11)
untilα<SLsi/L≤1(12)
end if(13)
DEPLOYANORMAL STATIONsi+1,G(14)
end procedureAlgorithm 3: Deploying a normal station.
(1)
procedureDEPLOYANORMAL STATIONsi+1,G(2)
SRsi+1←R(3)
Selectoi+1whereEDqi−1oi+1=SRsi+1MinLA∈CHG(4)
Locationsi+1←oi+1(5)
ifSLsi+1>Lthen(6)
repeat(7)
SRsi+1←ReductionSRsi+1(8)
Selectoi+1whereEDqi−1oi+1=SRsi+1MinLA∈CHG(9)
Locationsi+1←oi+1(10)
untilα<SLsi+1/L≤1(11)
end if(12)
ifEDqbqi−1>2∗SRsi+2then(13)
DEPLOYANEXT NORMAL STATIONsi+3,G(14)
else(15)
DEPLOY THE LAST STATIONsi+3,G(16)
end if(17)
end procedureStep 6: if all the internal spaces of the convex hullCHGk have been seamlessly covered, then the station deployment ends; otherwise, connect the intersections q0′,q1′,q2′,…,qi′,qi+1′,qb′ of the covering circle with a line segment to generate a new inner convex hull CHGk+1k=k+1, and jump to Step 1 to cover the newly generated convex hull. The following details how to generate the inner convex hull.This algorithm uses circles to cover the convex hull layer by layer. When the boundary of a convex hull is completely covered, the internal space may still not be covered. Therefore, the algorithm will continue to generate a new inner convex hull and cover the boundary of the convex hull again, according to the method of Step 5, until all the inner space of the convex hull is completely covered. Figure3 shows an example of generating the inner convex hull, but there will be a special case when generating the inner convex hull, as shown in Figure 3(a). If the radius of the last circle Cob,SRSb is greater than EDqb,oi−Roi/2, then the circle Cob,SRSb intersects the circle Coi,SRsi, which is not adjacent to it twice. This intersection will cause the inner convex hull to be unable to be generated normally and cause a waste of coverage, as shown in Figure 3(b). Therefore, this paper requires that the radius of the last circle Cob,SRSb be smaller than EDqb,oi−Roi/2, that is, SRb<EDqb,oi−Roi/2, so that no such intersection will occur in order to ensure that the inner convex hull can be generated normally.Figure 3
Schematic diagram of generating an inner convex hull.
(a)(b)Since the service radius of each covered circle is not the same, the polygon generated by sequentially connecting the intersectionsq0′,q1′,q2′,q3′,…,qb′ of two adjacent circles is not necessarily a convex polygon. The dashed line in the upper half of Figure 3 is the actual boundary of the polygon generated by connecting the intersection points, but it is not the corresponding convex hull boundary. Therefore, in Step 6 you first need to connect these intersections and then generate the convex hull from the formed polygon.When the station deployment is over, since the convex hull of the vertices of the graphGkVk,Ek is completely covered, then all the vertices in set Vk are covered by the circle centered on the swapping station, and SRsi≤R&a<SLsi/L≤1 is guaranteed. The algorithm for deploying the last swapping station is shown in Algorithm 4.Algorithm 4: Deploying the last station.
(1)
procedureDEPLOY THE LAST STATIONsi+2,G(2)
SRsb←MinEDqboi−SRi/2,R(3)
SelectobwhereEDqiob=SRsbEDqbob=SRsbMinLA∈CHG(4)
Locationsb←ob(5)
ifSLsb>Lthen(6)
repeat(7)
SRsb←ReductionSRsb(8)
ifEDqbqi−1>2∗SRsbthen(9)
SelectobwhereEDqiob=SRsbMinLA∈CHG(10)
Locationsb←ob(11)
else(12)
SelectobwhereEDqiob=SRsbEDqbob=SRsbMinLA∈CHG(13)
Locationsb←ob(14)
end if(15)
untilα<SLsb/L≤1(16)
end if(17)
ifqi+2≠qbthen(18)
DEPLOY THE LAST STATIONsi+3,G(19)
else(20)
COVERG1=q0′,q1′,…,qi+1′,qb′(21)
end if(22)
end procedure


## 3.1. Overview of the Convex Hull Covering Algorithm

Define the convex hullCHG as the smallest convex set covering the road section set V in graph GV,E. Its boundary BCHG is a convex polygon whose vertices are all from the road section set V, and the line connecting the two vertices with the farthest distance on the boundary is the convex hull diameter. The main idea of the convex hull coverage algorithm is to gradually cover all of the demand points by using a circle coverage method, where the center of each circle is the station location of the swapping station, and the coverage covered by the circle is the service coverage of the station. The deployment algorithm described in this paper is shown in Algorithm 1: (1) calculate the convex hull CHG and its diameter of the road network GV,E using taxi trajectory data and road network data; (2) calculate the smallest set of covering circles covering the boundary of the convex hull with the load of the service coverage as the constraint condition; (3) generate the inner convex hull according to the results of the circle coverage; (4) iteratively execute steps (2) and (3) until all areas are covered. The set of the center of the covering circle is the deployment scheme of the swapping station, the details of which are described as follows .Figure 1
Schematic diagram of selecting the next circle center.Figure 2
Schematic diagram of the convex hull coverage.Algorithm 1: The BSSP algorithm
(1)
procedure produce convex hull (G)(2)
based on traffic flow generate (GV,E)(3)
Graham algorithm (G)(4)
end procedure(5)
procedure cover (G)(6)
Preparata algorithm (CHG)(7)
deploy the first stations0,G(8)
cover (G1)(9)
end procedure


## 3.2. Convex Hull Covering Algorithm Steps

Step 1: initializek=0, GkVk,Ek=G, call the Graham algorithm [29] to obtain the convex hull CHGk, and find the convex hull boundary BCHGk.Step 2: call the Preparata algorithm [30] according to the convex hull CHGk to obtain the convex hull diameter vikvjk¯,vik,vjk∈V.Step 3: generate the first covering circle around the convex hull. Take the pointo0 on the convex hull diameter vikvjk¯ to maximize SRs0=EDvi,o0 and SLs0≤L.The service loadSLs0 is calculated using formulas (4) and (5), and SRs0 and o0 represent the service radius and candidate station of s0, respectively. According to the properties of the convex hull, it can be proven that there must be two intersection points between the circle Co0,SRs0 and the boundary BCHGk of the convex hull. Let these two intersection points be q0 and qb, respectively. The specific process of generating the first covering circle is as follows.Step 4: if the recently drawn circle isCoi+1,SRSi+1, point qi and point qi+1 are the intersection points of the circle and the convex hull boundary BCHGk, respectively. Draw the next circle by using the following iterative algorithm.Step 4.1: initializeSRsi+1=R.Step 4.2: calculate the positionoi+1 of the center of the next circle according to the service radius SRsi+1. The detailed process is as follows.Assuming that the newly drawn circle isCoi,SRSi, which intersects the convex hull boundary at qi and intersects the circle Coi+1,SRsi+1 to be drawn at qi and qi′; then the center discriminant equation (8) can be obtained. However, knowing only the coordinates of point qi and the service radius SRsi+1, an infinite number of solutions for the center oi+1 can be obtained.(8)EDqi′,oi+1=EDqi,oi+1=EDqi+1,oi+1=SRsi+1.In this paper, in order to obtain the unique solution of the circle center coordinates, the solution with the highest coverage rate is selected from the solutions of the countless circle centeroi+1. Considering the two situations that may be encountered when drawing the next circle, mark them as PART I and PART II, respectively, and mark the arched areas covered by each situation as I, II, III, and IV with shadows, as shown in Figure 1. For PART I, the total area of the shadow area can be obtained by separately calculating the area of each arch and then summing them. If x=2R/c is defined and the arc corresponding to the arch is a minor arc, R is the radius of the circle, and c is the chord length of the arch, then the area S of an arch is calculated, as shown in the following equation:(9)S=c2∗14∗x2∗arcsin1x−x2−14,x∈1,∞.If equation (9) is directly used to calculate the areas of the four arched regions, then the final equation expression is very complicated. By analyzing the calculation formula of the arch area, it can be determined that if the length of the chord is given, then the final area of the arch is determined by the parameter x. Since the new center of the circle must be in the convex hull, the shaded area IV in PART I must be less than half of the circle oi+1. Under this premise, when the chord length c is a given constant, then the final arch area S decreases monotonically with respect to x. To simplify the calculation, this paper approximates the problem of minimizing the shadow area into a problem of minimizing the sum of chord length LA. For PART I in Figure 1, the judgment formula of LA minimization is adopted:(10)LA=EDqi−1′qi′SRsi+EDqi′qiSRsi+EDqi′qiSRsi+1+EDqiqi+1SRsi+1.Among them, pointqi−1′ and point qi are fixed, and point qi′ and point qi+1 are determined by the next circle. This paper derives formula (10) for minimizing LA, and it can be seen that the three chord lengths in the shadow area can be calculated in the following way.Assuming that the center of the circleoi+1 is a,b, the radius is r1, the center of the circle oi is c,d, the radius is r2, and the linear equation of the chord is qiqi+1 (where the convex hull boundary intersecting the circle is Ax+By+C=0), obviously the linear equation of the common chord qi′qi formed by the intersection of the two circles is 2xc−a+2yd−b+a2−c2+b2−d2−r12+r22=0. Among them, the linear equation where the common chord is located is obtained by subtracting the circle oi+1 and the circle oi. Suppose the equation of the line where the chord qi−1′qi′ is located is A1x+B1y+C1=0. First, according to the distance formula from the point to the straight line, the distance from the center of the circle a,b to the line Ax+By+C=0 is d1, the distance from the center of the circle c,d to the common chord qi′qi is d2, and the distance from the center c,d to the distance of the straight line A1x+B1y+C1=0 is d3. Then, according to the Pythagorean theorem, the lengths of the three chord lengths in the four shaded areas can be calculated separately and represented by L1, L2, L3. Then, the total chord length LA=L1+L2+L3 is obtained. Finally, the obtained three chord lengths are substituted into the LA minimization judgment formula (10) to obtain the minimum chord length.Figure1 PART II demonstrates the situation where the next circle drawn hits the corner of the convex hull. On the one hand, it is assumed that the covering circle is much smaller than the convex hull and that there are very few corners in a circle; on the other hand, most of the corners in the convex hull are obtuse angles. Therefore, the shaded area IV in PART II in Figure 1 can be approximately calculated as the area of the arcuate area formed by the intersection of circle Coi+1,SRsi+1 and line segment qiqi+1¯. Although shadow area IV in PART II is no longer an arched area at this time, it can still be considered that its area is proportional to the distance between the line segments qiqi+1¯, so equation (11) is used to simplify this special case.In Figure1, if circle Coi+1,SRsi+1 is the second circle in the convex hull, then circle Coi−1,SRsi−1 and shadow area I do not exist. At this time, this paper cannot directly calculate LA through equation (9). However, for the second circle this paper can still minimize the other three shaded areas II, III, and IV. In this case, the formula for calculating LA is as follows:(11)LA=EDqi′qiSRsi+EDqi′qiSRsi+1+EDqiqi+1SRsi+1.When calculating the position of the next circle center, we must strictly avoid the intersection of two nonadjacent circles. In this paper, the Euclidean distance between the two circle centersoi+1 and oi−1 is greater than the sum of their radii to meet the following conditions:(12)EDοi−1οi+1>SRsi−1+SRsi+1.Step 4.3: if swapping stationsi+1 is deployed in oi+1, then its current service radius and service load are SRsi,SLsi, and the service radius SRsi+1 is updated through iteration (13), where L is the maximum load of the swapping station. Obviously, the service radius SRsi can be continuously updated according to the relationship between the current service load SLsi and L, but the update process may cross an excessively large step and continue to iterate at this time. If α<SRsi+1′/SRsi+1≤1 is satisfied, where α is a fixed constant close to 1, then swapping station si+1 is determined to be deployed in oi+1. Otherwise, SRsi+1=SRsi+1′, so jump to Step 4.2.(13)SRsi′=SRsi∗LSLsi,SRsi=SRsi′.CircleCoi+1,SRsi+1 intersects the convex hull boundary BCHGk at qi and qi+1 and intersects circle Coi,SRsi at qi and qi′. As shown in Figure 2, Step 4 is repeated until EDqi+1,qb≤2R, where qi+1 is the intersection of circle Coi+1,SRsi+1 and convex hull boundary BCHGk.CircleCoi+1,SRsi+1 intersects the convex hull boundary BCHGk at qi and qi+1 and intersects circle Coi,SRsi at qi and qi′. As shown in Figure 2, Step 4 is repeated until EDqi+1,qb≤2R, where qi+1 is the intersection of circle Coi+1,SRsi+1 and convex hull boundary BCHGk.Step 5: draw a circleCob,SRSb with a radius SRsb=R and pass through points qi+1 and qb. Suppose that a swapping station sb is deployed at ob and calculate the service load SLsb of the station. If SLsb≤L, then ob is a candidate location for the deployment of swapping station sb. Otherwise, call an iterative program similar to Step 4 to adjust the position of the service radius SRsb and the center ob. The specific process of generating the first and next covering circles is shown in Algorithms 2 and 3.Algorithm 2: Deploying the first station.
(1)
procedureDEPLOY THE FIRST STATIONsi,G(2)
SRsi←R(3)
vivj¯←DiameterCHG(4)
Selectoifromvivj¯where EDvioi=SRsi(5)
Locationsi←oi(6)
ifSLsi>Lthen(7)
repeat(8)
SRsi←ReductionSRsi(9)
Selectoifromvivj¯where EDvioi=SRsi(10)
Locationsi←oi(11)
untilα<SLsi/L≤1(12)
end if(13)
DEPLOYANORMAL STATIONsi+1,G(14)
end procedureAlgorithm 3: Deploying a normal station.
(1)
procedureDEPLOYANORMAL STATIONsi+1,G(2)
SRsi+1←R(3)
Selectoi+1whereEDqi−1oi+1=SRsi+1MinLA∈CHG(4)
Locationsi+1←oi+1(5)
ifSLsi+1>Lthen(6)
repeat(7)
SRsi+1←ReductionSRsi+1(8)
Selectoi+1whereEDqi−1oi+1=SRsi+1MinLA∈CHG(9)
Locationsi+1←oi+1(10)
untilα<SLsi+1/L≤1(11)
end if(12)
ifEDqbqi−1>2∗SRsi+2then(13)
DEPLOYANEXT NORMAL STATIONsi+3,G(14)
else(15)
DEPLOY THE LAST STATIONsi+3,G(16)
end if(17)
end procedureStep 6: if all the internal spaces of the convex hullCHGk have been seamlessly covered, then the station deployment ends; otherwise, connect the intersections q0′,q1′,q2′,…,qi′,qi+1′,qb′ of the covering circle with a line segment to generate a new inner convex hull CHGk+1k=k+1, and jump to Step 1 to cover the newly generated convex hull. The following details how to generate the inner convex hull.This algorithm uses circles to cover the convex hull layer by layer. When the boundary of a convex hull is completely covered, the internal space may still not be covered. Therefore, the algorithm will continue to generate a new inner convex hull and cover the boundary of the convex hull again, according to the method of Step 5, until all the inner space of the convex hull is completely covered. Figure3 shows an example of generating the inner convex hull, but there will be a special case when generating the inner convex hull, as shown in Figure 3(a). If the radius of the last circle Cob,SRSb is greater than EDqb,oi−Roi/2, then the circle Cob,SRSb intersects the circle Coi,SRsi, which is not adjacent to it twice. This intersection will cause the inner convex hull to be unable to be generated normally and cause a waste of coverage, as shown in Figure 3(b). Therefore, this paper requires that the radius of the last circle Cob,SRSb be smaller than EDqb,oi−Roi/2, that is, SRb<EDqb,oi−Roi/2, so that no such intersection will occur in order to ensure that the inner convex hull can be generated normally.Figure 3
Schematic diagram of generating an inner convex hull.
(a)(b)Since the service radius of each covered circle is not the same, the polygon generated by sequentially connecting the intersectionsq0′,q1′,q2′,q3′,…,qb′ of two adjacent circles is not necessarily a convex polygon. The dashed line in the upper half of Figure 3 is the actual boundary of the polygon generated by connecting the intersection points, but it is not the corresponding convex hull boundary. Therefore, in Step 6 you first need to connect these intersections and then generate the convex hull from the formed polygon.When the station deployment is over, since the convex hull of the vertices of the graphGkVk,Ek is completely covered, then all the vertices in set Vk are covered by the circle centered on the swapping station, and SRsi≤R&a<SLsi/L≤1 is guaranteed. The algorithm for deploying the last swapping station is shown in Algorithm 4.Algorithm 4: Deploying the last station.
(1)
procedureDEPLOY THE LAST STATIONsi+2,G(2)
SRsb←MinEDqboi−SRi/2,R(3)
SelectobwhereEDqiob=SRsbEDqbob=SRsbMinLA∈CHG(4)
Locationsb←ob(5)
ifSLsb>Lthen(6)
repeat(7)
SRsb←ReductionSRsb(8)
ifEDqbqi−1>2∗SRsbthen(9)
SelectobwhereEDqiob=SRsbMinLA∈CHG(10)
Locationsb←ob(11)
else(12)
SelectobwhereEDqiob=SRsbEDqbob=SRsbMinLA∈CHG(13)
Locationsb←ob(14)
end if(15)
untilα<SLsb/L≤1(16)
end if(17)
ifqi+2≠qbthen(18)
DEPLOY THE LAST STATIONsi+3,G(19)
else(20)
COVERG1=q0′,q1′,…,qi+1′,qb′(21)
end if(22)
end procedure


## 4. Improving the Deployment Scheme of Swapping Stations Based on Empty Urban Areas

### 4.1. Empty Area

The deployment scheme described in this paper is based on the assumption that all areas of the city need to be covered and that any node in the convex hull can be deployed for swapping stations. However, most cities have areas that are particularly large and cannot be deployed, such as lakes, rivers, large buildings, mountains, and forests. These areas are empty areas. The above deployment scheme may have the problem of deploying the swap station in an empty area and invalid coverage of the deployment area. Therefore, this paper improves the convex hull coverage algorithm to solve the above two special cases and ensure the rationality and effectiveness of the deployment scheme.To identify the candidate locations of which stations in the station set are located in the empty area and the invalid coverage area, it is necessary to accurately measure the boundaries of all the empty areas in the city. For most areas, a loop is bound to be formed by connecting a series of road nodes. This paper will ignore all the areas within the loop, that is, whether the area within the loop is truly located in an empty area. Therefore, this paper takes the smallest loop surrounding the empty area as the empty boundary. Since a part of the road nodes in the area may be outside the polygon, the polygon constructed by connecting the road nodes around an area does not correspond to a real loop. Therefore, we must first add intersection nodes to update the road network topology graphGV,E to draw a real loop around an empty area. Suppose that the set C=c1,c2,c3,…,cc represents all of the intersections in the road network, and then change GV,E to G″V,C,E″, where ∀v∈V&c∈C, and side e=v,c∈E″ indicates that one end of the road section v is the intersection c. In a real map, all the empty areas have been marked as a point. The problem now is how to calculate the minimum loop through the given network graph G″ and the punctuation of the empty area.Assuming that the punctuation of an empty area isu, initialize the loop set Bu=∅, and then search for the vertex closest to the punctuation u in set C or set V. Without loss of generality, assuming that the vertex is vi or cj, the vertex vi must be added to the loop set Bu, and then all the intersections adjacent to the vertex vi are searched to find the intersection ck that is closest to the punctuation u. If vertex vi is not unique among the adjacent road sections of ck, that is, if ck is not the end point, continue to add ck to loop set Bu and search again for all road sections adjacent to vertex ck. Otherwise, take the intersection that is not the end point and the closest to the loop as the next intersection of vi. For intersection ck, an intersection in adjacent road sections is the intersection that is closest to u and connects at least one non-end point (except ck) and it is added to Bu. In this way, road sections and intersections are added to the loop set one by one. If the search and addition work ends due to the selection of the intersection in the current intersection or road section, then the search and addition task will continue to be performed at the intersection of other feasible paths. Through this depth-first search algorithm to perform search and backtracking work, road sections and intersections are searched and added one by one, until the intersection or the vertices of the road section are in the set of adjacent road sections of the current intersection. In this case, Bu is the boundary of the empty area u in the loop. If the initial vertex is ci, add the intersections and road sections to Bu one by one, according to the same principle described above.


### 4.2. Adjusting the Deployment Scheme

When improving the deployment scheme, the deployed stations related to the empty area will be divided into three parts by the road boundary polygon. The first part is that the station’s candidate location and service radius are completely within the polygon. Since it is meaningless to deploy such stations, you only need to delete them directly from the candidate stations. The second part is that the candidate location of the station is in the polygon, but its service radius exceeds the boundary of the polygon. If it is deleted directly, then the area outside the polygon that needs to be covered cannot be served, but if it is not adjusted, then the candidate location is not available. Therefore, they must be deleted first, and then all the stations located outside the polygon must be adjusted or new stations should be redeployed at appropriate locations outside the polygon to cover the unserved area. The last part is the station located outside the polygon, but its service radius is within the polygon. Since this type of station itself is outside the polygon, then there is no need to adjust its deployment, so this paper mainly focuses on the station located in the second part.For stationsi in the second part, since its candidate position is in the polygon, this paper can ensure that the width of the part exceeding the polygon must be smaller than SRsi. According to the proposed algorithm, there must be candidate stations that are adjacent to si and available. However, because the coverage circle is tightly coupled, the candidate positions cannot be changed at will. At the same time, since any station in the station set has reached its service radius limit R or service load limit L, the service radius of the swapping station cannot be changed arbitrarily. Therefore, this paper cannot re-cover the uncovered area by adjusting the station at a reasonable location but it can rebuild a new station to cover the uncovered area.This paper describes the problem of updating the deployment scheme as follows: given an empty area and the loop boundary surrounding the empty area, a new area that has not been covered appears due to the deletion of candidate stations in the boundary polygon during the redeployment of new station coverage. As shown in Figure4, the inner polygon Bu is the loop boundary of area u, and the adjacency circle is the available station adjacent to but not intersecting with Bu. By connecting the intersection points between the adjacent circles, the outer boundary polygon of the shaded area BOu is generated. The middle area between Bu and BOu is a new area that needs to be covered again, and its width must be less than the service radius limit R.Figure 4
Schematic diagram of the site for redeploying the hollow area.Because some road sections in the loop boundary may not have the second part of the station, the new uncovered area may be a ring or a segmented, incomplete ring. As shown in Figure4(b), the two ends of the incomplete ring are a line segment formed by connecting the intersection of two circles and the intersection of the two circles and the road boundary; that is, when there is an intersection between Bu and BOu, the ring will be cut. Because the uncovered new area is sufficiently narrow, the second part of the station only requires one layer of the service circle, so this paper can only use one layer of the coverage circle to complete the coverage under the condition of meeting the deployment restrictions. Under this premise, this paper proposes an improved algorithm for the empty problem. As shown in Figure 5, pi and qi are two points in the boundaries Bu and BOu, respectively, and a station needs to be deployed through these two points to cover a part of the shadow area. First, draw two arcs with pi and qi as the center and SRsi as the radius. They intersect BOu at xi and yi and intersect each other at zi. If zi is in the shaded area, then zi can be considered a candidate location for the station. Otherwise, like the points pi+1 and qi+1 in Figure 5, select a point closer to pi+1 and qi+1 from xi+1 and yi+1 as the candidate position of the circle oi+1. Keep adjusting the coverage circle until all the shadow areas are completely covered. In the improved algorithm, the principle of reducing the service radius of the station is in the same form as the algorithm proposed above. For the left part of Figure 5, first generate the convex hull of the polygon BOu and its diameter, then select one end of the diameter as the initial point q1, and select the relatively close intersection point between the diameter and Bu as the initial point p1. For the right part of Figure 5, use two initial points p1 and q1 generated by the intersection of the circle and Bu and BOu at one end of the shaded area of the incomplete ring.Figure 5
Schematic diagram of deploying new stations.


## 4.1. Empty Area

The deployment scheme described in this paper is based on the assumption that all areas of the city need to be covered and that any node in the convex hull can be deployed for swapping stations. However, most cities have areas that are particularly large and cannot be deployed, such as lakes, rivers, large buildings, mountains, and forests. These areas are empty areas. The above deployment scheme may have the problem of deploying the swap station in an empty area and invalid coverage of the deployment area. Therefore, this paper improves the convex hull coverage algorithm to solve the above two special cases and ensure the rationality and effectiveness of the deployment scheme.To identify the candidate locations of which stations in the station set are located in the empty area and the invalid coverage area, it is necessary to accurately measure the boundaries of all the empty areas in the city. For most areas, a loop is bound to be formed by connecting a series of road nodes. This paper will ignore all the areas within the loop, that is, whether the area within the loop is truly located in an empty area. Therefore, this paper takes the smallest loop surrounding the empty area as the empty boundary. Since a part of the road nodes in the area may be outside the polygon, the polygon constructed by connecting the road nodes around an area does not correspond to a real loop. Therefore, we must first add intersection nodes to update the road network topology graphGV,E to draw a real loop around an empty area. Suppose that the set C=c1,c2,c3,…,cc represents all of the intersections in the road network, and then change GV,E to G″V,C,E″, where ∀v∈V&c∈C, and side e=v,c∈E″ indicates that one end of the road section v is the intersection c. In a real map, all the empty areas have been marked as a point. The problem now is how to calculate the minimum loop through the given network graph G″ and the punctuation of the empty area.Assuming that the punctuation of an empty area isu, initialize the loop set Bu=∅, and then search for the vertex closest to the punctuation u in set C or set V. Without loss of generality, assuming that the vertex is vi or cj, the vertex vi must be added to the loop set Bu, and then all the intersections adjacent to the vertex vi are searched to find the intersection ck that is closest to the punctuation u. If vertex vi is not unique among the adjacent road sections of ck, that is, if ck is not the end point, continue to add ck to loop set Bu and search again for all road sections adjacent to vertex ck. Otherwise, take the intersection that is not the end point and the closest to the loop as the next intersection of vi. For intersection ck, an intersection in adjacent road sections is the intersection that is closest to u and connects at least one non-end point (except ck) and it is added to Bu. In this way, road sections and intersections are added to the loop set one by one. If the search and addition work ends due to the selection of the intersection in the current intersection or road section, then the search and addition task will continue to be performed at the intersection of other feasible paths. Through this depth-first search algorithm to perform search and backtracking work, road sections and intersections are searched and added one by one, until the intersection or the vertices of the road section are in the set of adjacent road sections of the current intersection. In this case, Bu is the boundary of the empty area u in the loop. If the initial vertex is ci, add the intersections and road sections to Bu one by one, according to the same principle described above.


## 4.2. Adjusting the Deployment Scheme

When improving the deployment scheme, the deployed stations related to the empty area will be divided into three parts by the road boundary polygon. The first part is that the station’s candidate location and service radius are completely within the polygon. Since it is meaningless to deploy such stations, you only need to delete them directly from the candidate stations. The second part is that the candidate location of the station is in the polygon, but its service radius exceeds the boundary of the polygon. If it is deleted directly, then the area outside the polygon that needs to be covered cannot be served, but if it is not adjusted, then the candidate location is not available. Therefore, they must be deleted first, and then all the stations located outside the polygon must be adjusted or new stations should be redeployed at appropriate locations outside the polygon to cover the unserved area. The last part is the station located outside the polygon, but its service radius is within the polygon. Since this type of station itself is outside the polygon, then there is no need to adjust its deployment, so this paper mainly focuses on the station located in the second part.For stationsi in the second part, since its candidate position is in the polygon, this paper can ensure that the width of the part exceeding the polygon must be smaller than SRsi. According to the proposed algorithm, there must be candidate stations that are adjacent to si and available. However, because the coverage circle is tightly coupled, the candidate positions cannot be changed at will. At the same time, since any station in the station set has reached its service radius limit R or service load limit L, the service radius of the swapping station cannot be changed arbitrarily. Therefore, this paper cannot re-cover the uncovered area by adjusting the station at a reasonable location but it can rebuild a new station to cover the uncovered area.This paper describes the problem of updating the deployment scheme as follows: given an empty area and the loop boundary surrounding the empty area, a new area that has not been covered appears due to the deletion of candidate stations in the boundary polygon during the redeployment of new station coverage. As shown in Figure4, the inner polygon Bu is the loop boundary of area u, and the adjacency circle is the available station adjacent to but not intersecting with Bu. By connecting the intersection points between the adjacent circles, the outer boundary polygon of the shaded area BOu is generated. The middle area between Bu and BOu is a new area that needs to be covered again, and its width must be less than the service radius limit R.Figure 4
Schematic diagram of the site for redeploying the hollow area.Because some road sections in the loop boundary may not have the second part of the station, the new uncovered area may be a ring or a segmented, incomplete ring. As shown in Figure4(b), the two ends of the incomplete ring are a line segment formed by connecting the intersection of two circles and the intersection of the two circles and the road boundary; that is, when there is an intersection between Bu and BOu, the ring will be cut. Because the uncovered new area is sufficiently narrow, the second part of the station only requires one layer of the service circle, so this paper can only use one layer of the coverage circle to complete the coverage under the condition of meeting the deployment restrictions. Under this premise, this paper proposes an improved algorithm for the empty problem. As shown in Figure 5, pi and qi are two points in the boundaries Bu and BOu, respectively, and a station needs to be deployed through these two points to cover a part of the shadow area. First, draw two arcs with pi and qi as the center and SRsi as the radius. They intersect BOu at xi and yi and intersect each other at zi. If zi is in the shaded area, then zi can be considered a candidate location for the station. Otherwise, like the points pi+1 and qi+1 in Figure 5, select a point closer to pi+1 and qi+1 from xi+1 and yi+1 as the candidate position of the circle oi+1. Keep adjusting the coverage circle until all the shadow areas are completely covered. In the improved algorithm, the principle of reducing the service radius of the station is in the same form as the algorithm proposed above. For the left part of Figure 5, first generate the convex hull of the polygon BOu and its diameter, then select one end of the diameter as the initial point q1, and select the relatively close intersection point between the diameter and Bu as the initial point p1. For the right part of Figure 5, use two initial points p1 and q1 generated by the intersection of the circle and Bu and BOu at one end of the shaded area of the incomplete ring.Figure 5
Schematic diagram of deploying new stations.


## 5. Analysis of Algorithms

### 5.1. Correctness Analysis

This algorithm can completely cover the entire convex hull under predefined restrictions. According to the main steps of the algorithm, the next circle drawn must intersect the previous circle and the convex hull at a point. Therefore, the next circle covers a part of the convex hull boundary close to the intersection. For the convex hullCHG, all covering circles are constructed one after another along the convex hull boundary to ensure seamless boundary coverage. In addition, after covering the convex hull boundary BCHG and generating the next convex hull boundary BCHG1, these circles must cover the annular area between BCHG and BCHG1. The algorithm is continuously iterated for the next convex hull boundary BCHG1, and another series of circles must cover the annular area between BCHG1 and BCHG2. This algorithm is iterated in a similar manner in turn and the convex hull CHG,CHG1,…,CHGk and CHGCHG1CHG2⋯CHGk is obtained. This sequence is a decreasing sequence because the initial convex hull CHG is limited in size, and the sequence is also limited. Assuming that the sequence is CHG1,CHG2,CHG3,…,CHGk−1,CHGk, obviously the entire area to be deployed will be divided into ABCHG,BCHG1,ABCHG1,BCHG2,…,ABCHGk−1,BCHGk in turn, and elements ABCHGi,BCHGi+1 represent the circular area between BCHGi and BCHGi+1. Any vertex in the initial convex hull CHG will be completely located in a ring area or convex hull CHGk. Without a loss of generality, assuming that the vertex is in the ring area ABCHGi,BCHGi+1, the covering circle constructed along the convex hull boundary BCHGi must completely cover the ring area; that is, this vertex must be covered by at least one covering circle. In addition, any covering circle must satisfy the constraints predefined by the algorithm. This theorem has been proven.


### 5.2. Performance Analysis

Theorem 2.
The approximate ratio of the algorithm is less than33+2/2.Proof.
Because the radius of the coverage circle is not fixed, it is almost impossible to determine the best deployment scheme, but it is still possible to evaluate the algorithm in this paper by defining the infimum of the number of stations under certain circumstances. Assuming the infimum isK, the service range limit and service load limit of a station are obtained by the following methods:(14)L=TSLK,R=TAe∗π∗K,where TA represents the total area of the city and e represents the coverage rate of the circle. When a circle is used to cover an area, there will be repeated coverage areas, so 0<e<1. In a scene covered by a circle, the most effective method is hexagonal, so the most effective coefficient is 33/2π. By generating L and R in this way, the number of stations in the optimal scheme must be greater than K, so the infimum of the number of the stations is K. Based on the service coverage update mechanism in the algorithm, this paper can ensure that the service load of the small circle with a radius less than R is as close to the load limit L as possible. The service load calculated based on real taxi traffic information must be within a relatively reasonable coverage, and the total service load is distributed across all road sections. In this paper, we can consider regular circles with radius R and small circles with a radius less than R one by one. In this case, assume that ks disjoint small circles are deployed to cover areas where the service load is high and the service load is very close to L. Obviously, the number of small circles is less than or equal to the infimum value, that is, ks≤K, and it is ensured that there are no other small circles except ks small circles in the urban area. Although ks small circles play a certain role in coverage, kn regular circles are now needed to cover the entire area. To avoid the irregularity of the convex hull, first try to find the smallest rectangle surrounding the convex hull, and then draw kn regular circles to cover the rectangle in the most effective way. Assuming that the area of the rectangle is RA, then at most kn=RA/e′∗π∗R2 regular circles are needed to cover the rectangle. According to the construction of the rectangle, this paper can ensure that the area of the rectangle is not greater than twice the total area of the city TA; that is, RA≤2∗TA. Among them, e′ is the coverage rate of the rectangle, and it is easy to calculate the coverage rate e′=2/π, so we obtain kn:(15)kn=RAe′∗π∗R2≤2∗TAe′∗π∗R2=2∗TAe′∗π∗TA/e∗π∗K2=2ee′=233/2πR22/πR2K=332K.
Therefore, by addingks small circles and kn regular circles, the total number of covered circles is(16)K=ks+kn≤K+332K=33+22K.
Since the number of swapping stations in the optimal deployment scheme is greater thanK, the approximate ratio of this algorithm is less than 33+2/2.


### 5.3. Time Complexity Analysis

Assuming that there aren road sections in the road network, the number of vertices is also n. Then, the time complexity of generating the first convex hull is Onlogn [31], and the time complexity of calculating the diameter of the convex hull is On [30]. The deployment of road sections is relatively uniform, and the time complexity of adjusting and drawing a circle is constant. Assuming that the final deployment scheme requires k circles and k≪n, then the time complexity of drawing all circles is Ok. The new convex hulls generated by connecting the intersection points of all covering circles are CHG1,CHG2,…,CHGk. Then, the time complexity of the polygon formed by connecting these intersection points is Ok, and the time complexity of generating a new convex hull is Ok. The time complexity is also Ok. Therefore, the total time complexity of the algorithm involved in this paper is Onlogn+On+Ok+Ok+Ok+Ok. Since k≪n, the time complexity of the BSSP algorithm is Onlogn.


## 5.1. Correctness Analysis

This algorithm can completely cover the entire convex hull under predefined restrictions. According to the main steps of the algorithm, the next circle drawn must intersect the previous circle and the convex hull at a point. Therefore, the next circle covers a part of the convex hull boundary close to the intersection. For the convex hullCHG, all covering circles are constructed one after another along the convex hull boundary to ensure seamless boundary coverage. In addition, after covering the convex hull boundary BCHG and generating the next convex hull boundary BCHG1, these circles must cover the annular area between BCHG and BCHG1. The algorithm is continuously iterated for the next convex hull boundary BCHG1, and another series of circles must cover the annular area between BCHG1 and BCHG2. This algorithm is iterated in a similar manner in turn and the convex hull CHG,CHG1,…,CHGk and CHGCHG1CHG2⋯CHGk is obtained. This sequence is a decreasing sequence because the initial convex hull CHG is limited in size, and the sequence is also limited. Assuming that the sequence is CHG1,CHG2,CHG3,…,CHGk−1,CHGk, obviously the entire area to be deployed will be divided into ABCHG,BCHG1,ABCHG1,BCHG2,…,ABCHGk−1,BCHGk in turn, and elements ABCHGi,BCHGi+1 represent the circular area between BCHGi and BCHGi+1. Any vertex in the initial convex hull CHG will be completely located in a ring area or convex hull CHGk. Without a loss of generality, assuming that the vertex is in the ring area ABCHGi,BCHGi+1, the covering circle constructed along the convex hull boundary BCHGi must completely cover the ring area; that is, this vertex must be covered by at least one covering circle. In addition, any covering circle must satisfy the constraints predefined by the algorithm. This theorem has been proven.


## 5.2. Performance Analysis

Theorem 2.
The approximate ratio of the algorithm is less than33+2/2.Proof.
Because the radius of the coverage circle is not fixed, it is almost impossible to determine the best deployment scheme, but it is still possible to evaluate the algorithm in this paper by defining the infimum of the number of stations under certain circumstances. Assuming the infimum isK, the service range limit and service load limit of a station are obtained by the following methods:(14)L=TSLK,R=TAe∗π∗K,where TA represents the total area of the city and e represents the coverage rate of the circle. When a circle is used to cover an area, there will be repeated coverage areas, so 0<e<1. In a scene covered by a circle, the most effective method is hexagonal, so the most effective coefficient is 33/2π. By generating L and R in this way, the number of stations in the optimal scheme must be greater than K, so the infimum of the number of the stations is K. Based on the service coverage update mechanism in the algorithm, this paper can ensure that the service load of the small circle with a radius less than R is as close to the load limit L as possible. The service load calculated based on real taxi traffic information must be within a relatively reasonable coverage, and the total service load is distributed across all road sections. In this paper, we can consider regular circles with radius R and small circles with a radius less than R one by one. In this case, assume that ks disjoint small circles are deployed to cover areas where the service load is high and the service load is very close to L. Obviously, the number of small circles is less than or equal to the infimum value, that is, ks≤K, and it is ensured that there are no other small circles except ks small circles in the urban area. Although ks small circles play a certain role in coverage, kn regular circles are now needed to cover the entire area. To avoid the irregularity of the convex hull, first try to find the smallest rectangle surrounding the convex hull, and then draw kn regular circles to cover the rectangle in the most effective way. Assuming that the area of the rectangle is RA, then at most kn=RA/e′∗π∗R2 regular circles are needed to cover the rectangle. According to the construction of the rectangle, this paper can ensure that the area of the rectangle is not greater than twice the total area of the city TA; that is, RA≤2∗TA. Among them, e′ is the coverage rate of the rectangle, and it is easy to calculate the coverage rate e′=2/π, so we obtain kn:(15)kn=RAe′∗π∗R2≤2∗TAe′∗π∗R2=2∗TAe′∗π∗TA/e∗π∗K2=2ee′=233/2πR22/πR2K=332K.
Therefore, by addingks small circles and kn regular circles, the total number of covered circles is(16)K=ks+kn≤K+332K=33+22K.
Since the number of swapping stations in the optimal deployment scheme is greater thanK, the approximate ratio of this algorithm is less than 33+2/2.


## 5.3. Time Complexity Analysis

Assuming that there aren road sections in the road network, the number of vertices is also n. Then, the time complexity of generating the first convex hull is Onlogn [31], and the time complexity of calculating the diameter of the convex hull is On [30]. The deployment of road sections is relatively uniform, and the time complexity of adjusting and drawing a circle is constant. Assuming that the final deployment scheme requires k circles and k≪n, then the time complexity of drawing all circles is Ok. The new convex hulls generated by connecting the intersection points of all covering circles are CHG1,CHG2,…,CHGk. Then, the time complexity of the polygon formed by connecting these intersection points is Ok, and the time complexity of generating a new convex hull is Ok. The time complexity is also Ok. Therefore, the total time complexity of the algorithm involved in this paper is Onlogn+On+Ok+Ok+Ok+Ok. Since k≪n, the time complexity of the BSSP algorithm is Onlogn.


## 6. Experiment and Evaluation

### 6.1. Introduction to Trajectory Data

This paper analyzes the driving trajectory and the status of passengers on and off the 3997 taxis in the Suzhou road network from March 2012 to June 2012. Among them, each taxi uploads its location to the server every 2 minutes between 7 : 00 and 20 : 00 every day.The operation trajectory data of the taxis are shown in Table1. CAR_ID represents the license plate number of the taxi; X and Y are the latitude and longitude coordinates of the taxi to the current position; SPEED represents the current speed of the taxi; DIRECTION indicates the clockwise offset angle of the current direction of the vehicle; GET_TIME is the time to obtain the data, OPERATE is whether the vehicle is currently in operation, where 1 is in operation and 0 is out of operation; and LOCATION indicates the current passenger carrying situation of the taxi, where 1 is carrying passengers and 0 is an empty load.Table 1
Data example of taxi trajectory.
CAR_IDXYSPEEDDIRECTIONGET_TIMEOPERATELOCATION2108120.647931.27949421832012-3-3 10:24:11112108120.630431.2669929752012-3-3 10:32:01002239120.619331.3015421712012-3-3 11:37:38115229120.662431.28246351742012-3-3 11:37:4600……………………


### 6.2. Data Preprocessing

To deploy the swapping stations in the central area of the city, this paper screens the trajectory data between the latitude and longitude range of119∘55′∼121∘20′ east longitude and 30∘47′∼32∘02′ north latitude and maps them to the Suzhou road network, as shown in Figure 6. The characteristics of urban traffic usually take a weekly cycle, so this paper intercepts nearly 50 million valid taxi historical trajectory data for 4 weeks from 2012/3/5 to 2012/4/1.Figure 6
Taxi running track diagram.The algorithm in this paper calculates the convex hull of the road section multiple times where the trajectory point is located. The time complexity of a single convex hull calculation isOnlogn, where n is the number of trajectory points in the convex hull. To avoid the impact of the large-scale number of trajectory points on the calculation time of the convex hull, this paper first simplified the data, as shown in Figure 7. For all points located in the coordinate system, first use the method of rotating straight lines to take the horizontal axis in the coordinate system as the starting point, rotate the straight line along the boundary point counter clockwise by 360∘, and then connect the intersection of two adjacent straight lines to form a boundary polygon. Then, only some points within the polygon that are close to the boundary of the polygon are used as candidate points to generate the convex hull. Finally, these candidate points are substituted into the Graham algorithm [26] to find the convex hull. The convex hull calculated using simplified data is shown in Figure 8.Figure 7
Schematic diagram of the simplified data.Figure 8
Convex hull formed by taxi trajectory points in the core area.


### 6.3. Algorithm Result Analysis

Assuming that the maximum service radius of the swapping station is 5 kilometers (referring to the service radius of diesel locomotive gas stations), then the mileage of the electric taxi is set to 200 km. The preprocessed data are substituted into the graph-covering algorithm to find the position of the swapping station. First, the first layer of the convex hull surrounding all road section nodes is obtained according to the Graham algorithm [26]. The convex hull of the first layer can be calculated by the coverage algorithm where 10 coverage circles can be used to complete the coverage of the outermost convex hull; that is, 10 swapping stations need to be deployed in the outermost layer. The maximum service radius of the swapping stations is 5 kilometers, and the minimum service radius is 3.239 kilometers.After the first layer of the convex hull is covered, connect the intersection of the circles, calculate and generate the second layer of the convex hull, and then continue to use the covering circle to cover the new convex hull. From the calculation results, we need to deploy 7 services with a radius of 5 kilometers to complete the coverage of the second layer of convex hull. By analogy, as shown in Figure9 it can be calculated that the minimum number of swapping stations in the main urban area of Suzhou is 21, which can meet the battery swapping demand of all taxis in the road network of Suzhou city and ensure the load balance of the swapping stations.Figure 9
The covering circle completely covering the convex hull.As shown in Figure10, under the condition that the service coverage R=5 remains unchanged and that the load amount L is changed, the number of required coverage circles shows a downward trend as the load amount changes; that is, the greater the total load amount is, the fewer the swapping stations that are required.Figure 10
The relationship between the coverage circle and load.The following is a detailed analysis of the main factors that affect the deployment of swapping stations. Due to the limitation of the mileage of electric taxis, the swapping station built on a single road node may not be able to meet the battery swapping demand of the taxi to and from the daily operation road; that is, it cannot fully capture the traffic flow on the road section. With the continuous increase in the mileage of electric taxis, the demand for battery swapping generated by taxis will decrease.Figure11 shows that there is a positive correlation between the mileage of electric taxis and the number of swapping stations. The greater the mileage of an electric taxi, the more vehicles that can be serviced by the swap station. At this time, fewer swap stations can be deployed to meet the battery swapping demand in the road network. When the mileage of electric taxis increases from 100 km to 800 km, the number of swapping stations drops from 105 to 2, indicating that the mileage of taxis has a greater impact on the number of swapping stations required. Due to technical limitations, in the current situation the mileage of electric taxis is difficult to increase significantly, so this paper can be adjusted by increasing the swapping station density.Figure 11
The relationship between the range of electric taxis and the number of stations and traffic flow.Figure12(a) shows that the number of swapping stations in the early stage is positively correlated with the utilization rate of the swapping stations. When the number of swapping stations is increased to 80, the utilization rate of the swapping stations is close to 97%; that is, the swapping stations are fully utilized. However, as the number of swapping stations continues to increase in the later period, the utilization rate of the swapping stations gradually decreases, indicating that the number of deployed swapping stations is too redundant, which affects the utilization rate of the swapping stations. The size of the service coverage of the swapping station is related to whether the swapping station can be effectively utilized. As shown in Figure 12(b), the service coverage of the swapping station is close to 200 km, and the utilization rate of the swapping station is close to 96%. Therefore, the service coverage of the swapping station should be reasonably determined when deploying the station.Figure 12
The relationship between the utilization rate of the substation and the number of stations and the service scope R. (a) Number of stations. (b) Service scope of a single station (R).
(a)(b)To further prove the advantages of the deployment scheme of this paper, this paper selects a set of bit data to test the unified deployment scheme and the load balancing deployment scheme of this paper and then compares and analyzes the four aspects of the number of swap stations that need to be deployed, the load of the swap stations, the balance of the load, and the coverage of the service.By calculating the mean square error of the swapping station load in each deployment scheme and then analyzing it, through calculations the mean square error of the swapping station load in the load balancing deployment scheme is 460, and the mean square error of the swapping station load in the unified deployment scheme is 530. This shows that the deployed swapping station has a good overlap with its service load and that unified deployment scheme is decentralized from the deployment scheme adopted in this paper. Compared with the unified deployment scheme, the load mean square error is reduced by approximately 13%. This also shows that the load balancing deployment scheme described in this paper is more stable during the deployment of the swap station, so that the deployment result can ensure the consistency of the distribution of the swap station and its load while ensuring the load balance of a single swap station.It can be seen from Figure13(a) that when the service coverage R of the swapping station is fixed, with an increase in the load, the load balancing scheme is better than the unified deployment scheme. As can be seen from Figure 13(b), under the premise that the service load L of the swapping station is fixed, with an increase in the service coverage, the load balancing scheme is still better than the unified deployment scheme; under the premise of a small load variance, from Figures 13(c) and 13(d), it can be seen that, with an increase in the service coverage and service load, the number of swapping stations required by the load balancing deployment scheme is less than the number of swapping stations required by the unified deployment scheme.Figure 13
Comparative analysis of the two deployment schemes. (a) Relationship between upper limit of service load L and the mean square error. (b) Relationship between upper limit of service scope R and the mean square error. (c) Relationship between upper limit of service scope R and number of stations. (d) Relationship between upper limit of service load L and number of stations.
(a)(b)(c)(d)Under the premise of a certain service load level of the swapping station, the number of battery swaps is inversely proportional to the mileage of electric taxis. The greater the mileage of the taxi, the fewer the battery swaps that are required during driving. The number of battery swaps is directly proportional to the traffic flow in the road network. The greater the traffic flow in the road network, the greater the demand for battery swapping, which increases the number of battery swaps.Figure14 shows that, under the premise of a certain service load of the swapping station, the number of battery swaps of a taxi is closely related to the maximum mileage and traffic flow. In the case that it is difficult to increase mileage, the reasonable deployment of swapping stations will become increasingly important.Figure 14
The relationship between the number of battery exchanges, mileage, and traffic flow. (a) The relationship between the number of battery exchanges and the mileage. (b) The relationship between the number of battery exchanges and traffic flow.
(a)(b)


## 6.1. Introduction to Trajectory Data

This paper analyzes the driving trajectory and the status of passengers on and off the 3997 taxis in the Suzhou road network from March 2012 to June 2012. Among them, each taxi uploads its location to the server every 2 minutes between 7 : 00 and 20 : 00 every day.The operation trajectory data of the taxis are shown in Table1. CAR_ID represents the license plate number of the taxi; X and Y are the latitude and longitude coordinates of the taxi to the current position; SPEED represents the current speed of the taxi; DIRECTION indicates the clockwise offset angle of the current direction of the vehicle; GET_TIME is the time to obtain the data, OPERATE is whether the vehicle is currently in operation, where 1 is in operation and 0 is out of operation; and LOCATION indicates the current passenger carrying situation of the taxi, where 1 is carrying passengers and 0 is an empty load.Table 1
Data example of taxi trajectory.
CAR_IDXYSPEEDDIRECTIONGET_TIMEOPERATELOCATION2108120.647931.27949421832012-3-3 10:24:11112108120.630431.2669929752012-3-3 10:32:01002239120.619331.3015421712012-3-3 11:37:38115229120.662431.28246351742012-3-3 11:37:4600……………………


## 6.2. Data Preprocessing

To deploy the swapping stations in the central area of the city, this paper screens the trajectory data between the latitude and longitude range of119∘55′∼121∘20′ east longitude and 30∘47′∼32∘02′ north latitude and maps them to the Suzhou road network, as shown in Figure 6. The characteristics of urban traffic usually take a weekly cycle, so this paper intercepts nearly 50 million valid taxi historical trajectory data for 4 weeks from 2012/3/5 to 2012/4/1.Figure 6
Taxi running track diagram.The algorithm in this paper calculates the convex hull of the road section multiple times where the trajectory point is located. The time complexity of a single convex hull calculation isOnlogn, where n is the number of trajectory points in the convex hull. To avoid the impact of the large-scale number of trajectory points on the calculation time of the convex hull, this paper first simplified the data, as shown in Figure 7. For all points located in the coordinate system, first use the method of rotating straight lines to take the horizontal axis in the coordinate system as the starting point, rotate the straight line along the boundary point counter clockwise by 360∘, and then connect the intersection of two adjacent straight lines to form a boundary polygon. Then, only some points within the polygon that are close to the boundary of the polygon are used as candidate points to generate the convex hull. Finally, these candidate points are substituted into the Graham algorithm [26] to find the convex hull. The convex hull calculated using simplified data is shown in Figure 8.Figure 7
Schematic diagram of the simplified data.Figure 8
Convex hull formed by taxi trajectory points in the core area.


## 6.3. Algorithm Result Analysis

Assuming that the maximum service radius of the swapping station is 5 kilometers (referring to the service radius of diesel locomotive gas stations), then the mileage of the electric taxi is set to 200 km. The preprocessed data are substituted into the graph-covering algorithm to find the position of the swapping station. First, the first layer of the convex hull surrounding all road section nodes is obtained according to the Graham algorithm [26]. The convex hull of the first layer can be calculated by the coverage algorithm where 10 coverage circles can be used to complete the coverage of the outermost convex hull; that is, 10 swapping stations need to be deployed in the outermost layer. The maximum service radius of the swapping stations is 5 kilometers, and the minimum service radius is 3.239 kilometers.After the first layer of the convex hull is covered, connect the intersection of the circles, calculate and generate the second layer of the convex hull, and then continue to use the covering circle to cover the new convex hull. From the calculation results, we need to deploy 7 services with a radius of 5 kilometers to complete the coverage of the second layer of convex hull. By analogy, as shown in Figure9 it can be calculated that the minimum number of swapping stations in the main urban area of Suzhou is 21, which can meet the battery swapping demand of all taxis in the road network of Suzhou city and ensure the load balance of the swapping stations.Figure 9
The covering circle completely covering the convex hull.As shown in Figure10, under the condition that the service coverage R=5 remains unchanged and that the load amount L is changed, the number of required coverage circles shows a downward trend as the load amount changes; that is, the greater the total load amount is, the fewer the swapping stations that are required.Figure 10
The relationship between the coverage circle and load.The following is a detailed analysis of the main factors that affect the deployment of swapping stations. Due to the limitation of the mileage of electric taxis, the swapping station built on a single road node may not be able to meet the battery swapping demand of the taxi to and from the daily operation road; that is, it cannot fully capture the traffic flow on the road section. With the continuous increase in the mileage of electric taxis, the demand for battery swapping generated by taxis will decrease.Figure11 shows that there is a positive correlation between the mileage of electric taxis and the number of swapping stations. The greater the mileage of an electric taxi, the more vehicles that can be serviced by the swap station. At this time, fewer swap stations can be deployed to meet the battery swapping demand in the road network. When the mileage of electric taxis increases from 100 km to 800 km, the number of swapping stations drops from 105 to 2, indicating that the mileage of taxis has a greater impact on the number of swapping stations required. Due to technical limitations, in the current situation the mileage of electric taxis is difficult to increase significantly, so this paper can be adjusted by increasing the swapping station density.Figure 11
The relationship between the range of electric taxis and the number of stations and traffic flow.Figure12(a) shows that the number of swapping stations in the early stage is positively correlated with the utilization rate of the swapping stations. When the number of swapping stations is increased to 80, the utilization rate of the swapping stations is close to 97%; that is, the swapping stations are fully utilized. However, as the number of swapping stations continues to increase in the later period, the utilization rate of the swapping stations gradually decreases, indicating that the number of deployed swapping stations is too redundant, which affects the utilization rate of the swapping stations. The size of the service coverage of the swapping station is related to whether the swapping station can be effectively utilized. As shown in Figure 12(b), the service coverage of the swapping station is close to 200 km, and the utilization rate of the swapping station is close to 96%. Therefore, the service coverage of the swapping station should be reasonably determined when deploying the station.Figure 12
The relationship between the utilization rate of the substation and the number of stations and the service scope R. (a) Number of stations. (b) Service scope of a single station (R).
(a)(b)To further prove the advantages of the deployment scheme of this paper, this paper selects a set of bit data to test the unified deployment scheme and the load balancing deployment scheme of this paper and then compares and analyzes the four aspects of the number of swap stations that need to be deployed, the load of the swap stations, the balance of the load, and the coverage of the service.By calculating the mean square error of the swapping station load in each deployment scheme and then analyzing it, through calculations the mean square error of the swapping station load in the load balancing deployment scheme is 460, and the mean square error of the swapping station load in the unified deployment scheme is 530. This shows that the deployed swapping station has a good overlap with its service load and that unified deployment scheme is decentralized from the deployment scheme adopted in this paper. Compared with the unified deployment scheme, the load mean square error is reduced by approximately 13%. This also shows that the load balancing deployment scheme described in this paper is more stable during the deployment of the swap station, so that the deployment result can ensure the consistency of the distribution of the swap station and its load while ensuring the load balance of a single swap station.It can be seen from Figure13(a) that when the service coverage R of the swapping station is fixed, with an increase in the load, the load balancing scheme is better than the unified deployment scheme. As can be seen from Figure 13(b), under the premise that the service load L of the swapping station is fixed, with an increase in the service coverage, the load balancing scheme is still better than the unified deployment scheme; under the premise of a small load variance, from Figures 13(c) and 13(d), it can be seen that, with an increase in the service coverage and service load, the number of swapping stations required by the load balancing deployment scheme is less than the number of swapping stations required by the unified deployment scheme.Figure 13
Comparative analysis of the two deployment schemes. (a) Relationship between upper limit of service load L and the mean square error. (b) Relationship between upper limit of service scope R and the mean square error. (c) Relationship between upper limit of service scope R and number of stations. (d) Relationship between upper limit of service load L and number of stations.
(a)(b)(c)(d)Under the premise of a certain service load level of the swapping station, the number of battery swaps is inversely proportional to the mileage of electric taxis. The greater the mileage of the taxi, the fewer the battery swaps that are required during driving. The number of battery swaps is directly proportional to the traffic flow in the road network. The greater the traffic flow in the road network, the greater the demand for battery swapping, which increases the number of battery swaps.Figure14 shows that, under the premise of a certain service load of the swapping station, the number of battery swaps of a taxi is closely related to the maximum mileage and traffic flow. In the case that it is difficult to increase mileage, the reasonable deployment of swapping stations will become increasingly important.Figure 14
The relationship between the number of battery exchanges, mileage, and traffic flow. (a) The relationship between the number of battery exchanges and the mileage. (b) The relationship between the number of battery exchanges and traffic flow.
(a)(b)


## 7. Conclusion and Outlook

This paper took the deployment of electric taxi swap stations in the main urban areas of Suzhou as the research object, and the related theoretical principles, such as traffic flow theory and graph-covering theory, were comprehensively applied to battery studies. Considering the service coverage and service load of a single station, the influencing factors of traffic load were analyzed, where load balancing was the main optimization goal. A graph-covering algorithm based on computational geometry was designed to solve the model. According to the experimental results, the effects of service coverage, mileage, and traffic flow on the number of stations and station utilization were analyzed, and the final comparison of the unified deployment scheme showed that the algorithm has certain advantages in the deployment of electric taxi stations. We will use gas stations and other urban infrastructure to deploy charging stations to reduce deployment costs. Future research can deeply analyze the deployment problems of swapping stations in urban empty areas to adjust the deployment scheme and reduce the invalidity of the algorithm.



---

*Source: 4014837-2021-09-16.xml*