# Bike-Sharing Static Rebalancing by Considering the Collection of Bicycles in Need of Repair

Authors: Sheng Zhang; Guanhua Xiang; Zhongxiang Huang
Journal: Journal of Advanced Transportation
 (2018)

---

## Abstract

Bike-sharing systems, which are used in many cities worldwide, need to maintain a balance between the availability of bicycles and the availability of unoccupied bicycle slots. This paper presents an investigation of the net flow of each bike-sharing station in Jersey City. The data was recorded at 1-minute intervals. The sum of the initial bicycle number and the minimum net flow value was determined to be the demand for static rebalancing, and this led to the proposal of a bike-sharing demand prediction method based on autoregressive integrated moving average models. Considering that the existence of bicycles in a state of disrepair may adversely affect demand prediction and routine planning, we present an integer linear programming formulation to model bike-sharing static rebalancing. The proposed formulation takes into account the problem introduced by the need to collect bicycles in need of repair. A hybrid Discrete Particle Swarm Optimization (DPSO) algorithm was proposed to solve the model, which incorporates a reduced variable neighborhood search (RVNS) functionality together with DPSO to improve the global optimization performance. The effectiveness of the algorithms was verified by a detailed numerical example.


---

## Body

## 1. Introduction

Bike-sharing systems are designed to solve “the first and the last mile problem” and to provide a connection between multitransit modes. The most recent decade has witnessed the rapid worldwide development of bicycles as a transport mode driven by bike-sharing companies, such as Mobike and OFO. In this arrangement, the residents can rent or return bicycles at any public area, and this has greatly enhanced the convenience in the residents’ daily commutes. These systems have since spread all over the world. In contrast, the rapid development of bike-sharing has resulted in many problems, such as damaged bicycles, bike parking needs, and bicycle maintenance. In many cities, piles of bicycles are parked along public areas, and broken bicycles tend to be left carelessly. Among these violations, the parking problem is especially critical and has aroused people’s concern. These phenomena occur mainly because the operation management is not standardized and there is no effective rebalancing method.Scholars and policy makers are committed to determining the causes and solution of these problems. Reinforcing the standardization of the bike-sharing system management is essential; however, more importantly, advanced methods should be used to rebalance the stations to their optimum levels of occupation to avoid the inefficiency due to full and/or empty stations. Today’s competitive environment also requires bike-sharing companies to increase the operational efficiency by adopting such measures.Bike-sharing systems are composed of several bike stations located at different sites across the city, and each station operates a number of bike slots from which bicycles can be collected or to which they can be returned. In a balanced situation, each bike station is required to have a certain number of empty slots, to allow arrivals, and a certain number of full slots, to allow departures. A certain amount of time after the beginning of the service, the users will have moved the bicycles among the stations of the system and the state will have deviated from the balanced one. Thus, these companies engage in the redistribution of bicycles, an operation known as rebalancing. Rebalancing, which is performed by vehicles with the necessary capacity, is normally required at the end of the day, when the system is closed or when the use of the system can be considered negligible. In this case, the rebalancing is considered to be static. Some bike-sharing operators may require the rebalancing to be performed when the system is open and in operation, incurring a situation known as dynamic rebalancing. In this paper, we present the static case, which is referred to as the Static Bike-sharing Rebalancing Problem (SBRP).As discussed elsewhere [1–4], SBRP is a vehicle route problem, modeled on a complete digraph G=(V,A), where V={0,1,…,n} is the set of vertices, including the depot (vertex 0) andn stations (vertices 1,…,n), andA is the set of arcs between each pair of vertices. The traveling costcij is associated with the arc (i,j)∈A. For each vertexi∈V, a requestqi is given, withq0=0. Requests can be either positive or negative. Ifqi≥0, theni is a pickup station, from whichqibicycles should be removed; ifqi<0, theni is a delivery station, to whichqibicycles should be supplied, fori∈V\{0}. Bicycles supplied to delivery stations can either come from the depot or from pickup stations. Vehicles do not necessarily have to leave the depot empty. The objective is to operate a fleet ofm identical vehicles of capacityQ that need to be available at the depot to respond to requests and to minimize the total cost of the traversed arcs.However, the existence of bicycles in disrepair would disrupt the demand prediction and rebalancing efforts. In reality, after a long period of usage, some of the bicycles would require repairs. If the bicycles in need of repair continue to remain in a station indefinitely, they may be mistaken for operational bicycles by the system operator who would not deliver new bicycles to that station. Meanwhile, the customer’s demand at that station cannot be satisfied. Thus, the vehicles would need to collect these bicycles in disrepair and return them to the depot, thereby increasing the complexity of the static rebalancing problem. Contrary to the traditional vehicle routing problem with simultaneous pickup and delivery (VRPSPD), the number of bicycles in need of repair collected from some stations cannot be used to respond to requests of delivery stations. As these bicycles are expected to take up a particular amount of capacity of a vehicle, the path of the vehicle would need to be newly planned and rescheduled.This study makes the following substantial contributions to the bike-sharing static rebalancing problem:(1)
We conducted a pattern analysis of each bike-sharing station over the data set of Jersey City. On the basis of high-resolution net flow data analysis, we found that there is always a fixed input/output pattern in each bicycle station. Consequently, the bicycle stations can be divided into two types: demand-station and supply-station. We also found that the bicycle demand is determined by the minimum net flow value. For a demand-station, the bicycle demand is equal to the current number of bicycles plus the minimum net flow value; for supply-station, the supply of station is equal to the current number of bicycles plus the total difference (net flow value at the end of the day). Therefore, we established the new demand prediction equations based on bike-sharing station pattern analysis. These equations are easy to implement in practice.(2)
To the best of our knowledge, our work is the first to formulate the bike-sharing static rebalancing problem by considering the need to collect bicycles requiring repair. The existence of these bicycles increases the complexity of the static rebalancing problem. The number of bicycles requiring repair collected from some stations cannot be used to respond to requests of delivery stations. The bicycles in disrepair will take up a particular amount of capacity of a vehicle. We build new integer linear programming formulations for the bike-sharing static rebalancing problem by considering the need to collect bicycles requiring repair. The presented model is expected to ensure the specific constraint that the sum of requests of the visited stations and the bicycles in disrepair plus the initial load is never negative or greater than the capacity along the route traveled by a vehicle.The remainder of this paper is organized as follows. Previous related research is reviewed in Section2. In Section 3 we briefly present the real-world case study conducted in Jersey City, by analyzing the travel flows, and propose the station classification model based on state curve analysis. We then propose a demand prediction method based on the station classification method and autoregressive integrated moving average models. In Section 4, we consider the problem presented by the need to collect bicycles in disrepair and formulate a new bike-sharing static rebalancing model. In Section 5, we propose the hybrid algorithm of DPSO-VNS to solve the proposed model. In Section 6, a case study is presented to show the solution procedure. Conclusions and future work are presented in Section 7.


## 2. Literature Review

Lately, the bike-sharing rebalancing problem has received increased attention from the practitioners as well as researchers because of its great importance in the bike-sharing system operation. Many researchers viewed the SBRP as a special Traveling Salesman Problem (TSP) [3], multiple-TSP [5], Pickup and Delivery Traveling Salesman Problem (PDTSP) [1, 2, 6–11], and Vehicle Routing Problem with Simultaneous Pickup and Delivery (VRPSPD) [12], or built mixed integer linear programming (MILP) formulations with different objections or constrains [13–18].Table1 summarizes the static BRP (Bike Rebalancing Problem) publications according to the formulations, problem objectives, solutions, determination of demand, and whether the collection of damaged bicycles is considered in the formulations. Only a few studies have presented the bicycle demand determination method. To the best of our knowledge, only Alvarez-Valdes et al. [2] studied the return of damaged bicycles, but the impact of damaged bicycles to the vehicle’s capacity constraint has not been considered in the published literature. There are three main solution methods, e.g., exact method [1, 3, 9], heuristics method [2, 10, 12, 13, 15–17], and hybrids of these two methods [8, 14, 18].Table 1
Summary of the characteristics of the static BRP in the literature.
References Bicycle demand Formulations Objective Damaged bicycles Solution Benchimol et al.J. H. Lin et al. Given Single-vehicle 1-PDTSP Minimize total travel cost Without Consideration Approximation algorithm  Raviv et al. A decision variable Extends the 1-PDTSP Minimize the weighted sum of total travel time and penalty cost Without Consideration A two-stage heuristic  Chemla et al. Dell’Amico et al. Given Extends the 1-PDTSP Minimize total travel cost Without consideration Branch-and-cut algorithm  Ho and Szeto. A decision variable A single-vehicle MIP Minimize total penalty cost Without consideration Iterated tabu search heuristic  Angeloudis et al. Determined by the history net flow data m-TSP Minimize the total duration of travel Without Consideration. Approximation algorithm  Forma et al. A decision variable 3-step mathematical programming-based heuristic Minimize total cost, consisting of the sum of penalties and the operating costs. Without Consideration. 3-step math heuristic  Alvarez-Valdes et al. A decision variable Multi-commodity PDTSP Minimize the weighted sum of the total service time and the coefficient of variations of the duration of all routes Yes, but the impact of the damaged bicycles to the vehicle’scapacity constraint is not considered A heuristic algorithm based on minimum cost flow problem and an insertion algorithm  Szeto et al. Given Single-vehicle 1-PDTSP Minimize the weighted sum of unmet customer demand and operational time on the vehicle route. Without Consideration Enhanced Chemical Reaction Optimization  Li et al. Given MIP Minimize total cost Without Consideration. A combined hybrid genetic algorithm  Kadri et al. Given TSP Minimize the total waiting time of the stations Without consideration Branch and bound algorithm  Cruz et al. Given 1-PDTSP Minimize the traveling cost Without Consideration. Iterated Local Research-based heuristic.  Ho and Szeto. Given 1-VRPSPD Minimize the sum of penalty cost Without Consideration. A hybrid large neighborhood search combined with a tabu search.  Schuijbroek et al. Determine service level requirements at each station A Clustered MIP Minimize maximum tour length Without consideration Cluster-first route-second heuristic  Szeto and Shui. Given MIP Minimize the total demand dissatisfaction (TDD) and total service time Without consideration An enhanced artificial bee colony algorithm  Bulhões et al. Given MIP Minimize the total travel time Without consideration Branch-and-cut algorithm combined with local search metaheuristic  This study Present a predict method using the operating data Multi-commodity VRPSPD Minimize the traveling cost Yes, the vehicle’scapacity constraint to the damaged bicycles is considered A hybrid DPSO-VNSTo solve a problem such as the VRPSPD, the general method is to establish a mixed integer linear programming (MILP) formulation with different objectives or constrains. Our model for the bike-sharing static problem extends the multicommodity VRPSPD [19, 20]. In the proposed model, the vehicles are required not only to deliver goods to customer, but also to pick up goods at customer location [21].Currently, there is no general algorithm to find the global optimal solution for this kind of nonconvex and nonlinear programming problems. Heuristic algorithm can get high-quality solutions, with good computational efficiency. We have searched the literature which uses heuristic and metaheuristic approaches solving such problems, such as Ant Colony System algorithm, D-Ants, and local search algorithm. Obviously, the static rebalancing problem is also a combinational optimization problem, and the solution space is a discrete space made up of a limited number of integral points. Hence, Discrete Particle Swarm Optimization (DPSO) is applicable to solve this kind of problem. In this paper, we adopt a hybridization algorithm combining DPSO and Variable Neighborhood Search (VNS) to solve the VRPSPD, which is proposed by F. P. Goksal et al. [21] in order to enhance the global search ability of DPSO and overcome the possibility of entrapment in local optimum and make full use of the local search ability of variable neighborhood search.


## 3. Demand Prediction

### 3.1. Pattern Analysis of the State of a Station

Most existing methods to predict the short-time demand of the stations are based on the historical pattern of the bicycle number change at each station. Andreas et al. [22] detected the temporal and geographic mobility patterns based on the data from the operator’s website. These patterns were applied to the prediction of the available bicycle number at any station.After a fairly long period since a bicycle sharing system operated at an area, the resident will develop a stationary trip habit gradually. And the trip habits will tend to be stable.There will be a state periodic similarity at each bike-sharing station. That is to say, there is a fixed input/output pattern in bicycle station. We utilize operational data to analyze and find this input/output pattern. On the basis of the detailed data analysis we find that there are mainly two types of bicycle station: demand and supply-station, and the bicycle demand is determined by the minimum net flow value. For a demand-station and a demand-supply-station, the bicycle demand is equal to the current number of bicycles plus the minimum net flow value; for supply-station, the supply of station is equal to the current number of bicycles plus the total difference (net flow value at the end of the day).One of the advantages is that forecasting the bicycle demand through analyzing the pattern of each station is more precise than other methods, such as traffic mode split forecast. Currently, the operational data is adequate. With the accumulation of operational data, the demand pattern recognition of each bicycle station will be more accurate. At the same time, the ARIMA mode we utilized is a method that has been widely used as a famous linear time series model, which usually has the advantage of high accuracy, efficiency, and adapting ability.The bike-sharing system used for the case study is located in Jersey City, which has only 51 stations as shown in Figure1. This bike-sharing system operates during the daytime but remains closed at night, which is quite common in small/medium cities. The redistribution of bicycles is carried out during the night by means of a single vehicle that visits each station exactly once. The twelve-month historical operation data is used for the case study, which records the details of the trips including the departure and arrival time of each trip.Figure 1
Net flow of typical stations.
(a)
Station 3186 (January 2, 2016) (b)
Station 3192 (January 2, 2016) (c)
Station 3185 (January 2, 2016)The station state is characterized by the net flow on one-minute basis, which is calculated as the difference between the inflow and the outflow. From the net flow curve of each station, we can determine the difference between the bicycles available at the beginning of the day and those remaining at the end of the day. For convenience of observation, we set the initial net flow to zero. Figure1 shows the typical net flow curves of three stations, where the x- and y-axes represent the times series over one day and the difference between arrivals at and departures from a station per minute, respectively.Based on the state curve, the stations can be divided into three categories, namely, demand-station, supply-station, and demand-supply-station. These different types of stations play different roles in the bike-sharing system. For the demand-station, the number of the outgoing bicycles always exceeds the number of incoming bicycles, such as station 3186 shown in Figure1(a). For the supply-station, the number of outgoing bicycles is smaller than the number of incoming bicycles. Station 3129 shown in Figure 1(b) is an example of the supply-station, whose net flow remains positive, which indicates that the bicycles of this station can be delivered to other stations. At the third type of stations, demand-supply-station, the net flow can be both positive and negative, such as station 3185 shown in Figure 1(c).


### 3.2. The Demand Equation

The bicycle demand is the number of bikes required in stations, which can be defined based on the station state curve. As shown in Figure2, the real maximum demand of a day is determined by the minimum net flow value denoted by Ht, which is illustrated by the vertical line in Figure 2. We use TDt to denote the total difference between the bicycle number at the end of a day and that at the beginning of a day, which reflects the imbalance of the bicycle arrival and departure. Let Nt be the number of bicycles at the end of the previous day. We can easily find that the bicycle demand of a station qt is equal to Nt+Ht.Figure 2
State curve of station 3186 from January 2 to January 4, 2016.Based on similar analysis, we can define the bicycle demand of three types of the stations as follows.(  1) Demand-Station and Demand-Supply-Station.For a demand-station and a demand-supply-station, Hi,t is negative. The bicycle demand is calculated by(1)qi,t=Ni,t+Hi,twhereqi,t is the demand of station i of the tth day; Hi,t is the minimum demand of station i on the tth day (to be predicted); and Ni,t is the number of bicycles at station i at the end of the (t-1)th day.(  2) Supply-Station.For supply-station, Hj,t is zero and TDj,t is positive. The bicycle demand of the supply-station is defined as(2)qj,t=Nj,t+TDj,twhereqj,t is the supply of station j on the tth day; TDj,tis the total difference between the bicycle number at the end of the tth day and that at the beginning of the that day of station j; and Nj,t is the number of bicycles at station j at the end of the (t-1)th day.(  3) Shortage of the Bike-Sharing System.For a bike-sharing system, the total bicycle demand is equal to ∑iqi,t, and total supply is ∑jqj,t. If the total demand is greater than the total supply, then the bike-sharing system needs more bicycles. The number of bicycles that need to be supplied is equal to ∑iqi,t-∑jqj,t.If the total demand is smaller than the total supply, the bike-sharing system does not need new supplies, and the redundant bicycles do not need to be sent to the depot.


### 3.3. Maximum Daily Demand Prediction Model

In (1) and (2), Ni,t and Nj,t are given and Hi,t and TDj,t need to be predicted. There are many methods that can be used to handle the time series, such as multiple linear regression [23], exponential smoothing [24], artificial neural networks (ANN) [25], and methods based on a multimodel fusion algorithm [26]. Among the abovementioned methods, autoregressive integrated moving average (ARIMA) model is the most straightforward method to predict Hi,t and TDj,t value. ARIMA method has been widely used as a famous liner time series model, which usually has the advantage of high accuracy, efficiency, and adapting ability[27, 28].The following shows the prediction process ofHi,t. One can predict TDj,t by using the same method, so these will not be reviewed here.Figure3 shows the maximum demand (Hi,t) curve of station 3186 from Sep. 21, 2015, to Jan. 31, 2017.Figure 3
Maximum demand (Hi,t) curve of station 3186 from Sep. 21, 2015 to Jan. 31, 2017.Because the demand pattern is related to the day of a week, we divided them into categories based on the day of the week. Figure4 shows the prediction results of the Hi,t value of the Monday series. After calibration, we find that the ARIMA(3,1,3) model is the best for Monday series. The forecasting result is shown in Table 2. Due to the fact that the real demand value in station 3186 is small, the absolute percent errors actually indicate a small absolute difference, so this result is acceptable for a bike-sharing system.Table 2
Prediction error analysis of theHi,t value.
DateWeekday Feb. 1, 2017Sunday Feb. 2, 2017Monday Feb. 3, 2017Tuesday Feb. 4, 2017Wednesday Feb. 5, 2017Thursday Feb. 6, 2017Friday Feb. 7, 2017Saturday Predicted valuey -1 -29 -20 -27 -28 -42 -16  Measured valuex -2 -24 -19 -35 -33 -33 -18  Absolute Differencexi-yi 1 5 1 8 5 9 2  Absolute Percentage Errors 50% 20% 5% 22% 15% 27% 11%Figure 4
Prediction results for theHi,t value of Mondays by using ARIMA(3,1,3).With the proposed method, we have already the demand of each station. We can obtain the number of bicycles in need of repair at each station by checking and user report, which are used to guide the rebalancing work.


## 3.1. Pattern Analysis of the State of a Station

Most existing methods to predict the short-time demand of the stations are based on the historical pattern of the bicycle number change at each station. Andreas et al. [22] detected the temporal and geographic mobility patterns based on the data from the operator’s website. These patterns were applied to the prediction of the available bicycle number at any station.After a fairly long period since a bicycle sharing system operated at an area, the resident will develop a stationary trip habit gradually. And the trip habits will tend to be stable.There will be a state periodic similarity at each bike-sharing station. That is to say, there is a fixed input/output pattern in bicycle station. We utilize operational data to analyze and find this input/output pattern. On the basis of the detailed data analysis we find that there are mainly two types of bicycle station: demand and supply-station, and the bicycle demand is determined by the minimum net flow value. For a demand-station and a demand-supply-station, the bicycle demand is equal to the current number of bicycles plus the minimum net flow value; for supply-station, the supply of station is equal to the current number of bicycles plus the total difference (net flow value at the end of the day).One of the advantages is that forecasting the bicycle demand through analyzing the pattern of each station is more precise than other methods, such as traffic mode split forecast. Currently, the operational data is adequate. With the accumulation of operational data, the demand pattern recognition of each bicycle station will be more accurate. At the same time, the ARIMA mode we utilized is a method that has been widely used as a famous linear time series model, which usually has the advantage of high accuracy, efficiency, and adapting ability.The bike-sharing system used for the case study is located in Jersey City, which has only 51 stations as shown in Figure1. This bike-sharing system operates during the daytime but remains closed at night, which is quite common in small/medium cities. The redistribution of bicycles is carried out during the night by means of a single vehicle that visits each station exactly once. The twelve-month historical operation data is used for the case study, which records the details of the trips including the departure and arrival time of each trip.Figure 1
Net flow of typical stations.
(a)
Station 3186 (January 2, 2016) (b)
Station 3192 (January 2, 2016) (c)
Station 3185 (January 2, 2016)The station state is characterized by the net flow on one-minute basis, which is calculated as the difference between the inflow and the outflow. From the net flow curve of each station, we can determine the difference between the bicycles available at the beginning of the day and those remaining at the end of the day. For convenience of observation, we set the initial net flow to zero. Figure1 shows the typical net flow curves of three stations, where the x- and y-axes represent the times series over one day and the difference between arrivals at and departures from a station per minute, respectively.Based on the state curve, the stations can be divided into three categories, namely, demand-station, supply-station, and demand-supply-station. These different types of stations play different roles in the bike-sharing system. For the demand-station, the number of the outgoing bicycles always exceeds the number of incoming bicycles, such as station 3186 shown in Figure1(a). For the supply-station, the number of outgoing bicycles is smaller than the number of incoming bicycles. Station 3129 shown in Figure 1(b) is an example of the supply-station, whose net flow remains positive, which indicates that the bicycles of this station can be delivered to other stations. At the third type of stations, demand-supply-station, the net flow can be both positive and negative, such as station 3185 shown in Figure 1(c).


## 3.2. The Demand Equation

The bicycle demand is the number of bikes required in stations, which can be defined based on the station state curve. As shown in Figure2, the real maximum demand of a day is determined by the minimum net flow value denoted by Ht, which is illustrated by the vertical line in Figure 2. We use TDt to denote the total difference between the bicycle number at the end of a day and that at the beginning of a day, which reflects the imbalance of the bicycle arrival and departure. Let Nt be the number of bicycles at the end of the previous day. We can easily find that the bicycle demand of a station qt is equal to Nt+Ht.Figure 2
State curve of station 3186 from January 2 to January 4, 2016.Based on similar analysis, we can define the bicycle demand of three types of the stations as follows.(  1) Demand-Station and Demand-Supply-Station.For a demand-station and a demand-supply-station, Hi,t is negative. The bicycle demand is calculated by(1)qi,t=Ni,t+Hi,twhereqi,t is the demand of station i of the tth day; Hi,t is the minimum demand of station i on the tth day (to be predicted); and Ni,t is the number of bicycles at station i at the end of the (t-1)th day.(  2) Supply-Station.For supply-station, Hj,t is zero and TDj,t is positive. The bicycle demand of the supply-station is defined as(2)qj,t=Nj,t+TDj,twhereqj,t is the supply of station j on the tth day; TDj,tis the total difference between the bicycle number at the end of the tth day and that at the beginning of the that day of station j; and Nj,t is the number of bicycles at station j at the end of the (t-1)th day.(  3) Shortage of the Bike-Sharing System.For a bike-sharing system, the total bicycle demand is equal to ∑iqi,t, and total supply is ∑jqj,t. If the total demand is greater than the total supply, then the bike-sharing system needs more bicycles. The number of bicycles that need to be supplied is equal to ∑iqi,t-∑jqj,t.If the total demand is smaller than the total supply, the bike-sharing system does not need new supplies, and the redundant bicycles do not need to be sent to the depot.


## 3.3. Maximum Daily Demand Prediction Model

In (1) and (2), Ni,t and Nj,t are given and Hi,t and TDj,t need to be predicted. There are many methods that can be used to handle the time series, such as multiple linear regression [23], exponential smoothing [24], artificial neural networks (ANN) [25], and methods based on a multimodel fusion algorithm [26]. Among the abovementioned methods, autoregressive integrated moving average (ARIMA) model is the most straightforward method to predict Hi,t and TDj,t value. ARIMA method has been widely used as a famous liner time series model, which usually has the advantage of high accuracy, efficiency, and adapting ability[27, 28].The following shows the prediction process ofHi,t. One can predict TDj,t by using the same method, so these will not be reviewed here.Figure3 shows the maximum demand (Hi,t) curve of station 3186 from Sep. 21, 2015, to Jan. 31, 2017.Figure 3
Maximum demand (Hi,t) curve of station 3186 from Sep. 21, 2015 to Jan. 31, 2017.Because the demand pattern is related to the day of a week, we divided them into categories based on the day of the week. Figure4 shows the prediction results of the Hi,t value of the Monday series. After calibration, we find that the ARIMA(3,1,3) model is the best for Monday series. The forecasting result is shown in Table 2. Due to the fact that the real demand value in station 3186 is small, the absolute percent errors actually indicate a small absolute difference, so this result is acceptable for a bike-sharing system.Table 2
Prediction error analysis of theHi,t value.
DateWeekday Feb. 1, 2017Sunday Feb. 2, 2017Monday Feb. 3, 2017Tuesday Feb. 4, 2017Wednesday Feb. 5, 2017Thursday Feb. 6, 2017Friday Feb. 7, 2017Saturday Predicted valuey -1 -29 -20 -27 -28 -42 -16  Measured valuex -2 -24 -19 -35 -33 -33 -18  Absolute Differencexi-yi 1 5 1 8 5 9 2  Absolute Percentage Errors 50% 20% 5% 22% 15% 27% 11%Figure 4
Prediction results for theHi,t value of Mondays by using ARIMA(3,1,3).With the proposed method, we have already the demand of each station. We can obtain the number of bicycles in need of repair at each station by checking and user report, which are used to guide the rebalancing work.


## 4. Formulation of the VRPSPD Problem

To solve a problem like Vehicle Routing Problem with Simultaneous Picked and Delivery (VRPSPD), the common method is building a mixed integer linear programming (MILP) formulations with different objections or constrains. An integer linear programming formulation is proposed in this paper to model bike-sharing static rebalancing by considering the collection of bicycles in need of repair. The formulation is derived from Min (1989) [20] and Dethloff (2001) [19]. We enhanced the existing SBRP methodology in the following three ways: (a) the number of vehicles that perform the task may be greater than one; (b) the number of items picked up or delivered to each customer is predicted by a dedicated method derived from the net flow pattern analysis; and (c) there are two types of commodity in the model, ready bicycles, and damaged bicycles. The impact of damaged bicycles on the vehicle’s capacity constraint is considered in the proposed model. In the proposed model, the vehicles are required not only to deliver goods to customer, but also to pick up goods at customer location. This is imposed by constraints (11)-(16) in the paper.Specifically, the static rebalancing problem should involve determining how to drive at mostm vehicles through the graph, with the aim of minimizing the total cost and ensuring that the following constraints are not violated: (i) Each vehicle performs a route that starts and ends at the depot. (ii) Each vehicle starts from the depot empty or with some initial load (i.e., with a number of bicycles that vary from 0 to Q). (iii) Each station is visited exactly once and its request is completely fulfilled as a result of this visit. (iv) The number of bicycles in need of repair at each stationpi is given, and the quantities collected at pickup stations can be used to respond to the requests of delivery stations or can be returned directly to the depot; however, bicycles in disrepair cannot be used to meet the requests of delivery stations and need to be returned to the depot. (v) Finally, the sum of requests of the visited stations and the bicycles in disrepair plus the initial load is never negative or greater than Q along the route followed by a vehicle.Notations V: set of vertices, where V0⊆V\0,S≠∅ A: set of arcs n: number of stations p i: bicycles in disrepair at vertex i, pi≥0 K: set of vehicles k: the kth vehicle, where k=1,2,⋯,m m: number of vehicles Q: vehicle capacity q i: demand at vertex i c i j: cost of the arc i,j L 0 k: initial load of the kth vehicle L k0: final load of the kth vehicle when it returns to the depot L kj: load of the kth vehicle after vertex j, j∈V0 f i j: flow of arc i,j, that is the load of the vehicle passing arc i,j, if any, for (i,j)∈A Q t o t: the total demand g i j: the space remaining in the vehicle after the arc i,j, if any, for (i,j)∈A~ x i j k: a binary variable.The formulation for bike-sharing static rebalancing by considering the collection of bicycles in need of repair is given as follows:(3)min∑i∈V∑j∈V∑k∈Kcijxijksubject to(4)∑i∈V∑k∈Kxijk=1∀j∈V0(5)∑i∈Vxihk=∑j∈Vxhjk∀h∈V0;∀k∈K(6)∑j∈V0∑k∈Kx0jk≤m(7)∑j∈V0x0jk=∑i∈V0xi0k∀k∈K(8)Qtot=∑j∈V0qj(9)L0k=max⁡0,-∑i∈V∑j∈V0qjxijk∀k∈K(10)∑k∈KL0k≥max⁡0,-Qtot(11)Lk0=∑i∈V∑j∈V0pjxijk∀k∈K(12)∑k∈KLk0=∑j∈V0pj(13)Ljk≥L0k+qj+pj-M1-x0jk∀k∈K;∀j∈V0(14)Lkj≥θki+qj+pj-M1-∑k∈Kxijk∀i∈V0;∀j∈V0;i≠j(15)L0k≤Q,∀k∈K(16)Lkj≤Q,j∈V0,k∈K(17)∑i∈V0∑j∈V0∑k∈Kxijk≤S-max⁡0,-Qtot+∑i∈V0piQ,(18)xijk=1ifthearci,jwasusedbyvehiclek0otherwise.Objective function (3) minimizes the traveling cost. Constraints (4) require that every node except the depot is visited exactly once. Constraints (6) and (7) ensure, respectively, that at mostm vehicles leave the depot, and that all vehicles that are used return to the depot after completing their route. Constraints (9) impose the initial load of thekth vehicle. The total load leaving the initial depot should be nonnegative in any case; moreover, in caseQtot takes a negative value, it should not be lower than this value. This fact is imposed by constraints (10).Constraints (11) impose the final load of thekth vehicle. Constraints (12) state that the total load entering the depot is the number of all collected bicycles in need of repair. Constraints (13) impose the vehicle load after the first station. Constraints (14) impose the vehicle load after stationj. Constraints (15) and (16) specify the upper bounds on the loads.Constraints (17) state that, for each subsetS of vertices, the number of arcs with both their tail and head inS should not exceed the cardinality ofS minus the minimum number of vehicles required to serveS. An estimation of the minimum number of vehicles is simply obtained by computing the value of the sum of the demands and the sum of the bicycles in need of repair, dividing it by the vehicle capacity and then rounding up the result.


## 5. A Hybrid Algorithm Based on DPSO and VNS

It is intractable to use exact methods such as those shown in Table1 to solve large, realistic rebalancing problems, because the proposed problem is an NP-hard problem[4]. Previous studies [8, 13] have also illustrated this point by conducting numerical experiments. Hence, most of the existing studies focused on heuristic and metaheuristic approaches that can produce high-quality solutions within limited computational times [7, 8, 29, 30]. Ai and Kachitvichyanukul proposed Particle Swarm Optimization (PSO) algorithms for the VRPSPD[31]. Gajpal and Abad proposed an Ant Colony System (ACS) algorithm to solve the VRPSPD. However, the ACS is more complicated than other intelligence algorithms, such as PSO algorithms[32]. Berbeglia et al. suggested that the problem of VRPSPD falls under the category of a one-to-one pickup and delivery problem [29]. Reimann et al. proposed a decomposition approach (D-Ants) for the VRP model [33]. Li et al. discovered an Iterated Local Search algorithm, which is based on a metaheuristic. Their results showed that this algorithm is more accurate than others. However, the algorithm is computationally less efficient [34]. Liu et al. applied an algorithm for solving the Dynamic VRPSPD model based on POS of a Hopfield Neural Network. This algorithm combined the ability of the Hopfield Neural Network to reach the nearest point from the initial minimum point towards the descendent direction of gratitude and the advantage of the global search capability of the POS algorithm [35].Obviously, the static rebalancing problem is a combinational optimization problem, and the solution space is a discrete space made up of a limited number of integral points. Hence, Discrete Particle Swarm Optimization (DPSO) is applicable to solve this kind of problem. Compared with a Genetic algorithm, there is no mutation rate and cross rate parameters in DPSO; compared with the Simulated Annealing Method, there is no cumbersome process involving temperature drop; compared with the Ant Colony algorithm, there is no variety of complex judgment steps and parameter setting. In addition, the formula used to update the speed and position is used to ensure that all the particles converge to the global optimal value. Thus, the DPSO has a strong global searching capability.However, there are some drawbacks in the DPSO algorithm. First, setting the different parameters will have an obvious effect on the DPSO performance. Second, the DPSO algorithm may easily converge to a local optimum. Hence, to overcome these shortcomings while retaining the advantages, a new hybridization algorithm combining DPSO and Variable Neighborhood Search (VNS) was proposed by Goksal et al. [21]. This algorithm can enhance the global search ability of DPSO, overcome the possibility of entrapment in local optimum, and make full use of the local search ability of variable neighborhood search.The hybrid algorithm combined DPSO and VNS employed here includes four steps: representation of the solution, initializing the discrete particles, updating the positions of particles, and the strengthened local search based on the VNS algorithm.


### 5.1. Representation of the Solution

Finding a suitable representation of the solution is important. A good representation of the solution can enhance the performance of the DPSO algorithm and decrease the complexity of the algorithm. In this study we used the following procedures.There are N bicycle stations in the bicycle system, and the request is that the maximum load per route does not exceed the vehicle capacity (Q) and the remaining vehicle volume is sufficiently large to accommodate the bicycles from the next station after collection from or delivery to the last station. Each particle is a solution which contains several routes. The total number of particles is initialized withN={1,2,3,⋯,N0} positive integer numbers. Then each particle inserts zero (standing for the depot) under the constraints of the capacity of the vehicle and the space remaining in the vehicle. For instance, X=[1,3,6,5,0,8,11,2,4,0,7,10,9] is a solution, for which the particle has three routes. Route 1 is 0-1-3-6-5-0, route 2 is 0-8-11-2-4-0, and route 3 is 0-7-10-9-0. The cost of each particle is calculated by (19)cost=C0,Xi+∑k=i+1j-1CXk,Xk+1+CXj,0where Ci,j is the cost in arc (i,j), X(k) is the number of X particles representing the bicycle stations, and zero is the depot.


### 5.2. Initialization of the Discrete Particle Swarm

Commonly, the initial particle swarm is often a random product and we also adapt this format. Particles have different dimensions for different problems. Thus, the revised algorithm for gaining an initial particle swarm is provided in Algorithm1.Algorithm 1:The initial particle swarm algorithm.
Input:The number of bike stations, the cost matrix and the supply and demand per station Output:The initial solution (also known as the initial particle) Begin Initialize the discrete particle swarm based on the population parameter and the dimensional of particles, known asXt Loading the base data i.e., the cost matrix and supply and demand per station Whilen<Ndo Begin Load ←0 Route ←null r  ← 0 for  i=1,2,,…,N  do repeat Load = Load + Pi+Qi (Pi is the supply account, Qi is the demand account) Free Space = Capacity - Load (Capacity is the capacity of the vehicle) ifLoad > Capacity orFree Space < Pi O j =  i - 1 end for i = 1,2 , … , N  do repeat r←r+1; k←Oj for  m=k+1  to  j  do R o u t e(r) = Xi⋃Route(r) end end end n ← n + 1 end Return Route end


### 5.3. Position and Velocity Updating Rule for Particles

Generally, the velocity and position of the DPSO entail four basic operations, which are briefly mentioned here. The velocity and position updating rule of a particle is based on the following two equations:( 1 ) The velocity of a particle at iteration t is calculated by using four parts, the position minus, the position multiply, the velocity multiply, and the velocity and position, as shown below. (20)Vijt+1=wVijt⊕c1r1Pijt-Xijt⊕c2r2Gijt-Xijt( 2 ) The position of a particle at iteration t is updated by considering the operation of position plus velocity. The position updating rule is given by (21)Xijt+1=Xijt⊕Vijtwhere ⊕ is a newly defined operation.


### 5.4. Local Search Based on VNS Algorithm

As has been proposed in the literature review, there are many neighborhood structures for us to use. In our implementation, we utilize three neighborhood structures, i.e.,insert(i,j), swap(i,j), and2-Opt(i,j). These neighborhood structures are briefly introduced as follows:I n s e r t  ( i , j ) : This structure enables us to randomly produce two positive integer numbers (iandj) limited to the range of the particle, which indicate, respectively, the number of bicycles collected by customers and the position at which the particle is inserted and that both of them can belong to a different route; then the randomly produced number is collected and inserted into a randomly produced position.S w a p  ( i , j ) : This structure indicates that a customer i is removed from route 1 and a customer jis removed from route 2. Then these chosen customers cs1 and cs2 are reinserted into different routes, route 2 and route 1, respectively. Moreover, the route and the customer number cannot be adjacent but can be noncontiguous.2 - O p t  ( i , j ) : This structure reverses the direction of a path lying between the customers i-1 and j by replacing nonadjacent arcs(i-1i) and (j-1,j) belonging to the same route with (i-1,j-1) and (i,j), respectively. Equation (22) is used to update the personal best.(22)Pit=Pit,EiffPit,E<fPit,oPit,oiffxPit,E≥fPit,owhere Pit is the final personal best particle, Pit,E is the updated personal best particle operated by the three neighborhood structures, and Pit,o is the original personal best particle without any neighborhood structure operating. The cost function is expressed by (19).Figure5 illustrates the application of these neighborhood structures, and the VNS algorithm is shown in Algorithm 2. The hybridization of DPSO with VNS is shown in Algorithm 3.Algorithm 2:Local search algorithm (VNS).
Input:Initial vehicle routes Output:Improved vehicle routes by VNS algorithm Begin Initialize sequence of three neighborhood structures (1-insert, 2-swap, 3-2-opt) and randomly produce a positive integer of no more than 3 Choose thekth neighbor (S′) of S′∈N(S) from these three neighborhood structures and apply kthneighborhood to P and obtain(P′) if  FP′<F(P)  thenP←P′,F(P)←FP′ and apply the remainder of the neighborhood structures in sequence one by one to(P) and obtain(P′) if  FP′<F(P)  then  P←P′,F(P)←FP′ else delete the kth neighborhood structure and apply the remainder of the neighborhood structures one by one to (P) and obtain(P′) return  P and F(P) endAlgorithm 3:Hybridization of DPSO with VNS.
Input: Based data of bike system, DPSO parameters Output:The best optimization solution (routes of the bike system) Begin t←1 Initialize the population of the swarm, known as Pop and the relative parameters i.e., max iteration Apply the initial particle swarm algorithm to represent each particle in Xt Calculate the personal best particle inXt and obtain the personal best of each particle, stored in Pt Calculate the global best particle inXtand obtain the global best of each particle, stored in Gt Whileiisnotmorethanthemaxiterationorstoppingcriterionisnotsatisfied  do Begin Update particles using Eq. (21) and Eq. (22) Apply the initial particle swarm algorithmto evaluate the updating particles(Xt) Apply the Local search algorithm(VNS)to optimize the personal best particles(Pit) Update the Global best (Gt) endFigure 5
Three neighborhood structures.


## 5.1. Representation of the Solution

Finding a suitable representation of the solution is important. A good representation of the solution can enhance the performance of the DPSO algorithm and decrease the complexity of the algorithm. In this study we used the following procedures.There are N bicycle stations in the bicycle system, and the request is that the maximum load per route does not exceed the vehicle capacity (Q) and the remaining vehicle volume is sufficiently large to accommodate the bicycles from the next station after collection from or delivery to the last station. Each particle is a solution which contains several routes. The total number of particles is initialized withN={1,2,3,⋯,N0} positive integer numbers. Then each particle inserts zero (standing for the depot) under the constraints of the capacity of the vehicle and the space remaining in the vehicle. For instance, X=[1,3,6,5,0,8,11,2,4,0,7,10,9] is a solution, for which the particle has three routes. Route 1 is 0-1-3-6-5-0, route 2 is 0-8-11-2-4-0, and route 3 is 0-7-10-9-0. The cost of each particle is calculated by (19)cost=C0,Xi+∑k=i+1j-1CXk,Xk+1+CXj,0where Ci,j is the cost in arc (i,j), X(k) is the number of X particles representing the bicycle stations, and zero is the depot.


## 5.2. Initialization of the Discrete Particle Swarm

Commonly, the initial particle swarm is often a random product and we also adapt this format. Particles have different dimensions for different problems. Thus, the revised algorithm for gaining an initial particle swarm is provided in Algorithm1.Algorithm 1:The initial particle swarm algorithm.
Input:The number of bike stations, the cost matrix and the supply and demand per station Output:The initial solution (also known as the initial particle) Begin Initialize the discrete particle swarm based on the population parameter and the dimensional of particles, known asXt Loading the base data i.e., the cost matrix and supply and demand per station Whilen<Ndo Begin Load ←0 Route ←null r  ← 0 for  i=1,2,,…,N  do repeat Load = Load + Pi+Qi (Pi is the supply account, Qi is the demand account) Free Space = Capacity - Load (Capacity is the capacity of the vehicle) ifLoad > Capacity orFree Space < Pi O j =  i - 1 end for i = 1,2 , … , N  do repeat r←r+1; k←Oj for  m=k+1  to  j  do R o u t e(r) = Xi⋃Route(r) end end end n ← n + 1 end Return Route end


## 5.3. Position and Velocity Updating Rule for Particles

Generally, the velocity and position of the DPSO entail four basic operations, which are briefly mentioned here. The velocity and position updating rule of a particle is based on the following two equations:( 1 ) The velocity of a particle at iteration t is calculated by using four parts, the position minus, the position multiply, the velocity multiply, and the velocity and position, as shown below. (20)Vijt+1=wVijt⊕c1r1Pijt-Xijt⊕c2r2Gijt-Xijt( 2 ) The position of a particle at iteration t is updated by considering the operation of position plus velocity. The position updating rule is given by (21)Xijt+1=Xijt⊕Vijtwhere ⊕ is a newly defined operation.


## 5.4. Local Search Based on VNS Algorithm

As has been proposed in the literature review, there are many neighborhood structures for us to use. In our implementation, we utilize three neighborhood structures, i.e.,insert(i,j), swap(i,j), and2-Opt(i,j). These neighborhood structures are briefly introduced as follows:I n s e r t  ( i , j ) : This structure enables us to randomly produce two positive integer numbers (iandj) limited to the range of the particle, which indicate, respectively, the number of bicycles collected by customers and the position at which the particle is inserted and that both of them can belong to a different route; then the randomly produced number is collected and inserted into a randomly produced position.S w a p  ( i , j ) : This structure indicates that a customer i is removed from route 1 and a customer jis removed from route 2. Then these chosen customers cs1 and cs2 are reinserted into different routes, route 2 and route 1, respectively. Moreover, the route and the customer number cannot be adjacent but can be noncontiguous.2 - O p t  ( i , j ) : This structure reverses the direction of a path lying between the customers i-1 and j by replacing nonadjacent arcs(i-1i) and (j-1,j) belonging to the same route with (i-1,j-1) and (i,j), respectively. Equation (22) is used to update the personal best.(22)Pit=Pit,EiffPit,E<fPit,oPit,oiffxPit,E≥fPit,owhere Pit is the final personal best particle, Pit,E is the updated personal best particle operated by the three neighborhood structures, and Pit,o is the original personal best particle without any neighborhood structure operating. The cost function is expressed by (19).Figure5 illustrates the application of these neighborhood structures, and the VNS algorithm is shown in Algorithm 2. The hybridization of DPSO with VNS is shown in Algorithm 3.Algorithm 2:Local search algorithm (VNS).
Input:Initial vehicle routes Output:Improved vehicle routes by VNS algorithm Begin Initialize sequence of three neighborhood structures (1-insert, 2-swap, 3-2-opt) and randomly produce a positive integer of no more than 3 Choose thekth neighbor (S′) of S′∈N(S) from these three neighborhood structures and apply kthneighborhood to P and obtain(P′) if  FP′<F(P)  thenP←P′,F(P)←FP′ and apply the remainder of the neighborhood structures in sequence one by one to(P) and obtain(P′) if  FP′<F(P)  then  P←P′,F(P)←FP′ else delete the kth neighborhood structure and apply the remainder of the neighborhood structures one by one to (P) and obtain(P′) return  P and F(P) endAlgorithm 3:Hybridization of DPSO with VNS.
Input: Based data of bike system, DPSO parameters Output:The best optimization solution (routes of the bike system) Begin t←1 Initialize the population of the swarm, known as Pop and the relative parameters i.e., max iteration Apply the initial particle swarm algorithm to represent each particle in Xt Calculate the personal best particle inXt and obtain the personal best of each particle, stored in Pt Calculate the global best particle inXtand obtain the global best of each particle, stored in Gt Whileiisnotmorethanthemaxiterationorstoppingcriterionisnotsatisfied  do Begin Update particles using Eq. (21) and Eq. (22) Apply the initial particle swarm algorithmto evaluate the updating particles(Xt) Apply the Local search algorithm(VNS)to optimize the personal best particles(Pit) Update the Global best (Gt) endFigure 5
Three neighborhood structures.


## 6. Case Study

We assume a fictitious small-scale bicycle system as an illustrative example to show the efficiency of the proposed model and solving method. The system contains 14 bike stations. The daily supply and demand, the number of bicycles in need of repair, and the adjacency matrix for distance serving as cost matrix are given, separately, in Tables3, 4, and 6.Table 3
Daily supply and demand.
Station-ID 1 2 3 4 5 6 7 8 9 10 11 12 13 14  Demand and Supply 12 -2 8 5 -7 4 -2 -8 2 1 1 -3 3 10Table 4
Number of bicycles in disrepair.
Station-ID 1 2 3 4 5 6 7 8 9 10 11 12 13 14  Bicycles in disrepair 1 5 4 2 2 4 3 2 1 0 3 2 4 7The hybrid algorithm based on DPSO and VNS for rebalancing bike sharing has the following process.


### 6.1. Generating the Initial Discrete Particle Swarm

Normally, the initial discrete particle swarm is produced randomly, but choosing an appropriate initial swarm can expedite the convergence rate of the algorithm. In order to gain the favorable population, we initialize the discrete particle swarm randomly, belonging to the arrangementxij, xij∈1,20, i∈1,popsize, j∈1,length, and obtain the initial particle Xi=xij∣xij∈1,20,i∈1,popsize,j∈1,length.According to the size of this fictitious bike-sharing system, we set 20 as the solution scale. Divide the particles according to the constraint condition, and obtain the pattern of particles that was revealed in Table5. Numbers against asterisk mark denote the depot and the end of one route.Table 5
Initial Particles.
1st 6 3 11 7 0∗ 14 8 5 0∗ 1 2 4 13 0∗ 9 10 12  2nd 7 1 13 2 0∗ 14 6 10 0∗ 12 11 4 0∗ 8 3 9 5  3rd 2 10 4 0∗ 5 3 8 0∗ 12 11 7 1 0∗ 13 14 6 9  4th 10 9 0∗ 6 13 2 1 0∗ 11 14 3 8 7 0∗ 4 5 12  5th 8 9 x 0∗ 12 2 10 3 0∗ 7 4 1 13 11 0∗ 5 6  6th 2 3 0∗ 11 1 6 14 0∗ 5 9 8 0∗ 4 13 12 7 10  7th 9 2 1 0∗ 6 11 0∗ 13 10 8 14 3 0∗ 7 4 12 5  8th 14 7 10 8 0∗ 2 6 12 0∗ 1 3 4 5 0∗ 11 9 13  9th 3 7 1 0∗ 13 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2  10th 11 13 0∗ 1 5 14 12 0∗ 8 10 9 4 0∗ 3 7 6 2  11th 2 10 13 0∗ 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4  12th 9 10 0∗ 5 1 11 6 7 0∗ 4 14 12 0∗ 13 8 2 3  13th 10 6 4 0∗ 2 7 14 0∗ 13 1 8 5 0∗ 3 12 9 11  14th 1 10 0∗ 7 2 3 9 13 0∗ 14 6 0∗ 5 8 4 12 11  15th 8 14 6 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12  16th 3 11 10 0∗ 4 7 12 0∗ 13 5 14 2 0∗ 6 8 9 1  17th 4 1 14 9 0∗ 2 10 7 0∗ 13 12 8 6 0∗ 3 5 11  18th 6 10 14 0∗ 11 12 4 3 0∗ 1 2 13 0∗ 9 7 8 5  19th 11 1 0∗ 10 8 14 0∗ 3 13 2 9 0∗ 5 12 4 7 6  20th 9 11 6 0∗ 7 10 14 0∗ 5 1 2 0∗ 13 4 3 8 12Table 6
Cost Matrix.
Number 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 0.00 0.26 0.18 0.21 0.35 0.43 0.40 0.29 0.23 0.29 0.27 0.20 0.08 0.20 0.27  2 0.26 0.00 0.11 0.16 0.40 0.44 0.45 0.46 0.13 0.17 0.45 0.34 0.30 0.22 0.29  3 0.18 0.11 0.00 0.07 0.30 0.36 0.36 0.36 0.18 0.12 0.34 0.23 0.20 0.12 0.19  4 0.21 0.16 0.07 0.00 0.24 0.29 0.29 0.32 0.25 0.08 0.30 0.20 0.20 0.06 0.13  5 0.35 0.40 0.30 0.24 0.00 0.08 0.06 0.20 0.48 0.24 0.20 0.17 0.28 0.19 0.11  6 0.43 0.44 0.36 0.29 0.08 0.00 0.06 0.27 0.54 0.27 0.28 0.25 0.36 0.25 0.17  7 0.40 0.45 0.36 0.29 0.06 0.06 0.00 0.22 0.54 0.29 0.23 0.21 0.33 0.25 0.16  8 0.29 0.46 0.36 0.32 0.20 0.27 0.22 0.00 0.50 0.36 0.02 0.12 0.21 0.26 0.22  9 0.23 0.13 0.18 0.25 0.48 0.54 0.54 0.50 0.00 0.29 0.48 0.38 0.30 0.29 0.37  10 0.29 0.17 0.12 0.08 0.24 0.27 0.29 0.36 0.29 0.00 0.35 0.25 0.27 0.11 0.14  11 0.27 0.45 0.34 0.30 0.20 0.28 0.23 0.02 0.48 0.35 0.00 0.10 0.19 0.24 0.21  12 0.20 0.34 0.23 0.20 0.17 0.25 0.21 0.12 0.38 0.25 0.10 0.00 0.12 0.14 0.13  13 0.08 0.30 0.20 0.20 0.28 0.36 0.33 0.21 0.30 0.27 0.19 0.12 0.00 0.17 0.21  14 0.20 0.22 0.12 0.06 0.19 0.25 0.25 0.26 0.29 0.11 0.24 0.14 0.17 0.00 0.08  15 0.27 0.29 0.19 0.13 0.11 0.17 0.16 0.22 0.37 0.14 0.21 0.13 0.21 0.08 0.00 Note: Number 1 represents the depot.


### 6.2. Fitness Evaluation

The feasible solutions for this system are produced in Section6.1, so the calculation of these feasible solutions is the main topic of this section. We utilize the fitness value of each particle as an optimized target; the fitness function is a standard for screening particles; hence we choose the reciprocal of the objective function as the fitness function. The specific calculation formula is (23)Fitnesst=c0,Xti+∑k=i+1j-1cXtk,Xtk+1+cXtj,0whereXt(i) represents the value of the ith data of the tth particle,ci,j represents the vehicle cost (distance) from stationi toj,fitness is the fitness value of the particle.The following process of the first particle shows the steps to calculate the fitness value.1st631170∗14850∗124130∗91012Labeling method with list: Beginning from (0, 0) in the cost matrix, we obey these rules: Odd and even numbers of particles are indicated as rows and columns in the cost matrix, respectively. Further, when the adjacent number in the cost matrix is a cross, it is often the cost of the route. When all the numbers, exactly all the adjacent numbers, have been labeled, this process is complete. Last, the consequence is equal to the sum of all the labeled points.We present a representation of the cost matrix of the first particle in Table7.Table 7
Cost calculation progress of routes 1, 2, 3, and 4.
Number 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 0 0.00∗ 0.26③$ 0.18 0.21 0.35 0.43②& 0.40①# 0.29 0.23 0.29④@ 0.27 0.20 0.08④@ 0.20 0.27②&  1 0.26 0.00 0.11 0.16 0.40 0.44 0.45 0.46 0.13 0.17 0.45 0.34 0.30 0.22 0.29  2 0.18 0.11③$ 0.00 0.07 0.30③$ 0.36 0.36 0.36 0.18 0.12 0.34 0.23 0.20 0.12 0.19  3 0.21 0.16 0.07 0.00 0.24 0.29 0.29①# 0.32 0.25 0.08 0.30 0.20①# 0.20 0.06 0.13  4 0.35 0.40 0.30 0.24 0.00 0.08 0.06 0.20 0.48 0.24 0.20 0.17 0.28 0.19 0.11  5 0.43 0.44 0.36 0.29 0.08 0.00 0.06 0.27 0.54 0.27 0.28 0.25 0.36 0.25 0.17  6 0.40 0.45 0.36 0.29 0.06 0.06 0.00 0.22 0.54 0.29 0.23 0.21 0.33 0.25 0.16  7 0.29①# 0.46 0.36 0.32 0.20 0.27 0.22 0.00 0.50 0.36 0.02 0.12①# 0.21 0.26 0.22  8 0.23 0.13 0.18 0.25 0.48 0.54②& 0.54 0.50 0.00 0.29 0.48 0.38 0.30 0.29 0.37②&  9 0.29 0.17 0.12 0.08 0.24 0.27 0.29 0.36 0.29 0.00 0.35 0.25 0.27 0.11 0.14  10 0.27 0.45 0.34 0.30 0.20 0.28 0.23 0.02 0.48 0.35④@ 0.00 0.10 0.19④@ 0.24 0.21  11 0.20 0.34 0.23 0.20 0.17 0.25 0.21 0.12 0.38 0.25 0.10 0.00 0.12 0.14 0.13  12 0.08 0.30 0.20 0.20 0.28 0.36 0.33 0.21 0.30 0.27 0.19 0.12 0.00 0.17 0.21  13 0.20③$ 0.22 0.12 0.06 0.19③$ 0.25 0.25 0.26 0.29 0.11 0.24 0.14 0.17 0.00 0.08  14 0.27 0.29 0.19 0.13 0.11 0.17 0.16 0.22 0.37 0.14 0.21 0.13 0.21 0.08 0.00  ①# Route 1 ②& Route 2 ③$ Route 3 ④@ Route 4The fitness value is obtained by calculating the remaining particles in this way. The result is shown in Table8. Numbers against asterisk mark denote the depot and the end of one route.Table 8
Fitness value.
Initial Particles fitness 1st 6 3 11 7 0∗ 14 8 5 0∗ 1 2 4 13 0∗ 9 10 12 4.89  2nd 7 1 13 2 0∗ 14 6 10 0∗ 12 11 4 0∗ 8 3 9 5 4.70  3rd 2 10 4 0∗ 5 3 8 0∗ 12 11 7 1 0∗ 13 14 6 9 4.58  4th 10 9 0∗ 6 13 2 1 0∗ 11 14 3 8 7 0∗ 4 5 12 3.25  5th 8 9 x 0∗ 12 2 10 3 0∗ 7 4 1 13 11 0∗ 5 6 4.21  6th 2 3 0∗ 11 1 6 14 0∗ 5 9 8 0∗ 4 13 12 7 10 5.24  7th 9 2 1 0∗ 6 11 0∗ 13 10 8 14 3 0∗ 7 4 12 5 6.35  8th 14 7 10 8 0∗ 2 6 12 0∗ 1 3 4 5 0∗ 11 9 13 4.25  9th 3 7 1 0∗ 13 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2 4.68  10th 11 13 0∗ 1 5 14 12 0∗ 8 10 9 4 0∗ 3 7 6 2 4.78  11th 2 10 13 0∗ 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4 5.21  12th 9 10 0∗ 5 1 11 6 7 0∗ 4 14 12 0∗ 13 8 2 3 5.36  13th 10 6 4 0∗ 2 7 14 0∗ 13 1 8 5 0∗ 3 12 9 11 4.25  14th 1 10 0∗ 7 2 3 9 13 0∗ 14 6 0∗ 5 8 4 12 11 6.25  15th 8 14 6 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12 5.48  16th 3 11 10 0∗ 4 7 12 0∗ 13 5 14 2 0∗ 6 8 9 1 5.36  17th 4 1 14 9 0∗ 2 10 7 0∗ 13 12 8 6 0∗ 3 5 11 7.21  18th 6 10 14 0∗ 11 12 4 3 0∗ 1 2 13 0∗ 9 7 8 5 4.87  19th 11 1 0∗ 10 8 14 0∗ 3 13 2 9 0∗ 5 12 4 7 6 4.98  20th 9 11 6 0∗ 7 10 14 0∗ 5 1 2 0∗ 13 4 3 8 12 4.56


### 6.3. Calculation of the Historical and Global Optimal Values of the Particle

The historical and global optimal values are the particle’s cognitive and social messages, which have a significant impact on the probability of finding the preferable particle. Hence, the definition for the historical optimal particle is as follows: particles with the best historical fitness and global optimal particle; the best particles with the best historical fitness are mentioned. According to the definition, we obtain all the particles’ historical optimal values and the particle swarm’s global value at each iteration, as shown in Table9.Table 9
Historically optimal and globally optimal particles.
x 1,1 t x 1,2 t x 1,3 t ⋯ x 1 , l e n g t h - 2 t x 1 , l e n g t h - 1 t x 1 , l e n g t h t f 1 t 1st Particle’s Personal Best  x 2,1 k ′ x 2,2 k x 2,3 k ⋯ x 2 , l e n g t h - 2 k x 2 , l e n g t h - 1 k x 2 , l e n g t h k f 2 k 2nd Particle’s Personal Best  ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯  x 20,1 l ′ x 20,2 l x 20,3 l ⋯ x 20 , l e n g t h - 2 l x 20 , l e n g t h - 1 l x 20 , l e n g t h l f 20 l 20th Particle’s Personal Best  Global Best min ⁡ f i t ∣ i ∈ 1 , p o p s i z e , t ∈ 1 , M a x i t e r a t eIn the first iteration, we initialize the original fitness value of the particles as the personal best and initialize the minimum of the personal best as the global best. The consequence is represented in Table10.Table 10
Initial fitness value, personal best, and global best of the particles.
Initial Particles fitness Personal best 1st 6 3 11 7 0∗ 14 8 5 0∗ 1 2 4 13 0∗ 9 10 12 4.89 4.89  2nd 7 1 13 2 0∗ 14 6 10 0∗ 12 11 4 0∗ 8 3 9 5 4.70 4.70  3rd 2 10 4 0∗ 5 3 8 0∗ 12 11 7 1 0∗ 13 14 6 9 4.58 4.58  4th 10 9 0∗ 6 13 2 1 0∗ 11 14 3 8 7 0∗ 4 5 12 3.25 3.25  (Global best)  5th 8 9 x 0∗ 12 2 10 3 0∗ 7 4 1 13 11 0∗ 5 6 4.21 4.21  6th 2 3 0∗ 11 1 6 14 0∗ 5 9 8 0∗ 4 13 12 7 10 5.24 5.24  7th 9 2 1 0∗ 6 11 0∗ 13 10 8 14 3 0∗ 7 4 12 5 6.35 6.35  8th 14 7 10 8 0∗ 2 6 12 0∗ 1 3 4 5 0∗ 11 9 13 4.25 4.25  9th 3 7 1 0∗ 13 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2 4.68 4.68  10th 11 13 0∗ 1 5 14 12 0∗ 8 10 9 4 0∗ 3 7 6 2 4.78 4.78  11th 2 10 13 0∗ 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4 5.21 5.21  12th 9 10 0∗ 5 1 11 6 7 0∗ 4 14 12 0∗ 13 8 2 3 5.36 5.36  13th 10 6 4 0∗ 2 7 14 0∗ 13 1 8 5 0∗ 3 12 9 11 4.25 4.25  14th 1 10 0∗ 7 2 3 9 13 0∗ 14 6 0∗ 5 8 4 12 11 6.25 6.25  15th 8 14 6 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12 5.48 5.48  16th 3 11 10 0∗ 4 7 12 0∗ 13 5 14 2 0∗ 6 8 9 1 5.36 5.36  17th 4 1 14 9 0∗ 2 10 7 0∗ 13 12 8 6 0∗ 3 5 11 7.21 7.21  18th 6 10 14 0∗ 11 12 4 3 0∗ 1 2 13 0∗ 9 7 8 5 4.87 4.87  19th 11 1 0∗ 10 8 14 0∗ 3 13 2 9 0∗ 5 12 4 7 6 4.98 4.98  20th 9 11 6 0∗ 7 10 14 0∗ 5 1 2 0∗ 13 4 3 8 12 4.56 4.56 Note: Numbers against asterisk mark denote the depot and the end of one route.


### 6.4. Local Search through the Variable Neighborhood Operation

Our algorithm utilizes the variable neighborhood search algorithm given in Table11 to enhance the ability to choose the solution randomly, ensure diversity of the individual, and overcome the disadvantage that the solution is a local optimum. First, we updated the individual particles’ historical best and global best. Then, in our implementation, three variable neighborhood structures are applied to the particle’s personal best mentioned at Section 5.4.Table 11
Variable neighborhood operation.
x 1,1 t ′ x 1,2 t ′ x 1,3 t ′ ⋯ x 1 , l e n g t h - 2 t ′ x 1 , l e n g t h - 1 t ′ x 1 , l e n g t h t ′ f 1 t ′ 1st Particle’s Personal Best  x 2,1 k ′ x 2,2 k ′ x 2,3 k ′ ⋯ x 2 , l e n g t h - 2 k ′ x 2 , l e n g t h - 1 k ′ x 2 , l e n g t h k ′ f 2 k ′ 2nd Particle’s Personal Best  ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯  x 20,1 l ′ x 20,2 l ′ x 20,3 l ′ ⋯ x 20 , l e n g t h - 2 l ′ x 20 , l e n g t h - 1 l ′ x 20 , l e n g t h l ′ f 20 l ′ 20th Particle’s Personal Best  After variable neighborhood operation Global Best min ⁡ f i t ′ ∣ i ∈ 1 , p o p s i z e , t ′ ∈ 1 , M a x i t e r a t e Note:xi,jk′ is the result of the variable neighborhood operation.Based on the particular variable neighborhood structure, we used the personal best as an operation by taking the 2-opt operation structure as an example. The steps and consequence are the following.Procedure 1.
Produce a group of 2-opt structures randomly; i.e., in the first particle, there are three points, (3, 7), (14, 5), and (2, 10), respectively.Procedure 2.
Subject the particle to the 2-opt operation. An example using the1st particle is provided in Table 12.Table 12
Particles subjected to neighborhood structure operation (2-Opt).
1st 6 3 11 7 0∗ 14 8 5 0∗ 1 2 4 13 0∗ 9 10 12  After 2-opt-1st 6 3 7& 11& 0∗ 14 5& 8& 0∗ 1 2 10& 9& 0& 13& 4& 12  2nd 7 1 13 2 0∗ 14 6 10 0∗ 12 11 4 0∗ 8 3 9 5  3rd 2 10 4 0∗ 5 3 8 0∗ 12 11 7 1 0∗ 13 14 6 9  4th 10 9 0∗ 6 13 2 1 0∗ 11 14 3 8 7 0∗ 4 5 12  5th 8 9 14 0∗ 12 2 10 3 0∗ 7 4 1 13 11 0∗ 5 6  6th 2 3 0∗ 11 1 6 14 0∗ 5 9 8 0∗ 4 13 12 7 10  7th 9 2 1 0∗ 6 11 0∗ 13 10 8 14 3 0∗ 7 4 12 5  8th 14 7 10 8 0∗ 2 6 12 0∗ 1 3 4 5 0∗ 11 9 13  9th 3 7 1 0∗ 13 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2  10th 11 13 0∗ 1 5 14 12 0∗ 8 10 9 4 0∗ 3 7 6 2  11th 2 10 13 0∗ 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4  12th 9 10 0∗ 5 1 11 6 7 0∗ 4 14 12 0∗ 13 8 2 3  13th 10 6 4 0∗ 2 7 14 0∗ 13 1 8 5 0∗ 3 12 9 11  14th 1 10 0∗ 7 2 3 9 13 0∗ 14 6 0∗ 5 8 4 12 11  15th 8 14 6 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12  16th 3 11 10 0∗ 4 7 12 0∗ 13 5 14 2 0∗ 6 8 9 1  17th 4 1 14 9 0∗ 2 10 7 0∗ 13 12 8 6 0∗ 3 5 11  18th 6 10 14 0∗ 11 12 4 3 0∗ 1 2 13 0∗ 9 7 8 5  19th 11 1 0∗ 10 8 14 0∗ 3 13 2 9 0∗ 5 12 4 7 6  20th 9 11 6 0∗ 7 10 14 0∗ 5 1 2 0∗ 13 4 3 8 12 Note: Numbers against asterisk mark denote the depot and the end of one route. Numbers with & mark were subjected to operation by the 2-opt structure.The remaining particles were subjected to operation according to the form procedure, and the final result is presented in Table13.Table 13
Particles subjected to 2-Opt structure operation.
After 2-opt-1st 6 3 7& 11& 0∗ 14 5& 8& 0∗ 1 2 10& 9& 0& 13& 4& 12  After 2-opt-2nd 7 13& 1& 2 0∗ 14 4& 11& 12& 0& 10& 6& 0∗ 8 5& 9& 3&  After 2-opt-3rd 4& 10& 2& 0∗ 5 3 7& 11& 12& 0& 8& 1 0∗ 13 14 6 9  After 2-opt-4th 10 9 0∗ 6 1& 2& 13& 0∗ 11 14 3 8 5& 4& 0& 7& 12  After 2-opt-5th 8 14& 9& 0∗ 12 2 10 3 0∗ 7 4 11& 13& 1& 0∗ 5 6  After 2-opt-6th 2 3 0∗ 11 1 6 14 0∗ 5 9 4& 0& 8& 13 12 7 10  After 2-opt-7th 9 6& 0& 1& 2& 11 0∗ 13 10 8 14 4& 7& 0& 3& 12 5  After 2-opt-8th 14 7 10 2& 0& 8& 6 12 0∗ 5& 4& 3& 1& 0& 11 9 13  After 2-opt-9th 3 7 13& 0& 1& 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2  After 2-opt-10th 11 13 0∗ 1 14& 5& 12 0∗ 8 10 9 2& 6& 7& 3& 0& 4&  After 2-opt-11th 2 0& 13& 10& 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4  After 2-opt-12th 9 10 0∗ 1& 5& 11 6 4& 0& 7& 14 12 0∗ 13 8 2 3  After 2-opt-13th 10 2& 0& 4& 6& 7 14 0∗ 13 1 5& 8& 0∗ 3 12 9 11  After 2-opt-14th 1 10 7& 0& 2 3 9 13 14& 0& 6 0∗ 5 8 4 12 11  After 2-opt-15th 8 6& 14& 0∗ 1 7 9 4 0∗ 2 11 0& 13& 10 3 5 12  After 2-opt-16th 3 11 10 0∗ 4 7 12 0∗ 13 5 0& 2& 14& 6 8 9 1  After 2-opt-17th 4 9& 14& 1& 0∗ 2 10 7 0∗ 13 12 8 3& 0& 6& 5 11  After 2-opt-18th 6 10 11& 0& 14& 12 4 3 0∗ 1 2 13 7& 9& 0& 8 5  After 2-opt-19th 11 1 0∗ 10 8 9& 2& 13& 3& 0& 14& 0∗ 5 12 4 7 6  After 2-opt-20th 9 11 7& 0& 6& 10 14 0∗ 5 1 8& 3& 4& 13& 0& 2& 12 Note: Numbers against asterisk mark denote the depot and the end of one route. Numbers with & mark were subjected to operation by the 2-opt structure.Then, according to the constraint of the model, we check and recode particles. If the operated particles can satisfy the constraint, we accept these particles; otherwise, we recode the operated particles. The result of checking and recoding can be seen in Table14.Table 14
Recoding particles after the 2-Opt operation.
After 2-opt-1st 6 3 7 0∗ 11 14 5 8 0∗ 1 2 10 9 0∗ 13 4 12  After 2-opt-2nd 7 13 0∗ 1 2 14 4 0∗ 11 12 10 6 0∗ 8 5 9 3  After 2-opt-3rd 4 10 2 0∗ 5 3 7 11 12 0 8 1 0∗ 13 14 6 9  After 2-opt-4th 10 9 0∗ 6 1 2 13 0∗ 11 14 3 8 5 4 0∗ 7 12  After 2-opt-5th 8 14 9 0∗ 12 2 10 3 0∗ 7 4 11 13 1 0∗ 5 6  After 2-opt-6th 2 3 0∗ 11 1 6 14 0∗ 5 9 4 0∗ 8 13 12 7 10  After 2-opt-7th 9 6 1 2 0∗ 11 13 10 8 0∗ 14 4 7 0∗ 3 12 5  After 2-opt-8th 14 7 10 2 0∗ 8 6 12 0∗ 5 4 3 1 0∗ 11 9 13  After 2-opt-9th 3 7 13 0∗ 1 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2  After 2-opt-10th 11 13 0∗ 1 14 5 12 0∗ 8 10 9 2 6 7 3 0∗ 4  After 2-opt-11th 2 13 10 0∗ 7 9 12 14 0∗ 1 6 5 11 0∗ 8 3 4  After 2-opt-12th 9 10 0∗ 1 5 11 6 4 0∗ 7 14 12 0∗ 13 8 2 3  After 2-opt-13th 10 2 4 0∗ 6 7 14 0∗ 13 1 5 8 0∗ 3 12 9 11  After 2-opt-14th 1 10 7 0∗ 2 3 9 13 14 0∗ 6 5 8 0∗ 4 12 11  After 2-opt-15th 8 6 14 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12  After 2-opt-16th 3 11 0∗ 10 4 7 12 0∗ 13 5 2 0∗ 14 6 8 9 1  After 2-opt-17th 4 9 14 1 0 2 10 7 0∗ 13 12 8 3 0∗ 6 5 11  After 2-opt-18th 6 10 11 0∗ 14 12 4 3 0∗ 1 2 13 7 9 0∗ 8 5  After 2-opt-19th 11 1 0∗ 10 8 9 2 13 3 0∗ 14 5 12 0∗ 4 7 6  After 2-opt-20th 9 11 7 0∗ 6 10 14 0∗ 5 1 8 3 4 13 0∗ 2 12 Note: Numbers against asterisk mark denote the depot and the end of one route.


### 6.5. Particle Updating Rules and Algorithm Termination Condition

The velocity and position of a particle are updated according to (21) and (22), after which iterative operation occurs. The particles are updated according to the relative operation relating to the position and velocity. The results are provided in Table 15. Subsequently, the algorithm calculates the fitness value, personal best, and global best of the updated particles until fitnessn+1-fitness(n)/fitnessn+1≤0.05 or the iterative value exceeds the value of the original iterative step.Table 15
Updated particles.
1st 11 2 9 0∗ 3 5 8 12 13 0∗ 10 14 7 0∗ 1 6 4  2nd 7 8 0∗ 12 10 13 6 4 0∗ 14 3 11 9 0∗ 1 5 2  3rd 1 11 4 ∗ 7 9 10 6 0∗ 8 13 3 0∗ 5 14 2 12  4th 5 14 0∗ 10 9 13 4 6 0∗ 2 3 11 0∗ 12 1 7 8  5th 8 12 2 3 0∗ 4 5 7 6 0∗ 10 11 13 0∗ 9 14 1  6th 8 7 5 0∗ 11 10 12 6 13 4 0∗ 14 3 2 0∗ 9 1  7th 1 7 12 0∗ 8 10 5 6 14 0∗ 3 9 2 0∗ 4 11 13  8th 11 8 14 9 0∗ 3 13 10 4 0∗ 5 12 1 0∗ 6 7 2  9th 6 13 0∗ 10 1 14 12 0∗ 4 9 11 3 0∗ 8 7 2 5  10th 6 4 2 0∗ 9 12 8 3 10 0∗ 14 5 11 0∗ 1 7 13  11th 4 2 0∗ 9 8 13 1 0∗ 3 5 12 7 0∗ 6 10 14 11  12th 4 5 12 2 0∗ 1 8 14 7 0∗ 6 13 10 9 0∗ 3 11  13th 14 2 7 0∗ 8 4 1 10 9 0∗ 5 13 6 0∗ 3 12 11  14th 1 2 8 9 0∗ 6 11 14 10 12 0∗ 5 7 13 0∗ 3 4  15th 6 14 13 0∗ 2 7 12 8 0∗ 10 9 3 1 0∗ 11 5 4  16th 14 10 8 0∗ 2 5 12 11 3 0∗ 7 13 4 0∗ 6 1 9  17th 2 13 3 5 0∗ 1 8 11 6 0∗ 12 9 0∗ 10 4 14 7  18th 8 13 1 2 0∗ 7 14 6 0∗ 10 11 4 12 0∗ 9 5 3  19th 6 5 11 0∗ 3 13 4 14 0∗ 7 10 12 0∗ 8 2 9 1  20th 14 1 2 0∗ 3 12 4 0∗ 13 7 8 6 5 0∗ 11 10 9 Note: Numbers against asterisk mark denote the depot and the end of one route.When the conditions of convergence are met, the operation is complete. Otherwise, return to process in Section6.1 and implement the iterative operation again.


### 6.6. Computational Result and Analysis

The parameters of the hybridization algorithm of DPSO with VNS are included in Table16.Table 16
Parameters of the algorithm.
Parameters’ name of Hybridization algorithm of DPSO with VNS Value Population of particles (N) 20  The length of per particle (L) 14  Weight Value (C) 0.6  Iteration 50  Learning Value (CP) 1.5  Social Value (Cg) 2The dispatching plan is presented in Table17 and the dispatching routes are shown in Figure 6.Table 17
Dispatching plan.
Route-ID Routes Cost of per route R1 9-2-1-8 5.4 km  R2 3-13-14 4.7 km  R3 10-11-12 6.2 km  R4 4-6-5-7 9.33 km  Total Cost 25.63 kmFigure 6
Dispatching route of the hybridization algorithm of DPSO with VNS.Computation using the hybridization algorithm of DPSO with VNS produces four routes for dispatching, route 1=[9-2-1-8], route 2= [3-13-14], route 3= [10-11-12], and route 4 = [4-6-5-7]. Furthermore, the cost per route scaled by the distance corresponds to 5.4 km, 4.7 km, 6.2 km, and 9.33 km. The total cost is 25.63 km.In order to test the performance of the hybridization algorithm of DPSO with VNS, Goksal et al. made an overall comparison with Reactive Tabu Search (RTS), Large Neighborhood Search (LNS), Particle Swarm Optimization (PSO), Ant Colony System (ACS), Parallel Iterative Local Search (PILS), and Adaptive Memory Methodology (AMM). Since they are the most effective heuristics proposed in the literature. According to the indexes of average deviation to the best known solution, average computation time, scaled average computation time over 40 instances, and number of best known solutions found, the comprehensive performance of hybridization algorithm of DPSO with VNS is the best. For the computational results of the comparison, see Goksal et al. [21].In this study, we implement Discrete Particle Swarm Optimization (DPSO) and Genetic algorithm (GA) separately over the case study data set to get a more comprehensive view of the performance of the hybridization algorithm of DPSO with VNS. The results are provided in Figures7 and 8. DPSO is considered in the comparison with the purpose of presenting the performance of VNS. We can see from Figure 7 that the cost of the solution without VNS is 27.07 km. The result of the hybridization algorithm of DPSO with VNS (25.63 km) dropped by 5.3%. GA [15] is a population-based heuristic which can be used to solve the problem like SBRP. Y. Li et al. reported the performance in their study. Thus, in order to make a direct comparison, we run a GA program over the same data set of case study. From Figure 8, we can see that the total cost is 26.33 km, while the hybridization algorithm of DPSO with VNS brings a drop by 2.7%. This analysis indicates that the hybridization algorithm of DPSO with VNS outperforms the other existing heuristics in the SBRP literature.Figure 7
Result of DPSO algorithm.Figure 8
Result of GA.


## 6.1. Generating the Initial Discrete Particle Swarm

Normally, the initial discrete particle swarm is produced randomly, but choosing an appropriate initial swarm can expedite the convergence rate of the algorithm. In order to gain the favorable population, we initialize the discrete particle swarm randomly, belonging to the arrangementxij, xij∈1,20, i∈1,popsize, j∈1,length, and obtain the initial particle Xi=xij∣xij∈1,20,i∈1,popsize,j∈1,length.According to the size of this fictitious bike-sharing system, we set 20 as the solution scale. Divide the particles according to the constraint condition, and obtain the pattern of particles that was revealed in Table5. Numbers against asterisk mark denote the depot and the end of one route.Table 5
Initial Particles.
1st 6 3 11 7 0∗ 14 8 5 0∗ 1 2 4 13 0∗ 9 10 12  2nd 7 1 13 2 0∗ 14 6 10 0∗ 12 11 4 0∗ 8 3 9 5  3rd 2 10 4 0∗ 5 3 8 0∗ 12 11 7 1 0∗ 13 14 6 9  4th 10 9 0∗ 6 13 2 1 0∗ 11 14 3 8 7 0∗ 4 5 12  5th 8 9 x 0∗ 12 2 10 3 0∗ 7 4 1 13 11 0∗ 5 6  6th 2 3 0∗ 11 1 6 14 0∗ 5 9 8 0∗ 4 13 12 7 10  7th 9 2 1 0∗ 6 11 0∗ 13 10 8 14 3 0∗ 7 4 12 5  8th 14 7 10 8 0∗ 2 6 12 0∗ 1 3 4 5 0∗ 11 9 13  9th 3 7 1 0∗ 13 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2  10th 11 13 0∗ 1 5 14 12 0∗ 8 10 9 4 0∗ 3 7 6 2  11th 2 10 13 0∗ 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4  12th 9 10 0∗ 5 1 11 6 7 0∗ 4 14 12 0∗ 13 8 2 3  13th 10 6 4 0∗ 2 7 14 0∗ 13 1 8 5 0∗ 3 12 9 11  14th 1 10 0∗ 7 2 3 9 13 0∗ 14 6 0∗ 5 8 4 12 11  15th 8 14 6 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12  16th 3 11 10 0∗ 4 7 12 0∗ 13 5 14 2 0∗ 6 8 9 1  17th 4 1 14 9 0∗ 2 10 7 0∗ 13 12 8 6 0∗ 3 5 11  18th 6 10 14 0∗ 11 12 4 3 0∗ 1 2 13 0∗ 9 7 8 5  19th 11 1 0∗ 10 8 14 0∗ 3 13 2 9 0∗ 5 12 4 7 6  20th 9 11 6 0∗ 7 10 14 0∗ 5 1 2 0∗ 13 4 3 8 12Table 6
Cost Matrix.
Number 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 0.00 0.26 0.18 0.21 0.35 0.43 0.40 0.29 0.23 0.29 0.27 0.20 0.08 0.20 0.27  2 0.26 0.00 0.11 0.16 0.40 0.44 0.45 0.46 0.13 0.17 0.45 0.34 0.30 0.22 0.29  3 0.18 0.11 0.00 0.07 0.30 0.36 0.36 0.36 0.18 0.12 0.34 0.23 0.20 0.12 0.19  4 0.21 0.16 0.07 0.00 0.24 0.29 0.29 0.32 0.25 0.08 0.30 0.20 0.20 0.06 0.13  5 0.35 0.40 0.30 0.24 0.00 0.08 0.06 0.20 0.48 0.24 0.20 0.17 0.28 0.19 0.11  6 0.43 0.44 0.36 0.29 0.08 0.00 0.06 0.27 0.54 0.27 0.28 0.25 0.36 0.25 0.17  7 0.40 0.45 0.36 0.29 0.06 0.06 0.00 0.22 0.54 0.29 0.23 0.21 0.33 0.25 0.16  8 0.29 0.46 0.36 0.32 0.20 0.27 0.22 0.00 0.50 0.36 0.02 0.12 0.21 0.26 0.22  9 0.23 0.13 0.18 0.25 0.48 0.54 0.54 0.50 0.00 0.29 0.48 0.38 0.30 0.29 0.37  10 0.29 0.17 0.12 0.08 0.24 0.27 0.29 0.36 0.29 0.00 0.35 0.25 0.27 0.11 0.14  11 0.27 0.45 0.34 0.30 0.20 0.28 0.23 0.02 0.48 0.35 0.00 0.10 0.19 0.24 0.21  12 0.20 0.34 0.23 0.20 0.17 0.25 0.21 0.12 0.38 0.25 0.10 0.00 0.12 0.14 0.13  13 0.08 0.30 0.20 0.20 0.28 0.36 0.33 0.21 0.30 0.27 0.19 0.12 0.00 0.17 0.21  14 0.20 0.22 0.12 0.06 0.19 0.25 0.25 0.26 0.29 0.11 0.24 0.14 0.17 0.00 0.08  15 0.27 0.29 0.19 0.13 0.11 0.17 0.16 0.22 0.37 0.14 0.21 0.13 0.21 0.08 0.00 Note: Number 1 represents the depot.


## 6.2. Fitness Evaluation

The feasible solutions for this system are produced in Section6.1, so the calculation of these feasible solutions is the main topic of this section. We utilize the fitness value of each particle as an optimized target; the fitness function is a standard for screening particles; hence we choose the reciprocal of the objective function as the fitness function. The specific calculation formula is (23)Fitnesst=c0,Xti+∑k=i+1j-1cXtk,Xtk+1+cXtj,0whereXt(i) represents the value of the ith data of the tth particle,ci,j represents the vehicle cost (distance) from stationi toj,fitness is the fitness value of the particle.The following process of the first particle shows the steps to calculate the fitness value.1st631170∗14850∗124130∗91012Labeling method with list: Beginning from (0, 0) in the cost matrix, we obey these rules: Odd and even numbers of particles are indicated as rows and columns in the cost matrix, respectively. Further, when the adjacent number in the cost matrix is a cross, it is often the cost of the route. When all the numbers, exactly all the adjacent numbers, have been labeled, this process is complete. Last, the consequence is equal to the sum of all the labeled points.We present a representation of the cost matrix of the first particle in Table7.Table 7
Cost calculation progress of routes 1, 2, 3, and 4.
Number 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 0 0.00∗ 0.26③$ 0.18 0.21 0.35 0.43②& 0.40①# 0.29 0.23 0.29④@ 0.27 0.20 0.08④@ 0.20 0.27②&  1 0.26 0.00 0.11 0.16 0.40 0.44 0.45 0.46 0.13 0.17 0.45 0.34 0.30 0.22 0.29  2 0.18 0.11③$ 0.00 0.07 0.30③$ 0.36 0.36 0.36 0.18 0.12 0.34 0.23 0.20 0.12 0.19  3 0.21 0.16 0.07 0.00 0.24 0.29 0.29①# 0.32 0.25 0.08 0.30 0.20①# 0.20 0.06 0.13  4 0.35 0.40 0.30 0.24 0.00 0.08 0.06 0.20 0.48 0.24 0.20 0.17 0.28 0.19 0.11  5 0.43 0.44 0.36 0.29 0.08 0.00 0.06 0.27 0.54 0.27 0.28 0.25 0.36 0.25 0.17  6 0.40 0.45 0.36 0.29 0.06 0.06 0.00 0.22 0.54 0.29 0.23 0.21 0.33 0.25 0.16  7 0.29①# 0.46 0.36 0.32 0.20 0.27 0.22 0.00 0.50 0.36 0.02 0.12①# 0.21 0.26 0.22  8 0.23 0.13 0.18 0.25 0.48 0.54②& 0.54 0.50 0.00 0.29 0.48 0.38 0.30 0.29 0.37②&  9 0.29 0.17 0.12 0.08 0.24 0.27 0.29 0.36 0.29 0.00 0.35 0.25 0.27 0.11 0.14  10 0.27 0.45 0.34 0.30 0.20 0.28 0.23 0.02 0.48 0.35④@ 0.00 0.10 0.19④@ 0.24 0.21  11 0.20 0.34 0.23 0.20 0.17 0.25 0.21 0.12 0.38 0.25 0.10 0.00 0.12 0.14 0.13  12 0.08 0.30 0.20 0.20 0.28 0.36 0.33 0.21 0.30 0.27 0.19 0.12 0.00 0.17 0.21  13 0.20③$ 0.22 0.12 0.06 0.19③$ 0.25 0.25 0.26 0.29 0.11 0.24 0.14 0.17 0.00 0.08  14 0.27 0.29 0.19 0.13 0.11 0.17 0.16 0.22 0.37 0.14 0.21 0.13 0.21 0.08 0.00  ①# Route 1 ②& Route 2 ③$ Route 3 ④@ Route 4The fitness value is obtained by calculating the remaining particles in this way. The result is shown in Table8. Numbers against asterisk mark denote the depot and the end of one route.Table 8
Fitness value.
Initial Particles fitness 1st 6 3 11 7 0∗ 14 8 5 0∗ 1 2 4 13 0∗ 9 10 12 4.89  2nd 7 1 13 2 0∗ 14 6 10 0∗ 12 11 4 0∗ 8 3 9 5 4.70  3rd 2 10 4 0∗ 5 3 8 0∗ 12 11 7 1 0∗ 13 14 6 9 4.58  4th 10 9 0∗ 6 13 2 1 0∗ 11 14 3 8 7 0∗ 4 5 12 3.25  5th 8 9 x 0∗ 12 2 10 3 0∗ 7 4 1 13 11 0∗ 5 6 4.21  6th 2 3 0∗ 11 1 6 14 0∗ 5 9 8 0∗ 4 13 12 7 10 5.24  7th 9 2 1 0∗ 6 11 0∗ 13 10 8 14 3 0∗ 7 4 12 5 6.35  8th 14 7 10 8 0∗ 2 6 12 0∗ 1 3 4 5 0∗ 11 9 13 4.25  9th 3 7 1 0∗ 13 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2 4.68  10th 11 13 0∗ 1 5 14 12 0∗ 8 10 9 4 0∗ 3 7 6 2 4.78  11th 2 10 13 0∗ 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4 5.21  12th 9 10 0∗ 5 1 11 6 7 0∗ 4 14 12 0∗ 13 8 2 3 5.36  13th 10 6 4 0∗ 2 7 14 0∗ 13 1 8 5 0∗ 3 12 9 11 4.25  14th 1 10 0∗ 7 2 3 9 13 0∗ 14 6 0∗ 5 8 4 12 11 6.25  15th 8 14 6 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12 5.48  16th 3 11 10 0∗ 4 7 12 0∗ 13 5 14 2 0∗ 6 8 9 1 5.36  17th 4 1 14 9 0∗ 2 10 7 0∗ 13 12 8 6 0∗ 3 5 11 7.21  18th 6 10 14 0∗ 11 12 4 3 0∗ 1 2 13 0∗ 9 7 8 5 4.87  19th 11 1 0∗ 10 8 14 0∗ 3 13 2 9 0∗ 5 12 4 7 6 4.98  20th 9 11 6 0∗ 7 10 14 0∗ 5 1 2 0∗ 13 4 3 8 12 4.56


## 6.3. Calculation of the Historical and Global Optimal Values of the Particle

The historical and global optimal values are the particle’s cognitive and social messages, which have a significant impact on the probability of finding the preferable particle. Hence, the definition for the historical optimal particle is as follows: particles with the best historical fitness and global optimal particle; the best particles with the best historical fitness are mentioned. According to the definition, we obtain all the particles’ historical optimal values and the particle swarm’s global value at each iteration, as shown in Table9.Table 9
Historically optimal and globally optimal particles.
x 1,1 t x 1,2 t x 1,3 t ⋯ x 1 , l e n g t h - 2 t x 1 , l e n g t h - 1 t x 1 , l e n g t h t f 1 t 1st Particle’s Personal Best  x 2,1 k ′ x 2,2 k x 2,3 k ⋯ x 2 , l e n g t h - 2 k x 2 , l e n g t h - 1 k x 2 , l e n g t h k f 2 k 2nd Particle’s Personal Best  ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯  x 20,1 l ′ x 20,2 l x 20,3 l ⋯ x 20 , l e n g t h - 2 l x 20 , l e n g t h - 1 l x 20 , l e n g t h l f 20 l 20th Particle’s Personal Best  Global Best min ⁡ f i t ∣ i ∈ 1 , p o p s i z e , t ∈ 1 , M a x i t e r a t eIn the first iteration, we initialize the original fitness value of the particles as the personal best and initialize the minimum of the personal best as the global best. The consequence is represented in Table10.Table 10
Initial fitness value, personal best, and global best of the particles.
Initial Particles fitness Personal best 1st 6 3 11 7 0∗ 14 8 5 0∗ 1 2 4 13 0∗ 9 10 12 4.89 4.89  2nd 7 1 13 2 0∗ 14 6 10 0∗ 12 11 4 0∗ 8 3 9 5 4.70 4.70  3rd 2 10 4 0∗ 5 3 8 0∗ 12 11 7 1 0∗ 13 14 6 9 4.58 4.58  4th 10 9 0∗ 6 13 2 1 0∗ 11 14 3 8 7 0∗ 4 5 12 3.25 3.25  (Global best)  5th 8 9 x 0∗ 12 2 10 3 0∗ 7 4 1 13 11 0∗ 5 6 4.21 4.21  6th 2 3 0∗ 11 1 6 14 0∗ 5 9 8 0∗ 4 13 12 7 10 5.24 5.24  7th 9 2 1 0∗ 6 11 0∗ 13 10 8 14 3 0∗ 7 4 12 5 6.35 6.35  8th 14 7 10 8 0∗ 2 6 12 0∗ 1 3 4 5 0∗ 11 9 13 4.25 4.25  9th 3 7 1 0∗ 13 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2 4.68 4.68  10th 11 13 0∗ 1 5 14 12 0∗ 8 10 9 4 0∗ 3 7 6 2 4.78 4.78  11th 2 10 13 0∗ 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4 5.21 5.21  12th 9 10 0∗ 5 1 11 6 7 0∗ 4 14 12 0∗ 13 8 2 3 5.36 5.36  13th 10 6 4 0∗ 2 7 14 0∗ 13 1 8 5 0∗ 3 12 9 11 4.25 4.25  14th 1 10 0∗ 7 2 3 9 13 0∗ 14 6 0∗ 5 8 4 12 11 6.25 6.25  15th 8 14 6 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12 5.48 5.48  16th 3 11 10 0∗ 4 7 12 0∗ 13 5 14 2 0∗ 6 8 9 1 5.36 5.36  17th 4 1 14 9 0∗ 2 10 7 0∗ 13 12 8 6 0∗ 3 5 11 7.21 7.21  18th 6 10 14 0∗ 11 12 4 3 0∗ 1 2 13 0∗ 9 7 8 5 4.87 4.87  19th 11 1 0∗ 10 8 14 0∗ 3 13 2 9 0∗ 5 12 4 7 6 4.98 4.98  20th 9 11 6 0∗ 7 10 14 0∗ 5 1 2 0∗ 13 4 3 8 12 4.56 4.56 Note: Numbers against asterisk mark denote the depot and the end of one route.


## 6.4. Local Search through the Variable Neighborhood Operation

Our algorithm utilizes the variable neighborhood search algorithm given in Table11 to enhance the ability to choose the solution randomly, ensure diversity of the individual, and overcome the disadvantage that the solution is a local optimum. First, we updated the individual particles’ historical best and global best. Then, in our implementation, three variable neighborhood structures are applied to the particle’s personal best mentioned at Section 5.4.Table 11
Variable neighborhood operation.
x 1,1 t ′ x 1,2 t ′ x 1,3 t ′ ⋯ x 1 , l e n g t h - 2 t ′ x 1 , l e n g t h - 1 t ′ x 1 , l e n g t h t ′ f 1 t ′ 1st Particle’s Personal Best  x 2,1 k ′ x 2,2 k ′ x 2,3 k ′ ⋯ x 2 , l e n g t h - 2 k ′ x 2 , l e n g t h - 1 k ′ x 2 , l e n g t h k ′ f 2 k ′ 2nd Particle’s Personal Best  ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯  x 20,1 l ′ x 20,2 l ′ x 20,3 l ′ ⋯ x 20 , l e n g t h - 2 l ′ x 20 , l e n g t h - 1 l ′ x 20 , l e n g t h l ′ f 20 l ′ 20th Particle’s Personal Best  After variable neighborhood operation Global Best min ⁡ f i t ′ ∣ i ∈ 1 , p o p s i z e , t ′ ∈ 1 , M a x i t e r a t e Note:xi,jk′ is the result of the variable neighborhood operation.Based on the particular variable neighborhood structure, we used the personal best as an operation by taking the 2-opt operation structure as an example. The steps and consequence are the following.Procedure 1.
Produce a group of 2-opt structures randomly; i.e., in the first particle, there are three points, (3, 7), (14, 5), and (2, 10), respectively.Procedure 2.
Subject the particle to the 2-opt operation. An example using the1st particle is provided in Table 12.Table 12
Particles subjected to neighborhood structure operation (2-Opt).
1st 6 3 11 7 0∗ 14 8 5 0∗ 1 2 4 13 0∗ 9 10 12  After 2-opt-1st 6 3 7& 11& 0∗ 14 5& 8& 0∗ 1 2 10& 9& 0& 13& 4& 12  2nd 7 1 13 2 0∗ 14 6 10 0∗ 12 11 4 0∗ 8 3 9 5  3rd 2 10 4 0∗ 5 3 8 0∗ 12 11 7 1 0∗ 13 14 6 9  4th 10 9 0∗ 6 13 2 1 0∗ 11 14 3 8 7 0∗ 4 5 12  5th 8 9 14 0∗ 12 2 10 3 0∗ 7 4 1 13 11 0∗ 5 6  6th 2 3 0∗ 11 1 6 14 0∗ 5 9 8 0∗ 4 13 12 7 10  7th 9 2 1 0∗ 6 11 0∗ 13 10 8 14 3 0∗ 7 4 12 5  8th 14 7 10 8 0∗ 2 6 12 0∗ 1 3 4 5 0∗ 11 9 13  9th 3 7 1 0∗ 13 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2  10th 11 13 0∗ 1 5 14 12 0∗ 8 10 9 4 0∗ 3 7 6 2  11th 2 10 13 0∗ 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4  12th 9 10 0∗ 5 1 11 6 7 0∗ 4 14 12 0∗ 13 8 2 3  13th 10 6 4 0∗ 2 7 14 0∗ 13 1 8 5 0∗ 3 12 9 11  14th 1 10 0∗ 7 2 3 9 13 0∗ 14 6 0∗ 5 8 4 12 11  15th 8 14 6 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12  16th 3 11 10 0∗ 4 7 12 0∗ 13 5 14 2 0∗ 6 8 9 1  17th 4 1 14 9 0∗ 2 10 7 0∗ 13 12 8 6 0∗ 3 5 11  18th 6 10 14 0∗ 11 12 4 3 0∗ 1 2 13 0∗ 9 7 8 5  19th 11 1 0∗ 10 8 14 0∗ 3 13 2 9 0∗ 5 12 4 7 6  20th 9 11 6 0∗ 7 10 14 0∗ 5 1 2 0∗ 13 4 3 8 12 Note: Numbers against asterisk mark denote the depot and the end of one route. Numbers with & mark were subjected to operation by the 2-opt structure.The remaining particles were subjected to operation according to the form procedure, and the final result is presented in Table13.Table 13
Particles subjected to 2-Opt structure operation.
After 2-opt-1st 6 3 7& 11& 0∗ 14 5& 8& 0∗ 1 2 10& 9& 0& 13& 4& 12  After 2-opt-2nd 7 13& 1& 2 0∗ 14 4& 11& 12& 0& 10& 6& 0∗ 8 5& 9& 3&  After 2-opt-3rd 4& 10& 2& 0∗ 5 3 7& 11& 12& 0& 8& 1 0∗ 13 14 6 9  After 2-opt-4th 10 9 0∗ 6 1& 2& 13& 0∗ 11 14 3 8 5& 4& 0& 7& 12  After 2-opt-5th 8 14& 9& 0∗ 12 2 10 3 0∗ 7 4 11& 13& 1& 0∗ 5 6  After 2-opt-6th 2 3 0∗ 11 1 6 14 0∗ 5 9 4& 0& 8& 13 12 7 10  After 2-opt-7th 9 6& 0& 1& 2& 11 0∗ 13 10 8 14 4& 7& 0& 3& 12 5  After 2-opt-8th 14 7 10 2& 0& 8& 6 12 0∗ 5& 4& 3& 1& 0& 11 9 13  After 2-opt-9th 3 7 13& 0& 1& 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2  After 2-opt-10th 11 13 0∗ 1 14& 5& 12 0∗ 8 10 9 2& 6& 7& 3& 0& 4&  After 2-opt-11th 2 0& 13& 10& 7 9 0∗ 12 14 1 6 5 11 0∗ 8 3 4  After 2-opt-12th 9 10 0∗ 1& 5& 11 6 4& 0& 7& 14 12 0∗ 13 8 2 3  After 2-opt-13th 10 2& 0& 4& 6& 7 14 0∗ 13 1 5& 8& 0∗ 3 12 9 11  After 2-opt-14th 1 10 7& 0& 2 3 9 13 14& 0& 6 0∗ 5 8 4 12 11  After 2-opt-15th 8 6& 14& 0∗ 1 7 9 4 0∗ 2 11 0& 13& 10 3 5 12  After 2-opt-16th 3 11 10 0∗ 4 7 12 0∗ 13 5 0& 2& 14& 6 8 9 1  After 2-opt-17th 4 9& 14& 1& 0∗ 2 10 7 0∗ 13 12 8 3& 0& 6& 5 11  After 2-opt-18th 6 10 11& 0& 14& 12 4 3 0∗ 1 2 13 7& 9& 0& 8 5  After 2-opt-19th 11 1 0∗ 10 8 9& 2& 13& 3& 0& 14& 0∗ 5 12 4 7 6  After 2-opt-20th 9 11 7& 0& 6& 10 14 0∗ 5 1 8& 3& 4& 13& 0& 2& 12 Note: Numbers against asterisk mark denote the depot and the end of one route. Numbers with & mark were subjected to operation by the 2-opt structure.Then, according to the constraint of the model, we check and recode particles. If the operated particles can satisfy the constraint, we accept these particles; otherwise, we recode the operated particles. The result of checking and recoding can be seen in Table14.Table 14
Recoding particles after the 2-Opt operation.
After 2-opt-1st 6 3 7 0∗ 11 14 5 8 0∗ 1 2 10 9 0∗ 13 4 12  After 2-opt-2nd 7 13 0∗ 1 2 14 4 0∗ 11 12 10 6 0∗ 8 5 9 3  After 2-opt-3rd 4 10 2 0∗ 5 3 7 11 12 0 8 1 0∗ 13 14 6 9  After 2-opt-4th 10 9 0∗ 6 1 2 13 0∗ 11 14 3 8 5 4 0∗ 7 12  After 2-opt-5th 8 14 9 0∗ 12 2 10 3 0∗ 7 4 11 13 1 0∗ 5 6  After 2-opt-6th 2 3 0∗ 11 1 6 14 0∗ 5 9 4 0∗ 8 13 12 7 10  After 2-opt-7th 9 6 1 2 0∗ 11 13 10 8 0∗ 14 4 7 0∗ 3 12 5  After 2-opt-8th 14 7 10 2 0∗ 8 6 12 0∗ 5 4 3 1 0∗ 11 9 13  After 2-opt-9th 3 7 13 0∗ 1 9 14 8 0∗ 11 4 10 0∗ 5 6 12 2  After 2-opt-10th 11 13 0∗ 1 14 5 12 0∗ 8 10 9 2 6 7 3 0∗ 4  After 2-opt-11th 2 13 10 0∗ 7 9 12 14 0∗ 1 6 5 11 0∗ 8 3 4  After 2-opt-12th 9 10 0∗ 1 5 11 6 4 0∗ 7 14 12 0∗ 13 8 2 3  After 2-opt-13th 10 2 4 0∗ 6 7 14 0∗ 13 1 5 8 0∗ 3 12 9 11  After 2-opt-14th 1 10 7 0∗ 2 3 9 13 14 0∗ 6 5 8 0∗ 4 12 11  After 2-opt-15th 8 6 14 0∗ 1 7 9 4 0∗ 2 11 13 0∗ 10 3 5 12  After 2-opt-16th 3 11 0∗ 10 4 7 12 0∗ 13 5 2 0∗ 14 6 8 9 1  After 2-opt-17th 4 9 14 1 0 2 10 7 0∗ 13 12 8 3 0∗ 6 5 11  After 2-opt-18th 6 10 11 0∗ 14 12 4 3 0∗ 1 2 13 7 9 0∗ 8 5  After 2-opt-19th 11 1 0∗ 10 8 9 2 13 3 0∗ 14 5 12 0∗ 4 7 6  After 2-opt-20th 9 11 7 0∗ 6 10 14 0∗ 5 1 8 3 4 13 0∗ 2 12 Note: Numbers against asterisk mark denote the depot and the end of one route.


## 6.5. Particle Updating Rules and Algorithm Termination Condition

The velocity and position of a particle are updated according to (21) and (22), after which iterative operation occurs. The particles are updated according to the relative operation relating to the position and velocity. The results are provided in Table 15. Subsequently, the algorithm calculates the fitness value, personal best, and global best of the updated particles until fitnessn+1-fitness(n)/fitnessn+1≤0.05 or the iterative value exceeds the value of the original iterative step.Table 15
Updated particles.
1st 11 2 9 0∗ 3 5 8 12 13 0∗ 10 14 7 0∗ 1 6 4  2nd 7 8 0∗ 12 10 13 6 4 0∗ 14 3 11 9 0∗ 1 5 2  3rd 1 11 4 ∗ 7 9 10 6 0∗ 8 13 3 0∗ 5 14 2 12  4th 5 14 0∗ 10 9 13 4 6 0∗ 2 3 11 0∗ 12 1 7 8  5th 8 12 2 3 0∗ 4 5 7 6 0∗ 10 11 13 0∗ 9 14 1  6th 8 7 5 0∗ 11 10 12 6 13 4 0∗ 14 3 2 0∗ 9 1  7th 1 7 12 0∗ 8 10 5 6 14 0∗ 3 9 2 0∗ 4 11 13  8th 11 8 14 9 0∗ 3 13 10 4 0∗ 5 12 1 0∗ 6 7 2  9th 6 13 0∗ 10 1 14 12 0∗ 4 9 11 3 0∗ 8 7 2 5  10th 6 4 2 0∗ 9 12 8 3 10 0∗ 14 5 11 0∗ 1 7 13  11th 4 2 0∗ 9 8 13 1 0∗ 3 5 12 7 0∗ 6 10 14 11  12th 4 5 12 2 0∗ 1 8 14 7 0∗ 6 13 10 9 0∗ 3 11  13th 14 2 7 0∗ 8 4 1 10 9 0∗ 5 13 6 0∗ 3 12 11  14th 1 2 8 9 0∗ 6 11 14 10 12 0∗ 5 7 13 0∗ 3 4  15th 6 14 13 0∗ 2 7 12 8 0∗ 10 9 3 1 0∗ 11 5 4  16th 14 10 8 0∗ 2 5 12 11 3 0∗ 7 13 4 0∗ 6 1 9  17th 2 13 3 5 0∗ 1 8 11 6 0∗ 12 9 0∗ 10 4 14 7  18th 8 13 1 2 0∗ 7 14 6 0∗ 10 11 4 12 0∗ 9 5 3  19th 6 5 11 0∗ 3 13 4 14 0∗ 7 10 12 0∗ 8 2 9 1  20th 14 1 2 0∗ 3 12 4 0∗ 13 7 8 6 5 0∗ 11 10 9 Note: Numbers against asterisk mark denote the depot and the end of one route.When the conditions of convergence are met, the operation is complete. Otherwise, return to process in Section6.1 and implement the iterative operation again.


## 6.6. Computational Result and Analysis

The parameters of the hybridization algorithm of DPSO with VNS are included in Table16.Table 16
Parameters of the algorithm.
Parameters’ name of Hybridization algorithm of DPSO with VNS Value Population of particles (N) 20  The length of per particle (L) 14  Weight Value (C) 0.6  Iteration 50  Learning Value (CP) 1.5  Social Value (Cg) 2The dispatching plan is presented in Table17 and the dispatching routes are shown in Figure 6.Table 17
Dispatching plan.
Route-ID Routes Cost of per route R1 9-2-1-8 5.4 km  R2 3-13-14 4.7 km  R3 10-11-12 6.2 km  R4 4-6-5-7 9.33 km  Total Cost 25.63 kmFigure 6
Dispatching route of the hybridization algorithm of DPSO with VNS.Computation using the hybridization algorithm of DPSO with VNS produces four routes for dispatching, route 1=[9-2-1-8], route 2= [3-13-14], route 3= [10-11-12], and route 4 = [4-6-5-7]. Furthermore, the cost per route scaled by the distance corresponds to 5.4 km, 4.7 km, 6.2 km, and 9.33 km. The total cost is 25.63 km.In order to test the performance of the hybridization algorithm of DPSO with VNS, Goksal et al. made an overall comparison with Reactive Tabu Search (RTS), Large Neighborhood Search (LNS), Particle Swarm Optimization (PSO), Ant Colony System (ACS), Parallel Iterative Local Search (PILS), and Adaptive Memory Methodology (AMM). Since they are the most effective heuristics proposed in the literature. According to the indexes of average deviation to the best known solution, average computation time, scaled average computation time over 40 instances, and number of best known solutions found, the comprehensive performance of hybridization algorithm of DPSO with VNS is the best. For the computational results of the comparison, see Goksal et al. [21].In this study, we implement Discrete Particle Swarm Optimization (DPSO) and Genetic algorithm (GA) separately over the case study data set to get a more comprehensive view of the performance of the hybridization algorithm of DPSO with VNS. The results are provided in Figures7 and 8. DPSO is considered in the comparison with the purpose of presenting the performance of VNS. We can see from Figure 7 that the cost of the solution without VNS is 27.07 km. The result of the hybridization algorithm of DPSO with VNS (25.63 km) dropped by 5.3%. GA [15] is a population-based heuristic which can be used to solve the problem like SBRP. Y. Li et al. reported the performance in their study. Thus, in order to make a direct comparison, we run a GA program over the same data set of case study. From Figure 8, we can see that the total cost is 26.33 km, while the hybridization algorithm of DPSO with VNS brings a drop by 2.7%. This analysis indicates that the hybridization algorithm of DPSO with VNS outperforms the other existing heuristics in the SBRP literature.Figure 7
Result of DPSO algorithm.Figure 8
Result of GA.


## 7. Concluding Remarks

We quantify the station state using the bicycle net flow by calculating the difference between the inflow and the outflow. The maximum daily demand can be determined by the minimum value on the state curve. It should be noted that this minimum daily demand reveals the imbalance among the bicycle arrival and departure, rather than the real demand of the static rebalancing.An extended VRPSPD program is developed to model the bike-sharing static rebalancing problem considering the number of bicycles that need repair. We employed a hybrid DPSO-VNS algorithm to solve the proposed model. The computational result shows that the effectiveness of the algorithm is attributed to a combination of the following reasons. First, the idea of vehicle orientation ensures that each route only covers a restricted area. Second, the quality of the solution is improved by the cheapest insertion heuristic and 2-opt method, both of which are applied during the route construction. Third, the DPSO mechanism is capable of generating diverse solutions and maintaining the best solution found during the iterative process. Moreover, this approach is flexible, because it is observed that it can be applied to the related problems by modifying the initial solution and a feasibility check of neighborhood methods. The approach can be easily extended and adjusted by adding or replacing these neighborhood methods.Some aspects may improve the performance of the proposed algorithm, such as the parameter optimization and programming implementation. Although the PSO parameter set used in this study was obtained from a preliminary experiment, it may not be the best one. In addition, the implementation of the algorithm may be further improved. These aspects would be further investigated in future studies.



---

*Source: 8086378-2018-09-13.xml*