# Highway Traffic Flow Nonlinear Character Analysis and Prediction

Authors: Meng Hui; Lin Bai; YanBo Li; QiSheng Wu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

In order to meet the highway guidance demand, this work studies the short-term traffic flow prediction method of highway. The Yu-Wu highway which is the main road in Chongqing, China, traffic flow time series is taken as the study object. It uses phase space reconstruction theory and Lyapunov exponent to analyze the nonlinear character of traffic flow. A new Volterra prediction method based on model order reduction via quadratic-linear systems (QLMOR) is applied to predict the traffic flow. Compared with Taylor-expansion-based methods, these QLMOR-reduced Volterra models retain more information of the system and more accuracy. The simulation results using this new Volterra model to predict short time traffic flow reveal that the accuracy of chaotic traffic flow prediction is enough for highway guidance and could be a new reference for intelligent highway management.


---

## Body

## 1. Introduction

The prediction of highway traffic flow is a crucial part of transportation planning, traffic control, and intelligent transportation systems [1]. In particular, short-term traffic flow predicting can provide real time information for highway guidance. Also accurate predictions of traffic flow are essential part of efficient operations in advanced traffic management systems. In recent years, the short-term predicting technologies have gotten the attention of many researchers.Short-term traffic flow is a typical nonlinear time series in a time horizon of 15 min or less [2]. During the decades a variety of techniques have been applied to forecast the short-term traffic flow, such as fuzzy theory [3], neural networks [4], Kalman filter [5], and wavelet analysis [6] et al. But the models generated by these methods could not capture some strongly nonlinear characteristics of short-term traffic flow data. The prediction results may have big error. Chaos theory is an effective tool to study nonlinear system. This tool has been applied to forecast short-term traffic flow time series [7, 8]. The forecasting results are better than other methods applied before. However, the accuracy of chaotic forecasting result is not enough for real time highway guidance. There remains room for proposing new prediction technology to make the prediction more accurate.To use chaos theory we should judge whether the traffic flow data is a chaotic time series or not. The largest Lyapunov exponent is usually used to achieve this goal. Wolf algorithm [9], Jacobin algorithm [10], and small data method [11] are three algorithms often used to calculate Lyapunov exponent. If the traffic flow data is a chaotic time series, the phase space reconstruction theory could be applied to predict the traffic flow.In this paper, the traffic flow data of Yu-Wu highway in Chongqing, China, is taken as study object. We study the general day traffic flow character of highway and the holiday traffic flow model is not contained in this paper. Firstly, the approved small data method is used to judge whether the traffic flow data is a chaotic time series. Secondly, phase space reconstruction theory is used to reconstruct the traffic flow system. Thirdly, a new predicting model based on model order reduction via quadratic-linear systems is used to predict the chaotic traffic flow. This new prediction model has no truncation error and contains more information of original system. It is more accurate than other prediction models ever used before.


## 2. Traffic Flow Data Phase Space Reconstruction

### 2.1. Traffic Flow Data

All the traffic flow data are collected at Yu-Wu highway in Chongqing, China. Sampling time begins at 0:00 November 20, 2013, and ends at 0:00 November 27, 2013. Time interval is 5 minutes. Figure1 shows the traffic data.Figure 1
Traffic flow of Yu-Wu highway.


### 2.2. Mutual Information Method

Takens theorem [12] is the mathematical base of time series phase space reconstruction. Choosing optimal time delay τ and embedding dimension d is the key of phase space reconstruction. Mutual information method is an effective tool to calculate τ.Assuming chaotic time seriesr1,r2,… time delay is τ and embedding dimension is d. The reconstruction phase space Q of this time series is(1)Qt=q0t,q1t,…,qnt,where qn(t)=r(t+nτ). The average information of variable r is called information entropy which can be described as(2)HR=-∑iPrrilogPrri,where Pr(ri) is the probability of event r in i. Let r,u=[qt,qt+τ], if system (R,U) is a coupling system and r is given as ri the uncertainty of u is(3)HU∣ri=-∑jPu∣ruj∣rilog⁡Pu∣ruj∣ri=-∑iPruri,uiPrri×log⁡Pruri,uiPrri,where Pu∣r(uj∣ri) is the probability that a measurement of u yields uj, given that the measured value of r is ri. Given that q has been measured at time t the average uncertainty in a measurement of r at time t+τ is(4)HU∣R=∑iPrriHU∣ri=HR,U-HRHR,U=-∑i,jPruri,uilogPruri,ui,where H(R,U) is the uncertainty of u in isolation. H(U∣R) is the uncertainty of u in the measurement of r. The amount of a measurement of r reduces the uncertainty of u. The mutual information is(5)IU,R=HU-HU∣R=HQ+HR-HR,U=IR,U.The mutual information is a function of Pru which is the joint distribution probability of the event r and u. If a vector is a phase space reconstruction while the mutual information I first reaches the minimum the time delay could be the delay time of phase space reconstruction.Figure2 shows the traffic flow delay time calculated by the mutual information method. When τ=2 the system reaches the minimum so the optimal delay time τ=2.Figure 2
Mutual information method calculates delay time.


### 2.3. CAO Method

Time seriesy can be reconstructed in phase space as follows: Yi={y(i),y(i+τ),…,y(i+(d-1)τ)}, where τ is the delay time and d is the embedding dimension. The nearest neighbor point of i is YNN(i); the distance between two points is (6)Ld2i=Yi-YNNi2=∑j=1dyi+j-1τ-yNNi+j-1τ2.When the space dimension changes from d to d+1, the distance between two points is(7)Ld+12i=∑j=1d+1yi+j-1τ-yNNi+j-1τ2=Ld2i+yi+τd-yNNi+τd2.Define(8)a1i,d=yi+τd-yNNi+τdLdi.If a1(i,d)>Rτ,  Rτ∈[10,50], YNN(i) is the false nearest neighbor point. Substitute (6) into (8); we can get (9)a1i,d=yi+τd-yNNi+τdYi-YNN(i).CAO [13] made some changes to (9) yielding(10)a2i,d=Yd+1i-Yd+1NNiYdi-YdNNi.Here ‖·‖ is the Euclidian distance. Let E(d)=1/(N-dτ)∑i=1N-dτa2(i,d), E1d=Ed+1/E(d); E(d) is only dependent on the dimension d and delay time τ. It is found that E1(d) stops changing when d is greater than d0. d0 is the optimal embedding dimension of reconstruction space.For time series data from a random set of numbers,E1(d) will never attain a saturation value as d increases in principle. But in practical computations, it is difficult to resolve whether E1(d) is slowly increasing or has stopped changing if d is sufficiently large. In fact, since available observed data samples are limited, it may happen that E1(d) stops changing at some dimension although the time series is random. Therefore, we define(11)E′d=1N-dτ∑i=1N-dτyi+τd-yNNi+τd,E2d=E′d+1E′(d).For random series,E2(d)=1. For deterministic series, there must exist d making E2(d)≠1.Figure3 shows the embedding dimension calculated by CAO method. The optimal embedding dimension d=10.Figure 3
Embedding dimension calculated by CAO method.


## 2.1. Traffic Flow Data

All the traffic flow data are collected at Yu-Wu highway in Chongqing, China. Sampling time begins at 0:00 November 20, 2013, and ends at 0:00 November 27, 2013. Time interval is 5 minutes. Figure1 shows the traffic data.Figure 1
Traffic flow of Yu-Wu highway.


## 2.2. Mutual Information Method

Takens theorem [12] is the mathematical base of time series phase space reconstruction. Choosing optimal time delay τ and embedding dimension d is the key of phase space reconstruction. Mutual information method is an effective tool to calculate τ.Assuming chaotic time seriesr1,r2,… time delay is τ and embedding dimension is d. The reconstruction phase space Q of this time series is(1)Qt=q0t,q1t,…,qnt,where qn(t)=r(t+nτ). The average information of variable r is called information entropy which can be described as(2)HR=-∑iPrrilogPrri,where Pr(ri) is the probability of event r in i. Let r,u=[qt,qt+τ], if system (R,U) is a coupling system and r is given as ri the uncertainty of u is(3)HU∣ri=-∑jPu∣ruj∣rilog⁡Pu∣ruj∣ri=-∑iPruri,uiPrri×log⁡Pruri,uiPrri,where Pu∣r(uj∣ri) is the probability that a measurement of u yields uj, given that the measured value of r is ri. Given that q has been measured at time t the average uncertainty in a measurement of r at time t+τ is(4)HU∣R=∑iPrriHU∣ri=HR,U-HRHR,U=-∑i,jPruri,uilogPruri,ui,where H(R,U) is the uncertainty of u in isolation. H(U∣R) is the uncertainty of u in the measurement of r. The amount of a measurement of r reduces the uncertainty of u. The mutual information is(5)IU,R=HU-HU∣R=HQ+HR-HR,U=IR,U.The mutual information is a function of Pru which is the joint distribution probability of the event r and u. If a vector is a phase space reconstruction while the mutual information I first reaches the minimum the time delay could be the delay time of phase space reconstruction.Figure2 shows the traffic flow delay time calculated by the mutual information method. When τ=2 the system reaches the minimum so the optimal delay time τ=2.Figure 2
Mutual information method calculates delay time.


## 2.3. CAO Method

Time seriesy can be reconstructed in phase space as follows: Yi={y(i),y(i+τ),…,y(i+(d-1)τ)}, where τ is the delay time and d is the embedding dimension. The nearest neighbor point of i is YNN(i); the distance between two points is (6)Ld2i=Yi-YNNi2=∑j=1dyi+j-1τ-yNNi+j-1τ2.When the space dimension changes from d to d+1, the distance between two points is(7)Ld+12i=∑j=1d+1yi+j-1τ-yNNi+j-1τ2=Ld2i+yi+τd-yNNi+τd2.Define(8)a1i,d=yi+τd-yNNi+τdLdi.If a1(i,d)>Rτ,  Rτ∈[10,50], YNN(i) is the false nearest neighbor point. Substitute (6) into (8); we can get (9)a1i,d=yi+τd-yNNi+τdYi-YNN(i).CAO [13] made some changes to (9) yielding(10)a2i,d=Yd+1i-Yd+1NNiYdi-YdNNi.Here ‖·‖ is the Euclidian distance. Let E(d)=1/(N-dτ)∑i=1N-dτa2(i,d), E1d=Ed+1/E(d); E(d) is only dependent on the dimension d and delay time τ. It is found that E1(d) stops changing when d is greater than d0. d0 is the optimal embedding dimension of reconstruction space.For time series data from a random set of numbers,E1(d) will never attain a saturation value as d increases in principle. But in practical computations, it is difficult to resolve whether E1(d) is slowly increasing or has stopped changing if d is sufficiently large. In fact, since available observed data samples are limited, it may happen that E1(d) stops changing at some dimension although the time series is random. Therefore, we define(11)E′d=1N-dτ∑i=1N-dτyi+τd-yNNi+τd,E2d=E′d+1E′(d).For random series,E2(d)=1. For deterministic series, there must exist d making E2(d)≠1.Figure3 shows the embedding dimension calculated by CAO method. The optimal embedding dimension d=10.Figure 3
Embedding dimension calculated by CAO method.


## 3. Chaotic Character of Traffic Flow

Correlation dimension and largest Lyapunov exponent are two main indexes to distinguish chaotic system from others.


### 3.1. Correlation Dimension

G-P algorithm is the most popular method to calculate correlation dimension for its simplicity. Figure4 shows the lnC(r,m)-lnr curves. These curves got through increasing embedding dimension until the slope of the curve’s linear part is almost constant. Then, the correlation dimension estimation can be obtained using least square regression.Figure 4
Correlation dimension.Through Figure4 we can get the correlation dimension g(m)=5.4116; its fractal reveals that the system is a chaotic system.


### 3.2. The Largest Lyapunov Exponent of Traffic Flow

In the reconstructed phase space the distance between each pointXj and its nearest neighbor point Xj′ is expressed as(12)dj0=minXj⁡Xj-Xj′,j-j′>P,where P is the mean period. Based on the study in [14] the largest Lyapunov exponent is defined as (13)λ1i=1iΔt·1TN-i∑j=T0TN-ilndjidj0,where Δt is sampling period and dji is the distance between the jth pair of nearest neighbors after i discrete-time steps. Sato also gives an alternate form of (13)(14)λ1i,k=1kΔt·1TN-k∑j=T0TN-ilndji+kdji,where k is constant. Combining (13) and (14) yields(15)dji≈Cjeλ1iΔt,Cj=dj0.By taking the logarithm of both sides of above equation, we obtain(16)lndji=lnCj+λ1iΔt.The largest Lyapunov exponent equals the slope of (16).We can obtain the largest Lyapunov exponent of traffic flow dataλ=0.1760. Its a positive value which means the system exist chaos phenomena.


## 3.1. Correlation Dimension

G-P algorithm is the most popular method to calculate correlation dimension for its simplicity. Figure4 shows the lnC(r,m)-lnr curves. These curves got through increasing embedding dimension until the slope of the curve’s linear part is almost constant. Then, the correlation dimension estimation can be obtained using least square regression.Figure 4
Correlation dimension.Through Figure4 we can get the correlation dimension g(m)=5.4116; its fractal reveals that the system is a chaotic system.


## 3.2. The Largest Lyapunov Exponent of Traffic Flow

In the reconstructed phase space the distance between each pointXj and its nearest neighbor point Xj′ is expressed as(12)dj0=minXj⁡Xj-Xj′,j-j′>P,where P is the mean period. Based on the study in [14] the largest Lyapunov exponent is defined as (13)λ1i=1iΔt·1TN-i∑j=T0TN-ilndjidj0,where Δt is sampling period and dji is the distance between the jth pair of nearest neighbors after i discrete-time steps. Sato also gives an alternate form of (13)(14)λ1i,k=1kΔt·1TN-k∑j=T0TN-ilndji+kdji,where k is constant. Combining (13) and (14) yields(15)dji≈Cjeλ1iΔt,Cj=dj0.By taking the logarithm of both sides of above equation, we obtain(16)lndji=lnCj+λ1iΔt.The largest Lyapunov exponent equals the slope of (16).We can obtain the largest Lyapunov exponent of traffic flow dataλ=0.1760. Its a positive value which means the system exist chaos phenomena.


## 4. The Volterra Prediction Model

For single-input-multiple-output system in the form of(17)x˙=Gx+Bu,where x∈Rn are the state variables, u=acos⁡(ωt),  ω∈R is the input signal.Assume an autonomous system as follows:(18)v˙=Sv,S=0ω-ω0,v0=10.Let v=[v1,v2T,  v∈R2. The system matrix S could be decomposition(19)SQ=QΛ,Λ=jω00-jω,Q=121jj1.Q is unitary; the close form solution to (18) is(20)vt=cos⁡(ωt)sin⁡(ωt)-sin⁡(ωt)cos⁡(ωt)10=cos⁡ωt-sin⁡(ωt).Also the input u=acos⁡(ωt) can be taken as(21)u=Lv,L=a0.Combining (17)–(21) yields(22)x˙v˙=GBL0Sxv.If the eigenvalue λi(G)-λj(S)≠0,  (i=1,2,…,n,j=1,2), the Sylvester equation(23)GΠ+BL=ΠS,Π∈Rn×2has unique solution. Then let Ik denote the kth order identity matrix; we get(24)GBL0SInΠ0I2=InΠ0I2G00Sfrom which the solution to (22) is (25)xt-Πvt=eGtx0-Πv0.Equation (25) shows the system solution dependent on initial condition x0-Πv0 and steady state Πv(t). If x(0)=Πv(0) the system is driven into period steady state regardless of the fact that G is stable or not.The study [14] reveals that nonlinear ordinary differential equations can be rewritten equivalently in a special representation, quadratic-linear differential algebraic equations (QLDAEs). The QLDAEs take the form as(26)Cx˙=G1x+G2x⊗x+D1xu+D2x⊗xu+Bu.Similarly, any generalized polynomial system can be converted into generalized QLDAEs as follows:(27)ddtC1x+C2x⊗x=G1x+G2x⊗x+D1xu+D2x⊗xu+Bu.The system (27) is equivalent representation of the original system; this step makes no approximation of original system.Based on the Volterra theory that the response of a nonlinear system can be decomposed into responses of a series of homogeneous nonlinear systems, the solution of system can be written as follows:(28)xt=∑n=1∞xnt,where xn(t)=∫-∞+∞⋯∫-∞+∞hnσ1,…,σnu(t-σ1),…, and u(t-σn)dσ1,…,dσn. hn(σ1,…,σn) is nth order Volterra kernel. If the input of system (17) is αu the response of system is (29)xt=αx1t+α2x2t+α3x3t+⋯.Setting α=1 and substituting (29) into (26) we can get(30)Cx˙1=G1x1+Bu,(31)Cx˙2=G1x2+G2x1⊗x1+D1x1u,(32)Cx˙3=G1x3+G2x1⊗x2+x2⊗x1+D1x2u+D2x1⊗x1uand so on.Substituting (21) into (30) the steady state solution of system is x1(t)=Π1v(t), and  Π1∈Rn×2 solves the Sylvester equation(33)G1Π1+B1=Π1S,B1=BL.A unique Π1 always exists since the eigenvalues of G1 and the imaginary eigenvalues in S never add to zero. When the initial condition x1(0)=Π1v(0) and input u=acos⁡(ωt) the steady state solution of system is x1(t)=Π1v(t). Therefore, the time argument in the steady state solution can be omitted without ambiguity; namely, x1=Π1v.Setx⊗x=x②, because x1②=Π1②v② and u is a scalar so x1u=x1⊗u=Π1v⊗Lv=(Π1⊗L)v②. Equation (31) can be rewritten as(34)Cx˙2=G1x2+G2Π1②+D1Π1⊗Lv②≜G1x2+B2v②.Moreover, the derivative of v② is(35)v˙②=v˙⊗v+v⊗v˙=I2⊗S+S⊗I2v②≜S2v②.Equations (34) and (35) could construct another Sylvester equation like (23) with unknown Π2∈Rn×4. The new Sylvester equation is (36)G1Π2+B2=Π2S2.The steady state solution of (36) is x2=Π2v②. Also x1⊗x2=(Π1⊗Π2)v③, x2u=(Π2⊗L)v③ and so on. Setting B3=G2(Π1⊗Π2+Π2⊗Π1)+D1(Π2⊗L)+D2LΠ1②, B3=G2(Π1⊗Π2+Π2⊗Π1)+D1(Π2⊗L)+D2LΠ1②, the third Sylvester equation is(37)G1Π3+B3=Π3S3.The steady state solution of (37) is x3=Π3v③.Based on the analysis we could get the pattern of the steady state responsexk=ΠkvⓀ,  Πk∈Rn×2k,  k=1,2,…, and the steady state solution of system is(38)xt=∑n=1∞xnt=Π1vt+Π2v②t+Π3v③t+⋯.The initial condition x(0)=∑n=1∞xn(0). Equation (38) is exactly the same desired solution sought by SN algorithm.In this Volterra algorithm because‖vⓚ‖=1, ‖xkt‖=‖Πkvⓚ‖≤‖Πk‖, the number of Sylvester equations has limitation. This algorithm could converge below any preset tolerance.In order to simplify the algorithm programing, the traffic flow data is normalized as follows:(39)Xn=zn-z(n)-Maxzn-Maxzn,where X(n) is the normalization time series, z(n) is the original time series, z(n)- is average value of original time series, and Max[z(n)] and Min[z(n)] are the maximum and minimum value of original time series.


## 5. Simulation Results

The effectiveness of the Volterra prediction algorithm is verified by the traffic flow shown in Figure1. The 1500 data of 2500 traffic flow data are used to train the Volterra algorithm. The last 2 hours and 4 hours traffic data of November 27, 2013, are used to verify the effectiveness of the prediction algorithm mentioned in this paper. Also comparing with RBFNN prediction algorithm and old Volterra algorithm, the results show that the accuracy of new Volterra algorithm is better than RBFNN algorithm and old Volterra algorithm. The prediction error is shown in Table 1. The new Volterra algorithm prediction results are shown in Figure 5, and the RBFNN prediction results are shown in Figure 6; the old Volterra algorithm results are shown in Figure 7.Table 1
Error of two algorithms.
Model
2 h MAPE%
4 h MAPE%
New Volterra algorithm
3.98
7.2
RBFNN algorithm
5.21
10.04
Volterra algorithm
7.82
7.93Figure 5
Prediction results of new Volterra algorithm: (a) 2 hours, (b) 4 hours.
(a)
(b)Figure 6
Prediction results of RBFNN algorithm: (a) 2 hours, (b) 4 hours.
(a)
(b)Figure 7
Prediction results of Volterra algorithm: (a) 2 hours, (b) 4 hours.
(a)
(b)Through Figures5, 6, and 7 it can be found that the trend of prediction is similar with real traffic flow. But the new Volterra algorithm is more accurate than RBFNN and old Volterra algorithm.


## 6. Conclusion

This paper analyzes the nonlinear character of Yu-Wu highway traffic flow time series. Using the largest Lyapunov exponent and correlation dimension to show that the traffic flow data are a chaotic time series, the traffic flow is reconstructed by phase space theory and the optimal time delay and embedding dimension are confirmed. Finally a new Volterra algorithm is used to predict the short-term traffic flow. The prediction results show that this algorithm could predict the short-term traffic flow effectively. The accuracy could satisfy the demand of highway guidance.



---

*Source: 902191-2015-09-16.xml*