Wiener Process

Solution algorithms, including Kalman Filter and Gibbs Sampler, mentioned in section 3.2 are applied to each rolling stage. After the estimation, state vectors of the same time in each stage will be averaged to give an estimated value.

4.2.2 Wiener Process

Wiener process is a continuous-time stochastic process, often called Brownian motion. A Wiener process Z is characterized by _t

1. Z₀ =0.

2. Z is almost surely continuous. _t

3. Z has independent increments with distribution _t ^Zt −^Zs ~^N

(

0,^t−^s

)

. As in discrete form,

−1 drift term W can be estimated by least square method, that

44  variance-covariance matrix is

( )( )

^T

4.2.3 Solution Framework

A rolling horizon structure is incorporated in the solution framework. The solution procedure first takes observations from t=1~m+1 into account, and utilize the solution framework presented in section 3.2 to estimate the state vectors,

+1

Ψ_m , and transition matrix, F_m+2. Second, the procedure continues to use the observations from t=2~m+2 to estimate their corresponding variables. The solution procedure will keep rolling until it reaches its final stage. State vectors of the same time in each rolling stage, i.e. ψ₂ would exists in stage 1 and 2, will be averaged to give an estimated value.

After deriving F by rolling horizon method. The estimation of _t F_t+h can be computed by

( )

t

4.3 Estimation of Path Flow Considering Link Travel Time

In the previous section (section 4.2), we relax the assumption of fixed transition matrix by introducing time-varying coefficient dynamic system. However, travel time is not addressed in that section. The travel time issue can be easily taking into consideration in this section. Notations used in this section is the same as that of section 3.2, 3.3, and 4.2.

The time-varying coefficient state space model considering travel time effect is shown as follows.

46 

t t

t W F

F = + ₋₁+ρ (4.12)

n m

m m t

F_{t t t}

t = ₁Ψ ₁+ , = , +1, +2,...,

Ψ _{− −} σ ^(4.13)

n t

H

y_t = _tψ_t^lag +γ_t, =1,2,3,... (4.14)

Equation 4.12 that describe the relationship of time-dependent coefficient is the same as equation 4.2 in the previous section (section 4.2). The transition equation 4.13 is also the same as equation 4.3 in the previous; while equation 4.14 the describe the link observation considering time lag effect is introduced in section 3.3.

The model proposed in this section is developed from models proposed in previous chapters. Therefore, solution framework can be straight forward derived by combining solution frameworks discussed in the previous chapters. The solution framework is illustrated in figure 4.1.

Figure 4.1 The Solution framework of Time-varying coefficient state space model considering travel time effect.

Rolling Horizon Structure Gibbs Sampler

Kalman Filter

Path - Link flow Transformation

Link Dynamics Path Flows

Time- dependent

path-link incidence matrix

Chapter 5

Convergence Control of Gibbs Sampler and Parallel Implementation

In this chapter, the convergence control method of Gibbs Sampler and parallel implementation of solution framework are discussed. Begin with the convergence assessment for single chain in section 5.1, parallel chain method is then illustrated in section 5.2. Finally, Parallel computing technique and its implementation is addressed in section 5.3.

5.1 Single Chain Convergence Assessments

Gibbs sampler, an example of a Markov Chain Monte Carlo method, is an algorithm that generate sequences of samples from a joint probability distribution. It has been used to tackle a wide variety of statistical problems. If the method had been used naively without convergence control, it might result a misleading answer.

Properties of applying Gibbs Sampler to state space model can be referred to Frühwirth-Schnatter (1994), Carter and Kohn (1994), Gamerman (1998), and Jungbacker and Koopman (2007); while analytic convergence rates for Gibbs sampler applied to uni-variate state space model can be referred to Pitt and Shephard (1996).

Theoretical guarantee of convergence is different from the practical requirement.

It is thus necessary to develop a tool that determines whether the chain is converged or not. Monitoring the chain that Gibbs sampler produce is a reasonable method.

48 

However, monitoring the elements in a chain would only expose strong non-stationarities. Therefore, it is more relevant to consider the cumulated averages;

since they need to stabilize for convergence to be achieved. Let the number of Gibbs iterations be denoted as g , and the corresponding element x^(g). The difference between cumulated averages should be less than ε,

ε

− ≤

−

∑

∑

−1 ) ( )

(

1 1 1

g g g

g x

x g

g . (5.1)

Possible alternatives to empirical averages would be importance sampling, conditional expectations, or Riemann approximation techniques. However, the average may only correspond to the exploration of a single mode of the distribution by the chain; the single chain methods can never guarantee that the whole support of target distribution has been explored (Robert and Cellier, 1998).

5.2 Multiple Chain Convergence Assessments

The main idea of Gibbs sampler, an iterative simulation method, is to draw values of a random variable x from a sequence of distributions that converge to a desired target distribution of x. Single chain (sequence) methods hardly bring information on the regions of the space it does not visit. Parallel chain methods try to overcome such defect by generating parallel chains, aiming at eliminating the dependence on initial conditions. The convergence control is most often based on the comparison of the estimations of different quantities for the parallel chains. More precisely, the criterion is based on the difference between a weighted estimator of variance for each chain and the variance of the estimators on the different chains. The estimation method composed of two steps. First, create an estimate of the target

distribution, centered about its mode (or modes), and over-dispersed in the sense of being more variable than the target distribution. The approximate distribution is then used to start several independent chains of the iterative simulation. The second step is to analyze the multiple chains to form a distributional estimate of the target random variable.

The monitor of convergence of the iterative simulation is to estimate the factor by which the scale of the current distribution for the target distribution might be reduced if the simulations were continued to the limit n→∞. This potential scale reduction can be estimated by

) 2 within–chain variance, W. n is the count of iterations, S_m are the standard deviation of each chain, B is the between-chain variance and the M is the number of parallel chains.

The between- and within-chain variances, B and W, are defined as follows.

( )

50 

∑

∑

= =

=

= ^M

m m n

i i m

m n ₁ M ₁

)

( 1

1 ξ , ξ ξ

ξ (5.6)

where ξ_m⁽ⁱ⁾ is the i-th element in chain m.

Once Rˆ is near 1, it is conclude that each set of the simulated values is close to the target distribution. In practice, the potential scale reduction is chosen to be 10%

(Rˆ≤1.1). After it converges, we can calculate the desired sample value (the state vector) based on the empirical distribution of the n simulated iterates for each simulated chain.

Some Gibbs Sampler estimation examples are demonstrated in this section. In the following examples, four parallel chains with different random seeds are illustrated. Example in Figure 5.1 converged (Rˆ ≤1.1) after 1.92×10⁵ iterations, and the estimated value is 23.75; example in Figure 5.2 converged after 4.12×10⁵ iteration, and the estimated value is 0.11.

0 5 10 15 20 25 30 35 40 45

Estimated Value

Iteration Counts

Chain 1 Chain‐2 Chain‐3 Chain‐4

Figure 5.1 Example 1 of four parallel chains.

‐6

‐4

‐2 0 2 4 6 8 10 12

Est im ate d V alu e

Iteration Counts

Chain 1 Chain‐2 Chain‐3 Chain‐4

Figure 5.2 Example 2 of four parallel chains.

Another example of eight parallel chains is shown in Figure 5.3.

30 40 50 60 70 80 90 100 110

Est im ate d V alu e

Iteration Counts

Chain 1 Chain‐2 Chain‐3 Chain‐4 Chain 5 Chain 6 Chain 7 Chain 8

Figure 5.3 Example of eight parallel chains.

These chains converged (Rˆ≤1.1) after 5.5×10⁴ iterations, and the estimated value is 68.57.

52 

5.3 Parallel Implementation

Gibbs sampler, a particular type of Markov Chain Monte Carlo method, requires tremendous iterations during computation; normally, there would be tens of thousands iterations in the proposed algorithm. To achieve real-time information requirement, parallel computing technique is introduced to increase the performance. The solution algorithm should be modified to adopt the parallel implementation. Consider the Gibbs sampler algorithm in section 3.3.2, the algorithm is separated into several computing parts at the WHILE-LOOP. With different random seed, each computing part will lead to a different solution chain. The chain in each computing part will then be gathered to check the convergence. In this situation, communication between computing nodes is minimum, and computing power can be easily increased without communication bandwidth limitation. Figure 5.4 describes the parallel architecture.

The parallel environment of this PC-cluster consists of 16 computing nodes; each contains 2 processors equivalent to Intel XEON 3.2 GHz and 1GB memory. Nodes are connected with a Gigabits Ethernet switch for MPI protocol and a 100 Mbits PCI fast Ethernet switch for Network File System (NFS) and Network Information System (NIS).

Figure 5.4 The Parallel computation structure.

In the server node, parameters used in our algorithm are initialized, so does the necessary input data. When assign jobs, these input data are been sent to computing nodes existing in the cluster through TCP/IP base intranet with Message Passing Interface (MPI) Library through gigabit switch. The computational procedure for the parallel process consists of:

1. Load input data and parameters. Initialize MPI environment.

2. Count the computing nodes exists in the cluster environment. Send data to each computing nodes.

3. Each computing nodes generate its own ψ^ˆ and Fˆ′ sequences by given input data. These results were sending to the server for convergence check.

Server Node

2  ×  3.2GHz processors  2 Gb Memory  Task: 

1. Load Network Data. 

2. Assign Jobs to Computing Nodes. 

3. Check Convergence. 

4. Final Estimation of ψ and F from  multiple chains after convergent.

Computing Node 1 

2  ×  3.2GHz processors  1 Gb Memory  Task: 

1. Receive Job from Server Node. 

2. Calculate Single Sequence of ψ  and F. 

3. Send Local Computation Result  to Server Node. 

MPI  Data  Bus  (Gigabits  Switch)

NFS & NIS Data Bus (Fast Ethernet) 

…

Computing Node 16 

2  ×  3.2GHz processors  1 Gb Memory  Task: 

1. Receive Job from Server Node. 

2. Calculate Single Sequence of ψ and F. 

3. Send Local Computation Result  to Server Node. 

54 

4. The server check the convergence of each ψˆ and Fˆ′ . If Rˆ>1.1, the computing

nodes will continue step 3. As Rˆ≤1.1, the server estimate the global ψˆ and Fˆ′ by sequences generated by each computing nodes.

5. Stop MPI environment. Output data.

Figure 5.5 shows the speedups and efficiencies of the proposed algorithm, where the speedups is the ratio of the code execution time on a single processor to that on multiple processors and efficiency is defined as the speedup divided by the number of processors. As shown in Figure 5.5, the parallel scheme has efficiency of 72.5% in a 32 processors environment.

Figure 5.5 Speedups and efficiencies for the parallel computing.

Chapter 6

Numerical Examples

In this section, the proposed models are tested with real networks to discuss their performance. These networks include i) Taipei Mass Rapid Transit Network in section 6.1 and ii) Taiwan Freeway Network with Electronic Toll Collection information in section 6.2.

6.1 Examples 1 (Taipei Mass Rapid Transit Network)

A real network of Taipei MRT is tested in this section. The test network consists of 9 stations in Nangkang line, the topology of test network is demonstrated in figure 6.1. Real path-flow data is not likely to derive on a real road network, but it is possible in the MRT network. Therefore, MRT information is used here to give a measurement of model performance.

Figure 6.1 The MRT test network.

Real path-flow data are the numbers of passengers traverse between stations from 18:05~20:00 in five-minute interval; for example, there are 30 persons travel

City Hall Zhongxiao

Dunhua Zhongxiao

Fuxing

Zhongxiao Xinsheng Shandao

Temple

Taipei Main Station Ximen

Longshan Temple

A ^a  B C b c  D  d ^{E F G} e f g ^{H I} h 

Sun Yat-Sun Memorial Hall

56 

from Ximen to Taipei Main Station in the time-period of 18:00~18:35. Eight paths are included in this example defined as travelers toward Taipei Main Station from the rest of stations. These paths are demonstrated in Table 6.1.

Table 6.1 Path set of example MRT network.

Path Origin-Destination Link Set

( )

1

ψ A Æ C a , b

( )

²

ψ B Æ C b

( )

³

ψ D Æ C c

( )

4

ψ E Æ C d, c

( )

5

ψ F Æ C e, d , c

( )

⁶

ψ G Æ C f, e , d , c

( )

7

ψ H Æ C g, f , e , d ,c

( )

8

ψ I Æ C h, g , f , e , d , c

In the MRT test, travel time is not considered. Therefore, link flows can be easily obtained by the summation path flows on them, i.e., flow on link c is equal to

( ) ( ) ( ) ( ) ( ) ( )

3 ψ 4 ψ 5 ψ 6 ψ 7 ψ 8

ψ + + + + + , and the flow on link e is equal to

( ) ( ) ( ) ( )

⁵ ψ ⁶ ψ ⁷ ψ ⁸

ψ + + + . The time-dependent link flows are illustrated in Table 6.2.

Link flows on b and c are treated as the observation vector, y , in this test. Both _t time-invariant and time-varying coefficient dynamic model proposed in chapter 3 and 4 are conducted in this test.

Table 6.2 Time-dependent link flows of MRT network

6.1.1 Time-invariant coefficient model

In the time-invariant coefficient model, transition matrix is assumed to be fixed.

The estimated path flows compare to real path flows are demonstrated in Figure 6.2.

The black line indicates the real O-D data derived from the Taipei Mass Rapid Transit company, and the grey line indicates the estimation result of the proposed algorithm.

Time Interval

Link flow a

Link flow b

Link flow c

Link flow d

Link flow e

Link flow f

Link Flow g

Link Flow h

1 1 12 31 27 27 17 13 7

2 3 14 116 108 103 50 30 23

3 0 12 85 84 77 39 14 9

4 6 15 82 77 72 33 21 11

5 5 25 74 69 63 42 16 10

6 2 10 71 69 65 34 17 11

7 3 33 115 108 101 56 28 12

8 3 26 62 58 56 29 12 7

9 6 29 62 59 58 34 23 15

10 1 28 84 81 71 45 21 10

11 3 19 80 72 65 28 15 7

12 3 29 113 107 103 63 31 20

13 7 35 129 119 110 51 25 17

14 2 18 50 45 44 18 7 7

15 6 39 128 121 113 58 20 12

16 6 35 97 91 83 43 15 11

17 1 12 107 100 93 58 35 25

18 3 16 108 101 97 50 25 15

19 0 22 76 76 71 32 19 8

20 4 29 69 66 63 26 13 7

21 1 17 101 100 96 63 18 12

22 2 9 41 40 37 29 10 5

23 2 6 59 57 55 35 14 8

58 

Figure 6.2 The comparison of real and estimated data on MRT network by time-invariant coefficient model.

The correlation between real and estimated data is 0.816, which indicates a medium to strong association between real and estimated data. The mean absolute error (MAE) is 4.78, and the mean absolute percentage error (MAPE) is 49%.

To evaluate the prediction result of this model, observations at the last time-period is deducted from the observation vector. After applying the proposed model to the rest of observation vector, estimated transition matrix, Fˆ , and state vector at time-period 22, ψ , are obtained. The predicted path flows at time-period ˆ₂₂ 23 can be calculated by ψˆ₂₃ =Fˆψˆ₂₂. The mean absolute error of the prediction is 5.43; and the mean absolute percentage error is 51.3%. The comparisons of real and predicted path flows are demonstrated in figure 6.3. The predicted path flows are indicated as dash lines.

60 

Figure 6.3 The comparison of real and predicted data on MRT network by time-invariant coefficient model.

6.1.2Time-varying coefficient model

This subsection demonstrates the result of proposed time-varying coefficient model. The path flow estimation is based on a rolling horizon concept; in this example, there are nine states in a horizon period, is illustrated on table 6.3.

To evaluate the estimation and prediction result of this model, same procedures are conducted as that of time-invariant coefficient model. The correlation between real and estimated data is 0.969, which indicates a strong association between real and estimated data. The mean absolute error is 1.57, and the mean absolute percentage error is 15.1%. Six out of eight path flows, without path 5 and 8, pass the chi-square test of 95% confidence interval. From the solution framework, we know the first time-period of estimation tends to have larger error, because the ψ₀ is generated randomly. If the estimation of path flows in first time-period are deducted from the dataset, all of the path flows pass the chi-square test of 95% confidence interval.

The predicted path flows at time-period 23 can be calculated by ψˆ₂₃ =Fˆ₂₂ψˆ₂₂. The mean absolute error of the prediction is 3.54; and the mean absolute percentage error is 72.2%. The comparison of real, estimated and predicted path flows are demonstrated in figure 6.4.

Time 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 Real Path 11 11 12 9 20 8 30 23 23 27 16 26 28 16 33 29 11 13 22 25 16 7 4

Estimated Path

10.5 15.0 14.9 21.1 9.2 26.3 27.1 23.6 27.8 17.0 24.0 29.9 20.9 32.2 27.6 11.3 12.5 21.5 23.6 14.8 5.6 0

Period 1 10.5 13.4 9.2 20.5 8.9 31.1 24.8 23.7 28.2

Period 2 16.6 20.0 23.6 13.0 22.0 29.9 20.8 34.6 16.2

Period 3 15.5 19.2 9.9 23.3 26.2 18.0 31.7 13.2 28.0

Period 4 21.2 4.9 24.8 19.1 22.8 23.2 15.6 22.9 28.2

Period 5 9.2 26.7 34.5 19.0 28.1 18.5 17.3 31.3 29.5

Period 6 29.7 31.2 18.5 28.6 17.7 21.2 28.6 22.5 38.1

Period 7 23.9 32.9 24.4 16.6 23.7 29.8 23.0 28.5 25.9

Period 8 32.8 24.2 16.4 23.3 29.6 23.0 28.4 25.4 12.1

Period 9 27.3 20.2 26.7 37.2 26.3 34.1 27.1 14.7 12.5

Period 10 18.5 28.1 30.3 18.1 35.2 30.5 14.3 15.0 25.2

Period 11 24.9 26.9 14.7 31.9 27.6 10.5 11.8 21.2 23.9

Period 12 26.9 14.7 31.7 27.3 10.2 11.5 21.0 23.8 14.4

Period 13 16.8 33.4 29.8 11.0 13.6 22.2 25.4 16.6 7.5

Period 14 28.9 27.4 6.4 10.3 17.9 21.1 13.4 3.8 0

62 

Figure 6.4 The comparison of real, estimated, and predicted path flows on MRT network by time-varying coefficient model with 9 states in a horizon period.

6.1.3 Discussion of MRT network example

From the numerical example of MRT test network, we could discover a significant improvement comparing the time-varying coefficient model to the time-invariant coefficient mode. As for the estimation, the time-invariant model has a MAE of 4.78 and MAPE of 49%; while the time-varying model has a MAE 1.57 of and MAPE of 15.1%. As for the prediction, the time-varying coefficient model also has smaller mean absolute error. However, the mean absolute percentage error of time-varying coefficient model is larger.

6.2 Examples 2 (Taiwan Freeway Network)

Models considering travel time effect are demonstrated in this section. The test network is a real freeway network with Electronic Toll Collection (ETC) information.

ETC can provide where and when the particular vehicle appeared, this information can help us tracing a particular vehicle. The test network consists of northern part of Taiwan freeway network, including part of freeway No.1, No. 2, and No. 3. The network is illustrated as figure 6.5.

64 

Figure 6.5 The freeway test network.

The test network consists of 22 O-D pairs and 30 paths, demonstrated in Table 6.4. Path flow data are the number of vehicles depart from their origin toll station to their destination toll station, only vehicles with ETC equipment and travel longer than two toll stations will be count. Since ETC information provides when and where a particular vehicle appears, travel time between every two stations can be calculated.

Therefore, link flow on certain time can be obtained by the summation of corresponding path flows departs from their origin on a particular time. For example, flow on link e at time t can be calculated by the summation of path 2, 4, 6, 8, 10 at time t−2 and path 12, 14, 16, 18, 20 at time t−1.

Taishan Toll Station

Yangmei Toll Station

Zaociao Toll Station

Houli Toll Station Shulin Toll

Station

Longtan Toll Station

Houlong Toll Station

Dajia Toll Station

A

C

F  D B

H E 

G

a  c  b 

e  f 

g h 

i  j 

k  l 

I  J 

K d 

Table 6.4 Path set of example freeway network.

Path Origin-Destination Link Set

( )

1

66 

6.2.1 Time-invariant coefficient model considering travel time effect

Time-invariant coefficient model considering travel time effect is discussed in this subsection. In this subsection, transition matrix is assumed to be fixed, while path-link incidence matrix, H , is generated by real data. The observation vector is _t time-dependent link flows on the network, aggregated as 30-minute interval. The estimated path flows compare to real path flows are demonstrated in Figure 6.6 as follows.

0

68 

0

Figure 6.6 The comparison of real and estimated data on freeway network by time-invariant coefficient model considering travel time effect.

It is observed from the data, there exists nearly zero flows on path 4, 6, 8, 10, 13, 15, 17, and 19. These paths are excluded from the MAE and MAPE calculation, therefore, only 22 paths are considered in the calculation.

The correlation between real and estimated data is 0.982, which indicates a strong association between real and estimated data. The mean absolute error is 2.48, and the mean absolute percentage error is 31.6%. Eight out of 22 paths pass the paired-sample T test of 90% confidence interval; we further test the rest 14 paths, 10 out of them pass the chi-square test of 90% interval.

70 

6.2.2 Time-varying coefficient model considering travel time effect

This subsection demonstrated the result of proposed time-varying coefficient model considering travel time effect on the test freeway network. The estimation of path flow is also based on a rolling horizon concept as the MRT network, except of 33 states in a horizon period. The estimated path flows compare to real path flows are illustrated in figure 6.7.

72 

Figure 6.7 The comparison of real and estimated data on freeway network by time-varying coefficient model considering travel time effect.

Same as section 6.2.1, path 4, 6, 8, 10, 13, 15, 17, and 19 are excluded from the MAE and MAPE calculation with their nearly zero flows. The correlation between real and estimated data is 0.989, which indicates a strong association between real and estimated data. The mean absolute error is 1.66, and the mean absolute percentage error is 13.29%. Nine out of 22 paths pass the paired-sample T test of 90% confidence interval; we further test the rest 13 paths, 10 out of them pass the chi-square test of 90% confidence interval.

In the numerical example of freeway test network, the time-varying model also has a better performance compared to time-invariant model in estimation. However, the difference between these models is not as significant as that of MRT example.

Some of the paths in this numerical example have zero or very little flows during the

74 

whole estimation period, i.e. path 4, 6, 8, 10 and so on. Neither time-invariant nor time-varying model are able to describe those paths.

6.2.3 Comparison with Predetermined Transition Matrix Method

Existing researches concerning time-dependent O-D estimation usually assume the prior information of the O-D matrix (or transition matrix) is know (or at least

Existing researches concerning time-dependent O-D estimation usually assume the prior information of the O-D matrix (or transition matrix) is know (or at least

在文檔中 以貝氏方法構建與求解路徑基礎之時變交通指派模式 (頁 52-0)