Inverse Analysis of Origin-Destination matrix for Microscopic Traffic Simulator

K.Abe,H.Fujii and S.Yoshimura

1 Introduction

Microscopic traffic simulations are useful for solving various traffic-related problems,e.g.traffic jams and accidents,local and global environmental and energy problems,maintenance of mobility in aging societies,and evacuation planning for natural as well asman-made disasters.To use such microscopic traffic simulators,we need to input various types of traffic data.Data is typically input as an origindestination(OD)matrix,which describes demands between origin-destination pairs in a traffic network,and is particularly necessary in microscopic simulators.Since the OD matrix cannot beobserved directly,it has to beestimated in someway.

The approaches for OD matrix estimation can be roughly classified into two categories.The firstis based on the populationdistribution.This approach iscommonly used for traffic and civil planning using the four step model[Mc Nally(2008)].Since the population distribution isderived from traffic censusdata,wedo not have to measure the actual traffic fl ow.However,the resolution of the censusdataislow,and the accuracy of the estimated resultsisnot guaranteed.

The second approachisan inverseanalysis of the link traffic volumedata.Link traffic volumeis the traffic volumecounted atafixed location.This approachisusually more accurate.Here,the OD matrix is optimized by minimizing the distance between the observed and estimated link traffic volume.This can be accomplished with a bi-level programming approach[Bera and Rao(2011)].

The inverse analysis approaches can be classified into two categories depending on the method used to solve the direct problem involved in the estimation process.Here,the direct problem refers to the assignment from the OD matrix to the link traffic volume.The solution of the direct problem can be obtained analytically,or approximated with the equilibrium assignment algorithm[Larss on and Patriksson(1992)].Estimation of the OD matrix using this assignment algorithm was shown to be effective for large-scale road networks[Lundgren and Peterson(2008)].For dynamic estimates,Off-line time-sliced OD estimation based on dynamic equilibrium,which is similar to the static OD estimation method,has been proposed[Barceló and Montero(2015)].However,owing to the differencesin network handling,the results of the equilibriumal gorithmmay notbecompatible with the traffic simulator.The alternativeis to usethetraffic simulator to solvethedirect problem.However,this latter method requires a high computational cost.Few studies have been conducted on the use of the traffic simulator for solving the direct problem,so the stability and robustness of the subsequent resultsare unclear.

In this paper,we newly propose an OD matrix estimation method using a microscopic traffic simulator.We then examine the accuracy and stability of our results for the inverse analysis.

2 Method

2.1 Outline of method

An outline of our method isshown in Figure1.

The proposed method consists of the following steps:(1)calculating link traffic volume from the OD matrix,and(2)updating the optimal OD matrix solution.In step(1),a multi-agent based microscopic traffic simulator“ADVENTURE_Mates(MATES)”[Yoshimura(2006);Fujii,Yoshimura,and Seki(2010)]is used.In step(2),we use the Levenberg-Marquardt gradient method.After cal-culating link traffic volume,the residuals of the link traffic volume are calculated.The optimal estimated OD matrix isupdated using the seresiduals.The sestepsare iterated until the estimated OD matrix converges to a tolerance level.

Figure1:Flow diagram of theproposed method for estimating the OD matrix

2.2 Formulation

The OD matrix is estimated so as to minimize the residual norm of the link traffic volume between the estimated and the observed values.This process can be expressed asfollows:

a Tnhde traffic simulation is used to evaluate QQQˆcorresponds to the traffic simulation,where link traffic volumeiscalculated from the OD matrix xxx.

We now define the dimensions of the vectors. N denotes the dimension of xxx,i.e.the number of OD pairs,and M is the dimension of QQQ,i.e.the number of observation points.In general,N＞M.Therefore,in the OD matrix estimation,we have to estimate more variables from smaller datasets.

xxx : OD matrix(asaform of vector)

rrr : Residual(vector)of link traffic volume

In this paper,we use Euclid norms.If the variance and covariance of the observed link traffic volumeareknown and we use the Mahalanobis distance to calculate the norms,we can regularize the residuals,with consideration of the variance and covariance.However,this requires multiple observation of the traffic volume,which takes much observation cost.Therefore we assume that there is no variance and covariance in the observed link traffic volume data,and we employ Euclid norms.For simplicity,we assume that the traffic flow is in a steady state and the relation between the OD matrix and the link traffic volume is linear.The assumption of such a linear relation is appropriate if traffic congestion does not occur.Here,the linear relation means that we can describe the relation in Equation 3 using amatrix

In this study,we usetheobserved link traffic volumedata which satisfies Equation 3.The assumption of a steady state means that both the OD matrix and the link traffic volumearestationary in time.This assumption isreasonableif an appropriate time period is used,e.g.a short time period where a drastic change in traffic fl ow doesnot occur.

J in Equation 3 corresponds to the Jacobian matrix of F with respect to xxx,which is an M×N matrix.Since this system has no explicit constitutive equation,we have to approximate JJJ.The value of JJJ isestimated asin Equation 4 using the results of the traffic simulator.

Sinceachangein xxx corresponds to that in the routechoice,were-evaluate JJJ at each iteration.

2.3 Modification of simulator

MATES is a multi-agent microscopic traffic simulator,which models individual driver behavior.This isuseful for modeling individual vehicleson the microscopic scale and for extrapolating traffic flow on the macroscopic scale in a form of the emergence phenomena.Since MATESuses OD matrix as input data and outputs link traffic volume,we use this simulator to calculate link traffic volumefrom OD matrix.

In order to model individual driver behavior,stochastic elements are necessary,so MATES simulates random numbers using some random seeds.In this study,wefix the serandom seedsto make MATES deterministic for simplicity.

In addition,the original version of generates vehicles according to a Poissondistribution whose rate parameter corresponds to the OD traffic volume,i.e.the element of the OD matrix.However,when using a Poisson distribution,the number of generated vehicles in fixed interval of time is not always equal to that described in the OD matrix,which causes stochastic error in the OD estimation.Here,we modify the method to generate vehicles at fixed time intervals corresponding to the OD traffic volume.The number of generated vehicles thus is guaranteed to be exactly equal to the OD traffic volume.

2.4 Solution updating step

The OD matrix is updated using one of the gradient methods.This updatng step is given in the following general form:

where

其一，加大种植业保险保费补贴政策宣传力度，提高农户对种植业保险政策的认知度。互联网等信息技术在农村的应用与推广为种植业保险的宣传提供了新的载体，为宣传工作的进一步开展提供了便利。保险公司和政府相关部门应在原有传统宣传方式的基础上采用更多形式的宣传手段，如建立微信公众平台、与农业保险相关的手机移动客户端等，拓宽农户的认知渠道，让农户足不出户就可以了解最新的农业保险政策。另外，为了达到宣传的视听效果，相关机构还应制作拍摄与种植业保险相关的微视频、微广告进行宣传教育；为了发挥宣传的联动效应，在普遍宣传中要对目标群体特别是村干部和新型农业经营主体等有针对性地进行走访式宣传。

In this study,weuse the Levenberg-Marquardt method(LMM)[Levenberg(1944);Marquardt(1963)].The definition ofΔxxx in the LMMisgiven asfollows:

where

GGG : Gradient of F with respect to xxx;GGG=−JJJTrrr

HHH : Hessian matrix of F with respect to xxx.

The steepest decent method(SDM)is commonly used in optimization problems,because it only requires a gradient matrix as input.However,the rate of convergence is slower than other methods.On the other hand,the LMM is based on the Gauss-Newton method(GNM)and is applicable to nonlinear problems.In this method,the Hessian matrix is assumed to be the product of Jacobian matrices:HHH=JJJTJJJ.However the GNM can only beapplied when N＜M.In the LMM,the Hessian matrix isregularized asin Equation 7

where λ denotes a regularization parameter for LMM and DDD denotes a diagonal matrix,e.g.,the identity matrix.The addition of the diagonal matrix makes the Hessian matrix non-singular,so that its inverse can be calculated.Marquardt proposed initially settingλ to alargevalue and decreasingλ with successive iterations[Marquardt(1963)].Ifλbecomes zero,Equation 7 is equivalent to the updating scheme for the GNM.Ifλ is sufficiently large,the search direction approaches that of the SDM.According to Marquardt’sproposal,the initial iteration stepshave the robustness of SDM,and the later iteration ones have the high rate of convergence of the GNM.In this study,evaluating Fiscomputationally expensive.Thus,we decreaseλ from an initial valueλ0(＞0)at constant rateγ(0＜γ＜1)with successive iteration sas follows:

The convergence is judged comparing relative residual norm(RRN)‖rrr‖/‖¯QQQ‖and a to leranceε.Since RRN is regularized by observed link traffic volume¯QQQ,it can be compared among different networks.

2.5 Introduction of non-negativeconstraintsin to the estimation method

Non-negative constraint of traffic volume is necessary for OD estimation.When the OD matrix satisfies the constraint,the link traffic volume then satisfies it.Thus it is sufficient to apply the constraint only to OD matrix.Since the original LMM has no constraint for variables,we have to in troduceit into the LMM.We propose the following two kinds of methods to do so.

First,we describe the common assumption in both the methods.In this study,we makethenon-negativeconstraint stricter,i.e.for all elements i,xi＞δ,whereδisa small value and is larger than zero.The reason is that too small OD traffic volume causes some estimated link traffic volume to be zero.If it occurs,some elements of Jacobian matrix become zero,then the zero elements are propagated to some elements ofΔxxx.Consequently,the search direction getslimited.

Next,we describe the methods in sequence.The first,named method A,is designed to satisfy the constraint strictly.Here,the initial value ofαis set to 1.If＜δ for any elementi,αis modified to be so smaller as to satisfy=δfor all elements i.It does not modify the search direction from that given in the LMM.Nevertheless,it is expected to that residuals stop decreasing within insufficient large values becauseαbecomes smaller toward zero with successiveiterations.

The second,named method B,is the method using heuristics.Here,the search length coeffi cientαis fixed to be 1 strictly.After the solution updating step,the updated xiis forcibly modified to beδ individually when xi＜δ for any index i.Although it makes different search direction from that in the LMM,sinceαdoes not becomezero,iteration isexpected not to halt except the convergence.

We discuss these characteristics of the two methodsin Section 4.

3 Numerical experiments

3.1 Outline of experiments

Before using the proposed method we describe previously, we have to determine thefollowing two things; the optimal tolerance for convergence ε and which methodto be used for the non-negative constraint. To determine them, it is necessary to usesome evaluation index which can be discussed in terms of traffic engineering, e.g.the proportion of reproduced link traffic volume to observed one and the correlationcoefficient between reproduced and observed link traffic volume. Here, we applythe proposed method to some cases considering observed link traffic volume withand without noise. Using these results, we determine these two things and discussaccuracy, stability and characteristics of each case or method.

3.2 Networksused for experiments

In the numerical experiments, we use two types of road networks, i.e. grids or roadmaps. The first map is a regular grid with simple topology, shown in Figure 2(a). Inthe figure, each line corresponds to a link. Numerical instabilities may arise whensolving symmetric maps. Therefore, we slightly displace all grid nodes in a randommanner. In this map, each link has 2 lanes, and all intersections have traffic signals.

The second map is an actual road network topology,i.e.a 3 km x 3 km area in Tokyo.Themap isshown in Figure2(b).Thismap includes one way links,so that thenumber of ODpairs N islessthan n(n−1),where n denotesthenumber of OD nodes.

The properties of these networks are shown in Table 1. In the table, detectors are thepoints at which link traffic volume is measured in a simulation. As a general rule inthis study, detectors are located on all links, 5 m from the end point. Consequently,there are twice as many detectors as links.

Figure2:Networksused for experiments

Table1:Network properties

3.3 Initial conditionsfor experiments

In the numerical experiments, we use two types of road networks, i.e. grids or roadmaps. The first map is a regular grid with simple topology, shown in Figure 2(a). Inthe figure, each line corresponds to a link. Numerical instabilities may arise whensolving symmetric maps. Therefore, we slightly displace all grid nodes in a randommanner. In this map, each link has 2 lanes, and all intersections have traffic signals.

Table2:Initial conditionsin theexperiments

Next,we consider the observed link traffic volume.In this study,we useasi mulated dataset in advance in all experiments.To generate this link traffic volume,we run the traffic simulator using OD matrices which satisfy the assumption of linearity described in Section 2.2.The results are given in Table 1.Since the data has been simulated under the seconditions,the existence of the solution isguaranteed.

Using this data as a reference case, we can consider other cases where artificialnoise has been added to the data in Table 1, to examine the accuracy and stability of the proposal method. Noise is added to the networks as follows.

Here, δ is a noise vector, where each element i follows a uniform distribution whoseupper and lower limits are within ±10% of ¯Qi, For the experiments with each map,10 kinds of observed link traffic volume are prepared using 10 different randomseeds.

4 Results and Discussions

4.1 Experiments on the regular grid

First of all,we apply the proposed method to the regular grid.As mentioned in Section 3.1,we apply both of the two methods of the non-negative constraint,i.e.methods A and B.

We apply the methods to the cases with and with outnoise.RRN transi tionsin these cases are shown in Figures 3 and 4.The x-axis denotes the number of iteration counts,while the y-axis is the RRN on the logarithmic scale.Red and green lines denote the methods A and B,respectively.

Figure 3:RRN transition in the regular grid without noise

The sefigures show that RRN transitions in both the methods A and B decrease in the same rate until early iteration steps,i.e.6th iteration step without noise and 2nd iteration step with noise.The reason is that all elements of OD matrices are positive until those steps.Since the methods A and B consider both non-negativeconstraints,there is no difference unless the constraints are applied to the solution updating step.On the other hand,RRNs in method B shows different transition between the cases with and without noise.The RRN in the case without noise oscillates,while the RRN in the case with noise almost continuously increases after the RRN becomes minimum.This difference is caused by the difference of the number of the elements to which the non-negative constraint is applied in OD matrices.If the non-negative constraint is applied for an element in OD matrices,its value approachesδ.Thus,the number of the elements of xi=1 denotes the number of the elements to which the constraint is applied.Here,we show the resultsin Figure5.

Figure4:RRN transition in the regular grid with noise

This figure shows that there isadifferencein the number of the elements of xi=1 between the cases with and without noise.In the case without noise,the number oscillates in small,i.e.from 0 to 3,while in the case with noise,the number increases almost continuously.As a result,the case with noise has more elements to which the constraint is applied than the case without noise.This suggests that the situation of the regular grid without noise has much simplicity and stability in it erative solution search.

To discuss the reasonable tolerance for convergenceε and the characteristics in estimated link traffic volume,we confirm the estimated link traffic volume using the minimum RRN in the following cases:method A without noise,method B without noise and method B with noise.These are considered with the difference of scale of minimum RRN among the cases.All the results are shown in Figure6.The x-axis denotes observed link traffic volume¯QQQ,while the y-axis is the link traffic volume estimated from the OD matrixˆQQQ.The approximate straight line is fixed at the origin point.

Figure5:Transition of the number of theelementsof xi=1 in theregular grid

These figures of Estimated link traffic volume shows following two evaluation indexes:gradient of approximate straight line a and correlation coefficient R. Thegradient of approximate straight line a shows the scale of estimated link traffic volume.If a equals to 1, it suggests that the scales of observed and the estimated linktraffic volume are same. The scale of the estimated link traffic volume is higherthan the observed one if a is greater than 1, otherwise the scale is lower. The correlationcoefficient R shows the dispersion of accuracy of reproduction. If R equals to1, the following relation is satisfied: ˆQQ = a¯QQ. On the other hand, if some elementsof link traffic volume are plotted away from approximate straight line, R becomeslower.

In terms of these evaluation indexes, the values of a in all of these cases equal to 1in error by at most 0.006. The values of R in all cases also equal to 1 in error byat most 0.004. This high correlations arise from the characteristic that the LMMweights residuals equally. These results suggest the accuracy in these cases are thesame or higher than 99% in proportion of the estimated link traffic volume to theobserved one. Since these errors are sufficiently small to ignore and the highestRRN in these cases is 2.96×10−2 at that time, it is sufficient to set convergencecoefficient ε to 0.03.

Figure6:Estimated link traffic volume for the regular grid

4.2 Application to an actual road network

Next,we apply the proposed method to an actual road network,i.e.a3 km x 3 km areain Tokyo.We show the RRN transitionsin the cases with and without noisein Figures7 and 8.

In this map,RRN transitions in both the two cases are alike.Until early several iteration steps,RRNs in method B are smaller than those in method A,while after that,RRNs in method B become larger than those in method A.The increase of RRNs in method B is explained by the same reason mentioned in Section 4.1,i.e.the increase of the elements of xi=1 in OD matrices.Actually,in this map,the transitions of the number of the elements of xi=1 are similar to the green line in Figure5 in both of the cases with and without noise.We show it in Figure9.

Figure7:RRN transition for the3 km x 3 km areain Tokyo without noise

Figure8:RRN transition for the3 km x 3km areain Tokyo with noise

On the other hand,the seresults are different from those in the regular grid,where RRNs are the samein method A and Batearly several iteration steps.This suggests that in method A,the search lengthαapproaches zero at early few iteration steps and there fore the RRNsdecreaseinsuffi ciently.Actually,unlikein the regular grid,the non-negative constraint is applied to OD matrices at only the first iteration in Tokyo map.This is caused by the large number of variables in OD matrices in Tokyo map.

Figure 9:Transition of the number of the elements of xi=1 for the 3 km x 3 km areain Tokyo

For method B, the iteration steps when RRNs become minimum are alike as follows:3rd iteration steps in the case without noise and 2nd iteration steps in the casewith noise. These values are almost the same as that in the regular grid with noise,whose value is 2, and these are all small number. Therefore the results suggest thateven if not adding noise, solution search is more difficult in Tokyo map than in theregular grid, which is also caused by the large number of variables.

In terms of minimum RRNs, as similar as the regular grid, there is some differencebetween in the cases with and without noise. Here, we show the estimated linktraffic volume when RRN is minimum in the following cases, considering the scaleof minimum RRN: method A without noise, method B without noise, method Awith noise and method B with noise. All the results are shown in Figure 10.

These figures show that the gradient of approximate straight line a is distant from 1when method A is used, i.e. 0.852 without noise and 0.618 with noise. On the otherhand, when method B is used, a equals to 1 in error by at most 0.022 as like as inthe regular grid. The correlation coefficient R almost equals to 1 for all constraintmethods and cases about noise. These results suggest method B is more applicable.Then, the accuracy is the same or higher than 97% in proportion of estimated linktraffic volume to observed one. Even though this value is less than that in the regular grid, this is expected to be still sufficiently accurate. Consequently, it is sufficientto set the tolerance for convergence ε to 0.8, where only method B converges.

Figure 10:Estimated link traffic volume for the 3 km x 3 km areain Tokyo

5 Conclusions

We newly proposed an OD estimation method using a traffic simulator,whose result is designed to be suitable for use of the simulator directly.In addition,we introduced two kinds of approaches of applying non-negative constraints to the proposed method.Through numerical experiments,we demonstrated the validity of the method.Thees timated link traffic volumeisstrongly correlated with the ob-served link traffic volume.This is due to the iterative process in the LMM.When we introduce the non-negative constraint into our method,the one using heuristics is found to beeff ectiveand givessmaller RRN.Then,the tolerance for convergence for the RRN is found to set 0.08,when the estimated link traffic volume has 97%or more accuracy to the observed one.To improve the accuracy,especially even if the number of observation points becomes smaller,it is necessary to improve the method of non-negativeconstraint in terms of heuristics.Infuturework,weplan to examine the accuracy and stability of the proposed method with fewer dataset and to consider the stateof congestion.

Acknowledgement: This work wassupported by JSPSKAKENHI Grant Number 15H01785.

Barceló,J.;Montero,L.(2015): A Robust Framework for the Estimation of Dynamic OD Trip Matricesfor Reliable Traffi c Management.Transportation Research Procedia,vol.10,pp.134–144.

Bera,S.;Rao,K.(2011):Estimation of Origin-Destination Matrix from Traffi c Counts:the State of the Art.European Transport-Trasporti Europei,vol.49,pp.3–23.

Fujii,H.;Yoshimura,S.;Seki,K.(2010):Multi-agent Based Traffi c Simulation at Merging Section Using Coordinative Behavior Model.Computer Modeling in Engineering and Sciences,vol.63,no.3,pp.265–282.

Larsson,T.;Patriksson,M.(1992): Simplicial Decomposition with Disaggregated Representationfor the Traffic Assignment Problem.Transportation Science,vol.26,no.1,pp.4–17.

Levenberg,K.(1944):A Method for the Solution of Certain Non-Linear Problemsin Least Squares.The Quarterly of Applied Mathematics,vol.2,pp.164–168.

Lundgren,J.T.;Peterson,A.(2008): A Heuristic for the Bilevel Origin-Destination-Matrix Estimation Problem.Transportation Research Part B:Methodological,vol.42,pp.339–354.

Marquardt,D.W.(1963):An Algorithm for Least-Squares Estimation of Nonlinear Parameters.Journal of the Society for Industrial and Applied Mathematics,vol.11,no.2,pp.431–441.

McNally,M.G.(2008): The Four Step Model.Center for Activity Systems Analysis.

Yoshimura,S.(2006): MATES:Multi-Agent Based Traffi c and Environment Simulator-Theory,Implementationand Practical Application.Computer Modeling in Engineering and Sciences,vol.11,no.1,pp.17–25.