• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Efficient Origin-Destination Estimation Using Microscopic Traffic Simulation with Restricted Rerouting

    2023-02-26 10:17:00KazukiAbeHidekiFujiiandShinobuYoshimura

    Kazuki Abe,Hideki Fujii and Shinobu Yoshimura

    1Vector Research Institute,Inc.,Shibuya-ku,Tokyo,150-0002,Japan

    2School of Engineering,The University of Tokyo,Bunkyo-ku,Tokyo,113-8656,Japan

    ABSTRACT Traffic simulators are utilized to solve a variety of traffic-related problems.For such simulators,origin-destination(OD) traffic volumes as mobility demands are required to input,and we need to estimate them.The authors regard an OD estimation as a bi-level programming problem,and apply a microscopic traffic simulation model to it.However,the simulation trials can be computationally expensive if full dynamic rerouting is allowed,when employing multi-agent-based models in the estimation process.This paper proposes an efficient OD estimation method using a multi-agent-based simulator with restricted dynamic rerouting to reduce the computational load.Even though,in the case of large traffic demand,the restriction on dynamic rerouting can result in heavier congestion.The authors resolve this problem by introducing constraints of the bi-level programming problem depending on link congestion.Test results show that the accuracy of the link traffic volume reproduced with the proposed method is virtually identical to that of existing methods but that the proposed method is more computationally efficient in a wide-range or high-demand context.

    KEYWORDS OD estimation;microscopic traffic simulation;dynamic rerouting;bi-level programming;multi-agent based model

    1 Introduction

    Traffic simulation is used to solve various traffic problems and inform associated policies such as handling congestion or accidents,controlling CO2emissions,and guaranteeing mobility in aging societies.Many simulation models have been developed,and they can be roughly divided into macroscopic,mesoscopic,and microscopic models.The macroscopic models approximate aggregated vehicles as a continuous flow on links,where some form of traffic assignment models assign each link demand.The microscopic models treat each vehicle as a particle,making it possible to distinguish individual vehicles and thus to describe vehicles’precise behaviors.In contrast to macroscopic models,each vehicle has its own route along which it travels in some microscopic models.The mesoscopic models are just the middle of the macroscopic models and the microscopic models.While each vehicle is identified like microscopic models,the behavior is expressed in a macroscopic manner.

    The simulation model and input data must be reliable and realistic to reproduce a realistic traffic situation in a simulator.And in most cases,a single simulator is used for parametric studies by changing some parameters sets.In this research,the authors focus on the reliability of traffic demand as the input data given to the simulator.

    We can easily obtain road structure information(e.g.,intersection configuration)and observable and controllable information (e.g.,traffic light patterns) in the real world.However,traffic demand data,which directly affect link traffic volume,congestion,and service levels,are more difficult to obtain because of their unobservability and significant uncertainties.

    Origin-destination (OD) matrices describe traffic demand.An OD matrix consists of traffic demand(traffic volume)for each OD pair,i.e.,a pair of an origin and a destination point.Rows and columns of an OD matrix denote origins and destinations,respectively.For example,the(r,c)element in an OD matrix represents the number of vehicles per unit time that flow into the road network at origin pointOrand flow out at destinationDc.The unit of traffic demand is often given in vehicles per hour,veh/h in short.

    The authors regard OD estimation as a bi-level programming problem and target OD estimation problem for microscopic traffic simulation,especially.In this situation,the same traffic model should be employed in both lower and upper problems to make the input and the output for those problems consistent.

    It is concerned that the computational time increases due to the dynamic rerouting in each simulation trial and the increased number of trial iterations when using such a method.Dynamic rerouting is a route search operation required when a vehicle driver faces a high-cost link to pass(e.g.,a congested link)during driving along the route acquired by the initial search.Thus,dynamic rerouting is reasonable because it corresponds to each driver’s detour behavior.By introducing dynamic rerouting,vehicles detour to avoid congestion,and consequently,congestions are relaxed.However,each route search operation,including dynamic rerouting,is computationally expensive.In addition,each vehicle driver should know the traffic situation of the route candidates at a specific time point to reroute dynamically.It is concerned to be an unrealistic assumption beyond the vehicle drivers’capability.The authors propose a way to suppress the number of dynamic rerouting in OD estimation and deal with the problems resulting from this restriction.

    2 Literature Review

    Three general approaches are available for OD estimation:statistical approaches,traffic volume optimization approaches with demand assignment,and traffic volume optimization approaches with microscopic simulators.In the statistical approaches,the OD matrices are generated by aggregating survey data,which can come from personal interviews.A typical example is the four-step model[1],widely utilized for long-term demand prediction since the 1950s.

    The traffic volume optimization approaches with demand assignment optimize link traffic volume,i.e.,match the volume in the simulation to that in the real world by changing the solution candidates for the OD matrices.The assignment of traffic demand for each link is determined under the assumption of logit-based decision rules and user equilibrium (UE) theory.Most such methods handle traffic flow with stochastic turbulence based on minimum entropy[2–4],maximum likelihood[5,6],or generalized least squares optimization[7].For further details,see[8].

    The statistical and traffic volume optimization approaches with demand assignment handle wideranging and aggregated traffic flows.Therefore,these methods are often applied in macroscopic simulators because of the characteristics of aggregated flow.

    Microscopic simulators are also widely used.These simulators seek to capture more detailed traffic phenomena.A multi-agent system (MAS) is commonly used to implement such models.The MAS captures each agent’s autonomous recognition,decision-making,and action.When applied to traffic simulation,the MAS regards each driver or each vehicle as an agent.Entire phenomena such as congestion emerge from the interactions among the drivers.Since this type of modeling directly incorporates the individual behavior of actual drivers,the MAS is suitable for simulating traffic phenomena in the real world at a fine-scale resolution.The MAS employs route search algorithms(e.g.,Dijkstra’s algorithm [9] or the A?algorithm [10]) where each vehicle agent constructs its preferable route considering heterogeneity and the dynamic changes in the traffic situation.In recent years,many multi-agent-based microscopic simulators have been applied,including VISSIM [11],Aimsun [12],SUMO[13],ADVENTURE_Mates[14,15]and MATSim[16].

    Microscopic traffic simulations based on a MAS require OD matrices as traffic demand input data.In one traffic volume optimization approach using a MAS-based simulator for OD estimation,Nguyen et al.applied the DFROUTER tool to SUMO.This tool reconstructs vehicle routes based on detector measurements [17].Shafiei et al.proposed a gradient-based method using a linear approximation for the demand-link traffic volume relation [18].The simultaneous perturbation stochastic approximation model(SPSA)[19]has also been proposed for OD estimation in microscopic simulations,including MAS-based models.SPSA is one of the gradient-based models that minimize the absolute error in the simulated link traffic volume (i.e.,the difference between the simulated volume and the observed volume).Unlike the original gradient-based model,the gradient is calculated stochastically in the SPSA so that it is tolerant to noise.Yang et al.employed the SPSA for OD matrix optimization with VISSIM[20].Here,VISSIM is not only the simulation model itself;it is also a part of the OD matrix estimation model.Thus,the estimated OD matrices provide consistent input data for the VISSIM simulator.Djukic et al.applied the SPSA method to Aimsun,where both microscopic and macroscopic models were employed[21].There are also other OD estimation studies using SPSA,e.g.[22,23].

    Omrani and Lina proposed a genetic-algorithm (GA)-based approach for OD estimation [24].They used a genetic algorithm to calibrate parameters affecting route choice behavior in the simulation.Marzano et al.applied a Kalman-filter for OD estimation [25].Yousefikia et al.proposed the TflowFuzzy method,which is based on fuzzy set theory[26].These methods use bi-level programming;the update process for the OD matrix solution candidates and the simulation process are divided and conducted alternately.Theoretically,we can apply these methods based on bi-level programming to congested traffic situations in the real world[8].

    Because drivers’dynamic rerouting behaviors avoid the uneven distribution of vehicles in a network,bi-level-programming-based approaches employing MAS-based microscopic simulators within the model tend to work well.However,these approaches can involve significantly higher computational costs for the dynamic rerouting process,as each driver search for the shortest path frequently.For example,SPSA requires two simulation runs to determine the gradient[27].Thus,using a fine-modeled simulator for a large-scale network can be computationally quite expensive.

    Nowadays,several OD estimation methods using advanced deep learning models have been proposed[28,29].For example,Tang et al.proposed the method using the 3D convolution-based deep neural network (NN),Res3D,to estimate the OD matrices in Shanghai.Koca et al.proposed the graph NN to estimate the OD matrices in New York.These researches aim to estimate OD matrices by inputting numerous vehicle coordinate data into learned NNs.However,there are difficulties in employing advanced deep learning models to estimate OD matrices used in microscopic traffic simulations.For example,Koca et al.[29]reported that the coefficient of determination in link traffic volume is less than or equal to 0.3 in the cross-validation in which data obtained in an area were input into the model learned in another area.In other words,the versatility of NN-based methods is not sufficient.Consequently,as with other methods that do not incorporate microscopic simulators,the OD matrices estimated by the NN-based methods may be inconsistent with the simulator that intends to use them.

    3 Methodology

    3.1 Aims

    As mentioned in the previous section,this research aims to solve the OD estimation problem using a MAS-based simulator with restricted dynamic rerouting functions.Here,the simulator allows a vehicle agent to reroute only after finding itself in an unintended route,which can occur when a lane is blocked due to congestion.By restricting dynamic rerouting,it is expected that the computational cost will be reduced and that vehicle drivers’routing behavior will be more feasible.

    However,employing such a restricted MAS-based simulator faces several challenges in the OD estimation process.If high traffic demand is assigned to a low-throughput route in a certain iteration step,congestion may occur at intersections.This is a natural characteristic in real situations,and if full dynamic rerouting is allowed,the demand is adequately distributed to mitigate the congestion.But,when the restricted dynamic rerouting proposed in this paper is employed,a kind of positive feedback occurs so that heavier congestion arises at successive iteration steps.Once congestion incurs stuck of vehicles upstream of the observation point,the link traffic volume,i.e.,the number of vehicles that passed the link,is under-measured compared to the link traffic demand,i.e.,the number of vehicles that intend to pass the link.In order to resolve this shortage of link traffic volume,more vehicles are assigned to pass the link,regardless of whether it is the congested link or not.Consequently,the congestion becomes even heavier,and the traffic volume is further reduced unexpectedly.This traffic congestion generated by the restricted simulator functions is hereafter referred to as“fictitious congestion”.

    To resolve the fictitious congestion problem,the authors introduce upper-bound constraints for link demand reflecting a congested state in the OD estimation process.It is handled outside the simulator used in the OD estimation process.As a result,the proposed method is applicable to simulators even when dynamic rerouting is fully available.

    3.2 Overview

    The authors aim to solve the following two problems: (1) calculation cost for bi-levelprogramming-based OD estimation and(2)fictitious congestion when employing the simulator with restricted dynamic rerouting.Especially,the OD estimation method is not established in context of restricted dynamic rerouting.In this paper,the OD estimation problem is solved based on a framework that deals with OD traffic volume assignment to links by alternately using simulation and optimization of the solution candidate of the OD matrix via a numerical solver.Within this framework,a MASbased simulator with restricted dynamic rerouting is employed.By using this simplified framework,the computational cost can be substantially reduced.To demonstrate the practical performance of the proposed method,three cases were tested:a low-demand case,a high-demand case with link demand constraints,and a high-demand case without link demand constraints.

    A flowchart of the proposed OD estimation method is shown in Fig.1.This is based on the authors’previously proposed OD estimation method that embeds a MAS-based simulator[30,31].

    Figure 1:Flowchart of proposed OD estimation method

    Here,xdenotes the OD matrix in vector form.The dimension ofxisN,the number of OD pairs.anddenote the observed and simulated link traffic volumes,respectively.The observed link traffic volumes are given from the beginning;the simulated link traffic volumes are obtained from simulation results with the OD matrixx.(k)denotes the iteration number.Tdenotes transposition.rdenotes the residual of the observed and simulated link traffic volumes;that is,r=?.The dimensions ofr,andare allM.In general,Nis larger thanM.jdenotes the index of OD pairs.xudenotes the common upper bound of each OD traffic volume.This value is a constant regardless of the importance of the OD pair such as trunk or branch ones.

    The estimation process starts with an initial OD matrix in vector form,x(k)=x(0).In this form,each element of the vector denotes an OD traffic volume.Next,a simulation is conducted with the OD matrixx(k)to obtain the simulated link traffic volume in vector formand assignment matrixA(k),which is explained in Subsection 3.3.Convergence is then examined.If the root mean square error(RMSE),as defined in Eq.(1),in thek-th iteration step is smaller than a threshold valueεor the iteration number reaches the maximum number of iterationskmax,the solution is judged to have converged.εis often determined considering the total traffic demand or variant of that.

    If the solution does not converge,a new OD matrix candidate is calculated by solving a quadratic programming(QP)problem.When the QP is solved,some constraints corresponding to the congestion situation are considered.The Fig.1 indicates this operation as Congestion detection and Sign inversion.The detail is described in Section 3.3.These operations are iterated until a solution candidate converges.The flow of alternately conducting the simulation and solving the QP problem corresponds to a bi-level programming problem.To avoid an overly long computational time,a QP solver is used iteratively instead of the SPSA.Details of the formulation and variables are given in Subsection 3.3.

    To simplify the OD estimation problem,the authors assume that the traffic situation is nearly steady during a one-hour simulation.This assumption is based on the fact that link traffic volume is aggregated every hour in many data sources.Importantly,this OD estimation is performed offline,i.e.,all the link traffic volumes are determined before the estimation starts.By focusing on a particular onehour period to be estimated,the model can be considered static.There is no problem discussing the function of dynamic rerouting under the assumption of static conditions because the MAS itself shows a dynamic process.The behavior of an agent means that the environment is dynamically changing for other agents.

    The fictitious congestion appears in a simulation during the OD estimation process.In reality,even if congestion occurs,it will be resolved over time due to drivers’dynamic rerouting behavior.However,under the assumption of the static OD matrix,with the difficulty of detouring as explained in Section 2,the congestion will not be resolved since link demand is constant during the targeted hour,even if the link is congested.In such a situation,congestions must be avoided,i.e.,there should be no assignment of higher demand than the link traffic capacity.Here,the authors define fictitious congestion more concretely:It is the emergence of a vehicle queue that continuously exceeds a certain length in a simulation during the OD estimation.Then,whether or not fictitious congestion appears can be distinguished in terms of the simulation results.The authors seek to reduce the fictitious congestion by modifying the OD matrix rather than invoking dynamic rerouting in the simulation.

    3.3 Formulation

    In the proposed method,the residual between the observed and simulated link traffic volumes is minimized subject to constraints on OD traffic volume,as a kind of QP problem.The optimization problem with the objective function and its constraints is

    where

    To solve Eq.(2),the authors assume a relation between the OD matrix and link traffic volume as

    This assumption means that the product of matrixAand OD matrixxis equal to the simulated link traffic volume.MatrixAis an OD-link assignment matrix.ElementAijis defined as

    Here,yijthe traffic volume of vehicles that passed through linkiamong the vehicles traveled from the originOjto the destinationDjofODpairj=(Oj,Dj).Thus,elementAijis regarded as the fraction of the vehicles that selected linkifor entire vehicles inODpairj.The parameterδdenotes a minimum traffic volume threshold.Because the selection fraction of linkicannot be calculated if the OD traffic volumexjis too small,the relevant element in the assignment matrix is not updated in such case.

    Based on the relation in Eq.(4),a standard form of the QP problem is obtained as

    where

    λI is a regularization term that renders the square matrixATAnon-singular.In this study,the authors setλto 1.In Spiess’s original formulation,the objective function isf (x)=‖r‖2+λ‖?x‖2,whereis a target OD matrix that is typically obtained from a data source,such as a trip-based model[5].The authors assume no target OD matrix due to the difficulty of accessibility to the detailed OD matrix generated by another method.Then,parameterλin the proposed method is equivalent to that of Spiess,and Eq.(7)is regarded as a case where the target OD matrix is zero in Spiess’s formulation.

    The assumption in Eq.(4)is expected to be satisfied in the free-flow state.The QP problems are solved using the open-source solver qpOASES [32] in this paper.This solver employs the active set method,which is suitable for large-scale problems.

    3.4 Congestion Detection

    The model described in the previous subsection is a basic optimization model that works well only in a free-flow state for a network.As described in Section 2,the authors introduce congestion detection into the simulation during the OD estimation process and consider congestion effects.In the estimation method assuming static traffic demand,the traffic volume at a congested link becomes smaller than in a non-congested,free-flow state.In such a situation,adding traffic demand in order to increase link traffic volume leads to heavier congestion,and as a result,link traffic volume drops.

    The number of stuck vehiclessiis counted at each linkito detect congestion.Though we can also consider occupancy at each link for this purpose,the number of stuck vehicles is preferable because it is independent of the link length.The stuck vehicles are counted at a specific point in time,when sufficient time passed.A vehicle whose most recentTcminute mean velocity is lower thanvc[km/h]is judged to be stuck.These values are obtained from each vehicle’s trajectory data.The purpose of the detection is to reduce fictitious congestion in a simulation during the OD estimation process.The number of stuck vehicles corresponds to the congestion length,which can be measured by sensors at signalized intersections.

    The authors add a constraint regarding stuck vehicles to the optimization problem in Eq.(6).The new constraint is shown in Eq.(8).The constraint is applied only for links where observed real-world link traffic volume data are available.

    Here,Aidenotes an assignment column vector for linkiandAixdenotes the traffic demand calculated in the QP problem without simulation.The right term?sidenotes the acceptable traffic demand of linkiin the simulation.

    This constraint is given for a particular linkiwhere the number of stuck vehiclessiis greater than or equal to the critical value ?s.The constraint prevents the traffic volume from being suppressed unexpectedly at too many links.

    In addition,to reduce the fictitious congestion,the authors introduce sign inversion of the elements of matrixA.When linkiis judged to be congested,i.e.,siis greater than or equal to ?s,the sign ofAijfor alljis inverted.This sign inversion is a simplified strategy for handling the change in link traffic volume over time in static OD matrix contexts.

    The sign inversion is based on typical traffic flow characteristics [33].Originally,the relation between traffic volume and spatial density was configured as an inverse-lambda-shaped function,as shown in Fig.2,and the authors assume that the density reflects demand indirectly.If a link is congested,the throughput will decrease when the demand increases.This effect works as positive feedback in the OD matrix optimization and thus renders the solution candidates unstable.If the sign ofAijis inverted,it works as negative feedback,and the OD matrix solution candidates are likely to be less congested.

    Figure 2:Conceptual diagram of the relation between traffic volume and spatial density

    The BPR function,the empirical velocity-traffic volume relation,is often used for traffic assignment in traffic engineering conventionally [34].This function is also based on the relation of this diagram,i.e.,it indicates same velocity in lower traffic volume and indicates slower velocity over a critical traffic volume.The sign inversion is an analogy of such an empirical relationship.In addition,though the sign inversion is applied to QP,the basic idea of that is based on the gradient method.In the other words,since the sign which is to be inverted is an element of the gradient of the objective function,that element works as a sort of barriers in the QP.

    Since the shape of this diagram and the value of the critical density,where the relation of the volume and the density inverses,depends on the links,it is difficult to determine this relation directly.Therefore,it is important to employ congestion detection to estimate congestion situations.

    4 Numerical Experiments

    4.1 Simulator

    In the numerical experiments in this study,a MAS-based microscopic simulator,ADVENTURE_Mates[14,15],was employed.In the simulator,each vehicle travels from an origin node whose out-degree is 1 to a destination node whose in-degree is 1.Vehicles are modeled in detail: In the longitudinal direction,they proceed based on a car-following model,specifically the Intelligent Driver Model [35]; in the lateral direction,they change lanes considering the space to be occupied.They are generated stochastically at each origin node in accordance with a Poisson distribution.For this operation,the random seed is given exogenously.

    Vehicle searches the shortest travel time route,referring to the mean travel time of each link during the nearest 5 min.If a vehicle cannot keep its desired route because of some unexpected event,e.g.,a lane-change failure due to congestion,rerouting is invoked.

    The authors conducted 90-min simulations and measured link traffic volumes during the last 60-min for the evaluation.Since there were no vehicles in the initial state and the system was in a transient state for the first 30 min,the authors excluded these irregular states from consideration.

    4.2 Simulation Configuration

    The authors applied the method to a real city road network to test the proposed method.The network and observation point locations are shown in Fig.3.

    Figure 3:Target network and observation points(in red)

    The north-south and east-west ranges are both approximately 8 km in length.The network is in the central area of Fukui City,Japan.In this network,there are 5,852 possible OD pairs;however,in an effort to simplify,the authors extracted 1336 pairs(N=1336)for the estimation.The simplification method used to select the final OD pairs is explained below.

    First,the authors assigned a small OD traffic volume-2 veh/h (vehicles per hour)-to all 5,852 possible OD pairs to create the initial OD matrixx(0).This OD traffic volume was sufficient to calculate the matrixAand did not produce congestion at any of the links.Giving a small OD traffic volume as an initial condition makes it easy to calculate gradients in the solution space without congestions,where the relationship between the link traffic volumeand the OD traffic volumexis almost linear.Candidates in the solution space considering congestions,where the relationship betweenandxis nonlinear,are searched through iterations.Note that though the initial OD traffic volume is only 2 veh/h,the traffic volume of the link connected to an origin node is aboutveh/h.Because the number of possible destination nodes from one origin node is aboutNext,Eq.(6)was solved,not with a QP solver,but rather by using one of the gradient-based methods,the Levenberg-Marquardt method,for a single iteration.The gradient was calculated using

    A roughly-estimated OD matrix was then obtained by calculatingx(1)=x(0)+Δx,where negative traffic volumes were forced to zero to satisfy the constraint.

    Then,OD pairs with large OD traffic volumes were extracted.Fig.4 shows the cumulative relative frequency of the roughly-estimated OD traffic volumes obtained from this procedure.

    Figure 4:Cumulative relative frequency for OD traffic volume in roughly estimated OD matrix in range of 0–10 veh/h

    As seen in the figure,a jump in frequency occurs at 3.7–3.8 veh/h.The authors regarded 3.8 veh/h as a threshold and assumed that the OD pairs with a traffic volume smaller than the threshold were ineffective in terms of the overall traffic phenomena.Therefore,these pairs were excluded to reduce the total number of OD pairs.

    The network includes 118 observation points(M=118).The observed link traffic volume was aggregated from traffic census data for each of the observation points.The link traffic volumes were observed from 8 to 9 a.m.on Saturdays,2014,2016,and 2017.That is the peak congestion time on weekends.

    In the experiments,the criteria for the congestion detection are set as follows:Duration criterionTcis set to 5 min,which is based on the signal cycle.The velocity criterionvcis set to 1.0 km/h,which is chosen based on the numerical errors in the simulation.

    4.3 Evaluation Indices

    Among many evaluation indices,at least it is necessary to evaluate based on the objective function.Thus,the authors confirmed the reproducibility of the link traffic volumes by evaluating the correlation coefficientRand the regression coefficientabetween the observed and the simulated link traffic volumes.Here,the authors assumed that the y-intercept of the regression line was zero.The closer the values ofRandaare to 1,the better the estimated OD matrix reproduces the link traffic volume.

    In addition,the second index is used.This one indicates the degree of congestion measured as the number of stuck vehicles.Heat maps for this index are used to visually identify and evaluate the congestion range and magnitude.In this study,if the number of stuck vehicles is closer to zero,the result is recognized to be more effective in suppressing fictitious congestion.

    4.4 Experimental Cases

    The authors considered four experimental cases,as shown in Table 1.The cases are categorized into two groups: verification in free-flow conditions and performance in congested conditions.In the former,simulations were conducted where vehicle collisions were ignored,i.e.,in the country where vehicles are driven on the left,and right-turning vehicles do not yield even when there are oncoming vehicles.Implementing this modification increases traffic capacity and,as a result,congestions are suppressed.In such a situation,the problem assumption is equivalent to the authors’previously proposed methods[30,31].The authors regarded this case as a reference for comparisons with existing research.

    Table 1: Experimental cases

    The latter category is composed of three cases:Control,Weak and Strong.In all three cases,vehicle collisions should be avoided.In the Control case,congestion detection is not considered.If traffic volume is large,fictitious congestion is expected to occur.As the label suggests,this case is used as the reference case for the group.In contrast,congestion detection,as described in the previous section,is considered in both the Weak and Strong cases.The differences between the Weak and Strong cases are(1)the timing of the congestion measurement,i.e.,the threshold used to judge congestion,and(2)sign inversion.In the Strong case,the constraints are applied to more links than in the Weak case.

    5 Results

    5.1 Verification

    The proposed method was verified in the free-flow condition using five trials with different random seeds for the simulator.The results are shown in Table 2.Fig.5 shows the scatter plot for the link traffic volume.As expected,using random seeds varied the number of generated vehicles and the timing they flowed the network.The relative root mean square error(RRMSE)shown in Table 2 is calculated using Eq.(10).The numerator of the right hand denotes the RMSE.

    Table 2: Comparison of results between the proposed method and existing studies in free-flow condition

    Figure 5:Scatter plot for verification using regression line

    Table 2 also shows the estimation performance of the methods proposed by Yang et al.[20]and Omrani et al.[24].Like the method proposed in this paper,both of these existing approaches employ MAS-based microscopic traffic simulators; however,the three approaches differ in the upper-level optimization method and whether the dynamic rerouting functions are fully available.

    Under the free-flow assumption,the traffic situations reproduced by the three approaches showed no substantial differences regardless of whether the full dynamic rerouting functions were available,as congestion rarely occurred.Furthermore,all three methods produced small RRMSEs and large correlation coefficients.Namely,we can judge the proposed model has a similar estimation performance as the methods of Yang et al.and Omrani et al.

    5.2 Reproducibility of Link Traffic Volume in Congestion Conditions

    The link traffic volume reproducibilities under congestion conditions are shown in Table 3.The mean RMSEs of 10 trials with different random seeds are plotted in Fig.6 for the Control,Weak and Strong cases.The error bars indicate the 30–70th percentiles.In these tests,the authors set the convergence thresholdεto zero to investigate the transition of RMSE and examined 9 iterations for each trial.

    Table 3: Comparison of results between the proposed method and existing studies in congested condition

    Figure 6:Change in RMSE for control,weak and strong cases

    In all cases,the RMSEs were observed at the 3rd or 4th iteration step,and increased subsequently.This can be explained by the fact that congestion grew worse in later iterations.There were no significant differences among the cases in terms of their RMSE.The number of iteration steps needed to produce a feasible optimal result with the proposed method was nearly identical to that for Yang et al.’s method [20].Considering that SPSA requires simulations twice as lower-level optimizations to calculate the gradient[27]and that simulation runs occupy a large portion of the estimation time(see also Subsection 6.2),the proposed method is computationally more efficient than Yang et al.’s method.

    Figs.7a–7c show scatter plots obtained from a trial whose RMSE was the median at the 4th iteration.The correlation and regression coefficients are sufficiently high and similar to those in the verification case.Therefore,the authors concluded that,the link traffic volumes were reproduced well as a whole.However,from a different perspective,the fact that the RMSE remains at approximately 80 veh/h shows a limitation of the application in its current stage.Some large differences between the simulated and the observed link traffic volumes were seen at various observation points,with differences of 463 and 61 representing the most remarkable examples (Fig.7a).The reason for this is traffic volume drop due to congestion.Congestion detection leads to suppressing traffic volume on all links evenly.Thus,the RMSE can be improved only to a certain extent,as long as the entire reproducibility of the link traffic volumes across all links is limited.Further improvements require simulation parameters that more accurately describe driver behavior.

    Figure 7:Scatter plots for cases with/without congestion detection with regression lines

    Next,the results were compared to a study by Shafiei et al.[18] that used gradient-based optimization.Their method is also implemented offline and is nearly the same as the method described in Subsection 4.2 for preparing an initial OD matrix.

    Shafiei et al.pointed out that if the gradient of the objective function was simply used for updating the OD matrix solution candidate,it would be invalid under congested conditions.Though their approach is similar to the congestion detection method employed in this research,as given in Eq.(8),they dealt with the problem by considering the results of a sensitivity analysis of link traffic volume to determine the assignment matrixA.In brief,they attempted to substitute Eq.(4)with a first-order Taylor approximation,similar to Eq.(11).Note that the assignment matrixAchanges from a constant to a function of the OD vectorx.

    Here,subscriptjdenotes the index of the row,and the differential symbol denotes vector differentiation.x?is the current OD matrix candidate solution.

    Shafiei et al.calculated the second term of Eq.(11) in partial rather than total terms.This operation is based on introducing the initial OD matrix into the objective function.

    As a result of testing their approach on a small-scale network,Shafiei et al.obtained an RMSE of 20.79 in a heavy congestion scenario.Clearly,this is much smaller than the RMSE of 70–80 obtained in the present study.The authors consider that the difference arises from the dynamic rerouting feature implemented in the traffic simulator used by Shafiei et al.In terms of distributing demand,full dynamic rerouting affects each vehicle’s routing strategy dynamically,whereas the presently proposed method affects traffic volume statically outside the simulation.Therefore,the authors doubt that the larger RMSE for the proposed method arises from fictitious congestion occurring in the simulation.Subsection 6.1 clarifies the congestion behavior.

    6 Discussion

    6.1 Suppression Effect in Congestion

    The change of the total number of stuck vehicles,which is calculated simply by summing all values for all links,is shown in Fig.8.

    Figure 8:Change in total number of stuck vehicles

    Fig.8 shows that there were fewer stuck vehicles in the Strong case than in the Control case.Although there was no statistically significant difference,in comparing the mean values at the 4th iteration step,the number of stuck vehicles is 12%less in the Strong case than in the Control case.This result shows that the proposed congestion detection successfully reduces the number of such vehicles under the situation employing restricted dynamic rerouting.

    The results also show that the total number of stuck vehicles increases monotonically at each iteration in all cases.In iteration steps 1–4,since the total number of generated vehicles increases with each iteration,there is an increase in the total number of stuck vehicles.In later iteration steps,as noted in Subsection 5.2,there are many stuck vehicles and it is difficult to measure the true demand on the link because of the heavy congestion.In addition,the regression coefficient is always below 1.0,i.e.,the renewal of the OD matrices in the estimation loop always leads to increased demand in some OD pairs.This renewal causes positive feedback for congestion.It is also noteworthy that the number of stuck vehicles increases greatly in these later iteration steps.Considering above,these results suggest that the maximum number of iterationskmaxshould be 5.

    6.2 Computational Efficiency and Scalability

    The authors used a simulator with restricted dynamic rerouting to reduce the computational time in the OD estimation process.The most calculation-expensive operation in the simulation is the route search.In Dijkstra’s algorithm,the origin of the A?algorithm,the worst-case time complexity isO(E+VlogV),whereEandVare the numbers of links and nodes in the network,respectively.The A?algorithm is more efficient than Dijkstra’s method,but not to the extent that it significantly reduces the order of the time complexity.When invoked dynamic rerouting based on the A?algorithm,the computational cost is multiplied by the total traffic demand and the total number of rerouting operations per vehicle.Especially,the latter corresponds to the size of the network.Because driving on a larger network increases rerouting chance during driving.Though the computational complexity for this cannot be defined strictly,it can be formulated likeO(NE(E+VlogV))with some assumptions;The first assumption is that total traffic demand is proportional to the number of OD pairs.The second is that number of rerouting operations is proportional to the number of links.

    Here,considering additional assumptions,the formulation is expressed likeO(N2logN).The assumptions are like below: First,the number of links is proportional to the number of nodes:O(E)=O(V).There is no significant variation in the number of links that connect to one node.Second,the number of the origin or destination nodesPis proportional to the total number of nodes:O(P)=O(V).Third,the number of the OD pairs is the total number of combinations of origin and destination nodes:O(N1/2)=O(P).

    Quadratic programming was used in this study reduce the number of simulation runs in one iteration,and a majority of the computational time is consumed in the simulation portion of the procedure.For example,the mean values in one iteration of the Control case were 131.4 s in running simulation and 38.6 s in running QP solver.Even in the free-flow state,the computational time for the simulation predominates in the OD estimation.Our work supports that 75%of the computation time is spent in routing operation through a simulation run[36].By introducing dynamic rerouting into the simulator,the computational time for the simulation is expected to be several times longer than in the Control case.

    Here,Fig.9 shows the relationship between computational time for QP and simulation.The red points indicate these typical values in the examined cases.Purple and green lines are the approximated computation time discussed above.This figure shows that if the number of nodes is smaller than about 400–500,the simulation time is lower than QP.This result suggests that when computing the OD estimation in a typical medium-sized city in Japan like this experiment,since the computational cost for simulation is larger than upper-level optimization,the proposed method has an effect in reducing calculation cost.

    Figure 9:Estimated computational time for QP and simulation

    When dealing with extremely large networks,the QP solver may take longer computational time than the simulation.The worst-case time complexity for the QP algorithm is roughly approximated asO(N3) under the assumption that the linear equations are solved with LU decomposition and that upper or lower bounds are not considered.These rough estimates show that the QP solver could potentially cause a bottleneck in the estimation procedure.However,it is worthwhile to reduce the simulation time by using the proposed method even for large networks where QP is dominant.

    If we need to reduce the time for solving QP,reducing the number of variables is preferable.Since it is possible to conduct clustering for the results of a traffic simulation[37],a reduced-order estimation model can be developed using such a method.

    7 Conclusions

    In this paper,the OD estimation problem was formulated as a bi-level programming problem;the authors introduced congestion detection into the process of updating the OD matrices (upperlevel optimization) to combine with a MAS-based microscopic simulator (lower-level optimization)with restricted dynamic rerouting.Because full dynamic rerouting functions are computationally expensive,restricted dynamic rerouting was used to reduce the computational cost.At the same time,a congestion detection was introduced to suppress excess traffic demand that cannot be treated well by the restricted dynamic rerouting.Although the congestion detection was combined with restricted dynamic rerouting,the method can also be applied to a simulator with full dynamic rerouting.The proposed congestion detection could reduce fictitious congestions that result from restricting dynamic rerouting.Test results showed that stuck vehicles were reduced by 12%.The reduction of fictitious congestion can improve the numerical stability of the entire estimation method.

    Additionally,the reproducibility of link traffic volume using the proposed method proved to be nearly the same as that of existing methods.At the same time,the proposed method is computationally more efficient since the calculation cost of dynamic rerouting is reduced.This advantage is particularly evident when targeting wider-area or higher-demand conditions.However,the proposed method guarantees that reproducibility at the link without observation.It is required additional observation points or additional data source to fit the volume at such a link,velocities or other measurements.The authors regard it as for future task.

    Importantly,the proposed method can be applied to any other traffic simulators that can output link traffic volume and the number of stuck vehicles.In addition,the proposed approach of adding information to the analysis may contribute to enhancing the efficiency of NN-based methods.For example,the authors have shown in the previous work [38] that a similar approach improves the efficiency of traffic state estimation in the event of a traffic incident using a graph convolutional recurrent neural network(GCRNN).

    It would be prudent for future work to consider utilizing other observable traffic attributes as data sources for validation.For example,as described in Subsection 3.3,since the number of stuck vehicles corresponds to congestion length in the real world,it should be taken into account in the proposed method.

    Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author,K.Abe,upon reasonable request,except the observed traffic volume data.The observed traffic volume data are parts of the traffic survey result conducted by Fukui prefecture and the authors are permitted only to use them,not to distribute them.The simulator which is employed,ADVENTURE_Mates,is an open source software and is available at https://adventure.sys.t.u-tokyo.ac.jp/.Everyone can freely use this simulator for the purposes permitted under the simulator’s license and terms.

    Acknowledgement: The first author was supported by the H-UTokyo Lab.Energy Project.

    Funding Statement:This work was financially supported by JSPS KAKENHI(Grant Nos.15H01785 and 19H02377).

    Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

    美女大奶头视频| 亚洲一区二区三区色噜噜| 又紧又爽又黄一区二区| 网址你懂的国产日韩在线| 一二三四在线观看免费中文在| av天堂中文字幕网| 国产精品久久久久久久电影 | 18禁观看日本| 动漫黄色视频在线观看| 三级男女做爰猛烈吃奶摸视频| 一个人看的www免费观看视频| 亚洲欧美精品综合一区二区三区| 亚洲 欧美 日韩 在线 免费| 午夜福利在线观看免费完整高清在 | 人妻丰满熟妇av一区二区三区| 99热精品在线国产| 日韩欧美国产在线观看| 国产 一区 欧美 日韩| 在线播放国产精品三级| 国产又黄又爽又无遮挡在线| 精品久久久久久成人av| 午夜成年电影在线免费观看| 久久人妻av系列| 日日夜夜操网爽| 婷婷亚洲欧美| 国产精品久久久久久亚洲av鲁大| 欧美日本亚洲视频在线播放| 国产av在哪里看| 一级a爱片免费观看的视频| 国产又色又爽无遮挡免费看| 男女下面进入的视频免费午夜| 男女下面进入的视频免费午夜| 九九久久精品国产亚洲av麻豆 | 69av精品久久久久久| 岛国在线免费视频观看| 给我免费播放毛片高清在线观看| 国产麻豆成人av免费视频| 99久久精品国产亚洲精品| 成人三级做爰电影| 真实男女啪啪啪动态图| 在线观看舔阴道视频| h日本视频在线播放| 1000部很黄的大片| 亚洲精品在线观看二区| 亚洲乱码一区二区免费版| 国产一区在线观看成人免费| 免费av毛片视频| 丰满人妻一区二区三区视频av | 超碰成人久久| 欧美绝顶高潮抽搐喷水| 欧美日韩瑟瑟在线播放| 成人高潮视频无遮挡免费网站| 婷婷亚洲欧美| 中文亚洲av片在线观看爽| 国产精品亚洲一级av第二区| 精品一区二区三区视频在线观看免费| 国产精品一区二区三区四区免费观看 | 国产成年人精品一区二区| 人妻夜夜爽99麻豆av| 国产精品亚洲一级av第二区| 亚洲午夜理论影院| 少妇熟女aⅴ在线视频| 精品福利观看| 亚洲精品在线美女| 我的老师免费观看完整版| 午夜视频精品福利| 亚洲成人精品中文字幕电影| 亚洲国产精品999在线| 国产又色又爽无遮挡免费看| 国产精品久久久人人做人人爽| 中文字幕人妻丝袜一区二区| 久久精品综合一区二区三区| 免费大片18禁| 亚洲国产欧美一区二区综合| 国产精品久久久久久精品电影| 免费电影在线观看免费观看| 亚洲熟妇中文字幕五十中出| 久久久国产精品麻豆| 岛国在线免费视频观看| 啪啪无遮挡十八禁网站| 男人的好看免费观看在线视频| 亚洲色图 男人天堂 中文字幕| 国产亚洲av嫩草精品影院| 午夜激情福利司机影院| 久久久久国产一级毛片高清牌| 九色国产91popny在线| 成人精品一区二区免费| h日本视频在线播放| 亚洲 欧美 日韩 在线 免费| 日韩精品青青久久久久久| 国产野战对白在线观看| avwww免费| 深夜精品福利| 国产精品九九99| 一夜夜www| 欧美乱妇无乱码| av片东京热男人的天堂| 日本三级黄在线观看| 99热精品在线国产| 特级一级黄色大片| 变态另类成人亚洲欧美熟女| 1000部很黄的大片| 少妇人妻一区二区三区视频| av国产免费在线观看| 欧美黑人欧美精品刺激| 国产三级黄色录像| 精品国产乱码久久久久久男人| 欧美成人一区二区免费高清观看 | 人妻夜夜爽99麻豆av| 亚洲国产看品久久| 免费高清视频大片| 最近最新中文字幕大全免费视频| 香蕉国产在线看| 久久精品国产99精品国产亚洲性色| 男女那种视频在线观看| 欧美乱码精品一区二区三区| 亚洲精品在线美女| 97人妻精品一区二区三区麻豆| 黄色女人牲交| 九九热线精品视视频播放| 国产精品久久久久久人妻精品电影| 在线观看日韩欧美| 欧美中文综合在线视频| 国产伦精品一区二区三区四那| 国产乱人伦免费视频| av中文乱码字幕在线| 黄片小视频在线播放| 少妇裸体淫交视频免费看高清| 欧美xxxx黑人xx丫x性爽| 国产美女午夜福利| 欧美3d第一页| 亚洲在线观看片| 激情在线观看视频在线高清| 麻豆久久精品国产亚洲av| 在线视频色国产色| 丰满人妻一区二区三区视频av | 精品国产亚洲在线| 午夜免费成人在线视频| 久久久成人免费电影| 视频区欧美日本亚洲| 丝袜人妻中文字幕| 美女高潮的动态| 啦啦啦韩国在线观看视频| 一二三四在线观看免费中文在| 亚洲自拍偷在线| www.999成人在线观看| 91在线精品国自产拍蜜月 | bbb黄色大片| 亚洲天堂国产精品一区在线| 久久久久久久久中文| 久久香蕉精品热| 男女午夜视频在线观看| 波多野结衣高清作品| 亚洲精华国产精华精| 欧美丝袜亚洲另类 | 熟妇人妻久久中文字幕3abv| 999久久久精品免费观看国产| 美女免费视频网站| 久久精品综合一区二区三区| 美女cb高潮喷水在线观看 | 精品国产三级普通话版| 国产精品自产拍在线观看55亚洲| 一本久久中文字幕| 亚洲avbb在线观看| 亚洲一区二区三区不卡视频| 精品国产亚洲在线| 午夜福利欧美成人| 国产精品自产拍在线观看55亚洲| 99在线视频只有这里精品首页| 日韩欧美国产在线观看| 日韩国内少妇激情av| 国产三级中文精品| 国产私拍福利视频在线观看| 国产精品女同一区二区软件 | 国产人伦9x9x在线观看| 天天一区二区日本电影三级| 亚洲 国产 在线| 九色成人免费人妻av| 亚洲五月天丁香| 久久草成人影院| 精品熟女少妇八av免费久了| e午夜精品久久久久久久| 99热精品在线国产| 亚洲 欧美 日韩 在线 免费| 国产精品久久久人人做人人爽| 国产高清有码在线观看视频| 两人在一起打扑克的视频| 啦啦啦韩国在线观看视频| 日本黄色视频三级网站网址| 国内揄拍国产精品人妻在线| 国产黄色小视频在线观看| 欧美黄色淫秽网站| 精品久久久久久久人妻蜜臀av| 99久久久亚洲精品蜜臀av| 国产成年人精品一区二区| 俺也久久电影网| 精品久久久久久久久久久久久| 欧美不卡视频在线免费观看| 日韩成人在线观看一区二区三区| 国产精品九九99| 给我免费播放毛片高清在线观看| 在线播放国产精品三级| 国产精品98久久久久久宅男小说| 97人妻精品一区二区三区麻豆| 成在线人永久免费视频| 三级毛片av免费| 精品国内亚洲2022精品成人| 三级男女做爰猛烈吃奶摸视频| 色噜噜av男人的天堂激情| 欧美日韩乱码在线| xxxwww97欧美| 看片在线看免费视频| 悠悠久久av| 99热只有精品国产| 久久久久国内视频| 中文字幕最新亚洲高清| 免费无遮挡裸体视频| 五月伊人婷婷丁香| 51午夜福利影视在线观看| 国产亚洲欧美98| 韩国av一区二区三区四区| 午夜免费成人在线视频| 国产精品九九99| 午夜激情福利司机影院| 一区二区三区国产精品乱码| 亚洲专区中文字幕在线| 床上黄色一级片| 亚洲欧美日韩高清专用| 九色国产91popny在线| 久久香蕉精品热| 一进一出好大好爽视频| 观看美女的网站| 日本撒尿小便嘘嘘汇集6| 日日摸夜夜添夜夜添小说| 国内少妇人妻偷人精品xxx网站 | 最新美女视频免费是黄的| 国产亚洲av高清不卡| 国产精品香港三级国产av潘金莲| 黄频高清免费视频| 中出人妻视频一区二区| 国产免费男女视频| 亚洲美女视频黄频| 嫁个100分男人电影在线观看| 亚洲第一电影网av| 激情在线观看视频在线高清| 一a级毛片在线观看| 99国产精品一区二区蜜桃av| 法律面前人人平等表现在哪些方面| 黄色成人免费大全| 国产精品久久久人人做人人爽| 高清毛片免费观看视频网站| 国产 一区 欧美 日韩| 国产亚洲精品久久久久久毛片| 18禁国产床啪视频网站| 免费看十八禁软件| 精品久久久久久久久久免费视频| 全区人妻精品视频| 亚洲精品美女久久久久99蜜臀| 久久午夜亚洲精品久久| 高潮久久久久久久久久久不卡| 国产午夜福利久久久久久| 十八禁网站免费在线| 亚洲黑人精品在线| 国产精华一区二区三区| 国产亚洲欧美在线一区二区| 在线观看一区二区三区| 亚洲国产高清在线一区二区三| 国产 一区 欧美 日韩| 亚洲色图 男人天堂 中文字幕| 国产伦人伦偷精品视频| 精品熟女少妇八av免费久了| 国产成人av激情在线播放| 色哟哟哟哟哟哟| 综合色av麻豆| 亚洲欧美激情综合另类| 色综合婷婷激情| 村上凉子中文字幕在线| 51午夜福利影视在线观看| 校园春色视频在线观看| 精品福利观看| 国产蜜桃级精品一区二区三区| 嫁个100分男人电影在线观看| 99热这里只有是精品50| 国产亚洲av高清不卡| 亚洲av免费在线观看| 嫁个100分男人电影在线观看| 天堂网av新在线| 99久久久亚洲精品蜜臀av| 亚洲精品粉嫩美女一区| 欧美日韩一级在线毛片| 亚洲成av人片在线播放无| 国产野战对白在线观看| 精品免费久久久久久久清纯| 亚洲欧美精品综合久久99| 亚洲国产精品合色在线| 国产蜜桃级精品一区二区三区| 俺也久久电影网| 国内精品久久久久精免费| 老汉色∧v一级毛片| 麻豆国产97在线/欧美| 日本撒尿小便嘘嘘汇集6| 白带黄色成豆腐渣| 国产精品免费一区二区三区在线| 欧美乱妇无乱码| 高潮久久久久久久久久久不卡| 两个人的视频大全免费| 国产三级黄色录像| 无限看片的www在线观看| 久久久久国产精品人妻aⅴ院| 美女午夜性视频免费| 国产精品 国内视频| 最近最新中文字幕大全电影3| 国产爱豆传媒在线观看| 国产精品女同一区二区软件 | 在线免费观看不下载黄p国产 | 动漫黄色视频在线观看| 亚洲av免费在线观看| 青草久久国产| 精品久久久久久久久久免费视频| 999精品在线视频| 亚洲精品国产精品久久久不卡| 三级毛片av免费| 欧美乱码精品一区二区三区| 少妇人妻一区二区三区视频| 欧美xxxx黑人xx丫x性爽| 91老司机精品| 好男人电影高清在线观看| 国产精品爽爽va在线观看网站| 在线免费观看不下载黄p国产 | 午夜福利在线观看免费完整高清在 | 久久久久亚洲av毛片大全| 久久久久精品国产欧美久久久| 久久久成人免费电影| 欧美性猛交黑人性爽| 美女高潮的动态| 变态另类成人亚洲欧美熟女| 中文字幕人妻丝袜一区二区| 无人区码免费观看不卡| 欧美绝顶高潮抽搐喷水| 久久午夜亚洲精品久久| 久久婷婷人人爽人人干人人爱| 琪琪午夜伦伦电影理论片6080| 国产成人av教育| 日韩国内少妇激情av| 日本 av在线| 日韩高清综合在线| 淫妇啪啪啪对白视频| 男人舔女人的私密视频| 91av网站免费观看| 午夜成年电影在线免费观看| 精品国产超薄肉色丝袜足j| 午夜福利免费观看在线| 国产精品亚洲av一区麻豆| 曰老女人黄片| 免费av毛片视频| 国产成人影院久久av| 99精品在免费线老司机午夜| 国产av在哪里看| 国产成人精品久久二区二区91| 国产精品 欧美亚洲| 免费在线观看影片大全网站| 长腿黑丝高跟| 亚洲国产色片| 国产在线精品亚洲第一网站| 色老头精品视频在线观看| 免费av毛片视频| a级毛片在线看网站| 人妻久久中文字幕网| 国产伦在线观看视频一区| 18禁观看日本| 老汉色av国产亚洲站长工具| 亚洲五月天丁香| 日本一二三区视频观看| 中文在线观看免费www的网站| 九色国产91popny在线| 欧美日韩精品网址| 国产精品av视频在线免费观看| 少妇裸体淫交视频免费看高清| 成人性生交大片免费视频hd| 看黄色毛片网站| 人人妻,人人澡人人爽秒播| 亚洲av日韩精品久久久久久密| 亚洲色图 男人天堂 中文字幕| 中文字幕久久专区| www.自偷自拍.com| 国产一区二区在线av高清观看| 久久午夜亚洲精品久久| 免费观看精品视频网站| 男人舔女人的私密视频| 欧美zozozo另类| 中文字幕最新亚洲高清| 成人特级黄色片久久久久久久| 美女午夜性视频免费| 亚洲五月天丁香| 美女扒开内裤让男人捅视频| 色老头精品视频在线观看| 国产一区二区在线观看日韩 | 中文字幕av在线有码专区| 女人高潮潮喷娇喘18禁视频| 嫩草影院精品99| 97超级碰碰碰精品色视频在线观看| 欧美乱码精品一区二区三区| 日本与韩国留学比较| cao死你这个sao货| 国产又黄又爽又无遮挡在线| 免费在线观看影片大全网站| 亚洲无线在线观看| 中文字幕精品亚洲无线码一区| 99久久99久久久精品蜜桃| 亚洲专区中文字幕在线| 国产亚洲欧美在线一区二区| 熟女人妻精品中文字幕| 久久久精品大字幕| 成人国产一区最新在线观看| 久久香蕉国产精品| x7x7x7水蜜桃| 好看av亚洲va欧美ⅴa在| 国产高清激情床上av| 高清毛片免费观看视频网站| 黄色丝袜av网址大全| 我的老师免费观看完整版| 国产精品精品国产色婷婷| 91麻豆av在线| 九九在线视频观看精品| 免费电影在线观看免费观看| 俺也久久电影网| 久久午夜综合久久蜜桃| 日本黄色视频三级网站网址| 亚洲国产欧美人成| www国产在线视频色| 一级黄色大片毛片| 国产成年人精品一区二区| 美女cb高潮喷水在线观看 | 草草在线视频免费看| 精品欧美国产一区二区三| 性欧美人与动物交配| 一区二区三区国产精品乱码| 成人三级黄色视频| 黄片大片在线免费观看| 99久久精品一区二区三区| 亚洲国产精品999在线| 久久精品91无色码中文字幕| 国内精品久久久久精免费| www.自偷自拍.com| 国产高清有码在线观看视频| 国产精品 欧美亚洲| 日本黄色视频三级网站网址| 视频区欧美日本亚洲| 欧美一级毛片孕妇| 国产精品久久久久久久电影 | 国产精品久久久久久精品电影| 两性夫妻黄色片| 欧美高清成人免费视频www| 一区福利在线观看| 好男人在线观看高清免费视频| 国产亚洲精品av在线| 婷婷六月久久综合丁香| 欧美色视频一区免费| 久99久视频精品免费| 久久精品亚洲精品国产色婷小说| 久久精品国产清高在天天线| 在线看三级毛片| 手机成人av网站| 亚洲精品美女久久av网站| 精品无人区乱码1区二区| 特级一级黄色大片| 97超视频在线观看视频| 欧美丝袜亚洲另类 | 国产激情偷乱视频一区二区| 午夜精品久久久久久毛片777| 黑人操中国人逼视频| 久久婷婷人人爽人人干人人爱| 亚洲色图 男人天堂 中文字幕| 欧美乱码精品一区二区三区| 国产高清激情床上av| 欧美日韩一级在线毛片| 亚洲欧美日韩东京热| 十八禁网站免费在线| 国产亚洲欧美98| 女同久久另类99精品国产91| 视频区欧美日本亚洲| 国产激情久久老熟女| 两个人的视频大全免费| 岛国在线免费视频观看| 国产成人av激情在线播放| 色在线成人网| 在线观看午夜福利视频| 一级毛片高清免费大全| 久久亚洲精品不卡| 99久久国产精品久久久| 一区二区三区激情视频| 精品久久久久久久末码| 日韩精品中文字幕看吧| 久久中文字幕一级| 色综合婷婷激情| 久久香蕉国产精品| 久久久精品大字幕| 国内揄拍国产精品人妻在线| 国产97色在线日韩免费| 好男人电影高清在线观看| 禁无遮挡网站| 在线观看66精品国产| cao死你这个sao货| 国产免费男女视频| 露出奶头的视频| 天堂av国产一区二区熟女人妻| 国产欧美日韩精品亚洲av| 99热这里只有是精品50| 美女黄网站色视频| 国产av一区在线观看免费| 亚洲美女视频黄频| 特大巨黑吊av在线直播| 亚洲国产精品sss在线观看| 国产69精品久久久久777片 | 国产精品美女特级片免费视频播放器 | 男人舔女人的私密视频| 久久久国产成人免费| 全区人妻精品视频| 日本 欧美在线| 婷婷精品国产亚洲av| 亚洲aⅴ乱码一区二区在线播放| 一本综合久久免费| 男女之事视频高清在线观看| 搡老熟女国产l中国老女人| 一个人免费在线观看的高清视频| 91老司机精品| 男女床上黄色一级片免费看| 日本一本二区三区精品| 亚洲精品色激情综合| 真实男女啪啪啪动态图| 波多野结衣高清作品| 91老司机精品| 中文字幕人妻丝袜一区二区| 亚洲av成人精品一区久久| 精品乱码久久久久久99久播| 久久午夜综合久久蜜桃| 亚洲av免费在线观看| 国产精品一区二区免费欧美| 18禁国产床啪视频网站| 久久久久精品国产欧美久久久| 午夜免费观看网址| 国产高清videossex| 久久久成人免费电影| 亚洲欧洲精品一区二区精品久久久| 免费av毛片视频| 亚洲天堂国产精品一区在线| 国产精品一及| 丰满人妻熟妇乱又伦精品不卡| 久久精品国产清高在天天线| 久久人妻av系列| 波多野结衣高清作品| 成人一区二区视频在线观看| 十八禁网站免费在线| 欧美最黄视频在线播放免费| 九色成人免费人妻av| 一区二区三区国产精品乱码| 俄罗斯特黄特色一大片| 亚洲欧美日韩东京热| 日本 欧美在线| 又紧又爽又黄一区二区| 禁无遮挡网站| 人妻夜夜爽99麻豆av| 熟女人妻精品中文字幕| 国产精品久久久久久亚洲av鲁大| 老司机午夜十八禁免费视频| 少妇人妻一区二区三区视频| 国产单亲对白刺激| 人人妻人人澡欧美一区二区| 老司机午夜福利在线观看视频| 亚洲片人在线观看| 中文在线观看免费www的网站| 亚洲精品色激情综合| 窝窝影院91人妻| 听说在线观看完整版免费高清| 搞女人的毛片| 国产成人精品久久二区二区91| 在线a可以看的网站| 两性夫妻黄色片| 在线观看66精品国产| 国产精华一区二区三区| 精品电影一区二区在线| 两人在一起打扑克的视频| 亚洲国产精品999在线| 欧美成人免费av一区二区三区| 久久久成人免费电影| 夜夜躁狠狠躁天天躁| 19禁男女啪啪无遮挡网站| 国产熟女xx| 国产成人精品无人区| 国产精品亚洲一级av第二区| 最近在线观看免费完整版| 人人妻,人人澡人人爽秒播| 欧美一级毛片孕妇| 国产欧美日韩一区二区三| 久久久精品大字幕| 亚洲无线在线观看| 超碰成人久久| 最近最新免费中文字幕在线| 三级国产精品欧美在线观看 | 久久精品91蜜桃| bbb黄色大片| 一级a爱片免费观看的视频| 伦理电影免费视频| 亚洲精品美女久久久久99蜜臀| 可以在线观看毛片的网站| 国内毛片毛片毛片毛片毛片| 老司机在亚洲福利影院| 美女cb高潮喷水在线观看 | 1024手机看黄色片| 好男人电影高清在线观看| 亚洲精品456在线播放app | 91在线精品国自产拍蜜月 | 一个人看视频在线观看www免费 | 香蕉国产在线看| www.自偷自拍.com|