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.
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. Section 4 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 graph to complete station deployment. First, the service load capacity of any swapping station should be close to, but not exceed, the maximum load limit ; 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 . The service load of swapping station comes from the total battery swapping demand of taxis traveling on all road sections within its service radius . 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 graph , find the smallest station set so that the service coverage of the stations in the set can completely cover the convex hull of the vertex set . 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, Theorem 1 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 an lattice in the traffic network graph , where the distance between adjacent points is , and the lattice just covers the convex hull of . 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 , the service load is limited to , and the traffic load of section is . The undirected graph can then be obtained through the following steps: Step 1: initialize . For each vertex in , find the nearest point in the lattice. If , add to . Step 2: for and , if the Euclidean distance and the total traffic load of , then edge is added to edge set . Step 3: for , set the vertex in to , for , if , add to the edge set . Step 4: for , add to in turn.Then, the covering problem is transformed based on the undirected graph , finding the minimum vertex subset from , where any edge of has at least one vertex in . 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 section in the road network on day as . By using the Kolmogorov-Smirnov test (K-S test), reference [28] proposes that the time statistics of vehicles arriving at road section are similar to the exponential distribution; that is, the traffic load statistics of road section in a fixed time period obey a Poisson distribution. Therefore, the probability distribution function of iswhere represents the average traffic load statistics of road section in a day and represents the incidence of taxis randomly arriving at road section in a day. is calculated using the trajectory data of all taxis in Suzhou in 366 days in 2012:
Then, calculate and count the average traffic load of road section by the following formula:
2.2.2. Service Load Analysis of the Swapping Station
As mentioned earlier, the service load of swapping station is determined by the total traffic load of all sections within the service radius, namely,
Although the total traffic load cannot be directly used as the service load of the swapping station, because they are still linearly proportional, this paper defines the service load of the swapping station aswhere 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 is the number of taxis. Assuming that each taxi needs to change electricity times a day on average, can be calculated as
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 section is abstracted as node to obtain road section set . If and intersect, then connect and to generate edge , thereby obtaining the road network topology graph , which is located in the latitude and longitude coordinate system.
Define the taxi set as , where represents the -th electric taxi; the trajectory of taxi on day is , and represents a certain road section in the road network; represents the frequency of taxi passing through section on day , and the traffic load of section on day is . The objective function and constraint conditions based on the above conditions are as follows:where is the total number of deployed swapping stations; is the variance function; indicates that the service radius of swapping station is less than or equal to the maximum service radius of the station; indicates that the service load of swapping station is less than or equal to the maximum load of the station; is the total service load of all swapping stations; and is the Euclidean distance between swapping station and section . 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 hull as the smallest convex set covering the road section set in graph . Its boundary is a convex polygon whose vertices are all from the road section set , 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 and its diameter of the road network 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 .
|
3.2. Convex Hull Covering Algorithm Steps
Step 1: initialize , , call the Graham algorithm [29] to obtain the convex hull , and find the convex hull boundary . Step 2: call the Preparata algorithm [30] according to the convex hull to obtain the convex hull diameter . Step 3: generate the first covering circle around the convex hull. Take the point on the convex hull diameter to maximize and . The service load is calculated using formulas (4) and (5), and and represent the service radius and candidate station of , respectively. According to the properties of the convex hull, it can be proven that there must be two intersection points between the circle and the boundary of the convex hull. Let these two intersection points be and , respectively. The specific process of generating the first covering circle is as follows. Step 4: if the recently drawn circle is , point and point are the intersection points of the circle and the convex hull boundary , respectively. Draw the next circle by using the following iterative algorithm. Step 4.1: initialize . Step 4.2: calculate the position of the center of the next circle according to the service radius . The detailed process is as follows. Assuming that the newly drawn circle is , which intersects the convex hull boundary at and intersects the circle to be drawn at and ; then the center discriminant equation (8) can be obtained. However, knowing only the coordinates of point and the service radius , an infinite number of solutions for the center can be obtained. 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 center . 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 is defined and the arc corresponding to the arch is a minor arc, is the radius of the circle, and is the chord length of the arch, then the area of an arch is calculated, as shown in the following equation: 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 . 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 . Under this premise, when the chord length is a given constant, then the final arch area decreases monotonically with respect to . To simplify the calculation, this paper approximates the problem of minimizing the shadow area into a problem of minimizing the sum of chord length . For PART I in Figure 1, the judgment formula of minimization is adopted: Among them, point and point are fixed, and point and point are determined by the next circle. This paper derives formula (10) for minimizing , 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 circle is , the radius is , the center of the circle is , the radius is , and the linear equation of the chord is (where the convex hull boundary intersecting the circle is ), obviously the linear equation of the common chord formed by the intersection of the two circles is . Among them, the linear equation where the common chord is located is obtained by subtracting the circle and the circle . Suppose the equation of the line where the chord is located is . First, according to the distance formula from the point to the straight line, the distance from the center of the circle to the line is , the distance from the center of the circle to the common chord is , and the distance from the center to the distance of the straight line is . 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 , , . Then, the total chord length is obtained. Finally, the obtained three chord lengths are substituted into the LA minimization judgment formula (10) to obtain the minimum chord length. Figure 1 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 and line segment . 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 , so equation (11) is used to simplify this special case. In Figure 1, if circle is the second circle in the convex hull, then circle and shadow area I do not exist. At this time, this paper cannot directly calculate 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 is as follows: 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 centers and is greater than the sum of their radii to meet the following conditions: Step 4.3: if swapping station is deployed in , then its current service radius and service load are , and the service radius is updated through iteration (13), where is the maximum load of the swapping station. Obviously, the service radius can be continuously updated according to the relationship between the current service load and , but the update process may cross an excessively large step and continue to iterate at this time. If is satisfied, where is a fixed constant close to 1, then swapping station is determined to be deployed in . Otherwise, , so jump to Step 4.2. Circle intersects the convex hull boundary at and and intersects circle at and . As shown in Figure 2, Step 4 is repeated until , where is the intersection of circle and convex hull boundary . Circle intersects the convex hull boundary at and and intersects circle at and . As shown in Figure 2, Step 4 is repeated until , where is the intersection of circle and convex hull boundary . Step 5: draw a circle with a radius and pass through points and . Suppose that a swapping station is deployed at and calculate the service load of the station. If , then is a candidate location for the deployment of swapping station . Otherwise, call an iterative program similar to Step 4 to adjust the position of the service radius and the center . The specific process of generating the first and next covering circles is shown in Algorithms 2 and 3.
|
|
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. Figure 3 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 is greater than , then the circle intersects the circle , 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 be smaller than , that is, , so that no such intersection will occur in order to ensure that the inner convex hull can be generated normally.
(a)
(b)
Since the service radius of each covered circle is not the same, the polygon generated by sequentially connecting the intersections 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 graph is completely covered, then all the vertices in set are covered by the circle centered on the swapping station, and is guaranteed. The algorithm for deploying the last swapping station is shown in Algorithm 4.
|
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 graph to draw a real loop around an empty area. Suppose that the set represents all of the intersections in the road network, and then change to , where , and side indicates that one end of the road section is the intersection . 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 and the punctuation of the empty area.
Assuming that the punctuation of an empty area is , initialize the loop set , and then search for the vertex closest to the punctuation in set or set . Without loss of generality, assuming that the vertex is or , the vertex must be added to the loop set , and then all the intersections adjacent to the vertex are searched to find the intersection that is closest to the punctuation . If vertex is not unique among the adjacent road sections of , that is, if is not the end point, continue to add to loop set and search again for all road sections adjacent to vertex . Otherwise, take the intersection that is not the end point and the closest to the loop as the next intersection of . For intersection , an intersection in adjacent road sections is the intersection that is closest to and connects at least one non-end point (except ) and it is added to . 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, is the boundary of the empty area in the loop. If the initial vertex is , add the intersections and road sections to 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 station 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 . According to the proposed algorithm, there must be candidate stations that are adjacent to 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 or service load limit , 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 Figure 4, the inner polygon is the loop boundary of area , and the adjacency circle is the available station adjacent to but not intersecting with . By connecting the intersection points between the adjacent circles, the outer boundary polygon of the shaded area is generated. The middle area between and is a new area that needs to be covered again, and its width must be less than the service radius limit .
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 Figure 4(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 and , 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, and are two points in the boundaries and , 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 and as the center and as the radius. They intersect at and and intersect each other at . If is in the shaded area, then can be considered a candidate location for the station. Otherwise, like the points and in Figure 5, select a point closer to and from and as the candidate position of the circle . 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 and its diameter, then select one end of the diameter as the initial point , and select the relatively close intersection point between the diameter and as the initial point . For the right part of Figure 5, use two initial points and generated by the intersection of the circle and and at one end of the shaded area of the incomplete ring.
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 hull , 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 and generating the next convex hull boundary , these circles must cover the annular area between and . The algorithm is continuously iterated for the next convex hull boundary , and another series of circles must cover the annular area between and . This algorithm is iterated in a similar manner in turn and the convex hull and is obtained. This sequence is a decreasing sequence because the initial convex hull is limited in size, and the sequence is also limited. Assuming that the sequence is , obviously the entire area to be deployed will be divided into in turn, and elements represent the circular area between and . Any vertex in the initial convex hull will be completely located in a ring area or convex hull . Without a loss of generality, assuming that the vertex is in the ring area , the covering circle constructed along the convex hull boundary 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 than .
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 is , the service range limit and service load limit of a station are obtained by the following methods:where represents the total area of the city and represents the coverage rate of the circle. When a circle is used to cover an area, there will be repeated coverage areas, so . In a scene covered by a circle, the most effective method is hexagonal, so the most effective coefficient is . By generating and in this way, the number of stations in the optimal scheme must be greater than , so the infimum of the number of the stations is . 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 is as close to the load limit 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 and small circles with a radius less than one by one. In this case, assume that disjoint small circles are deployed to cover areas where the service load is high and the service load is very close to . Obviously, the number of small circles is less than or equal to the infimum value, that is, , and it is ensured that there are no other small circles except small circles in the urban area. Although small circles play a certain role in coverage, 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 regular circles to cover the rectangle in the most effective way. Assuming that the area of the rectangle is , then at most 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 ; that is, . Among them, is the coverage rate of the rectangle, and it is easy to calculate the coverage rate , so we obtain :Therefore, by adding small circles and regular circles, the total number of covered circles isSince the number of swapping stations in the optimal deployment scheme is greater than , the approximate ratio of this algorithm is less than .
5.3. Time Complexity Analysis
Assuming that there are road sections in the road network, the number of vertices is also . Then, the time complexity of generating the first convex hull is [31], and the time complexity of calculating the diameter of the convex hull is [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 circles and , then the time complexity of drawing all circles is . The new convex hulls generated by connecting the intersection points of all covering circles are . Then, the time complexity of the polygon formed by connecting these intersection points is , and the time complexity of generating a new convex hull is . The time complexity is also . Therefore, the total time complexity of the algorithm involved in this paper is . Since , the time complexity of the BSSP algorithm is .
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 Table 1. 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.
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 of east longitude and 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.
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 is , where 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 , 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.
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 Figure 9 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.
As shown in Figure 10, under the condition that the service coverage remains unchanged and that the load amount 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.
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.
Figure 11 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 12(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.
(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 Figure 13(a) that when the service coverage 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 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.
(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.
Figure 14 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.
(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.
Data Availability
Suzhou Taxi trace data cannot be used publicly. Part of the research data can be acquired from the corresponding author through e-mail: [email protected].
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This project was supported by the National Social Science Fund of China (16BTQ086).