# Adaptive Cultural Algorithm-Based Cuckoo Search for Time-Dependent Vehicle Routing Problem with Stochastic Customers Using Adaptive Fractional Kalman Speed Prediction

Authors: H. Xue
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

For the Time-Dependent Vehicle Routing Problem with Stochastic Customers (TDVRPSC), an adaptive Cultural Algorithm-Based Cuckoo Search (CACS) has been proposed in this paper. The convergence of the new algorithm is proved. An adaptive fractional Kalman filter (AFKF) for traffic speed prediction is proposed. An adaptive mechanism for choosing the covariance of state noise is designed. Its mathematical process is proved. Several benchmark instances with different scales are tested, and new solutions are discovered, which are better than the published solutions. The effects of the parameters on the convergence and the results are studied. According to cargo weight of customers to be delivered, the customers can be divided into large, small, and retail customers. The algorithm is tested with fixed demand probability and also different customer types with stochastic demand. The traffic speeds in different business districts in Xiamen at different times are predicted by AFKF. The results show that AFKF has smaller prediction error and better prediction accuracy than fractional Kalman filter and Kalman filter. The effect of different fractional orders on prediction error is compared. The performance of the new algorithm is compared with that of the cultural algorithm and the Cuckoo Search. The result shows that the new algorithm can efficiently and effectively solve DTVRPSC and improve the accuracy of vehicle routing planning of time-varying actual urban traffic road.


---

## Body

## 1. Introduction

With the rapid development of the logistics industry, Vehicle Routing Problem (VRP) gets more and more attention. The traditional VRP uses static road network model and cannot accurately estimate the travelling time based on the actual changing traffic conditions. In the actual traffic applications, there are many uncertainties, such as random customer demand, random number of customers, random traffic speed, random conditions of road traffic, and rush hours. Therefore, Time-Dependent Vehicle Routing Problem with Stochastic Customers (TDVRPSC) is a stochastic optimization problem with extensive practical significance. It is important to study new intelligent algorithms to realize logistics automation under dynamic traffic information. The results of this paper can be applied to vehicle routing, artificial intelligence, and data filtering.Guo studied dynamic multiobjective Vehicle Routing Problem with a corresponding carbon emission model and it was set as an optimization objective [1]. After finding optimal robust virtual routes for all customers by adopting multiobjective particle swarm optimization in the first phase, static vehicle routes for static customers are formed by removing all dynamic customers from robust virtual routes in next phase. The dynamically appearing customers append to be served according to their service time and the vehicles’ statues. Global vehicle routing optimization is triggered only when no suitable locations can be found for dynamic customers. A metric measuring the algorithms robustness is given. Guo proposed an adaptive immune clonal selection cultural algorithm [2]. Dual structure of cultural algorithm was adopted, and a hybrid selection strategy integrating selection and clonal selection was put forward. The proportion of population influenced by each selection method was adaptively adjusted according to implicit knowledge extracted from the evolution process. Guo proposed a novel multipopulation cultural algorithm adopting knowledge migration derived from human cultural interaction [3]. Knowledge extracted from the evolution process contained more effective information. By migrating knowledge among subpopulations at constant intervals, the algorithm realized more effective interaction with less communication cost. Geng presented a modified centralized algorithm based on particle swarm optimization (MCPSO) [4]. Three distributed algorithms were compared against three centralized algorithms. Results showed that MCPSO could be used as a benchmark for distributed algorithms for determining the performance gap.Random VRP was first proposed by Stewart and Golden and has gotten more and more attention [5]. Bertsimas considered a natural probabilistic variation of the classical vehicle routing problem and used stochastic demands [6]. Bianchi et al. analyzed the performance of metaheuristics on the Vehicle Routing Problem with stochastic demands [7]. Bozorgi-Amiri et al. developed a multiple objective robust stochastic programming approach for disaster relief logistics under uncertainty [8]. Sungur et al. introduced a robust optimization approach to solve the VRP with demand uncertainty [9]. Tan et al. considered VRP with stochastic demand, where actual demand was revealed only when the vehicle arrived at the customer [10]. Ibarra-Rojas et al. studied a VRP to optimize accessibility based on six indicators [11]. Wollmer et al. used a cutting plane algorithm for solving the stochastic programming and proved the algorithm to be finite [12]. Kenyon and Morton used Monte Carlo sampling based solution procedure for stochastic VRP with large sample spaces [13]. The solution space of the classical NP-hard problem is of large scale and complex. Large-scale networks may be found in many fields and a crucial step is to estimate them from the data. Sparse reciprocal graphical models have been considered for learning large-scale networks using efficient algorithms. Alpago et al. proposed an identification procedure of a sparse graphical model associated with a Gaussian stationary stochastic process [14].At present, perfect solutions can only be reached as close as possible within a reasonable run time. Therefore, new intelligent methods are required. In 2009, Yang and Deb proposed Cuckoo Search (CS) to solve optimization problems [15]. Mareli and Twala developed three Cuckoo Search algorithms based on dynamically changing switch parameters [16]. Ishak Boushaki et al. proposed a quantum chaotic Cuckoo Search algorithm for data clustering [17]. Ayoubi et al. proposed a nonhomogeneous Cuckoo Search algorithm for allocation of passive filters [18].To accelerate the convergence speed, this paper adopts the double mechanism in cultural algorithm (CA) [19] and proposed a Cultural Algorithm Based Cuckoo Search (CACS). CACS can overcome the defects of slow convergence rate of CS, as well as the blind guidance of CA group space, which makes the whole population algorithm not easy to converge. Its adaptive mechanism can overcome the dependence of parameters on algorithm performance. In CACS, the individuals with high fitness form the corresponding culture of the group in the Belief space by accepting the rules. The representative individuals in the population, utilizing the outstanding individuals after adjustment, will further affect the evolution of next-generation population.Kalman filter was proposed by Kalman in 1960 [20]. In recent years, it has been widely applied to adaptive control and system recognition [21]. With the presence of uncertainties, the robust Kalman filter based on the tau-divergence which accounts for uncertainties could be used. San Gultekin et al. considered the nonlinear Kalman filtering problem using α-divergence measures as optimization criteria [22]. In order to extend the application of Kalman filter from integer-order to fractional-order system, Solís-Pérez et al. presented a fractional-order Kalman filter according to the Riemann–Liouville definition [23]. Sun et al. used the extended Kalman filtering for fractional-order nonlinear discrete system with Lévy noises [24]. Liu et al. proposed a generalized fractional central difference Kalman filter for nonlinear discrete fractional dynamic systems [25].In this paper, an adaptive fractional Kalman filter for traffic speed prediction is proposed to update and estimate the time-varying vehicle speed. An adaptive mechanism is adopted in computing the covariance of prediction error. Its mathematical process is proved. A solution of classic VRP cases that is better than those optimal solutions published before is found. Using AFKF and real-time urban traffic data, vehicle speed is predicted to adjust vehicle routing in time to avoid congested roads and improve the accuracy of VRP. The contributions of the paper are summarized as follows:(1)
An adaptive Cultural Algorithm Based Cuckoo Search is proposed(2)
An adaptive fractional Kalman filter for traffic speed prediction is proposed(3)
The convergence of CACS and the mathematical process of AFKF are proved(4)
New solutions of VRP are discovered, which are better than the published solutions(5)
With real-time urban traffic data, vehicle speed is dynamically predicted to adjust the solution scheme of TDVRPSC in time


## 2. Mathematical Model of TDVRPSC

### 2.1. Mathematical Model

Stochastic time-dependent road network can be represented by a directed networkGC,E,vijk. C represents a set of customer nodes, and E represents a set of road segments. V=N.vijk represents the travel speed from customer node i to j during time period k.The goal is to obtain the shortest expected travel time and service time.p0<p<1 denotes the customer demand probability. The total number of vehicles is K. tij denotes the travelling time from customer i to j. x¯ijl denotes whether vehicle l is responsible for the route from customer j to i. uil denotes whether customer i is serviced by vehicle l. τil denotes the service time of customer i by vehicle l. qi is the demand of the customer i.Wl is the total load of the vehicle l. The model can be described as follows:(1)minEmax1≤l≤K∑1≤i,j≤N,j≠itijlx¯ijl+∑1≤i≤Npτiluil,(2)s.t.∑1≤i≤N,i≠j∑1≤l≤Kxijl=1,1≤j≤N,(3)∑1≤j≤N,i≠j∑1≤l≤Kxijl=1,1≤i≤N,(4)∑1≤l≤Kuil=1,1≤i≤N,(5)uil≤∑1≤j≤Nx¯ijl=1,1≤i≤N,1≤l≤K,(6)∑1≤i≤Nqiuil≤W,1≤l≤Kl,(7)x¯ijl=1,lselects arci,j,1≤i,j≤N,1≤l≤K,0,else,(8)uil=1,iserviced byl,1≤i≤N,1≤l≤K,0,else.Hypothesis is stated as follows:(1)
Each customer is serviced exactly once. The loading and unloading services of each customer are completed only once. Equations (2) and (3) ensure that each customer is served exactly once.(2)
Each customer is loaded immediately after being visited by a vehicle. The total time includes travelling time and service time with no other times wasted.(3)
The total weight of goods demanded by customers cannot exceed the vehicle load servicing them. Equation (6) ensures the total weight of goods demanded by customers cannot exceed the vehicle load servicing them.(4)
Each customer is serviced by a vehicle only if he is visited on the route of the vehicle. Equation (5) ensures a customer is serviced by a vehicle only if he is visited on the route of the vehicle.(5)
Each customer can only be serviced by exactly one vehicle. Equation (4) ensures each customer is serviced exactly one vehicle.(6)
The positions of all customers are known.(7)
The position of the depot is known. The source of vehicles is determinant.(8)
The distances from the depot to each customer are known.(9)
The customer demand probability is a stochastic real number between zero and one.(10)
The travelling time between two customers is stochastic time-dependent due to the time-varying speed and real-time traffic conditions on roads.The existence of at least one solution of the above problem is proved in literature [12, 13]. A cutting plane algorithm for solving the above stochastic programming was proved to be finite [8]. The above stochastic objective function (1) was proved to be a convex function on the convex hull of the set of decision vectors satisfying constraints (2)–(8) [9]. Since the local minimum of convex optimization problem is the global minimum, the local minimum solution of the above problem will also be the global minimum.


### 2.2. Adaptive Fractional-Order Kalman Filter for Speed Prediction

The distance between customeri and j is represented by dij. Δt denotes the sampling time interval. The calculation procedure for travelling time tij is shown in Figure 1. vijk is updated and estimated by Kalman filter algorithm.Figure 1
Calculation procedure for travelling time.The fractional-order difference is formulated as(9)Δαxk=1hα∑j=0k−1jγjxk−j,(10)γj=αj=1,j=0,αα−1⋯α−j+1j!,j>0,where α is the fractional order. h is the sampling interval. k is the number of samples. The fractional-order model of speed prediction is given by(11)Δαxk+1=Akxk+wk,(12)xk+1=Δαxk+1−∑j=1k+1−1jγjxk+1−j,(13)zk=Gkxk+εk,where xk is the state variable of the traffic speed. zk is the measurement output. wk is the system state noise. εk is the measurement noise. x˜k is calculated as(14)x˜k+1=Exk+1|zk∗,where zk∗ is the measurement sequence containing measurement output z1,z2,⋯,zk. P˜k denotes the covariance of prediction error at time instant k, which is calculated as(15)P˜k+1=Exk+1−x˜k+1xk+1−x˜k+1T.Qk denotes the covariance of system noise at the time instant k, which is calculated as(16)Qk=Ewk−1wk−1T.Rk denotes the covariance of measurement noise at the time instant k, which is calculated as(17)Rk=Eεk−1εk−1T.x^k denotes the state estimation at time instant k, which is calculated as(18)x^k+1=Exk+1|zk+1∗.P^k is the covariance of estimation error at the time instant k, which is calculated as(19)P^k+1=Exk+1−x^k+1xk+1−x^k+1T.The covariance of state noise is designed as follows:(20)Qk=λEwk−1wk−1T.where λ>0 is a gain coefficient.When the covariance of state noise is small, the convergence speed is slow and the convergence accuracy is high. If the state changes too fast in a short time, the slow convergence rate will easily lead to tracking failure. When the covariance of state noise is large, the convergence speed is fast and the accuracy is low. Before reaching a stable state, choosing a larger covariance of state noise can accelerate the convergence speed. After reaching a stable state, the convergence accuracy can be improved by choosing smaller covariance of state noise.In standard Kalman filter,λ is a constant. However, it is difficult to determine an optimal gain in the whole process. When the covariance of state noise is small, the number of estimations will update slowly if λ is too small. When the error changes sharply, the estimation will oscillate and even the convergence process will be unstable if λ is too large. To accelerate the convergence process, an adaptive gain λ is used to make the estimation process adjusted reasonably and adapt to different data characteristics. λ increases when the covariance of state noise is small and decreases when the covariance of state noise is large. According to the noise covariance information, λ is calculated dynamically. Therefore, an adaptive mechanism for choosing λ is designed as follows:(21)λ=λmax1+P^k+12.The Kalman gainKk is calculated as(22)Kk+1=P˜k+1Gk+1TGk+1P˜k+1Gk+1T+Rk+1−1.The fractional Kalman filter is formulated as(23)x˜k+1=Akx^k−∑j=1k+1−1jγjx^k+1−j,(24)Rk+1=zk+1−Gk+1x˜k+1zk+1−Gk+1x˜k+1T−Gk+1P˜k+1Gk+1T,(25)x^k+1=x˜k+1+Kk+1zk+1−Gk+1x˜k+1,(26)P^k+1=1−Kk+1Gk+1P˜k+11−Kk+1Gk+1T+Kk+1Rk+1Kk+1T,(27)P˜k+1=Ak+γ1P^kAk+γ1T+Qkλ+∑j=2k+1γjP^k+1−jγjT.Assumption 1.
Ewk=0. wk is not related to zk.Assumption 2.
Eεk=0. εk is not related to xk and x˜k.Theorem 1.
For the fractional-order system with system state noise and measurement noise described by formulas (11)–(13), AFKF is given by formulas (23)–(27).Proof.
Substituting (11) and (12) into (14), one can obtain(28)x˜k+1=EAkxk+wk−∑j=1k+1−1jγjxk+1−j|zk∗=Akx^k+AkExk|zk∗−x^k+Ewk|zk∗−∑j=1k+1−1jγjExk+1−j|zk∗.
From Assumption1, one can obtain(29)Ewk|zk∗=∑i=1∞wkPw=wk|z=zk∗=∑i=1∞wkPw=wk,z=zk∗Pz=zk∗=∑i=1∞wkPw=wkPz=zk∗Pz=zk∗=∑i=1∞wkPw=wk=Ewk=0.
Substituting (29) into (28), one can obtain(30)x˜k+1=Akx^k+AkExk|zk∗−x^k−∑j=1k+1−1jγjExk+1−j|zk∗.
Substituting (18) into (29), one can obtain(31)x˜k+1=Akx^k+Akx^k−x^k−∑j=1k+1−1jγjExk+1−j|zk∗=x^k−∑j=1k+1−1jγjExk+1−j|zk∗=x^k−∑j=1k+1−1jγjExk+1−j|zk+1−j∗.
Substituting (18) into (31), one can obtain(32)x˜k+1=Akx^k−∑j=1k+1−1jγjx^k+1−j.
Formula (23) is obtained.
Substituting (11) and (12) into (15), one can obtain(33)P˜k+1=Exk+1−x˜k+1xk+1−x˜k+1T=EAkxk+wk−∑j=1k+1−1jγjxk+1−j−x˜k+1Akxk+wk−∑j=1k+1−1jγjxk+1−j−x˜k+1T.
Substituting (32) into (33), one can obtain(34)P˜k+1=EAkxk+wk−∑j=1k+1−1jγjxk+1−j−Akx^k+∑j=1k+1−1jγjx^k+1−jAkxk+wk−∑j=1k+1−1jγjxk+1−j−Akx^k+∑j=1k+1−1jγjx^k+1−jT=Ak+γ1Exk−x^kxk−x^kTAk+γ1T+∑j=2k+1γjExk+1−j−x^k+1−jxk+1−j−x^k+1−jTγjT+Ewk−1wk−1T.
Substituting (19) and (20) into (34), one can obtain(35)P˜k+1=Ak+γ1P^kAk+γ1T+Qkλ+∑j=2k+1γjP^k+1−jγjT.
Formula (27) is obtained.
The cost function is defined by(36)x^k+1=argminxk+1x˜k+1−xk+1P˜k+1−1x˜k+1−xk+1T+zk+1−Gk+1x˜k+1−Gk+1xk+1−x˜k+1Rk+1−1⋅zk+1−Gk+1x˜k+1−Gk+1xk+1−x˜k+1T.
Differentiating (36) with respect to xk+1 and making the result equal to zero, one can obtain(37)P˜k+1−1x^k+1−x˜k+1=Gk+1TRk+1−1zk+1−Gk+1x˜k+1−Gk+1x^k+1−x˜k+1.
From (37), one can obtain(38)x^k+1=Rk+1Gk+1−1P˜k+1−1+Gk+1−1⋅Rk+1Gk+1−1P˜k+1−1x˜k+1+zk+1−Gk+1x˜k+1+Gk+1x˜k+1=x˜k+1+Rk+1Gk+1−1P˜k+1−1+Gk+1−1zk+1−Gk+1x˜k+1=x˜k+1+P˜k+1Gk+1TGk+1P˜k+1Gk+1T+Rk+1−1zk+1−Gk+1x˜k+1.
Substituting (22) into (38) yields(39)x^k+1=x˜k+1+Kk+1zk+1−Gk+1x˜k+1.
Formula (25) is obtained.
Substituting (39) into (19) yields(40)P^k+1=Exk+1−x˜k+1−Kk+1zk+1−Gk+1x˜k+1xk+1−x˜k+1−Kk+1zk+1−Gk+1x˜k+1T.
Substituting (13) into (40) yields(41)P^k+1=Exk+1−x˜k+1−Kk+1Gk+1xk+1+εk+1−Gk+1x˜k+1xk+1−x˜k+1−Kk+1Gk+1xk+1+εk+1−Gk+1x˜k+1T=EI−Kk+1Gk+1xk+1−x˜k+1−Kk+1εk+1I−Kk+1Gk+1xk+1−x˜k+1−Kk+1εk+1T.
Substituting (15) and (17) into (41) yields(42)P^k+1=1−Kk+1Gk+1P˜k+11−Kk+1Gk+1T+Kk+1Rk+1Kk+1T−2EI−Kk+1Gk+1xk+1−x˜k+1εk+1TKk+1T.
From Assumption2, one can obtain(43)EI−Kk+1Gk+1xk+1−x˜k+1εk+1TKk+1T=I−Kk+1Gk+1Exk+1εk+1T−Ex˜k+1εk+1TKk+1T=I−Kk+1Gk+1Exk+1Eεk+1T−Ex˜k+1Eεk+1TKk+1T=0.
Substituting (43) into (42) yields(44)P^k+1=1−Kk+1Gk+1P˜k+11−Kk+1Gk+1T+Kk+1Rk+1Kk+1T.
Formula (26) is obtained.
Substituting (13) into (17) yields(45)Rk+1=Ezk+1−Gk+1xk+1zk+1−Gk+1xk+1T=Ezk+1−Gk+1x˜k+1−Gk+1xk+1−x˜k+1zk+1−Gk+1x˜k+1−Gk+1xk+1−x˜k+1T.
Substituting (15) into (45) yields(46)Rk+1=zk+1−Gk+1x˜k+1zk+1−Gk+1x˜k+1T+Gk+1P˜k+1Gk+1T−2EGk+1xk+1−x˜k+1+εk+1xk+1−x˜k+1TGk+1T=zk+1−Gk+1x˜k+1zk+1−Gk+1x˜k+1T−Gk+1P˜k+1Gk+1T−2Eεk+1xk+1−x˜k+1TGk+1T.
From Assumption2, one can obtain(47)Eεk+1xk+1−x˜k+1TGk+1T=Eεk+1xk+1TGk+1T−Eεk+1x˜k+1TGk+1T=0.
Substituting (15) into (45) yields(48)Rk+1=zk+1−Gk+1x˜k+1zk+1−Gk+1x˜k+1T−Gk+1P˜k+1Gk+1T.
Formula (24) is obtained.


## 2.1. Mathematical Model

Stochastic time-dependent road network can be represented by a directed networkGC,E,vijk. C represents a set of customer nodes, and E represents a set of road segments. V=N.vijk represents the travel speed from customer node i to j during time period k.The goal is to obtain the shortest expected travel time and service time.p0<p<1 denotes the customer demand probability. The total number of vehicles is K. tij denotes the travelling time from customer i to j. x¯ijl denotes whether vehicle l is responsible for the route from customer j to i. uil denotes whether customer i is serviced by vehicle l. τil denotes the service time of customer i by vehicle l. qi is the demand of the customer i.Wl is the total load of the vehicle l. The model can be described as follows:(1)minEmax1≤l≤K∑1≤i,j≤N,j≠itijlx¯ijl+∑1≤i≤Npτiluil,(2)s.t.∑1≤i≤N,i≠j∑1≤l≤Kxijl=1,1≤j≤N,(3)∑1≤j≤N,i≠j∑1≤l≤Kxijl=1,1≤i≤N,(4)∑1≤l≤Kuil=1,1≤i≤N,(5)uil≤∑1≤j≤Nx¯ijl=1,1≤i≤N,1≤l≤K,(6)∑1≤i≤Nqiuil≤W,1≤l≤Kl,(7)x¯ijl=1,lselects arci,j,1≤i,j≤N,1≤l≤K,0,else,(8)uil=1,iserviced byl,1≤i≤N,1≤l≤K,0,else.Hypothesis is stated as follows:(1)
Each customer is serviced exactly once. The loading and unloading services of each customer are completed only once. Equations (2) and (3) ensure that each customer is served exactly once.(2)
Each customer is loaded immediately after being visited by a vehicle. The total time includes travelling time and service time with no other times wasted.(3)
The total weight of goods demanded by customers cannot exceed the vehicle load servicing them. Equation (6) ensures the total weight of goods demanded by customers cannot exceed the vehicle load servicing them.(4)
Each customer is serviced by a vehicle only if he is visited on the route of the vehicle. Equation (5) ensures a customer is serviced by a vehicle only if he is visited on the route of the vehicle.(5)
Each customer can only be serviced by exactly one vehicle. Equation (4) ensures each customer is serviced exactly one vehicle.(6)
The positions of all customers are known.(7)
The position of the depot is known. The source of vehicles is determinant.(8)
The distances from the depot to each customer are known.(9)
The customer demand probability is a stochastic real number between zero and one.(10)
The travelling time between two customers is stochastic time-dependent due to the time-varying speed and real-time traffic conditions on roads.The existence of at least one solution of the above problem is proved in literature [12, 13]. A cutting plane algorithm for solving the above stochastic programming was proved to be finite [8]. The above stochastic objective function (1) was proved to be a convex function on the convex hull of the set of decision vectors satisfying constraints (2)–(8) [9]. Since the local minimum of convex optimization problem is the global minimum, the local minimum solution of the above problem will also be the global minimum.


## 2.2. Adaptive Fractional-Order Kalman Filter for Speed Prediction

The distance between customeri and j is represented by dij. Δt denotes the sampling time interval. The calculation procedure for travelling time tij is shown in Figure 1. vijk is updated and estimated by Kalman filter algorithm.Figure 1
Calculation procedure for travelling time.The fractional-order difference is formulated as(9)Δαxk=1hα∑j=0k−1jγjxk−j,(10)γj=αj=1,j=0,αα−1⋯α−j+1j!,j>0,where α is the fractional order. h is the sampling interval. k is the number of samples. The fractional-order model of speed prediction is given by(11)Δαxk+1=Akxk+wk,(12)xk+1=Δαxk+1−∑j=1k+1−1jγjxk+1−j,(13)zk=Gkxk+εk,where xk is the state variable of the traffic speed. zk is the measurement output. wk is the system state noise. εk is the measurement noise. x˜k is calculated as(14)x˜k+1=Exk+1|zk∗,where zk∗ is the measurement sequence containing measurement output z1,z2,⋯,zk. P˜k denotes the covariance of prediction error at time instant k, which is calculated as(15)P˜k+1=Exk+1−x˜k+1xk+1−x˜k+1T.Qk denotes the covariance of system noise at the time instant k, which is calculated as(16)Qk=Ewk−1wk−1T.Rk denotes the covariance of measurement noise at the time instant k, which is calculated as(17)Rk=Eεk−1εk−1T.x^k denotes the state estimation at time instant k, which is calculated as(18)x^k+1=Exk+1|zk+1∗.P^k is the covariance of estimation error at the time instant k, which is calculated as(19)P^k+1=Exk+1−x^k+1xk+1−x^k+1T.The covariance of state noise is designed as follows:(20)Qk=λEwk−1wk−1T.where λ>0 is a gain coefficient.When the covariance of state noise is small, the convergence speed is slow and the convergence accuracy is high. If the state changes too fast in a short time, the slow convergence rate will easily lead to tracking failure. When the covariance of state noise is large, the convergence speed is fast and the accuracy is low. Before reaching a stable state, choosing a larger covariance of state noise can accelerate the convergence speed. After reaching a stable state, the convergence accuracy can be improved by choosing smaller covariance of state noise.In standard Kalman filter,λ is a constant. However, it is difficult to determine an optimal gain in the whole process. When the covariance of state noise is small, the number of estimations will update slowly if λ is too small. When the error changes sharply, the estimation will oscillate and even the convergence process will be unstable if λ is too large. To accelerate the convergence process, an adaptive gain λ is used to make the estimation process adjusted reasonably and adapt to different data characteristics. λ increases when the covariance of state noise is small and decreases when the covariance of state noise is large. According to the noise covariance information, λ is calculated dynamically. Therefore, an adaptive mechanism for choosing λ is designed as follows:(21)λ=λmax1+P^k+12.The Kalman gainKk is calculated as(22)Kk+1=P˜k+1Gk+1TGk+1P˜k+1Gk+1T+Rk+1−1.The fractional Kalman filter is formulated as(23)x˜k+1=Akx^k−∑j=1k+1−1jγjx^k+1−j,(24)Rk+1=zk+1−Gk+1x˜k+1zk+1−Gk+1x˜k+1T−Gk+1P˜k+1Gk+1T,(25)x^k+1=x˜k+1+Kk+1zk+1−Gk+1x˜k+1,(26)P^k+1=1−Kk+1Gk+1P˜k+11−Kk+1Gk+1T+Kk+1Rk+1Kk+1T,(27)P˜k+1=Ak+γ1P^kAk+γ1T+Qkλ+∑j=2k+1γjP^k+1−jγjT.Assumption 1.
Ewk=0. wk is not related to zk.Assumption 2.
Eεk=0. εk is not related to xk and x˜k.Theorem 1.
For the fractional-order system with system state noise and measurement noise described by formulas (11)–(13), AFKF is given by formulas (23)–(27).Proof.
Substituting (11) and (12) into (14), one can obtain(28)x˜k+1=EAkxk+wk−∑j=1k+1−1jγjxk+1−j|zk∗=Akx^k+AkExk|zk∗−x^k+Ewk|zk∗−∑j=1k+1−1jγjExk+1−j|zk∗.
From Assumption1, one can obtain(29)Ewk|zk∗=∑i=1∞wkPw=wk|z=zk∗=∑i=1∞wkPw=wk,z=zk∗Pz=zk∗=∑i=1∞wkPw=wkPz=zk∗Pz=zk∗=∑i=1∞wkPw=wk=Ewk=0.
Substituting (29) into (28), one can obtain(30)x˜k+1=Akx^k+AkExk|zk∗−x^k−∑j=1k+1−1jγjExk+1−j|zk∗.
Substituting (18) into (29), one can obtain(31)x˜k+1=Akx^k+Akx^k−x^k−∑j=1k+1−1jγjExk+1−j|zk∗=x^k−∑j=1k+1−1jγjExk+1−j|zk∗=x^k−∑j=1k+1−1jγjExk+1−j|zk+1−j∗.
Substituting (18) into (31), one can obtain(32)x˜k+1=Akx^k−∑j=1k+1−1jγjx^k+1−j.
Formula (23) is obtained.
Substituting (11) and (12) into (15), one can obtain(33)P˜k+1=Exk+1−x˜k+1xk+1−x˜k+1T=EAkxk+wk−∑j=1k+1−1jγjxk+1−j−x˜k+1Akxk+wk−∑j=1k+1−1jγjxk+1−j−x˜k+1T.
Substituting (32) into (33), one can obtain(34)P˜k+1=EAkxk+wk−∑j=1k+1−1jγjxk+1−j−Akx^k+∑j=1k+1−1jγjx^k+1−jAkxk+wk−∑j=1k+1−1jγjxk+1−j−Akx^k+∑j=1k+1−1jγjx^k+1−jT=Ak+γ1Exk−x^kxk−x^kTAk+γ1T+∑j=2k+1γjExk+1−j−x^k+1−jxk+1−j−x^k+1−jTγjT+Ewk−1wk−1T.
Substituting (19) and (20) into (34), one can obtain(35)P˜k+1=Ak+γ1P^kAk+γ1T+Qkλ+∑j=2k+1γjP^k+1−jγjT.
Formula (27) is obtained.
The cost function is defined by(36)x^k+1=argminxk+1x˜k+1−xk+1P˜k+1−1x˜k+1−xk+1T+zk+1−Gk+1x˜k+1−Gk+1xk+1−x˜k+1Rk+1−1⋅zk+1−Gk+1x˜k+1−Gk+1xk+1−x˜k+1T.
Differentiating (36) with respect to xk+1 and making the result equal to zero, one can obtain(37)P˜k+1−1x^k+1−x˜k+1=Gk+1TRk+1−1zk+1−Gk+1x˜k+1−Gk+1x^k+1−x˜k+1.
From (37), one can obtain(38)x^k+1=Rk+1Gk+1−1P˜k+1−1+Gk+1−1⋅Rk+1Gk+1−1P˜k+1−1x˜k+1+zk+1−Gk+1x˜k+1+Gk+1x˜k+1=x˜k+1+Rk+1Gk+1−1P˜k+1−1+Gk+1−1zk+1−Gk+1x˜k+1=x˜k+1+P˜k+1Gk+1TGk+1P˜k+1Gk+1T+Rk+1−1zk+1−Gk+1x˜k+1.
Substituting (22) into (38) yields(39)x^k+1=x˜k+1+Kk+1zk+1−Gk+1x˜k+1.
Formula (25) is obtained.
Substituting (39) into (19) yields(40)P^k+1=Exk+1−x˜k+1−Kk+1zk+1−Gk+1x˜k+1xk+1−x˜k+1−Kk+1zk+1−Gk+1x˜k+1T.
Substituting (13) into (40) yields(41)P^k+1=Exk+1−x˜k+1−Kk+1Gk+1xk+1+εk+1−Gk+1x˜k+1xk+1−x˜k+1−Kk+1Gk+1xk+1+εk+1−Gk+1x˜k+1T=EI−Kk+1Gk+1xk+1−x˜k+1−Kk+1εk+1I−Kk+1Gk+1xk+1−x˜k+1−Kk+1εk+1T.
Substituting (15) and (17) into (41) yields(42)P^k+1=1−Kk+1Gk+1P˜k+11−Kk+1Gk+1T+Kk+1Rk+1Kk+1T−2EI−Kk+1Gk+1xk+1−x˜k+1εk+1TKk+1T.
From Assumption2, one can obtain(43)EI−Kk+1Gk+1xk+1−x˜k+1εk+1TKk+1T=I−Kk+1Gk+1Exk+1εk+1T−Ex˜k+1εk+1TKk+1T=I−Kk+1Gk+1Exk+1Eεk+1T−Ex˜k+1Eεk+1TKk+1T=0.
Substituting (43) into (42) yields(44)P^k+1=1−Kk+1Gk+1P˜k+11−Kk+1Gk+1T+Kk+1Rk+1Kk+1T.
Formula (26) is obtained.
Substituting (13) into (17) yields(45)Rk+1=Ezk+1−Gk+1xk+1zk+1−Gk+1xk+1T=Ezk+1−Gk+1x˜k+1−Gk+1xk+1−x˜k+1zk+1−Gk+1x˜k+1−Gk+1xk+1−x˜k+1T.
Substituting (15) into (45) yields(46)Rk+1=zk+1−Gk+1x˜k+1zk+1−Gk+1x˜k+1T+Gk+1P˜k+1Gk+1T−2EGk+1xk+1−x˜k+1+εk+1xk+1−x˜k+1TGk+1T=zk+1−Gk+1x˜k+1zk+1−Gk+1x˜k+1T−Gk+1P˜k+1Gk+1T−2Eεk+1xk+1−x˜k+1TGk+1T.
From Assumption2, one can obtain(47)Eεk+1xk+1−x˜k+1TGk+1T=Eεk+1xk+1TGk+1T−Eεk+1x˜k+1TGk+1T=0.
Substituting (15) into (45) yields(48)Rk+1=zk+1−Gk+1x˜k+1zk+1−Gk+1x˜k+1T−Gk+1P˜k+1Gk+1T.
Formula (24) is obtained.


## 3. Cultural Algorithm-Based Cuckoo Search

### 3.1. Encoding Method

Construct the mapping relationship from the individual to the solution decision variabley. y is composed of randomly generated (N + K − 1)-dimensional real numbers. Take Table 1 for example. Suppose N is 4 and K is 3. N + K − 1 = 6.Table 1
Mapping from individual to decision variable.
jyjπj1−0.99221.80533.0144−0.7215−1.20662.163y4 is the largest dimension, so π4=1, y1 is the second largest dimension, so π1=2, and so on. Finally, the whole mapped individual is π=2,5,4,1,6,3.Then,πjN+1≤πj≤N+K−1 are regarded as the boundary lines of different vehicles. Customers within the same boundary lines are assigned to the same vehicle in turn. So π2=5 and π5=6 are regarded as boundary lines. So the mapped solution decision variable can be written as πit=2|4,1|3. This means that the second customer is assigned to the first vehicle, the fourth and the first customers are assigned to the second vehicle, and the third customer is assigned to the third vehicle.


### 3.2. Cuckoo Search

In CS, the position value of each cuckoo nest corresponds to each solution to VRP. Its goal is to search for the best cuckoo nest.Lévy flight is a random walk strategy, whose step length can be determined by the Lévy distribution. When the newyk+1 is discovered to be solution-suitable, a Lévy flight can be used as follows:(49)yik+1=yik+σ⊕levyβ,where i represents the ith cuckoo. σ represents the step size related to the scales of the problem. β is the index coefficient. An adaptive method can also be used:(50)σ=σ0yjt−yit,where σ0 is a constant. The product ⊕ means entry-wise multiplications. The step length of a random walk by the Lévy flight is computed from a Lévy distribution:(51)levyβ=t−1−β,0<β≤2.The steps essentially form a random walk process with a power-law step-length distribution with a heavy tail. Some new solutions should be generated by Lévy walk around the best solution obtained so far to speed up the local search. The locations of a new solution may be far enough from the current best solution; this will make sure the system will not be trapped in a local optimum. When the length of random step is determined by the maximum step length in the Lévy distribution, the Lévy flight can effectively provide a random walk model.To avoid parameters dependence, an adaptive mechanism is introduced. The step size of Lévy flight and discovery probability are computed as follows:(52)σi=σmin+σmax−σminfmax−fifmax−fmin,(53)pa=pa,min+pa,max−pa,min1−itIT,where σmax denotes the maximum step size of Lévy flight and σmin denotes the minimum step size of Lévy flight. fmax denotes the maximum fitness of all nests, and fmin denotes the minimum fitness of all nests. pa,max denotes the maximum discovery probability, and pa,min denotes the minimum discovery probability. it denotes the current iteration, and IT denotes the maximum iterations number. As the fitness of a nest decreases, its step size of Lévy flight will increase. As the time of iteration increases, the discovery probability of nests will decrease.


### 3.3. Cultural Algorithm

CA is a dual evolution mechanism proposed by Reynolds in 1994. For VRP, situation knowledge represents the optimization goal, and normative knowledge represents constraint conditions. The individuals gained in the evolution process are passed to Belief space by Accept function. Then the Belief space is updated through Update function. The Influence function uses the useful information of Belief space to guide the evolution process of population space. The main structure of the cultural algorithm is shown in Figure2.Figure 2
Structure of the cultural algorithm.In the Belief space,u denotes the upper bound of the individual, and l denotes the lower bound of the individual. Lj is lower limit of the fitness value of lj, and Uj is upper limit of the fitness value of uj. The updating rules of situation knowledge are as follows:(54)Sk+1=ybestk,fybestk<fSk,St,else,where Sk denotes the optimal individual in k generation and f denotes the fitness function. The updating rules of standard knowledge are(55)ljt+1=yi,jt,yi,jt≤ljtorfyit<Ljt,ljt,else,(56)Ljt+1=fyit,yi,jt≤ljtorfyit<Ljt,Ljt,else,(57)ujt+1=yk,jt,yk,jt≥ujtorfykt<Ujt,ujt,else,(58)Ujt+1=fykt,yk,jt≥ujtorfykt<Ujt,Ujt,else.


### 3.4. CACS Process

To accelerate the convergence speed of CS, this paper adopts the double mechanism of the CA. The individuals with high fitness form the culture in the Belief space through acceptance rules. The representative individuals further influence the evolution of next population through the outstanding individuals adjusted. The process of CACS with vehicle speed estimation is shown in Figure3.Figure 3
Pseudocode of CACS.In this paper, the adaptive factor strategy is introduced. Through the dynamic adjustment of the factor to the optimization step size, the algorithm can search with adaptive step size in the early stage of evolution and improve the convergence speed of the algorithm. In the late stage of evolution, it can quickly tend to the optimal solution and further improve the convergence speed of the algorithm. This reduces the computational cost of the algorithm.


### 3.5. Convergence of CACS

The solution generated in thekth iteration is denoted as yk. The optimal group state set is denoted as H. The optimal state set is denoted as R. All the states of all the cuckoos form a state space for the group denoted as Q.Lemma 1.
(see [26]). Assume that the objective function f is a measurable function. The feasible solution space A is a measurable subset on ℜn. Yε is a set of global optimal solutions. Pyk∈Yε denotes the probability that the solution yk belongs to the optimal solutions set Yε. For any Borei subset B, its Lebesgue measure νB>0. μkB is the probability of obtaining B by measuring μk.∏k=0∞1−μkB=0. Thus,(59)limk⟶+∞Pyk∈Rε=1.Lemma 2.
The finite homogeneous Markov chain starting from any nonrecurrent states will reach its recurrent state with probability of 1.Lemma 3.
Suppose Markov chain has a nonempty setC and there is not any nonempty closed set D satisfying C∩D=ϕ. Then when j∉C, limn⟶+∞Pxn=j=0.Theorem 2.
The Markov chain represented by CS is finite and homogeneous.Proof.
The search space of the optimization algorithm is limited, and the cuckoo nest state is limited. So the state space of nest position is limited. Since the population state of nest positions is composed ofm nest positions and m is a finite positive integer, the population state of nest position is also limited.
From the process of CS, the transfer probabilityPTyk−1=yk of any nest position is only related to the state of the previous state yk−1. T⋅ denotes the transfer operation. Therefore, the population states of nest position have Markov property.
The transfer probabilityPTyk−1=yk of any nest position is not related to the time instant k. That is to say, the population states sequence of nest location is homogeneous. Therefore, the population states sequence of cuckoo nest location is a finite homogeneous Markov chain.Theorem 3.
When the number of iterations approaches becoming sufficiently large, the group state sequence will converge to the optimal state solution setH.Proof.
∀yi∈R,∀yj∉R, fyi<fyj. From the process of CS, it is obtained that PTyj=yi=0. For the state sequence {y (k); k > 0}, the state set R of optimal nest position is a closed set on the state space Y.
SinceR is a closed set on Y, according to the definition of closed set, the state set H of the optimal nest position group is a closed set on the group state space Q.
Assume that there exists a closed setB in the group space Q, such that B∩H=ϕ. ∀yi∈R,∀yj∉R, fyi<fyj. Then, ∀yi∈B∩H⊆H, fyj≤fyi. Thus, it is obtained that PTyj=yi≠0. So B is not a closed set, which contradicts the hypothesis. Therefore, there is no nonempty closed set outside H in the nest position state space Q such that B∩H=ϕ.
From Lemma3, when the number of iterations approaches becoming sufficiently large, the group state sequence will converge to the optimal state solution set H.Each iteration of CS algorithm preserves the optimal position of the population and guarantees the nonincremental fitness. According to Theorem3, the population sequence of cuckoo nest position will enter the optimal state after continuous and infinite iterations. Therefore, the possibility that the global optimal solution cannot be searched continuously and infinitely is 0. So the CS algorithm converges to the global optimum.Theorem 4.
The random variation of the state in Belief space of CACS belongs to finite homogeneous Markov chain.Proof.
In CACS, the Belief space is used to record the evolutionary knowledge in the search process, and CACS searches for the optimal value in discrete and limited space. So the Belief space used for recording the evolutionary knowledge is also limited. Meanwhile, the state transition probability in the Belief space is only related to its current state and has nothing to do with time instantk. Therefore, the random variation of the state in the Belief space of CACS belongs to the finite homogeneous Markov chain.Theorem 5.
CACS has global convergence with probability of 1.Proof.
Denote the global knowledge state space in the Belief space asΩg. Denote the local knowledge state space as Ωl. The whole state space of Belief space is Π=Ig,Is, in which Ig denotes the global knowledge set and Is denotes the local set of population space.
Π is divided into two subspaces, namely, recurrent state space Π1 and nonrecurrent state space Π2. The recurrent state space can be expressed as follows:(60)Π1=Ig,Is|Ig∈Ωg,Is∈Ωl,and thus Π2∩Π2=ϕ and Π2∪Π2⊆Π.
If there is a knowledge statey∈Π1, according to the definition of global knowledge, it will not transfer to Π2. So it is a closed set. The state in Π1 is interlinked. So the state in Π1 is a recurrent state. If y∈Π2, each individual will move throughout the whole solution space due to evolution. So the convergence probability of the evolutionary population to the global optimal solution is greater than 0. That means the transition probability of y from Π2 to Π1 is greater than 0. Therefore, the state in Π2 is a nonrecurrent state. According to Lemma 2, the evolution of knowledge in Belief space from any nonrecurrent state will always transfer to the recurrent state with probability of 1, that is, converge to the global optimal solution according to probability of 1.


## 3.1. Encoding Method

Construct the mapping relationship from the individual to the solution decision variabley. y is composed of randomly generated (N + K − 1)-dimensional real numbers. Take Table 1 for example. Suppose N is 4 and K is 3. N + K − 1 = 6.Table 1
Mapping from individual to decision variable.
jyjπj1−0.99221.80533.0144−0.7215−1.20662.163y4 is the largest dimension, so π4=1, y1 is the second largest dimension, so π1=2, and so on. Finally, the whole mapped individual is π=2,5,4,1,6,3.Then,πjN+1≤πj≤N+K−1 are regarded as the boundary lines of different vehicles. Customers within the same boundary lines are assigned to the same vehicle in turn. So π2=5 and π5=6 are regarded as boundary lines. So the mapped solution decision variable can be written as πit=2|4,1|3. This means that the second customer is assigned to the first vehicle, the fourth and the first customers are assigned to the second vehicle, and the third customer is assigned to the third vehicle.


## 3.2. Cuckoo Search

In CS, the position value of each cuckoo nest corresponds to each solution to VRP. Its goal is to search for the best cuckoo nest.Lévy flight is a random walk strategy, whose step length can be determined by the Lévy distribution. When the newyk+1 is discovered to be solution-suitable, a Lévy flight can be used as follows:(49)yik+1=yik+σ⊕levyβ,where i represents the ith cuckoo. σ represents the step size related to the scales of the problem. β is the index coefficient. An adaptive method can also be used:(50)σ=σ0yjt−yit,where σ0 is a constant. The product ⊕ means entry-wise multiplications. The step length of a random walk by the Lévy flight is computed from a Lévy distribution:(51)levyβ=t−1−β,0<β≤2.The steps essentially form a random walk process with a power-law step-length distribution with a heavy tail. Some new solutions should be generated by Lévy walk around the best solution obtained so far to speed up the local search. The locations of a new solution may be far enough from the current best solution; this will make sure the system will not be trapped in a local optimum. When the length of random step is determined by the maximum step length in the Lévy distribution, the Lévy flight can effectively provide a random walk model.To avoid parameters dependence, an adaptive mechanism is introduced. The step size of Lévy flight and discovery probability are computed as follows:(52)σi=σmin+σmax−σminfmax−fifmax−fmin,(53)pa=pa,min+pa,max−pa,min1−itIT,where σmax denotes the maximum step size of Lévy flight and σmin denotes the minimum step size of Lévy flight. fmax denotes the maximum fitness of all nests, and fmin denotes the minimum fitness of all nests. pa,max denotes the maximum discovery probability, and pa,min denotes the minimum discovery probability. it denotes the current iteration, and IT denotes the maximum iterations number. As the fitness of a nest decreases, its step size of Lévy flight will increase. As the time of iteration increases, the discovery probability of nests will decrease.


## 3.3. Cultural Algorithm

CA is a dual evolution mechanism proposed by Reynolds in 1994. For VRP, situation knowledge represents the optimization goal, and normative knowledge represents constraint conditions. The individuals gained in the evolution process are passed to Belief space by Accept function. Then the Belief space is updated through Update function. The Influence function uses the useful information of Belief space to guide the evolution process of population space. The main structure of the cultural algorithm is shown in Figure2.Figure 2
Structure of the cultural algorithm.In the Belief space,u denotes the upper bound of the individual, and l denotes the lower bound of the individual. Lj is lower limit of the fitness value of lj, and Uj is upper limit of the fitness value of uj. The updating rules of situation knowledge are as follows:(54)Sk+1=ybestk,fybestk<fSk,St,else,where Sk denotes the optimal individual in k generation and f denotes the fitness function. The updating rules of standard knowledge are(55)ljt+1=yi,jt,yi,jt≤ljtorfyit<Ljt,ljt,else,(56)Ljt+1=fyit,yi,jt≤ljtorfyit<Ljt,Ljt,else,(57)ujt+1=yk,jt,yk,jt≥ujtorfykt<Ujt,ujt,else,(58)Ujt+1=fykt,yk,jt≥ujtorfykt<Ujt,Ujt,else.


## 3.4. CACS Process

To accelerate the convergence speed of CS, this paper adopts the double mechanism of the CA. The individuals with high fitness form the culture in the Belief space through acceptance rules. The representative individuals further influence the evolution of next population through the outstanding individuals adjusted. The process of CACS with vehicle speed estimation is shown in Figure3.Figure 3
Pseudocode of CACS.In this paper, the adaptive factor strategy is introduced. Through the dynamic adjustment of the factor to the optimization step size, the algorithm can search with adaptive step size in the early stage of evolution and improve the convergence speed of the algorithm. In the late stage of evolution, it can quickly tend to the optimal solution and further improve the convergence speed of the algorithm. This reduces the computational cost of the algorithm.


## 3.5. Convergence of CACS

The solution generated in thekth iteration is denoted as yk. The optimal group state set is denoted as H. The optimal state set is denoted as R. All the states of all the cuckoos form a state space for the group denoted as Q.Lemma 1.
(see [26]). Assume that the objective function f is a measurable function. The feasible solution space A is a measurable subset on ℜn. Yε is a set of global optimal solutions. Pyk∈Yε denotes the probability that the solution yk belongs to the optimal solutions set Yε. For any Borei subset B, its Lebesgue measure νB>0. μkB is the probability of obtaining B by measuring μk.∏k=0∞1−μkB=0. Thus,(59)limk⟶+∞Pyk∈Rε=1.Lemma 2.
The finite homogeneous Markov chain starting from any nonrecurrent states will reach its recurrent state with probability of 1.Lemma 3.
Suppose Markov chain has a nonempty setC and there is not any nonempty closed set D satisfying C∩D=ϕ. Then when j∉C, limn⟶+∞Pxn=j=0.Theorem 2.
The Markov chain represented by CS is finite and homogeneous.Proof.
The search space of the optimization algorithm is limited, and the cuckoo nest state is limited. So the state space of nest position is limited. Since the population state of nest positions is composed ofm nest positions and m is a finite positive integer, the population state of nest position is also limited.
From the process of CS, the transfer probabilityPTyk−1=yk of any nest position is only related to the state of the previous state yk−1. T⋅ denotes the transfer operation. Therefore, the population states of nest position have Markov property.
The transfer probabilityPTyk−1=yk of any nest position is not related to the time instant k. That is to say, the population states sequence of nest location is homogeneous. Therefore, the population states sequence of cuckoo nest location is a finite homogeneous Markov chain.Theorem 3.
When the number of iterations approaches becoming sufficiently large, the group state sequence will converge to the optimal state solution setH.Proof.
∀yi∈R,∀yj∉R, fyi<fyj. From the process of CS, it is obtained that PTyj=yi=0. For the state sequence {y (k); k > 0}, the state set R of optimal nest position is a closed set on the state space Y.
SinceR is a closed set on Y, according to the definition of closed set, the state set H of the optimal nest position group is a closed set on the group state space Q.
Assume that there exists a closed setB in the group space Q, such that B∩H=ϕ. ∀yi∈R,∀yj∉R, fyi<fyj. Then, ∀yi∈B∩H⊆H, fyj≤fyi. Thus, it is obtained that PTyj=yi≠0. So B is not a closed set, which contradicts the hypothesis. Therefore, there is no nonempty closed set outside H in the nest position state space Q such that B∩H=ϕ.
From Lemma3, when the number of iterations approaches becoming sufficiently large, the group state sequence will converge to the optimal state solution set H.Each iteration of CS algorithm preserves the optimal position of the population and guarantees the nonincremental fitness. According to Theorem3, the population sequence of cuckoo nest position will enter the optimal state after continuous and infinite iterations. Therefore, the possibility that the global optimal solution cannot be searched continuously and infinitely is 0. So the CS algorithm converges to the global optimum.Theorem 4.
The random variation of the state in Belief space of CACS belongs to finite homogeneous Markov chain.Proof.
In CACS, the Belief space is used to record the evolutionary knowledge in the search process, and CACS searches for the optimal value in discrete and limited space. So the Belief space used for recording the evolutionary knowledge is also limited. Meanwhile, the state transition probability in the Belief space is only related to its current state and has nothing to do with time instantk. Therefore, the random variation of the state in the Belief space of CACS belongs to the finite homogeneous Markov chain.Theorem 5.
CACS has global convergence with probability of 1.Proof.
Denote the global knowledge state space in the Belief space asΩg. Denote the local knowledge state space as Ωl. The whole state space of Belief space is Π=Ig,Is, in which Ig denotes the global knowledge set and Is denotes the local set of population space.
Π is divided into two subspaces, namely, recurrent state space Π1 and nonrecurrent state space Π2. The recurrent state space can be expressed as follows:(60)Π1=Ig,Is|Ig∈Ωg,Is∈Ωl,and thus Π2∩Π2=ϕ and Π2∪Π2⊆Π.
If there is a knowledge statey∈Π1, according to the definition of global knowledge, it will not transfer to Π2. So it is a closed set. The state in Π1 is interlinked. So the state in Π1 is a recurrent state. If y∈Π2, each individual will move throughout the whole solution space due to evolution. So the convergence probability of the evolutionary population to the global optimal solution is greater than 0. That means the transition probability of y from Π2 to Π1 is greater than 0. Therefore, the state in Π2 is a nonrecurrent state. According to Lemma 2, the evolution of knowledge in Belief space from any nonrecurrent state will always transfer to the recurrent state with probability of 1, that is, converge to the global optimal solution according to probability of 1.


## 4. Example and Analysis

### 4.1. Experimental Results with Deterministic Demand Probability

The CACS algorithm is tested on Intel® Core™ i3-4150T CPU @3.00 GHz, 64-bit operating system with memory of 4.00 GB, and x64-based processor. The VRP library used is published on the Networking and Emerging Optimization homepage [27]. The capacitated VRP instance library was selected. The first layer index is the first letter of word in the instance library. The second layer index is the number of customers, and the third layer index is the number of vehicles. Take A, B, P, and E sets as examples. Set A has 27 VRP instances from A-n32-k5 to A-n80-k10 proposed by Augerat et al. Set B has 23 VRP instances from B-n31-k5 to B-n78-k10 proposed by Augerat et al. Set P has 24 VRP instances from P-n16-k8 to P-n101-k4 proposed by Augerat et al. Set E has 15 VRP instances from E-n13-k4 to E-n101-k14 proposed by Augerat et al.Consider the case when the arrival probability is 1. Table2 gives the shortest path length, average iteration times, and average CPU computing time of each instance.Table 2
Calculation result of different instances when the demand probability is 1.
InstanceA-n32-k5A-n39-k6B-n31-k5P-n16-k8P-n19-k2N3239311619K56582Optimal solution (m)784.612831.2696662.3275451.9471212.6569Average iteration number427648329193308Average CPU computing time (s)18.968827.15636.81253.29696.7188Published optimal solution (m)784831672450212Improvement ratio (%)001.4393600For instance, B-n31-k5 algorithm has found a better solution than the published optimal solution, whose scheme is shown in Table3.Table 3
Comparison of new solution and published optimal solution for B-n31-k5.
RoutePublished optimal solutionNew solutionTravelling sequenceDistanceCapacityTravelling sequenceDistanceCapacity130 ⟶ 23 ⟶ 8 ⟶ 12 ⟶ 28 ⟶ 26109.6811977 ⟶ 30 ⟶ 5 ⟶ 25 ⟶ 18 ⟶ 16 ⟶ 21119.237289397221 ⟶ 16 ⟶ 18 ⟶ 25 ⟶ 5 ⟶ 4 ⟶ 29112.99158612 ⟶ 24 ⟶ 19 ⟶ 1 ⟶ 3 ⟶ 6 ⟶ 9102.36598029937 ⟶ 17 ⟶ 13 ⟶ 6 ⟶ 9 ⟶ 22110.82879620 ⟶ 27 ⟶ 10 ⟶ 2 ⟶ 14 ⟶ 15 ⟶ 11 ⟶ 28 ⟶ 8240.508851895420 ⟶ 27 ⟶ 10 ⟶ 2183.5752384 ⟶ 29 ⟶ 22 ⟶ 13 ⟶ 17 ⟶ 23108.715642299514 ⟶ 15 ⟶ 11 ⟶ 24 ⟶ 19 ⟶ 1 ⟶ 3155.4298952698.0612053822Total672.506412668.889412Figure4 shows the shortest path of 70 nodes when demand probability is 1. The horizontal axis represents the horizontal axis of a city map with the unit of meter. The vertical axis represents the longitudinal axis of a city map with the unit of meter.Figure 4
Optimal path of 70 nodes when demand probability is 1.Figure5 shows the objective function convergence curves of 70 nodes when demand probability is 1. The horizontal axis represents iteration times. The vertical axis represents the best individual in each iteration.Figure 5
Objective function convergence curves of 70 nodes when demand probability is 1.


### 4.2. Parameter Influence Analysis

Figure6 shows the average results of E-n22-k4, B-n31-k5, and B-n64-k9 after 20 times calculation under different demand probability p with IT 100. The updating rate of situational knowledge α is 0.3. The updating rate of standard knowledge β is 0.5. The population size m is 50. τil is 1.Figure 6
Calculation results of differentp with IT 100. (a) E-n22-k4, (b) B-n31-k5, and (c) B-n64-k9.
(a)(b)(c)Figure7 shows average results of 20 times calculation under different α with IT 100. p is 0.45. β is 0.5 and m is 50.Figure 7
Calculation results of differentα with IT 100. (a) E-n22-k4, (b) B-n31-k5, and (c) B-n64-k9.
(a)(b)(c)Figure8 shows average results of 20 times calculation under different β with IT 100. p is 0.45. α is 0.45 and m is 50.Figure 8
Calculation results of differentβ with IT 100. (a) E-n22-k4, (b) B-n31-k5, and (c) B-n64-k9.
(a)(b)(c)Figure9 shows average results of 20 times calculation under different m with IT 100. p is 0.45. α is 0.3 and β is 0.5.Figure 9
Calculation results of different m with IT 100. (a) E-n22-k4, (b) B-n31-k5, and (c) B-n64-k9.
(a)(b)(c)From the above figures, whenm is larger, the optimal solution is better. Increasing m is advantageous to searching for better solutions. However, when m reaches a certain value, the optimal solution is almost unchanged as m increases. In this instances, the optimal parameters are p=0.45, α=0.3, β=0.5, and m=50.The sensitivity analysis of iterations times is shown in Table4. In the case of E-n22-k4, when IT is too small, the calculation results still need to be improved and Δf/ΔNmaxis large. When IT reaches 130, Δf/ΔNmax remains close to zero as IT increases. The results are almost stable.Table 4
Calculation results of E-n22-k4 with different iteration times.
Nmax100110120130140Mean379.3766377.4962375.8675375.1023375.1023Std.0.76250.21660.10260.01250.0063ΔNmax1010101010ΔMean1.88041.62870.765200ΔMean/ΔNmax0.188040.112870.0765200We have also tried to vary the number of host nests and the discovery probabilitypa in Cuckoo Search algorithm. We have used n = 5, 10, 15, 20, 50, 100, 150, 250, and 500 and pa = 0.01, 0.05, 0.1, 0.15, 0.2, 0.25, 0.4, and 0.5. We found that n = 15 and pa = 0.25 are sufficient for most optimization problems. Results and analysis also imply that the convergence rate of Cuckoo Search algorithm is not sensitive to the parameters used.Apart from the number of host nests, there is one parameterpa. This also means that we do not have to fine-tune these parameters for a specific problem. Subsequently, Cuckoo Search algorithm is more generic and robust for many optimization problems. There are fewer parameters to be fine-tuned in Cuckoo Search than in other metaheuristic algorithms.


### 4.3. VRPSC Experiments with Fixed Demand Probability

p is set as different real numbers from 0.1 to 0.9. The parameters are set as follows: α=0.3, β=0.5, m=50, and τil=1. The initial vehicle speeds v are set as stochastic real numbers within the range of [1, 6]. Table 5 gives the average results of shortest expected total travelling time and servicing time for P-n16-k8 case with different IT and different p under 20 times calculation.Table 5
Results with different IT under differentp.
p1002003004005000.147.01664747.01664747.01664747.01664747.0166470.265.35669665.35668065.35668065.35668065.3566810.374.81729074.81697874.81697574.81697774.8169840.481.35662781.35402281.35401081.35401781.3540640.586.76676686.75332886.75322086.75316486.7533280.691.63513591.58497891.58353291.58332291.5833180.796.18311196.04008496.02895496.02695096.0269810.8100.469282100.163211100.110735100.106147100.1049060.9104.306277103.887785103.744337103.739400103.737179Table5 shows that, with the increasing of p, the objective function also increases. With the increasing of IT, the objective function reduces. However, when IT reaches 300, even if IT continues to increase, the shortest expected total travelling time and servicing time will not be reduced and will remain almost the same.


### 4.4. VRPSC Experiments with Distinguishing Customer Characteristics

Customers can be divided into large customers, small and medium customers, and retail customers according to the weight of the goods they order. Consider the influence of customer property in these three cases. In the case of B-n31-k5, the customers who order goods weighing in the range of [19, 25] belong to large customers. The customers who order goods weighing in the range of [9, 18] belong to small and medium customers. The customers who order goods weighing in the range of [2, 8] belong to small and retail customers. The parameters are set as follows:p=0.45, α=0.3, β=0.5, m=50, and τil=1. The initial vehicle speeds v are set as stochastic real numbers within the range of [1, 6].Table6 shows the optimal solution, average solution, and standard deviation in three cases of different customer types after 20 times calculation, when the demand probability is 0.8 or 0.9. In Table 6, 1st represents the result for the first calculation, 2nd denotes the results for the second calculation, and so on.Table 6
Objective function with different customer types and different demand probability.
Probability0.90.8CustomerLargeSmallRetailLargeSmallRetail1st596.155638.2462656.2919585.8801590.4435638.04122nd594.4813637.9749651.5232585.6425594.5347636.50343rd590.7583635.9657653.6638585.8801594.5347637.89914th596.155638.7321659.7856587.4265594.5347636.50345th590.7583639.7165653.4783582.6757590.4435634.85296th590.7583637.8215655.5595588.2996590.4435636.50347th590.7583637.9749657.6003582.6757590.4435634.85298th596.155637.8215656.7167588.2996594.5347636.86929th590.7583637.9749651.375582.6757590.4435636.503410th596.155638.6415654.8217588.2996594.5347636.5034Ave.593.2893638.087655.0816585.7755592.4891636.5032Std.2.7141530.9476772.6643242.3735792.1562521.04854Table6 shows when the large customers have least objective function. The small and medium customers have second shortest objective function. The retail customers have the longest objective function. This is because when the vehicle does not need to serve large customers, it can freely and flexibly arrange its routes. When vehicles need to serve large customers, with the limit of maximum carrying capacity, the choice of driving routes has more constraints. With the increase of p, the objective function also increases. The objective function with p of 0.9 is greater than that with p of 0.8.


### 4.5. Comparative Analysis

Table6 gives performance of CACS, CA, and CS with different p. Since CACS uses CA and CS, it is more complex than these two algorithms, while its time cost and storage cost are slightly higher compared to the latter two. However, the rapid development of computer hardware technology makes this difference for solving problems of general scale not apparent. From the experimental results, CACS performs better than CA, CS, and ACO.To avoid sensitivity to parameters selection, which is one of the main limitations in swarm intelligence algorithms, adaptive mechanism is adopted into CACS. ACO [28] used for comparison also proposed machine learning strategies to control the parameter adaptation. The objective function by different p and different algorithms is listed in Table 7.Table 7
Objective function with differentp and different algorithms.
pCACSCACSACOResultsRate (%)ResultsRate (%)ResultsRate (%)0.147.01664747.1056540.1889547.1341360.2492747.0910610.158020.265.35668165.3745870.0273965.3943740.0576465.3647340.012320.374.81698474.8348670.0239074.9212100.1391174.9014230.112730.481.35406481.5048640.1850281.4314570.0950481.4017630.058600.586.75332886.8912110.1586886.7931340.0458686.8014750.055470.691.58331891.7163010.1449991.6137060.0331791.6546010.077770.796.02698196.1053460.0815496.0713410.0461796.0578200.032100.8100.104906100.1159630.01104100.1100360.00512100.1134670.008550.9103.737179103.9213470.17722103.8641360.12223103.9163160.12223


### 4.6. TDVRPSC Using AFKF Speed Prediction

Based on the technology of AFKF, the highway data analysis platform is established. After extracting the information and constructing appropriate traffic prediction model, the future vehicle speed in different business districts can be effectively estimated. When the speed of a road is detected to be below the threshold, it can be inferred that this road is congested due to accidents or other reasons. Rapid identification of congested roads is completed. This is also useful for emergencies. When traffic congestion is detected in a certain road, decision support and vehicle dispatching optimization can adjust VRP scheme in time and avoid congested roads, so as to reduce transportation time and cost. In the field of real-time traffic forecasting, the fast information processing and mining ability of large data have very high reliability for real-time VRP scheduling. The accuracy of VRP scheduling is improved for the time-varying actual urban traffic road.Table8 lists the traffic speed and congestion delay index in different business districts of Xiamen at 16:25 on February 17, 2019 [29]. Congestion delay index is used to evaluate the degree of urban congestion. It is the ratio of the actual travel time to the free travel time of a trip.Table 8
Traffic speed and congestion delay index in Xiamen (2019.2.17, 16:25).
RankingBusiness districtCongestion delay indexSpeed (km/h)1Siming1.5827.212Huli1.3336.463Haicang1.2544.824Jimei1.2441.805Tongan1.1347.74Table9 lists the traffic speed in different business districts of Xiamen from 16:30 to 20:00 on February 17, 2019. The unit of speed is km/h.Table 9
Traffic speed in Xiamen during 2019.2.17, 16:30–20:00, in km/h.
TimeSimingHuliHaicangJimeiTongan16:3027.2536.8845.1342.0346.8916:4027.3336.3842.0146.1946.7816:5027.8536.2545.3942.3446.0517:0028.0736.5345.7242.5245.5517:1028.8336.6942.3845.9945.8917:2028.4636.4641.7045.7844.9017:3028.2736.7140.9342.9245.6417:4027.5736.6340.9242.8244.8317:5026.8841.8136.3341.1843.0118:0026.4140.5636.0942.4641.5818:1026.6840.1335.6842.1441.3818:2026.8140.2535.4140.8444.5618:3026.5935.3741.6140.7645.9518:4027.0441.4336.1441.1346.1818:5027.5436.5243.2941.6646.6519:0028.4337.0645.6042.4747.0719:1029.0737.0746.2343.4747.9219:2029.9437.7846.9343.547.9319:3030.2138.0646.6743.8348.1819:4031.0438.2147.0144.2248.7819:5031.2238.4647.0243.9548.5020:0031.6938.9447.4044.0249.10Figure10 shows Xiamen city map at 16:25 on February 17, 2019. Different business districts are marked with their congestion delay indexes.Figure 10
Xiamen map for real-time experiment (2019.2.17, 16:25).Several branches of Yuantong Express in Xiamen city distributed in different business districts are used. Table9 lists their business district, latitude, longitude, and the Mercator coordination converted. The No. 0 branch is the head office. VRP task is to distribute goods to all other branches.


### 4.7. Accuracy of Prediction Model

Figure11 shows the traffic speed prediction in Huli district using AFKF prediction from 16:30 to 20:00 on February 17, 2019. The horizontal axis represents time. The vertical axis represents the traffic speed in km/h. The star marking connected with solid wire represents the real speed. The joint spider connected with dotted lines represents the predicted speed.Figure 11
Traffic speed prediction in Huli using adaptive fractional Kalman prediction.Table10 shows the prediction mean square error (MSE) and mean absolute percentage error (MAPE) by AFKF, fractional Kalman filter (FKF), and Kalman filter (KF). Table 11 shows AFKF has smaller prediction error and better prediction accuracy than FKF and KF.Table 10
Prediction error with different algorithms (km/h).
Business districtAFKFFKFKFMSEMAPEMSEMAPEMSEMAPESiming1.44433.87731.98884.59952.62755.2535Huli2.76323.49292.78233.61893.10273.6632Haicang2.29123.42452.33343.43222.56533.4551Jimei1.75562.29081.76632.49312.66642.736Tongan1.49512.07321.94252.21981.92842.4390Table 11
Geographical location of Yuantong Express.
No.AreaLongitude (degree)Latitude (degree)Vertical coordinate (m)Horizontal ordinate (m)0Siming118.07654824.4496640061477627399761Siming118.08607424.4839640061630527456702Siming118.13810824.4774680062529827449303Siming118.1701624.4824230062981827452494Huli118.14282224.5146570062543727515055Huli118.17888724.5041870063222127495196Huli118.14729724.4889770062673727472517Huli118.14467624.5219740062595927524768Jimei118.10840724.6100440062012227684779Jimei118.1046424.60853600619048276918810Tongan118.15615524.72603900627708279065711Tongan118.14663524.69745200626189278561412Tongan118.13816224.70864600624935278773813Haicang118.01761924.53105700604864275386914Haicang118.04275824.49852800608609274878915Haicang118.01589924.530205006043892753588Figure12 shows that the prediction errors vary with different fractional orders of Kalman filter. The horizontal axis represents fractional order. The vertical axis represents the prediction error of traffic speed in km/h. The five-pointed star represents the mean relative error. The asterisk represents the mean square error. For certain business districts, there exists a most suitable fractional order for reducing its prediction error. Compared with integer-order filters, fractional-order Kalman filters have wider parameter selection range and higher prediction accuracy.Figure 12
Prediction error of traffic speed with different fractional orders.If the prediction model is not accurate, wrong route selection may be made. For example, the traffic flow is busy in Siming district and its vehicle speed is relatively slow. If the vehicle selects Siming district with no accurate speed prediction, it will take more travelling time, even if the distance is short.


### 4.8. Vehicle Routing Solution

With four vehicles, TDVRPSC schemes including the route number, customer visiting sequence, total route length, travelling starting time, travelling ending time, and total travelling time are shown in Table12. The parameter settings are the same as those in Section 4.5.Table 12
Vehicle routing solution schemes for TDVRPSC.
RouteCustomer visiting sequenceLength (km)StartEndTime (min)10 ⟶ 8 ⟶ 11 ⟶ 10 ⟶ 12 ⟶ 9 ⟶ 14 ⟶ 0109.603316:3018:4013020 ⟶ 6 ⟶ 4 ⟶ 7 ⟶ 5 ⟶ 3 ⟶ 047.314616:3017:578730 ⟶ 2 ⟶ 1 ⟶ 026.549016:3017:265640 ⟶ 13 ⟶ 15 ⟶ 034.740716:3017:2757Table12 gives the minimal expected total travelling time with different p, different departure time, and different algorithms. Its performance is compared with CA, CS, and ACO. The parameter settings of the algorithms are the same as those in Section 4.5.Table13 shows CACS is better than CA, CS, and ACO algorithms in solving TDVRPSC. With the increase of p, the minimal total expected time also increases. The total time also depends on the departure time. When the departure time is not the rush hour with traffic congestion, the total time is lower. When the departure time is the rush hour with traffic congestion, the minimal total expected time will increase.Table 13
Results with differentp, different departure time, and different algorithms (min).
p (%)Departure timeCACSCACSACO2516:3022322822522617:302312372342365016:3027628227928017:302862812772797516:3031231831431717:3032432932632810016:3033033733333617:30343349345347Compared with other kinds of cultural algorithms, such as adaptive immune clonal selection cultural algorithm, the convergence rate of Cuckoo Search algorithm is not sensitive to the parameters used. Therefore, adaptive Cultural Algorithm Based Cuckoo Search has fewer parameters to be fine-tuned in Cuckoo Search and is more robust for optimization problems.


## 4.1. Experimental Results with Deterministic Demand Probability

The CACS algorithm is tested on Intel® Core™ i3-4150T CPU @3.00 GHz, 64-bit operating system with memory of 4.00 GB, and x64-based processor. The VRP library used is published on the Networking and Emerging Optimization homepage [27]. The capacitated VRP instance library was selected. The first layer index is the first letter of word in the instance library. The second layer index is the number of customers, and the third layer index is the number of vehicles. Take A, B, P, and E sets as examples. Set A has 27 VRP instances from A-n32-k5 to A-n80-k10 proposed by Augerat et al. Set B has 23 VRP instances from B-n31-k5 to B-n78-k10 proposed by Augerat et al. Set P has 24 VRP instances from P-n16-k8 to P-n101-k4 proposed by Augerat et al. Set E has 15 VRP instances from E-n13-k4 to E-n101-k14 proposed by Augerat et al.Consider the case when the arrival probability is 1. Table2 gives the shortest path length, average iteration times, and average CPU computing time of each instance.Table 2
Calculation result of different instances when the demand probability is 1.
InstanceA-n32-k5A-n39-k6B-n31-k5P-n16-k8P-n19-k2N3239311619K56582Optimal solution (m)784.612831.2696662.3275451.9471212.6569Average iteration number427648329193308Average CPU computing time (s)18.968827.15636.81253.29696.7188Published optimal solution (m)784831672450212Improvement ratio (%)001.4393600For instance, B-n31-k5 algorithm has found a better solution than the published optimal solution, whose scheme is shown in Table3.Table 3
Comparison of new solution and published optimal solution for B-n31-k5.
RoutePublished optimal solutionNew solutionTravelling sequenceDistanceCapacityTravelling sequenceDistanceCapacity130 ⟶ 23 ⟶ 8 ⟶ 12 ⟶ 28 ⟶ 26109.6811977 ⟶ 30 ⟶ 5 ⟶ 25 ⟶ 18 ⟶ 16 ⟶ 21119.237289397221 ⟶ 16 ⟶ 18 ⟶ 25 ⟶ 5 ⟶ 4 ⟶ 29112.99158612 ⟶ 24 ⟶ 19 ⟶ 1 ⟶ 3 ⟶ 6 ⟶ 9102.36598029937 ⟶ 17 ⟶ 13 ⟶ 6 ⟶ 9 ⟶ 22110.82879620 ⟶ 27 ⟶ 10 ⟶ 2 ⟶ 14 ⟶ 15 ⟶ 11 ⟶ 28 ⟶ 8240.508851895420 ⟶ 27 ⟶ 10 ⟶ 2183.5752384 ⟶ 29 ⟶ 22 ⟶ 13 ⟶ 17 ⟶ 23108.715642299514 ⟶ 15 ⟶ 11 ⟶ 24 ⟶ 19 ⟶ 1 ⟶ 3155.4298952698.0612053822Total672.506412668.889412Figure4 shows the shortest path of 70 nodes when demand probability is 1. The horizontal axis represents the horizontal axis of a city map with the unit of meter. The vertical axis represents the longitudinal axis of a city map with the unit of meter.Figure 4
Optimal path of 70 nodes when demand probability is 1.Figure5 shows the objective function convergence curves of 70 nodes when demand probability is 1. The horizontal axis represents iteration times. The vertical axis represents the best individual in each iteration.Figure 5
Objective function convergence curves of 70 nodes when demand probability is 1.


## 4.2. Parameter Influence Analysis

Figure6 shows the average results of E-n22-k4, B-n31-k5, and B-n64-k9 after 20 times calculation under different demand probability p with IT 100. The updating rate of situational knowledge α is 0.3. The updating rate of standard knowledge β is 0.5. The population size m is 50. τil is 1.Figure 6
Calculation results of differentp with IT 100. (a) E-n22-k4, (b) B-n31-k5, and (c) B-n64-k9.
(a)(b)(c)Figure7 shows average results of 20 times calculation under different α with IT 100. p is 0.45. β is 0.5 and m is 50.Figure 7
Calculation results of differentα with IT 100. (a) E-n22-k4, (b) B-n31-k5, and (c) B-n64-k9.
(a)(b)(c)Figure8 shows average results of 20 times calculation under different β with IT 100. p is 0.45. α is 0.45 and m is 50.Figure 8
Calculation results of differentβ with IT 100. (a) E-n22-k4, (b) B-n31-k5, and (c) B-n64-k9.
(a)(b)(c)Figure9 shows average results of 20 times calculation under different m with IT 100. p is 0.45. α is 0.3 and β is 0.5.Figure 9
Calculation results of different m with IT 100. (a) E-n22-k4, (b) B-n31-k5, and (c) B-n64-k9.
(a)(b)(c)From the above figures, whenm is larger, the optimal solution is better. Increasing m is advantageous to searching for better solutions. However, when m reaches a certain value, the optimal solution is almost unchanged as m increases. In this instances, the optimal parameters are p=0.45, α=0.3, β=0.5, and m=50.The sensitivity analysis of iterations times is shown in Table4. In the case of E-n22-k4, when IT is too small, the calculation results still need to be improved and Δf/ΔNmaxis large. When IT reaches 130, Δf/ΔNmax remains close to zero as IT increases. The results are almost stable.Table 4
Calculation results of E-n22-k4 with different iteration times.
Nmax100110120130140Mean379.3766377.4962375.8675375.1023375.1023Std.0.76250.21660.10260.01250.0063ΔNmax1010101010ΔMean1.88041.62870.765200ΔMean/ΔNmax0.188040.112870.0765200We have also tried to vary the number of host nests and the discovery probabilitypa in Cuckoo Search algorithm. We have used n = 5, 10, 15, 20, 50, 100, 150, 250, and 500 and pa = 0.01, 0.05, 0.1, 0.15, 0.2, 0.25, 0.4, and 0.5. We found that n = 15 and pa = 0.25 are sufficient for most optimization problems. Results and analysis also imply that the convergence rate of Cuckoo Search algorithm is not sensitive to the parameters used.Apart from the number of host nests, there is one parameterpa. This also means that we do not have to fine-tune these parameters for a specific problem. Subsequently, Cuckoo Search algorithm is more generic and robust for many optimization problems. There are fewer parameters to be fine-tuned in Cuckoo Search than in other metaheuristic algorithms.


## 4.3. VRPSC Experiments with Fixed Demand Probability

p is set as different real numbers from 0.1 to 0.9. The parameters are set as follows: α=0.3, β=0.5, m=50, and τil=1. The initial vehicle speeds v are set as stochastic real numbers within the range of [1, 6]. Table 5 gives the average results of shortest expected total travelling time and servicing time for P-n16-k8 case with different IT and different p under 20 times calculation.Table 5
Results with different IT under differentp.
p1002003004005000.147.01664747.01664747.01664747.01664747.0166470.265.35669665.35668065.35668065.35668065.3566810.374.81729074.81697874.81697574.81697774.8169840.481.35662781.35402281.35401081.35401781.3540640.586.76676686.75332886.75322086.75316486.7533280.691.63513591.58497891.58353291.58332291.5833180.796.18311196.04008496.02895496.02695096.0269810.8100.469282100.163211100.110735100.106147100.1049060.9104.306277103.887785103.744337103.739400103.737179Table5 shows that, with the increasing of p, the objective function also increases. With the increasing of IT, the objective function reduces. However, when IT reaches 300, even if IT continues to increase, the shortest expected total travelling time and servicing time will not be reduced and will remain almost the same.


## 4.4. VRPSC Experiments with Distinguishing Customer Characteristics

Customers can be divided into large customers, small and medium customers, and retail customers according to the weight of the goods they order. Consider the influence of customer property in these three cases. In the case of B-n31-k5, the customers who order goods weighing in the range of [19, 25] belong to large customers. The customers who order goods weighing in the range of [9, 18] belong to small and medium customers. The customers who order goods weighing in the range of [2, 8] belong to small and retail customers. The parameters are set as follows:p=0.45, α=0.3, β=0.5, m=50, and τil=1. The initial vehicle speeds v are set as stochastic real numbers within the range of [1, 6].Table6 shows the optimal solution, average solution, and standard deviation in three cases of different customer types after 20 times calculation, when the demand probability is 0.8 or 0.9. In Table 6, 1st represents the result for the first calculation, 2nd denotes the results for the second calculation, and so on.Table 6
Objective function with different customer types and different demand probability.
Probability0.90.8CustomerLargeSmallRetailLargeSmallRetail1st596.155638.2462656.2919585.8801590.4435638.04122nd594.4813637.9749651.5232585.6425594.5347636.50343rd590.7583635.9657653.6638585.8801594.5347637.89914th596.155638.7321659.7856587.4265594.5347636.50345th590.7583639.7165653.4783582.6757590.4435634.85296th590.7583637.8215655.5595588.2996590.4435636.50347th590.7583637.9749657.6003582.6757590.4435634.85298th596.155637.8215656.7167588.2996594.5347636.86929th590.7583637.9749651.375582.6757590.4435636.503410th596.155638.6415654.8217588.2996594.5347636.5034Ave.593.2893638.087655.0816585.7755592.4891636.5032Std.2.7141530.9476772.6643242.3735792.1562521.04854Table6 shows when the large customers have least objective function. The small and medium customers have second shortest objective function. The retail customers have the longest objective function. This is because when the vehicle does not need to serve large customers, it can freely and flexibly arrange its routes. When vehicles need to serve large customers, with the limit of maximum carrying capacity, the choice of driving routes has more constraints. With the increase of p, the objective function also increases. The objective function with p of 0.9 is greater than that with p of 0.8.


## 4.5. Comparative Analysis

Table6 gives performance of CACS, CA, and CS with different p. Since CACS uses CA and CS, it is more complex than these two algorithms, while its time cost and storage cost are slightly higher compared to the latter two. However, the rapid development of computer hardware technology makes this difference for solving problems of general scale not apparent. From the experimental results, CACS performs better than CA, CS, and ACO.To avoid sensitivity to parameters selection, which is one of the main limitations in swarm intelligence algorithms, adaptive mechanism is adopted into CACS. ACO [28] used for comparison also proposed machine learning strategies to control the parameter adaptation. The objective function by different p and different algorithms is listed in Table 7.Table 7
Objective function with differentp and different algorithms.
pCACSCACSACOResultsRate (%)ResultsRate (%)ResultsRate (%)0.147.01664747.1056540.1889547.1341360.2492747.0910610.158020.265.35668165.3745870.0273965.3943740.0576465.3647340.012320.374.81698474.8348670.0239074.9212100.1391174.9014230.112730.481.35406481.5048640.1850281.4314570.0950481.4017630.058600.586.75332886.8912110.1586886.7931340.0458686.8014750.055470.691.58331891.7163010.1449991.6137060.0331791.6546010.077770.796.02698196.1053460.0815496.0713410.0461796.0578200.032100.8100.104906100.1159630.01104100.1100360.00512100.1134670.008550.9103.737179103.9213470.17722103.8641360.12223103.9163160.12223


## 4.6. TDVRPSC Using AFKF Speed Prediction

Based on the technology of AFKF, the highway data analysis platform is established. After extracting the information and constructing appropriate traffic prediction model, the future vehicle speed in different business districts can be effectively estimated. When the speed of a road is detected to be below the threshold, it can be inferred that this road is congested due to accidents or other reasons. Rapid identification of congested roads is completed. This is also useful for emergencies. When traffic congestion is detected in a certain road, decision support and vehicle dispatching optimization can adjust VRP scheme in time and avoid congested roads, so as to reduce transportation time and cost. In the field of real-time traffic forecasting, the fast information processing and mining ability of large data have very high reliability for real-time VRP scheduling. The accuracy of VRP scheduling is improved for the time-varying actual urban traffic road.Table8 lists the traffic speed and congestion delay index in different business districts of Xiamen at 16:25 on February 17, 2019 [29]. Congestion delay index is used to evaluate the degree of urban congestion. It is the ratio of the actual travel time to the free travel time of a trip.Table 8
Traffic speed and congestion delay index in Xiamen (2019.2.17, 16:25).
RankingBusiness districtCongestion delay indexSpeed (km/h)1Siming1.5827.212Huli1.3336.463Haicang1.2544.824Jimei1.2441.805Tongan1.1347.74Table9 lists the traffic speed in different business districts of Xiamen from 16:30 to 20:00 on February 17, 2019. The unit of speed is km/h.Table 9
Traffic speed in Xiamen during 2019.2.17, 16:30–20:00, in km/h.
TimeSimingHuliHaicangJimeiTongan16:3027.2536.8845.1342.0346.8916:4027.3336.3842.0146.1946.7816:5027.8536.2545.3942.3446.0517:0028.0736.5345.7242.5245.5517:1028.8336.6942.3845.9945.8917:2028.4636.4641.7045.7844.9017:3028.2736.7140.9342.9245.6417:4027.5736.6340.9242.8244.8317:5026.8841.8136.3341.1843.0118:0026.4140.5636.0942.4641.5818:1026.6840.1335.6842.1441.3818:2026.8140.2535.4140.8444.5618:3026.5935.3741.6140.7645.9518:4027.0441.4336.1441.1346.1818:5027.5436.5243.2941.6646.6519:0028.4337.0645.6042.4747.0719:1029.0737.0746.2343.4747.9219:2029.9437.7846.9343.547.9319:3030.2138.0646.6743.8348.1819:4031.0438.2147.0144.2248.7819:5031.2238.4647.0243.9548.5020:0031.6938.9447.4044.0249.10Figure10 shows Xiamen city map at 16:25 on February 17, 2019. Different business districts are marked with their congestion delay indexes.Figure 10
Xiamen map for real-time experiment (2019.2.17, 16:25).Several branches of Yuantong Express in Xiamen city distributed in different business districts are used. Table9 lists their business district, latitude, longitude, and the Mercator coordination converted. The No. 0 branch is the head office. VRP task is to distribute goods to all other branches.


## 4.7. Accuracy of Prediction Model

Figure11 shows the traffic speed prediction in Huli district using AFKF prediction from 16:30 to 20:00 on February 17, 2019. The horizontal axis represents time. The vertical axis represents the traffic speed in km/h. The star marking connected with solid wire represents the real speed. The joint spider connected with dotted lines represents the predicted speed.Figure 11
Traffic speed prediction in Huli using adaptive fractional Kalman prediction.Table10 shows the prediction mean square error (MSE) and mean absolute percentage error (MAPE) by AFKF, fractional Kalman filter (FKF), and Kalman filter (KF). Table 11 shows AFKF has smaller prediction error and better prediction accuracy than FKF and KF.Table 10
Prediction error with different algorithms (km/h).
Business districtAFKFFKFKFMSEMAPEMSEMAPEMSEMAPESiming1.44433.87731.98884.59952.62755.2535Huli2.76323.49292.78233.61893.10273.6632Haicang2.29123.42452.33343.43222.56533.4551Jimei1.75562.29081.76632.49312.66642.736Tongan1.49512.07321.94252.21981.92842.4390Table 11
Geographical location of Yuantong Express.
No.AreaLongitude (degree)Latitude (degree)Vertical coordinate (m)Horizontal ordinate (m)0Siming118.07654824.4496640061477627399761Siming118.08607424.4839640061630527456702Siming118.13810824.4774680062529827449303Siming118.1701624.4824230062981827452494Huli118.14282224.5146570062543727515055Huli118.17888724.5041870063222127495196Huli118.14729724.4889770062673727472517Huli118.14467624.5219740062595927524768Jimei118.10840724.6100440062012227684779Jimei118.1046424.60853600619048276918810Tongan118.15615524.72603900627708279065711Tongan118.14663524.69745200626189278561412Tongan118.13816224.70864600624935278773813Haicang118.01761924.53105700604864275386914Haicang118.04275824.49852800608609274878915Haicang118.01589924.530205006043892753588Figure12 shows that the prediction errors vary with different fractional orders of Kalman filter. The horizontal axis represents fractional order. The vertical axis represents the prediction error of traffic speed in km/h. The five-pointed star represents the mean relative error. The asterisk represents the mean square error. For certain business districts, there exists a most suitable fractional order for reducing its prediction error. Compared with integer-order filters, fractional-order Kalman filters have wider parameter selection range and higher prediction accuracy.Figure 12
Prediction error of traffic speed with different fractional orders.If the prediction model is not accurate, wrong route selection may be made. For example, the traffic flow is busy in Siming district and its vehicle speed is relatively slow. If the vehicle selects Siming district with no accurate speed prediction, it will take more travelling time, even if the distance is short.


## 4.8. Vehicle Routing Solution

With four vehicles, TDVRPSC schemes including the route number, customer visiting sequence, total route length, travelling starting time, travelling ending time, and total travelling time are shown in Table12. The parameter settings are the same as those in Section 4.5.Table 12
Vehicle routing solution schemes for TDVRPSC.
RouteCustomer visiting sequenceLength (km)StartEndTime (min)10 ⟶ 8 ⟶ 11 ⟶ 10 ⟶ 12 ⟶ 9 ⟶ 14 ⟶ 0109.603316:3018:4013020 ⟶ 6 ⟶ 4 ⟶ 7 ⟶ 5 ⟶ 3 ⟶ 047.314616:3017:578730 ⟶ 2 ⟶ 1 ⟶ 026.549016:3017:265640 ⟶ 13 ⟶ 15 ⟶ 034.740716:3017:2757Table12 gives the minimal expected total travelling time with different p, different departure time, and different algorithms. Its performance is compared with CA, CS, and ACO. The parameter settings of the algorithms are the same as those in Section 4.5.Table13 shows CACS is better than CA, CS, and ACO algorithms in solving TDVRPSC. With the increase of p, the minimal total expected time also increases. The total time also depends on the departure time. When the departure time is not the rush hour with traffic congestion, the total time is lower. When the departure time is the rush hour with traffic congestion, the minimal total expected time will increase.Table 13
Results with differentp, different departure time, and different algorithms (min).
p (%)Departure timeCACSCACSACO2516:3022322822522617:302312372342365016:3027628227928017:302862812772797516:3031231831431717:3032432932632810016:3033033733333617:30343349345347Compared with other kinds of cultural algorithms, such as adaptive immune clonal selection cultural algorithm, the convergence rate of Cuckoo Search algorithm is not sensitive to the parameters used. Therefore, adaptive Cultural Algorithm Based Cuckoo Search has fewer parameters to be fine-tuned in Cuckoo Search and is more robust for optimization problems.


## 5. Conclusion

The main contributions of this paper are summarized as follows:(1)
An adaptive CACS has been proposed and its convergence is proved.(2)
AFKF for traffic speed prediction is proposed. An adaptive mechanism is adopted with the covariance of prediction error. Its mathematical process is proved.(3)
New solutions of VRP are discovered, which are better than the published solutions.(4)
CACS and AFKF are used for TDVRPSC with real-time urban vehicle data. Traffic speed is predicted and vehicle routing can be adjusted in time to avoid congested roads to reduce transportation time and cost.The next step is to study better TDVRPSC models and design new methods to deal with uncertainties.



---

*Source: 7258780-2020-07-24.xml*