Table Of Content

IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 1 Accelerated Evaluation of Automated Vehicles Using Piecewise Mixture Models Zhiyuan Huang1 , Ding Zhao2, Henry Lam1, David J. LeBlanc2 Abstract—The process to certify highly Automated Vehicles both Google and Tesla update their software throughout the hasnotyetbeendefinedbyanycountryintheworld.Currently, process, which may have improved safety, but the newest companies test Automated Vehicles on public roads, which is version of the AV has not accumulated that many miles as time-consuming and inefficient. We proposed the Accelerated they have claimed. In summary, today’s best practice adopted 7 Evaluation concept, which uses a modified statistics of the 1 surrounding vehicles and the Importance Sampling theory to by the industry is time-consuming and inefficient. A better 0 reduce the evaluation time by several orders of magnitude, approach is needed. 2 while ensuring the evaluation results are statistically accurate. In this paper, we further improve the accelerated evaluation n concept by using PiecewiseMixture Distribution models, instead a A. Related Researches of Single Parametric Distribution models. We developed and J applied this idea to forward collision control system reacting BesidestheN-FOT,thetestmatrixapproach[6],[7]andthe 1 to vehicles making cut-in lane changes. The behavior of the worst-casescenariosapproach[8],[9],[10]aretwoalternative 3 cut-in vehicles was modeled based on more than 403,581 lane methods for vehicle evaluation. Our approach follows the changes collected by the University of Michigan Safety Pilot ] ModelDeploymentProgram.Simulationresultsconfirmthatthe AcceleratedEvaluationconceptweproposed[11]toprovidea Y accuracy and efficiency of the Piecewise Mixture Distribution brand-new alternative. The basic concept is that as high-level S method outperformed single parametric distribution methods in AVs just began to penetrate the market, they mainly interact . accuracy and efficiency, and accelerated the evaluation process s withhuman-controlledvehicles(HVs).Thereforewefocuson c by almost four orders of magnitude. modeling the interaction between the AV and the HV around [ it. The evaluation procedure involves four steps: 1 Index Terms—automated vehicles, testing, evaluation, safety v • Modelthebehaviorsoftheprimaryothervehicles(POVs) 5 represented by f(x) as the major disturbance to the AV I. INTRODUCTION 1 using large-scale naturalistic driving data 9 IT is critical to thoroughly and rigorously test and evaluate • Skew the disturbance statistics from f(x) to modified 8 an Automated Vehicle (AV) before its release. Recent statistics f∗(x) (accelerated distribution) to generate 0 crashes involving a Google self-driving car [1] and a Tesla more frequent and intense interactions between AVs and . 1 Autopilotvehicle[2]attractedthepublic’sattentiontoAVtest- POVs 0 ing and evaluation. While these AVs are generally considered • Conduct accelerated tests with f∗(x) 7 as industrial leaders, because they use public road for testing, 1 • Use the Importance Sampling (IS) theory to skew back statistically they have not yet accumulated enough miles. The : the results to understand real-world behavior and safety v TeslaAutopilot,inparticular,wascriticizedforbeingreleased benefits i X too early in the hands of the general public [3]. This approach has been successfully applied to evaluate Currently, there are no standards or protocols to test AVs r AVs in the frontal crash with a cut-in vehicle [11] and also a at automation level 2 or higher. Many companies adopt the frontalcrashwithaleadvehicle[12],[13].Thisapproachwas Naturalistic Field Operational Tests (N-FOT) approach [4]. confirmed to significantly reduce the evaluation time while However, this method is inefficient because safety critical accurately preserving the statistical behavior of the AV-HV scenarios rarely happen in daily driving. The Google Self- interaction. In the previous studies, the evaluation time was driving cars accumulated 1.9 million driving. This distance, reduced by two to five orders of magnitudes - the accelerated although sounds a lot, provides limited exposure to critical rate depends on the test scenarios, where rarer events achieve events, given that U.S. drivers encounter a police reported higher accelerated rate. The non-accelerated models and the crash every five hundred thousand miles on average and fatal accelerated models were built based on signal component crash every one hundred million miles [5]. In the meantime, distributions. While this method does benefit from its simple *This work was funded by the Mobility Transformation Center at the mathematical form, it has a few drawbacks as illustrated in University of Michigan with grant no. N021552. Z. Huang and D. Zhao Fig. 1. i) The fitting of the rare events (usually the tail part contributedequallytotheresearch. 1Zhiyuan Huang ([email protected]) and Henry Lam of the statistical distributions) would be dominated by the ([email protected])areintheDepartmentofIndustrialandOperations fitting of the normal driving behaviors (the majority part of EngineeringattheUniversityofMichigan. the distributions), which may induce large errors. ii) The full 2Ding Zhao (corresponding author: [email protected]) and potential in higher accelerated rate is not achieved due to the DavidLeBlanc([email protected])areintheUniversityofMichi- ganTransportationResearchInstitute. lack of flexibility of the modified accelerated models. IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 2 Fig. 1. Acceleration evaluation based on single parametric distribution and PiecewiseMixtureDistribution. Fig.2. LanechangedatacollectedbySPMDvehicle. B. Contribution In this paper, we proposed a more general framework for independently with regard to the value of v . By comparing L the Accelerated Evaluation method to overcome the afore- among 17 types of commonly used distribution templates, we mentioned limitations based on Piecewise Mixture Distribu- selected the Pareto distribution to model R−1 and used the L tion Models as illustrated in Fig. 1 b). In this paper, we exponential distribution for TTC−1 segments. L implemented the Accelerated Evaluation method under the Using the empirical distribution of v and parametric dis- L new framework. Comparing to our previous work [14], we tributions of R and TTC , we drew values from these L L thoroughly discuss the Cross Entropy method with proposed distributions as inputs to simulate the AV-HV interaction. framework in this paper. We present practical tips to over- The outcome from the simulation can be considered as an come numerical issues and reduce computational efforts. We event indicator function I (x) that returns {1,0} depending ε demonstratethismethodbyevaluatingthelongitudinalcontrol on the event of interest. Given the stochastic distribution of system reacting to vehicles making cut-in lane changes. the variables and the event indicator function, we obtained the optimal exponential distribution for Importance Sampling C. Paper Structure by implementing the Cross Entropy method [16]. As we have shown in Fig. 1 a), we used only single parametric distribu- Section II will introduce the lane change model based on tions.Inthenextsection,weintroduceournewapproachusing single parametric distributions. In Section III, we present the Piecewise Mixture Distributions. new lane change model with Piecewise Mixture Distributions. We establish the Accelerated Evaluation in Section IV and discusstheCrossEntropymethodwithPiecewiseMixtureDis- III. LANECHANGEMODELWITHPIECEWISEMIXTURE tributionmodelsinSectionV.Simulationresultsarediscussed DISTRIBUTIONS in Section VI. Section VII concludes this paper. Although many commonly used parametric distributions have concise and elegant forms, they do not always describe II. ACCELERATEDEVALUATIONWITHSINGLE the data distribution well. Instead, a better fitting can be PARAMETRICDISTRIBUTIONS achieved by dividing the dataset into several subsets. We estimatethemodelparametersusingtheMaximumLikelihood The lane change events were extracted from the Safety Estimation (MLE) [17] in each subset. The general process of Pilot Model Deployment (SPMD) database [15]. With over MLE is as follow. 2 million miles of vehicle driving data collected from 98 Assume we have a family of distribution with Cumula- cars over 3 years, we identify 403,581 lane change events. tive Distribution Function (CDF) F(x|θ), where θ is the Previously [11], we used 173,692 events with a negative parametervectorofF.ThecorrespondingProbabilityDensity range rate to build a statistical model focusing on three key Function (PDF) of F is f(x|θ). Assuming that data D = variables that captured the effects of gap acceptance of the {X ,X ,...,X } is independently and identically distributed lanechangingvehicle:velocityoftheleadvehicle(v ),range 1 2 N L and the distribution is in the family of F(x|θ), we want to totheleadvehicle(RL)andtimetocollision(TTCL).TTCL find the most “likely” parameter θˆ. was defined as: R We define the likelihood function [18] as TTC =− L, (1) L R˙L L(θ|D)=P(D|θ)=ΠNn=1f(Xn|θ). (2) where R˙L is the relative speed. We call the estimation of θˆ that maximizes the likelihood The modeling of these three variables was hard to handle function the mostly likely estimation MLE. because of dependency, so we simplified it based on a crucial For computation convenience, we introduce the log- observation. Although TTC is dependent on v generally, L L likelihood function we split the data into 3 segments: v at 5 to 15 m/s, 15 to 25 L N m/sand25to35m/s.Withineachsegment,R isindependent (cid:88) L L(θ|D)=lnL(θ|D)= lnf(X |θ). (3) withv andTTC .ThisallowedustomodelTTC andR n L L L L n=1 IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 3 Since the logarithm is monotone, the log-likelihood function 1) MLEforboundedexponentialdistribution: Thebounded preservestheoptimizeroftheoriginalfunction.[19]Theopti- exponential distribution with rate θ has the form mizeroflog-likelihoodfunction,θˆ,istheMLEofdistribution θe−θx family F. We have the MLE as f(x|γ ≤x<γ )= (10) 1 2 e−θγ1 −e−θγ2 θˆ=argmax L(θ|D). (4) for γ ≤x<γ . θ 1 2 In the following, we describe the Piecewise Mixture Dis- For dataset D ={X1,...,XN}, the log-likelihood function tribution fitting concept based on MLE and we present the is bounded distribution fitting results. All optimization problems N (cid:88) presented in this section are tractable and can be solved by L(D|θ)= lnθ−θX −ln(e−θγ1 −e−θγ2), (11) n fminunc in MATLAB. n=1 where L is concave over θ. Although we cannot solve the A. General Framework of the Piecewise Mixture Distribution maximization analytically, it is solvable through numerical Lane Change Model methods. We define Piecewise Mixture Distribution to be distribution Therefore, the MLE of θ is given by the optimization with CDF in the form of N (cid:88) (cid:88)k max Nlnθ−Nln(e−θγ1 −e−θγ2)− θXn. (12) F(x)= π F (x|γ ≤x<γ ), (5) θ i θi i−1 i n=1 i=1 2) MLE for bounded normal distribution: Consider a where k is the number of truncation, (cid:80)ki=1πi = 1, and bounded normal distribution with mean 0 and variance θ2 Fi(x|γi−1 ≤x<γi)istheconditionalcumulativedistribution conditional on 0≤γ1 ≤x<γ2. The PDF is function, meaning that F (γ |γ ≤ x < γ ) = 0 and i i−1 i−1 i 1φ(x) Fi(γi|γi−1 ≤ x < γi) = 1 for i = 1,...,k. θi denotes the f(x|γ ≤x<γ )= θ θ . (13) parameter(s) for F . We can consider that π = P(γ ≤ 1 2 Φ(γ2)−Φ(γ1) i i i−1 θ θ x < γ ) and when x ≥ 0, we have γ = 0 and γ = ∞. By i 0 k The MLE of the bounded normal distribution is given by this definition, the PDF of the Piecewise Mixture Distribution is k (cid:80)N X2 γ γ f(x)=(cid:88)π f (x|γ ≤x<γ ). (6) max− n=1 n −Nlnθ−Nln(Φ( 2)−Φ( 1)). (14) i θi i−1 i θ 2θ2 θ θ i=1 3) Fitting mixture model with EM algorithm: Compared to In our case, θ = {π ,...,π ,θ ,...,θ }. Splitting D into 1 k 1 k single parametric distributions, mixture distribution combines pieces regarding the truncation points {γ ,...,γ }, gives 1 k−1 several classes of distribution and thus is more flexible. We data index sets S = {j|γ ≤ X < γ } for i = 1,...,k. i i−1 j i consider the fitting problem of mixture bounded normal dis- We can write the log-likelihood function as tribution. L(θ|D)=(cid:80)k (cid:80) lnπ The PDF of mixture of m bounded normal distribution can +(cid:80)ki=1(cid:80) n∈Silnf i(X |γ ≤x<γ ). (7) be written as i=1 n∈Si θi n i−1 i m We obtain the MLE of θ can be obtained by maximizing (cid:88) f(x|γ ≤x<γ )= p f (x|γ ≤x<γ ) (15) L(θ|D) over θ. Since L is concave over π , we take 1 2 j j 1 2 i j=1 ∂L =0 (8) where f is bounded Gaussian distribution with mean ∂π j i 0 and variance σ2. The parameters here are θ = j and get {p ,...,p ,σ2,...,σ2 }. We want to find MLE of p and σ2 πˆ =|S |/N. (9) 1 m 1 m j j i i for j =1,...,m. For parameters θi in Fi, it is known (7) to be the same The log-likelihood function for data D ={Xn}Nn=1 is as computing the MLE of θ with corresponding dataset i N m Di = {X|γi−1 ≤ X < γi and X ∈ D}. Since we L(θ|D)= (cid:88)ln(cid:88)pjfj(Xn|γ1 ≤x<γ2). (16) use bounded distribution for each F , below we explain the i n=1 j=1 estimationofparametersforthethreedistributionsweapplied We note that this is hard to solve directly, because there in later sections. is a sum within the log function. Therefore, we apply the TosamplefromaPiecewiseMixtureDistribution,wecould Expectation-Maximization (EM) [20] algorithm to find the use the inverse function approach. See Appendix A for the optimizer, i.e. MLE, for the parameters. details. We define Zj to denote whether or not the random number n X comesfrommixturedistributionj,j =1,...,m,andZj = B. Bounded Distribution n n {0,1}. We also introduce the expectation Wedevelopthreeboundeddistributionsandusetheminthe lane change model. E[Zj|X ]:=τj. (17) n n n IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 4 The EM algorithm starts with initial parameters {p ,σ }, with smaller variance. This is known as Importance Sampling j j j = 1,...,m. For data D = {X }N , we set com- [21] and F∗ is the IS distribution. n n=1 plete data as D = {X ,Z }N . The EM algorithm op- For estimating P(X ∈ ε),we note that an optimal IS c n n n=1 timizes E[L(θ|D )|D] in every step. The E step updates distribution c E[L(θ|D )|D], and the M step optimizes this function. The P(X ≤x, ε) c F∗∗(x)=F(x|ε)= (25) algorithm iterates E step and M step until reaching the P(x∈ε) convergence criterion. couldreducethevarianceofISestimationto0,buttheoptimal In our case, requires the knowledge of P(X ∈ ε). However, it guides the E[L(θ|D )|D]=(cid:80)N (cid:80)m τj(lnp +lnf (X )). selection of the IS distribution. c n=1 j=1 n j j n (18) B. Exponential Change of Measure SinceobjectiveE[l (θ|D )|D]intheMstepisconcaveover c c p and σ , we could maximize the objective function through Exponential change of measure is commonly used to con- j j an analytic approach for p : structF∗.Althoughtheexponentialchangeofmeasurecannot j guarantee convergence to optimal distribution, it is easy to (cid:80)N τj p = n=1 n. (19) implement and the new distribution generally stays within the j N same class of distribution. For σ , we can solve the following maximization problem Exponential change of measure distribution takes the form j through numerical approach. of f (x)=exp(θx−κ(θ))f(x), (26) (cid:18) (cid:19) θ X σj =argmσijn−τnjlnσj +τnjlnφ σjn − where θ is the change of measure parameter and κ(θ) is (cid:18) (cid:19) the log-moment generating function of original distribution f. γ γ τjln Φ( 2)−Φ( 1) . (20) When θ =0, we have f (x)=f(x). n σ σ θ j j For a bounded exponential distribution, the exponential See Appendix B for the full EM algorithm. change of measure distribution is (λ−θ)e−(λ−θ)x IV. ACCELERATEDEVALUATIONWITHIMPORTANCE fθ(x|γ1 ≤x<γ2)= e−(λ−θ)γ1 −e−(λ−θ)γ2, (27) SAMPLING whereλistheparameterforexponentialdistribution.Wenote Importance Sampling (IS) is thus used to accelerate the that f is still a bounded exponential distribution and λ = θ θ evaluationprocess,becausecrudeMonteCarlosimulationsfor λ−θ. rare events can be time-consuming. Below we describe the IS For a bounded normal distribution, the exponential change method. of measure distribution is 1φ(x−σ2θ) A. Important Sampling and Optimal IS distribution fθ(x|γ1 ≤x<γ2)= Φ(γ2−θσσ2)−σΦ(γ1−θσ2), (28) Let x be a random variable generated from distribution F, σ σ where the original distribution truncated from a normal distri- and ε ⊂ Ω where ε is the rare event of interest and Ω is the bution with parameters µ=0 and σ. We note that the change sample space. Our objective is to estimate of measure distribution is still a bounded normal distribution (cid:90) P(X ∈ε)=E[I (X)]= I (x)dF (21) with µ=θσ2 and σ. ε ε V. CROSSENTROPYMETHODANDIMPLEMENTATION where (cid:40) 1 x∈ε, Section IV discussed optimal IS distribution F∗∗ providing I (x)= (22) ε 0 variance estimation to the value of interest, whereas this 0 otherwise. section describes the Cross Entropy method used to estimate We can write the evaluation of rare events as the sample the “optimal” parameters θ, which minimizes the “distance” mean of Iε(x) between a parametric distribution Fθ and F∗∗ without know- ing F∗∗. The description below is based on the Piecewise N Pˆ(X ∈ε)= 1 (cid:88)I (X ), (23) Mixture Distribution structure. N ε n n=1 A. Introduction where X ’s are drawn from distribution F. i Since we have TheCrossEntropy,whichisalsoknownasKullback-Leibler (cid:90) (cid:90) dF distance [22], measures the similarity between distributions. E[Iε(X)]= Iε(x)dF = Iε(x)dF∗dF∗, (24) We define the Cross Entropy between function g and h as g(X) (cid:90) dwiestrcibauntiocnomFpu∗t,ewthheichsamhapslethmeeasnamoef sIuεp(Xpo)rtddFFw∗ithoveFr,thtoe D(g,h)=Eg[lnf(X)]= g(x)lng(x)dx− (cid:90) obtain an unbiased estimation of P(X ∈ε). By appropriately g(x)lnh(x)dx. (29) selecting F∗, the evaluation procedure obtains an estimation IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 5 From (25), we know that the PDF of the optimal IS We note that if we have two independent variables where distribution F∗∗ is f(x,y) = f(x)f(y), we can take a parametric distribution f∗∗(x)= Iε(x)f(x). (30) fΘor=ea{cθh,vθar}i.abTlheeaonbdjehcativveeffuΘn(cxt,ioyn)c=orrfeθs1p(oxn)dfiθn2g(yt)o,(w33h)eries P(x∈ε) 1 2 Since P(x ∈ ε) is generally unavailable, we use a para- N mWeetrwicadnitsttroibfiuntidonthFeθptaoraamppertoearcθh∗ththeaotpmtiminaiml IiSzedsistthreibuCtrioosns. mθax N1 n(cid:88)=1Iε(Xn,Yn)ffΘ(sX(Xnn,Y,Ynn))(lnfθ1(Xn)+ Entropy [23] between f∗∗ and fθ. We denote θ∗ as the lnfθ2(Yn)), (34) optimal parameter for the parametric distribution. Then the which can be decoupled into two optimization problem over minimization problem θ and θ respectively and I (X ,Y ) f(Xn,Yn) is a known minD(fθ,f∗∗) (31) co1nstant2given {Xn,Yn}. ε n n fΘs(Xn,Yn) θ We implement the Cross Entropy on the Piecewise Mixture is equivalent to Distribution with one variable. We note that we can apply the f(X) results to the lane change model, since the Cross Entropy ob- max E [I (X) lnf (X)], (32) θ θs ε fθs(X) θ jective function of independent variables can be implemented in (34). where f denotes the sampling distribution with parameters θs θ . We note that this is a generalized setting, since we can s use any sampling distribution f as long as it has the same B. OptimizationFunctionforPiecewiseMixtureDistributions θs support with f. This is the baseline for iterations in the Cross We propose a parametric family of IS distribution for Entropy method. We use the same form as fθ because in the Piecewise Mixture Distribution following sections, we use a sampling distribution which is in k the same family as the parametric distribution. (cid:88) f (x)= π˜ exp(θ x−κ(θ ))f (x|γ ≤x<γ ), (35) Weestimateθ∗ bysolvingthestochasticcounterpartof(32) θ i i i i i−1 i i=1 max 1 (cid:88)N I (X ) f(Xn) lnf (X ), (33) where we use exponential change of measure for each piece θ N n=1 ε n fθs(Xi) θ n poafrdaimstertiebrutiisonθa=nd{θad,j.u.s.,tθth,eπ˜pr,o..p.o,πr˜tio}n. parameter to π˜i. The 1 k 1 k wdihsterriebustiaomnpfleθss.{X1,...,XN} are drawn from the sampling datIan,(s3o3)w,ecnsi=mpIlεif(yXtnh)effθfsu(X(nXncnt)i)onisaasknownconstantgiventhe We note that if I (X ) = 0 for all n = 1,..,N in (33), ε n N the objective equals to 0 constantly. To avoid this situation, 1 (cid:88) max c lnf (X ). (36) we select a sampling distribution which emphasizes the rarer θ N n θ n n=1 events. Fig. 3 shows the iteration procedure of the Cross Entropy We split the samples into index sets Si = {j|γi−1 ≤ method. The core part of the Cross Entropy method is to use Xj < γi} for i = 1,...,k for each bounded segment. Since theoptimizeroftheobjectivefunction(33)intheithiteration, fi(Xn|γi−1 ≤ x < γi) (cid:54)= 0 only if n ∈ Si, for each θi and θi∗, as the parameters for the sampling distribution in the next π˜i, the optimization function is equivalent to iteration. The underlying idea is that the IS distribution in 1 (cid:88) distribution family f should better approach the optimal IS max c ln(π˜ exp(θ X −κ(θ )) distribution.Thereforθe,asweiterate,weobtainmore“critical” θi,π˜i N n∈Si n i i n i rareeventsandhaveabetterestimationoftheoptimizerwhich f (X |x<γ ≤x<γ )). (37) i n i−1 i leads to even more “critical” rare events in the next iteration. We define the stopping criterion regarding the parameter or We can further rewrite the optimization function regarding the objective value. In practice, we want to start with an θi and π˜i respectively. For π˜i, we have appropriate sampling distribution to get a good solution with 1 (cid:88) max c lnπ˜, (38) laesssamiteprlaintigond.isSteriebuseticotnio.n V-C1 for a discussion of initializing π˜i N n∈Si n i which obtains an analytical form for the optimizer (cid:80) c 1{n∈S } π˜ = n∈Si n i . (39) i (cid:80) c n∈Si n For θ , we have i 1 (cid:88) max c lnexp(θ X −κ(θ )) θi N n∈Si n i n i f (X |γ ≤x<γ ), (40) Fig.3. IterationsofCrossEntropy. i n i−1 i IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 6 which is an exponential change of measure with D only. 2) Adjusting sample size N: The choice of sample size N i We note that we can simplify this optimization function by should not only depend on the total number of rare events rewriting the log term as obtained in each iteration. For each parameter of interest, we need sufficient non-zero c ’s to guarantee the qualification of 1 (cid:88) n max c (lnexp(θ X −κ(θ ))+ theestimation.Wenotethattheparametersestimationdepend θi N n∈Si n i n i onlyontherareeventinthecorrespondingpiece,soweadjust lnfi(Xn|γi−1 ≤x<γi)), (41) sample size N to ensure that each piece with large portion π˜i contains enough rare event samples. which is equivalent to 3) Setting a lower bound for π˜ : When we update π˜ in i i 1 (cid:88) mθaixN n∈Sicn(θiXn−κ(θi)), (42) e(3v9e)n,tisfacmnp=le 0infothrealplience∈, wSie,hmaevaeniπ˜nig=tha0t.tWhehreeniswneohraavree π˜ = 0, the support of the IS distribution will differ from since the latter term does not depend on θ . i i the original distribution. We note that it might cause bias in For a bounded exponential distribution with parameter λ, our simulation analysis. On the other hand, once π˜ hits 0, it the Cross Entropy iteration solves i will be 0 in the following iterations. Therefore, we need to 1 (cid:88) keep π˜ > 0. Setting a low bound for π˜ , for example, 0.01, max c (θ X − i i θi N n∈Si n i n when there is no rare event for piece i, gives an efficient IS distribution while avoiding the problems. e−(λ−θi)γi−1 −e−(λ−θi)γi ln ). (43) 4) Updating parameter θ : The absence of rare event λ−θ i i samples also leads to failures in updating θ . In this case, i For a bounded normal distribution with parameters µ = 0 we use either the value of θ in the last iteration, or we set i andσ,theoptimizationfunctionfortheCrossEntropyiteration it to 0, i.e. reset the distribution as the real distribution. We is note that we can tolerant some inaccurate estimation if π˜ is i max (cid:88) c X θ −((cid:88) c )(σ2θi2+ small, since a small π˜i indicates that this piece might not be θi n∈Si n n i n∈Si n 2 im5p)orCtahnatntgointhgetrraurnecaevtieonntsγ. : The truncations of the Piece- i Φ(γi−θiσ2)−Φ(γi−1−θiσ2) wise Mixture Distribution are fixed throughout the Cross En- ln σ σ ). (44) Φ(γi)−Φ(γi−1) tropymethod.Thus,ifthereisabadselectionoftruncationin σ σ ouroriginaldistributionmodel,theCrossEntropycannotgive C. Discussion on Numerical Implementation an efficient IS distribution. The changing of truncation points WehavepresentedtheoptimizationfunctionsforCrossEn- is hard to implement by optimization, so we use a heuristic tropy iterations, but we cannot reliably apply these equations approach for adjusting the truncation points to emphasize the in practice without considering some of the problematical tail part of the Piecewise IS distribution. numerical details. In this section, we discuss methods to In any iteration, if the number of rare events is not enough overcome these numerical issues. to properly update the parameters, we check π˜ of the current i 1) Initializing Cross Entropy Iterations for Rare Events: sampling distribution. If the π˜ of the tail piece is the largest k Since rare events occur with small probability, using the possible value, we increase the value of the all truncation originaldistributionassamplingdistributiontostarttheCross points except γ with a certain value. Shifting the truncation 0 Entropy iterations it becomes computationally burdensome gives more weight to the tail part. Then by sampling from to sample a single rare event. One possible approach is to the adjusted distribution, we check if the number of events of initialize with guess of sampling distribution. When we have interest is sufficient. We repeat these actions until we obtain some rough knowledge about the optimal IS distribution, enough rare events in the iteration. we can use the knowledge to construct a proper sampling We propose this heuristic approach, since the flexibility of distribution. the Piecewise Mixture Distribution is not fully exploited if Forcaseswherewehavelittleknowledgeabouttheoptimal we cannot change the truncation points. We note that finding IS distribution, we construct adaptive events that gradually a more systematic procedure to locate the knots remains an reduce the rarity. For rare events denoted by ε, we define open question. the sequence of events to be ε ⊃ ε ⊃ ... ⊃ ε ⊃ ε, 1 2 n where ε is not rare for our initializing sampling density. For 1 VI. SIMULATIONANALYSIS each iteration t, we gradually reduce the rare event set ε and t A. Automated Vehicle Model use ε to replace ε in the objective function. Since ε is a t t subsetofε ,theISdistributionforε alsoprovidesmore First,wepresentourPiecewiseMixtureModelsforR−1and t−1 t−1 chances for samples from ε . We use the optimal solution TTC−1 andthencomparetheresultswiththesingleparamet- t in (t−1)th iteration θ∗ as the sampling parameter θ for ric distribution model used in [11]. For both approaches, we t−1 t the next iteration and choose ε to have a relatively larger divide the data of TTC−1 into three segments regarding the t probability to occur under f . Since ε gradually approaches rangeofv.Sincethethreesegmentsaresimilarindistribution, θt t ε as we iterate, eventually we obtain the optimal parameters we only show the results of the segment for v in the range of for ε. 15 to 25 m/s. IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 7 Fig.4. PiecewiseMixtureDistributionfittingforR−1. Fig.6. ComparisonoffittingforR−1. 1) Piecewise mixture models for R−1 and TTC−1: Fig. 4 shows the fitting of R−1 using two bounded exponential dis- tributionsandthreeboundedexponentialdistributions.Adding one more truncation point provides a better fitting to the body part of distribution while having the same fitting of tail. In Fig. 5, we truncated the data into two parts. For the tail part, we use the exponential distribution. For the body part, the mixture of two normal distributions gives a better fit. The Piecewise Mixture Models enable us to use different distributions for the body part and the tail part. 2) Comparison with single parametric distribution models: Fig.6 and Fig.7 compare the new model and the previous Fig.7. ComparisonoffittingforTTC−1 givenvL between15and25m/s. model. We note that Piecewise Mixture Models provide more flexibility in data fitting. lane change model determines whether a crash happens by checking to see if the value of R, the range between two B. Cross Entropy Results vehicles,reaches0.Meanwhile,theTTC alsogoesto0when Here,weusethelanechangemodeltoexemplifytheCross a crash happens. To construct events less rare than a crash, Entropy method. For the three variables R,TTC,v, the dis- we relax the criterion for crash to be either R hits t > 0 R tribution is f(R,TTC,v)=f(v)f(R)f(TTC|v) where f(v) or TTC hits t > 0. By changing these two thresholds, TTC is the empirical distribution. Since we have three conditional t and t as shown in Fig. 8, we construct the adaptive R TTC distributions of TTC regarding the value of v, we find the rare events sequence for the Cross Entropy iterations. We use IS distributions independently for each case. We present the sample size N =1000 for each iteration. results for v from 5 to 15 m/s. Fig. 9 and 10 show the parameters present in each of the Weassumethatwehavelessinformationabouttherelation iterations. We observe that the parameters stabilize gradually. between the distribution of variables and the rare events. Fig. 11 shows how the distribution changes gradually from Our objective is to construct adaptive rare events to help the original distribution to the IS distribution. We note that us approach the IS distribution. We recall that our original the density moves toward the tail part as we iterate. C. Simulation Results In our simulation experiments, we set the convergence criterion as the relative half-width of 100(1−α)% confidence interval drops below β. In this case, we use α=0.2 and β = 0.2 to study the number of samples needed for convergence. OurgoalistocomparetheefficiencyofthePiecewiseMixture Distribution and single exponential distribution models. Fig. 12 shows that both models give a similar estima- tion as the number of experiments grows large, and that the Piecewise Mixture Distribution model converges slightly faster than the single parametric model. The circles show that the relative half-width of the Piecewise Mixture Distribution Fig.5. PiecewiseMixtureDistributionfittingforTTC−1givenvLbetween 15and25m/s. model reaches the target confidence value after 7800 samples, IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 8 Fig.8. CrossEntropyiterationswithsequenceofeventswiththresholdsfor crash.Weleaveiteration1blanktokeepthex-axisconsistentwithFig.9and 10. whereas the single parametric model needs about 13800 sam- Fig.10. CrossEntropyiterationswithsequenceofeventsofTTC−1 forv from5to15m/s. ples.UsingthePiecewiseMixtureDistributionmodelreduced the sample size by 44%. To reduce stochastic uncertainty, we repeat the tests 10 times and calculate the average. It takes 7840 samples on average to obtain a converged estimation using the Piecewise Mixture Distribution model, whereas it takes 12320 samples on average using the single accelerated distribution model to converge. Table I compares the two models with the crude MonteCarlomethod[24].Weestimatethenumberneededfor convergence of crude Monte Carlo by using the fact that the numberofeventsofinterestoccurringisBinomialdistributed. We compute the standard deviation of the crude Monte Carlo estimation Pˆ(x∈ε) by (cid:115) Fig.11. DistributionchangethroughCrossEntropyiterationswithsequence Pˆ(x∈ε)(1−Pˆ(x∈ε)) ofeventsofTTC−1 forv from5to15m/s. std(Pˆ(x∈ε))= , (45) n which allows us to estimate where zα/2 is the (1−α/2) quantile of normal distribution. WecalculatetherequiredsamplesizeN ofcrudeMonteCarlo z2 (1−Pˆ(x∈ε)) in Table I from an estimation Pˆ(x ∈ ε) = 7.4×10−7 with Nˆ = α/2 , (46) β2Pˆ(x∈ε) 80% confidence interval (7.0×10−7,7.8×10−7). Finally, we apply the heuristic approach in Section V-C5 to thedatasegmentwithv from5to15m/s.Werunsimulations with this segment and compare the results with the standard Fig.9. CrossEntropyiterationswithsequenceofeventsofR−1 forvfrom Fig.12. Estimationofcrashprobabilityforonelanechangeusingpiecewise 5to15m/s. andsingleaccelerateddistributions. IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 9 TABLEI TABLEII NUMBEROFSAMPLES(N)NEEDEDTOCONVERGE. COMPARISONOFTHECOMPUTATIONTIMEBETWEENSINGLE PARAMETRICMODELANDPIECEWISEMODEL. Piecewise Single Crude N 7840 12320 5.5×107 Stages Crude Single Piecewise RatiotoPiecewise 1 1.57 7×103 Fitting - 4parametersto 18parametersto estimate estimate 30,000simulations 24,000simulations CrossEntropy - 4parameters 18parameters 5.5×107 12,320 7840 Simulation simulations simulations simulations the same confidence level compared to the single parametric Model. APPENDIXA INVERSECDFOFPIECEWISEMIXTUREDISTRIBUTIONS We can sample from Piecewise Mixture Distribution by the inverse CDF approach. Here, we derive the inverse CDF for Piecewise Mixture Distribution. Fig. 13. Relative half-width of crash probability estimation for one lane changewithleadingvehicle’sspeedinrangeof5to15m/s,comparingsingle, The CDF of Piecewise Mixture Distribution (5) can split piecewiseandheuristicaccelerateddistributions. into  ...  approach for the ewise Mixture Distribution and single para- F(x)= (cid:80)i−1π +π F (x|γ ≤x<γ ) γ ≤x<γ . j=1 j i i i−1 i i−1 i metric distribution models. Fig. 13 shows the convergence of ... confidence half-width. We note that the relative half-width of (47) the heuristic, which is smaller than the standard approach for Therefore the inverse function can be written as the Piecewise Mixture Distribution model, indicates that the  ... latter model’s performance can be further improved. F−1(y)=F−1(y−(cid:80)ij−=11πj|γ ≤x<γ) (cid:80)i−1π ≤y<(cid:80)i π . ...i πi i−1 i j=1 j j=1 j VII. CONCLUSIONS (48) Thispaperproposedanewmodelforacceleratedevaluation where Fi−1 is the inverse conditional CDF of Fi. Below, we of AVs. The Piecewise Mixture Distribution Models provide give two example of inverse conditional CDF. more accurate fitting to the surrounding human-controlled FortheinverseCDFofconditionalexponentialdistribution, vehicle behaviors than the single parametric model used we have in the literature. The proposed model was more efficient F−1(y|F (γ )≤y <F (γ ))= and reduced the evaluation time by almost half than single θ θ 1 θ 2 parametric model. The Cross Entropy procedure described in Fθ−1((Fθ(γ2)−Fθ(γ1))y+Fθ(γ1)), (49) this paper effectively worked in this scenario analysis. We whereF andF−1aretheCDFandinverseCDFofexponential provided practical solutions to deal with the numerical issues distribution. which occurred while calculating the optimal parameters. The For conditional normal distribution, the inverse CDF is heuristic approach exploited the flexibility of the Piecewise Mixture Distribution structure. Testing the proposed model F−1(y|F (γ )≤y <F (γ ))= θ θ 1 θ 2 on a large dataset of cut-in crashes caused by improper lane γ −θσ2 γ −θσ2 changes, the Piecewise Mixture Distribution model reduced σΦ−1((Φ( 2 )−Φ( 1 ))y+ σ σ the simulation cases by about 33% compared with the single γ −θσ2 parametric model under the same convergence requirement. Φ( 1 ))+θσ2. (50) σ Moreover, the proposed model was 7000 times faster than the Crude Monte Carlo method. Table II summarizes the comparison of the computation efforts between the models. We note that using the Piecewise APPENDIXB Mixture Distribution model increases the number of parame- EMALGORITHMFORMIXTUREBOUNDEDNORMAL ters estimated, where the estimation of parameters is almost DISTRIBUTION instant.Inthe CrossEntropystage,the numberofsimulations required for the Piecewise model is not significantly less than Here, we present a numerical MLE algorithm with mixture thesingleparametricmodel,becauseweassumenoknowledge bounded normal distribution. The steps are as follows. about the optimal IS distribution for the Piecewise model. ALGORITHM: Overall,thePiecewisemodelneedsfewersimulationstoreach 1) Initialize {p ,σ }, j =1,...,m. j j IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 10 2) E step: update [16] R. Y. Rubinstein, “Rare Event Simulation via Cross-entropy and Im- portance Sampling,” in Second International Workshop on Rare Event p f (X |σ ) τj = j j n j . (51) Simulation,1999,pp.1–17. n (cid:80)m p f (X |σ ) [17] J. Aldrich, “RA Fisher and the Making of Maximum j=1 j j n j Likelihood 1912-1922,” Statistical Science, 1997. [Online]. Available: 3) M step: update http://projecteuclid.org/euclid.ss/1030037906 (cid:80)N τj [18] D.CoxandD.Hinkley,TheoreticalStatistics,1979. p = n=1 n (52) [19] S.BoydandL.Vandenberghe,ConvexOptimization,2004. j N [20] A. Dempster, N. Laird, and D. Rubin, “Maximum Likelihood from IncompleteDataviatheEMAlgorithm,”Journaloftheroyalstatistical and society,1977. σ (Φ(γ2)−Φ(γ1)) [21] J. Bucklew, Introduction to Rare Event Simulation. Springer Science σ =argmin−τjln j σj σj . (53) &BusinessMedia,2004. j σj n φ(Xσjn) [[2223]] VD..VPa.pKnirko,esSet,atRis.ticYa.lRLuebairnnsitneginT,haenodryP,.1W99.8G.lynn, The Cross-Entropy 4) Repeat 2 and 3 until L(θ|D) converges. MethodforEstimation. ElsevierB.V.,2013,vol.31. [24] S. Asmussen and P. Glynn, Stochastic Simulation: Algorithms and Analysis. Springer,2007. APPENDIXC [25] P.-T. d. Boer, “A Tutorial on the Cross-Entropy Method,” pp. 19–67, 2005. VANILLACROSSENTROPYMETHOD ALGORITHM:[25] 1) Initialize θ . s 2) Sample {X ,...,X } from f and update 1 N θs θ =argmax 1 (cid:88)N I (X ) f(Xi) lnf (X ). (54) Zprhei-ycuanadnidHatueanPgh.DZ.hisytuuadnenHtuinangIndisusatrisaelcoannddyOepa-r θ N i=1 ε i fθs(Xi) θ i erationsEngineeringattheUniversityofMichigan, AnnArbor.Hisresearchinterestsincludesimulation 3) Update andstochasticoptimization. θ =θ. (55) s 4) Repeat 2 and 3 until θ “converges”. REFERENCES [1] GoogleAutoLLC,“Monthlyreports?GoogleSelf-DrivingCarProject.” [Online].Available:https://www.google.com/selfdrivingcar/reports/ [2] “PreliminaryReport,?HighwayHWY16FH018.” [3] A.Evan,“FatalTeslaSelf-DrivingCarCrashRemindsUsThatRobots Aren’tPerfect,”IEEESpectrum,2016. DingZhaoDingZhaoreceivedthePh.D.degreein [4] FESTA-Consortium,“FESTAHandbookVersion2DeliverableT6.4of 2016 from theUniversity of Michigan, AnnArbor. theFieldopErationalteStsupporTAction,”FESTA,Tech.Rep.,2008. HeiscurrentlyaResearchFellowintheUniversity [5] NHTSA,“TrafficSafetyFacts2014,”DOT,Tech.Rep.,2014.[Online]. of Michigan Transportation Research Institute. His Available:http://www-nrd.nhtsa.dot.gov/Pubs/811620.pdf\nhttp://www- research interest includes evaluation of connected nrd.nhtsa.dot.gov/Pubs/809778.pdf and automated vehicles, vehicle dynamic control, [6] H.PengandD.Leblanc,“EvaluationofthePerformanceandSafetyof driverbehaviorsmodeling,andbigdataanalysis. AutomatedVehicles,”WhitePap.NSFTransp.CPSWork,2012. [7] M.Aust,“EvaluationProcessforActiveSafetyFunctions:Addressing KeyChallengesinFunctional,FormativeEvaluationofAdvancedDriver AssistanceSystems,”2012. [8] Y. Kou, “Development and Evaluation of Integrated Chassis Control Systems,”2010. [9] A.UngorenandH.Peng,“AnAdaptiveLateralPreviewDriverModel,” VehicleSystemDynamics,vol.43,no.4,pp.245–259,42005. [10] W. Ma and H. Peng, “A Worst-case Evaluation Method for Dynamic Systems,”Journalofdynamicsystems,,1999. HenryLamHenryLamreceivedtheB.S.degreein [11] D.Zhao,H.Lam,H.Peng,S.Bao,D.J.LeBlanc,K.Nobukawa,and actuarialsciencefromtheUniversityofHongKong C. S. Pan, “Accelerated Evaluation of Automated Vehicles Safety in in2005,andtheA.M.andPh.D.degreesinstatistics Lane-Change Scenarios Based on Importance Sampling Techniques,” from Harvard University, Cambridge, in 2006 and IEEETransactionsonIntelligentTransportationSystems,2016. 2011. [12] D.Zhao,H.Peng,S.Bao,K.Nobukawa,D.J.LeBlanc,andC.S.Pan, From2011to2014,hewasanAssistantProfessor “Accelerated evaluation of automated vehicles using extracted natural- in the Department of Mathematics and Statistics istic driving data,” in Proceeding for 24th International Symposium of at Boston University. Since 2015, he has been an VehiclesonRoadandTracks,2015. Assistant Professor in the Department of Indus- [13] D.Zhao,X.Huang,H.Peng,H.Lam,andD.J.Leblanc,“Accelerated trial and Operations Engineering at the University EvaluationofAutomatedVehiclesinCar-FollowingManeuvers,”ArXiv, of Michigan, Ann Arbor. His research focuses on p.12,2016. stochastic simulation, risk analysis, and simulation optimization. Dr. Lam’s [14] Z.Huang,D.Zhao,H.Lam,D.J.LeBlanc,andH.Peng,“Evaluation works have been funded by National Science Foundation and National of Automated Vehicles in the Frontal Cut-in Scenario - an Enhanced Security Agency. He has also received an Honorable Mention Prize in the Approach using Piecewise Mixture Models,” 10 2016. [Online]. Institute for Operations Research and Management Sciences (INFORMS) Available:http://arxiv.org/abs/1610.09450 GeorgeNicholsonBestStudentPaperAward,andFinalistinINFORMSJunior [15] D. Bezzina and J. R. Sayer, “Safety Pilot: Model Deployment Test FacultyInterestGroupBestPaperCompetition. Conductor Team Report,” NHTSA, Tech. Rep. June, 2014. [Online]. Available:http://safetypilot.umtri.umich.edu/

ebook img

Accelerated Evaluation of Automated Vehicles Using Piecewise Mixture Models PDF

1.4 MB·
by  Zhiyuan Huang
#journals #arxiv
Save to my drive
Quick download
Download
Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview Accelerated Evaluation of Automated Vehicles Using Piecewise Mixture Models

IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 1 Accelerated Evaluation of Automated Vehicles Using Piecewise Mixture Models Zhiyuan Huang1 , Ding Zhao2, Henry Lam1, David J. LeBlanc2 Abstract—The process to certify highly Automated Vehicles both Google and Tesla update their software throughout the hasnotyetbeendefinedbyanycountryintheworld.Currently, process, which may have improved safety, but the newest companies test Automated Vehicles on public roads, which is version of the AV has not accumulated that many miles as time-consuming and inefficient. We proposed the Accelerated they have claimed. In summary, today’s best practice adopted 7 Evaluation concept, which uses a modified statistics of the 1 surrounding vehicles and the Importance Sampling theory to by the industry is time-consuming and inefficient. A better 0 reduce the evaluation time by several orders of magnitude, approach is needed. 2 while ensuring the evaluation results are statistically accurate. In this paper, we further improve the accelerated evaluation n concept by using PiecewiseMixture Distribution models, instead a A. Related Researches of Single Parametric Distribution models. We developed and J applied this idea to forward collision control system reacting BesidestheN-FOT,thetestmatrixapproach[6],[7]andthe 1 to vehicles making cut-in lane changes. The behavior of the worst-casescenariosapproach[8],[9],[10]aretwoalternative 3 cut-in vehicles was modeled based on more than 403,581 lane methods for vehicle evaluation. Our approach follows the changes collected by the University of Michigan Safety Pilot ] ModelDeploymentProgram.Simulationresultsconfirmthatthe AcceleratedEvaluationconceptweproposed[11]toprovidea Y accuracy and efficiency of the Piecewise Mixture Distribution brand-new alternative. The basic concept is that as high-level S method outperformed single parametric distribution methods in AVs just began to penetrate the market, they mainly interact . accuracy and efficiency, and accelerated the evaluation process s withhuman-controlledvehicles(HVs).Thereforewefocuson c by almost four orders of magnitude. modeling the interaction between the AV and the HV around [ it. The evaluation procedure involves four steps: 1 Index Terms—automated vehicles, testing, evaluation, safety v • Modelthebehaviorsoftheprimaryothervehicles(POVs) 5 represented by f(x) as the major disturbance to the AV I. INTRODUCTION 1 using large-scale naturalistic driving data 9 IT is critical to thoroughly and rigorously test and evaluate • Skew the disturbance statistics from f(x) to modified 8 an Automated Vehicle (AV) before its release. Recent statistics f∗(x) (accelerated distribution) to generate 0 crashes involving a Google self-driving car [1] and a Tesla more frequent and intense interactions between AVs and . 1 Autopilotvehicle[2]attractedthepublic’sattentiontoAVtest- POVs 0 ing and evaluation. While these AVs are generally considered • Conduct accelerated tests with f∗(x) 7 as industrial leaders, because they use public road for testing, 1 • Use the Importance Sampling (IS) theory to skew back statistically they have not yet accumulated enough miles. The : the results to understand real-world behavior and safety v TeslaAutopilot,inparticular,wascriticizedforbeingreleased benefits i X too early in the hands of the general public [3]. This approach has been successfully applied to evaluate Currently, there are no standards or protocols to test AVs r AVs in the frontal crash with a cut-in vehicle [11] and also a at automation level 2 or higher. Many companies adopt the frontalcrashwithaleadvehicle[12],[13].Thisapproachwas Naturalistic Field Operational Tests (N-FOT) approach [4]. confirmed to significantly reduce the evaluation time while However, this method is inefficient because safety critical accurately preserving the statistical behavior of the AV-HV scenarios rarely happen in daily driving. The Google Self- interaction. In the previous studies, the evaluation time was driving cars accumulated 1.9 million driving. This distance, reduced by two to five orders of magnitudes - the accelerated although sounds a lot, provides limited exposure to critical rate depends on the test scenarios, where rarer events achieve events, given that U.S. drivers encounter a police reported higher accelerated rate. The non-accelerated models and the crash every five hundred thousand miles on average and fatal accelerated models were built based on signal component crash every one hundred million miles [5]. In the meantime, distributions. While this method does benefit from its simple *This work was funded by the Mobility Transformation Center at the mathematical form, it has a few drawbacks as illustrated in University of Michigan with grant no. N021552. Z. Huang and D. Zhao Fig. 1. i) The fitting of the rare events (usually the tail part contributedequallytotheresearch. 1Zhiyuan Huang ([email protected]) and Henry Lam of the statistical distributions) would be dominated by the ([email protected])areintheDepartmentofIndustrialandOperations fitting of the normal driving behaviors (the majority part of EngineeringattheUniversityofMichigan. the distributions), which may induce large errors. ii) The full 2Ding Zhao (corresponding author: [email protected]) and potential in higher accelerated rate is not achieved due to the DavidLeBlanc([email protected])areintheUniversityofMichi- ganTransportationResearchInstitute. lack of flexibility of the modified accelerated models. IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 2 Fig. 1. Acceleration evaluation based on single parametric distribution and PiecewiseMixtureDistribution. Fig.2. LanechangedatacollectedbySPMDvehicle. B. Contribution In this paper, we proposed a more general framework for independently with regard to the value of v . By comparing L the Accelerated Evaluation method to overcome the afore- among 17 types of commonly used distribution templates, we mentioned limitations based on Piecewise Mixture Distribu- selected the Pareto distribution to model R−1 and used the L tion Models as illustrated in Fig. 1 b). In this paper, we exponential distribution for TTC−1 segments. L implemented the Accelerated Evaluation method under the Using the empirical distribution of v and parametric dis- L new framework. Comparing to our previous work [14], we tributions of R and TTC , we drew values from these L L thoroughly discuss the Cross Entropy method with proposed distributions as inputs to simulate the AV-HV interaction. framework in this paper. We present practical tips to over- The outcome from the simulation can be considered as an come numerical issues and reduce computational efforts. We event indicator function I (x) that returns {1,0} depending ε demonstratethismethodbyevaluatingthelongitudinalcontrol on the event of interest. Given the stochastic distribution of system reacting to vehicles making cut-in lane changes. the variables and the event indicator function, we obtained the optimal exponential distribution for Importance Sampling C. Paper Structure by implementing the Cross Entropy method [16]. As we have shown in Fig. 1 a), we used only single parametric distribu- Section II will introduce the lane change model based on tions.Inthenextsection,weintroduceournewapproachusing single parametric distributions. In Section III, we present the Piecewise Mixture Distributions. new lane change model with Piecewise Mixture Distributions. We establish the Accelerated Evaluation in Section IV and discusstheCrossEntropymethodwithPiecewiseMixtureDis- III. LANECHANGEMODELWITHPIECEWISEMIXTURE tributionmodelsinSectionV.Simulationresultsarediscussed DISTRIBUTIONS in Section VI. Section VII concludes this paper. Although many commonly used parametric distributions have concise and elegant forms, they do not always describe II. ACCELERATEDEVALUATIONWITHSINGLE the data distribution well. Instead, a better fitting can be PARAMETRICDISTRIBUTIONS achieved by dividing the dataset into several subsets. We estimatethemodelparametersusingtheMaximumLikelihood The lane change events were extracted from the Safety Estimation (MLE) [17] in each subset. The general process of Pilot Model Deployment (SPMD) database [15]. With over MLE is as follow. 2 million miles of vehicle driving data collected from 98 Assume we have a family of distribution with Cumula- cars over 3 years, we identify 403,581 lane change events. tive Distribution Function (CDF) F(x|θ), where θ is the Previously [11], we used 173,692 events with a negative parametervectorofF.ThecorrespondingProbabilityDensity range rate to build a statistical model focusing on three key Function (PDF) of F is f(x|θ). Assuming that data D = variables that captured the effects of gap acceptance of the {X ,X ,...,X } is independently and identically distributed lanechangingvehicle:velocityoftheleadvehicle(v ),range 1 2 N L and the distribution is in the family of F(x|θ), we want to totheleadvehicle(RL)andtimetocollision(TTCL).TTCL find the most “likely” parameter θˆ. was defined as: R We define the likelihood function [18] as TTC =− L, (1) L R˙L L(θ|D)=P(D|θ)=ΠNn=1f(Xn|θ). (2) where R˙L is the relative speed. We call the estimation of θˆ that maximizes the likelihood The modeling of these three variables was hard to handle function the mostly likely estimation MLE. because of dependency, so we simplified it based on a crucial For computation convenience, we introduce the log- observation. Although TTC is dependent on v generally, L L likelihood function we split the data into 3 segments: v at 5 to 15 m/s, 15 to 25 L N m/sand25to35m/s.Withineachsegment,R isindependent (cid:88) L L(θ|D)=lnL(θ|D)= lnf(X |θ). (3) withv andTTC .ThisallowedustomodelTTC andR n L L L L n=1 IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 3 Since the logarithm is monotone, the log-likelihood function 1) MLEforboundedexponentialdistribution: Thebounded preservestheoptimizeroftheoriginalfunction.[19]Theopti- exponential distribution with rate θ has the form mizeroflog-likelihoodfunction,θˆ,istheMLEofdistribution θe−θx family F. We have the MLE as f(x|γ ≤x<γ )= (10) 1 2 e−θγ1 −e−θγ2 θˆ=argmax L(θ|D). (4) for γ ≤x<γ . θ 1 2 In the following, we describe the Piecewise Mixture Dis- For dataset D ={X1,...,XN}, the log-likelihood function tribution fitting concept based on MLE and we present the is bounded distribution fitting results. All optimization problems N (cid:88) presented in this section are tractable and can be solved by L(D|θ)= lnθ−θX −ln(e−θγ1 −e−θγ2), (11) n fminunc in MATLAB. n=1 where L is concave over θ. Although we cannot solve the A. General Framework of the Piecewise Mixture Distribution maximization analytically, it is solvable through numerical Lane Change Model methods. We define Piecewise Mixture Distribution to be distribution Therefore, the MLE of θ is given by the optimization with CDF in the form of N (cid:88) (cid:88)k max Nlnθ−Nln(e−θγ1 −e−θγ2)− θXn. (12) F(x)= π F (x|γ ≤x<γ ), (5) θ i θi i−1 i n=1 i=1 2) MLE for bounded normal distribution: Consider a where k is the number of truncation, (cid:80)ki=1πi = 1, and bounded normal distribution with mean 0 and variance θ2 Fi(x|γi−1 ≤x<γi)istheconditionalcumulativedistribution conditional on 0≤γ1 ≤x<γ2. The PDF is function, meaning that F (γ |γ ≤ x < γ ) = 0 and i i−1 i−1 i 1φ(x) Fi(γi|γi−1 ≤ x < γi) = 1 for i = 1,...,k. θi denotes the f(x|γ ≤x<γ )= θ θ . (13) parameter(s) for F . We can consider that π = P(γ ≤ 1 2 Φ(γ2)−Φ(γ1) i i i−1 θ θ x < γ ) and when x ≥ 0, we have γ = 0 and γ = ∞. By i 0 k The MLE of the bounded normal distribution is given by this definition, the PDF of the Piecewise Mixture Distribution is k (cid:80)N X2 γ γ f(x)=(cid:88)π f (x|γ ≤x<γ ). (6) max− n=1 n −Nlnθ−Nln(Φ( 2)−Φ( 1)). (14) i θi i−1 i θ 2θ2 θ θ i=1 3) Fitting mixture model with EM algorithm: Compared to In our case, θ = {π ,...,π ,θ ,...,θ }. Splitting D into 1 k 1 k single parametric distributions, mixture distribution combines pieces regarding the truncation points {γ ,...,γ }, gives 1 k−1 several classes of distribution and thus is more flexible. We data index sets S = {j|γ ≤ X < γ } for i = 1,...,k. i i−1 j i consider the fitting problem of mixture bounded normal dis- We can write the log-likelihood function as tribution. L(θ|D)=(cid:80)k (cid:80) lnπ The PDF of mixture of m bounded normal distribution can +(cid:80)ki=1(cid:80) n∈Silnf i(X |γ ≤x<γ ). (7) be written as i=1 n∈Si θi n i−1 i m We obtain the MLE of θ can be obtained by maximizing (cid:88) f(x|γ ≤x<γ )= p f (x|γ ≤x<γ ) (15) L(θ|D) over θ. Since L is concave over π , we take 1 2 j j 1 2 i j=1 ∂L =0 (8) where f is bounded Gaussian distribution with mean ∂π j i 0 and variance σ2. The parameters here are θ = j and get {p ,...,p ,σ2,...,σ2 }. We want to find MLE of p and σ2 πˆ =|S |/N. (9) 1 m 1 m j j i i for j =1,...,m. For parameters θi in Fi, it is known (7) to be the same The log-likelihood function for data D ={Xn}Nn=1 is as computing the MLE of θ with corresponding dataset i N m Di = {X|γi−1 ≤ X < γi and X ∈ D}. Since we L(θ|D)= (cid:88)ln(cid:88)pjfj(Xn|γ1 ≤x<γ2). (16) use bounded distribution for each F , below we explain the i n=1 j=1 estimationofparametersforthethreedistributionsweapplied We note that this is hard to solve directly, because there in later sections. is a sum within the log function. Therefore, we apply the TosamplefromaPiecewiseMixtureDistribution,wecould Expectation-Maximization (EM) [20] algorithm to find the use the inverse function approach. See Appendix A for the optimizer, i.e. MLE, for the parameters. details. We define Zj to denote whether or not the random number n X comesfrommixturedistributionj,j =1,...,m,andZj = B. Bounded Distribution n n {0,1}. We also introduce the expectation Wedevelopthreeboundeddistributionsandusetheminthe lane change model. E[Zj|X ]:=τj. (17) n n n IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 4 The EM algorithm starts with initial parameters {p ,σ }, with smaller variance. This is known as Importance Sampling j j j = 1,...,m. For data D = {X }N , we set com- [21] and F∗ is the IS distribution. n n=1 plete data as D = {X ,Z }N . The EM algorithm op- For estimating P(X ∈ ε),we note that an optimal IS c n n n=1 timizes E[L(θ|D )|D] in every step. The E step updates distribution c E[L(θ|D )|D], and the M step optimizes this function. The P(X ≤x, ε) c F∗∗(x)=F(x|ε)= (25) algorithm iterates E step and M step until reaching the P(x∈ε) convergence criterion. couldreducethevarianceofISestimationto0,buttheoptimal In our case, requires the knowledge of P(X ∈ ε). However, it guides the E[L(θ|D )|D]=(cid:80)N (cid:80)m τj(lnp +lnf (X )). selection of the IS distribution. c n=1 j=1 n j j n (18) B. Exponential Change of Measure SinceobjectiveE[l (θ|D )|D]intheMstepisconcaveover c c p and σ , we could maximize the objective function through Exponential change of measure is commonly used to con- j j an analytic approach for p : structF∗.Althoughtheexponentialchangeofmeasurecannot j guarantee convergence to optimal distribution, it is easy to (cid:80)N τj p = n=1 n. (19) implement and the new distribution generally stays within the j N same class of distribution. For σ , we can solve the following maximization problem Exponential change of measure distribution takes the form j through numerical approach. of f (x)=exp(θx−κ(θ))f(x), (26) (cid:18) (cid:19) θ X σj =argmσijn−τnjlnσj +τnjlnφ σjn − where θ is the change of measure parameter and κ(θ) is (cid:18) (cid:19) the log-moment generating function of original distribution f. γ γ τjln Φ( 2)−Φ( 1) . (20) When θ =0, we have f (x)=f(x). n σ σ θ j j For a bounded exponential distribution, the exponential See Appendix B for the full EM algorithm. change of measure distribution is (λ−θ)e−(λ−θ)x IV. ACCELERATEDEVALUATIONWITHIMPORTANCE fθ(x|γ1 ≤x<γ2)= e−(λ−θ)γ1 −e−(λ−θ)γ2, (27) SAMPLING whereλistheparameterforexponentialdistribution.Wenote Importance Sampling (IS) is thus used to accelerate the that f is still a bounded exponential distribution and λ = θ θ evaluationprocess,becausecrudeMonteCarlosimulationsfor λ−θ. rare events can be time-consuming. Below we describe the IS For a bounded normal distribution, the exponential change method. of measure distribution is 1φ(x−σ2θ) A. Important Sampling and Optimal IS distribution fθ(x|γ1 ≤x<γ2)= Φ(γ2−θσσ2)−σΦ(γ1−θσ2), (28) Let x be a random variable generated from distribution F, σ σ where the original distribution truncated from a normal distri- and ε ⊂ Ω where ε is the rare event of interest and Ω is the bution with parameters µ=0 and σ. We note that the change sample space. Our objective is to estimate of measure distribution is still a bounded normal distribution (cid:90) P(X ∈ε)=E[I (X)]= I (x)dF (21) with µ=θσ2 and σ. ε ε V. CROSSENTROPYMETHODANDIMPLEMENTATION where (cid:40) 1 x∈ε, Section IV discussed optimal IS distribution F∗∗ providing I (x)= (22) ε 0 variance estimation to the value of interest, whereas this 0 otherwise. section describes the Cross Entropy method used to estimate We can write the evaluation of rare events as the sample the “optimal” parameters θ, which minimizes the “distance” mean of Iε(x) between a parametric distribution Fθ and F∗∗ without know- ing F∗∗. The description below is based on the Piecewise N Pˆ(X ∈ε)= 1 (cid:88)I (X ), (23) Mixture Distribution structure. N ε n n=1 A. Introduction where X ’s are drawn from distribution F. i Since we have TheCrossEntropy,whichisalsoknownasKullback-Leibler (cid:90) (cid:90) dF distance [22], measures the similarity between distributions. E[Iε(X)]= Iε(x)dF = Iε(x)dF∗dF∗, (24) We define the Cross Entropy between function g and h as g(X) (cid:90) dwiestrcibauntiocnomFpu∗t,ewthheichsamhapslethmeeasnamoef sIuεp(Xpo)rtddFFw∗ithoveFr,thtoe D(g,h)=Eg[lnf(X)]= g(x)lng(x)dx− (cid:90) obtain an unbiased estimation of P(X ∈ε). By appropriately g(x)lnh(x)dx. (29) selecting F∗, the evaluation procedure obtains an estimation IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 5 From (25), we know that the PDF of the optimal IS We note that if we have two independent variables where distribution F∗∗ is f(x,y) = f(x)f(y), we can take a parametric distribution f∗∗(x)= Iε(x)f(x). (30) fΘor=ea{cθh,vθar}i.abTlheeaonbdjehcativveeffuΘn(cxt,ioyn)c=orrfeθs1p(oxn)dfiθn2g(yt)o,(w33h)eries P(x∈ε) 1 2 Since P(x ∈ ε) is generally unavailable, we use a para- N mWeetrwicadnitsttroibfiuntidonthFeθptaoraamppertoearcθh∗ththeaotpmtiminaiml IiSzedsistthreibuCtrioosns. mθax N1 n(cid:88)=1Iε(Xn,Yn)ffΘ(sX(Xnn,Y,Ynn))(lnfθ1(Xn)+ Entropy [23] between f∗∗ and fθ. We denote θ∗ as the lnfθ2(Yn)), (34) optimal parameter for the parametric distribution. Then the which can be decoupled into two optimization problem over minimization problem θ and θ respectively and I (X ,Y ) f(Xn,Yn) is a known minD(fθ,f∗∗) (31) co1nstant2given {Xn,Yn}. ε n n fΘs(Xn,Yn) θ We implement the Cross Entropy on the Piecewise Mixture is equivalent to Distribution with one variable. We note that we can apply the f(X) results to the lane change model, since the Cross Entropy ob- max E [I (X) lnf (X)], (32) θ θs ε fθs(X) θ jective function of independent variables can be implemented in (34). where f denotes the sampling distribution with parameters θs θ . We note that this is a generalized setting, since we can s use any sampling distribution f as long as it has the same B. OptimizationFunctionforPiecewiseMixtureDistributions θs support with f. This is the baseline for iterations in the Cross We propose a parametric family of IS distribution for Entropy method. We use the same form as fθ because in the Piecewise Mixture Distribution following sections, we use a sampling distribution which is in k the same family as the parametric distribution. (cid:88) f (x)= π˜ exp(θ x−κ(θ ))f (x|γ ≤x<γ ), (35) Weestimateθ∗ bysolvingthestochasticcounterpartof(32) θ i i i i i−1 i i=1 max 1 (cid:88)N I (X ) f(Xn) lnf (X ), (33) where we use exponential change of measure for each piece θ N n=1 ε n fθs(Xi) θ n poafrdaimstertiebrutiisonθa=nd{θad,j.u.s.,tθth,eπ˜pr,o..p.o,πr˜tio}n. parameter to π˜i. The 1 k 1 k wdihsterriebustiaomnpfleθss.{X1,...,XN} are drawn from the sampling datIan,(s3o3)w,ecnsi=mpIlεif(yXtnh)effθfsu(X(nXncnt)i)onisaasknownconstantgiventhe We note that if I (X ) = 0 for all n = 1,..,N in (33), ε n N the objective equals to 0 constantly. To avoid this situation, 1 (cid:88) max c lnf (X ). (36) we select a sampling distribution which emphasizes the rarer θ N n θ n n=1 events. Fig. 3 shows the iteration procedure of the Cross Entropy We split the samples into index sets Si = {j|γi−1 ≤ method. The core part of the Cross Entropy method is to use Xj < γi} for i = 1,...,k for each bounded segment. Since theoptimizeroftheobjectivefunction(33)intheithiteration, fi(Xn|γi−1 ≤ x < γi) (cid:54)= 0 only if n ∈ Si, for each θi and θi∗, as the parameters for the sampling distribution in the next π˜i, the optimization function is equivalent to iteration. The underlying idea is that the IS distribution in 1 (cid:88) distribution family f should better approach the optimal IS max c ln(π˜ exp(θ X −κ(θ )) distribution.Thereforθe,asweiterate,weobtainmore“critical” θi,π˜i N n∈Si n i i n i rareeventsandhaveabetterestimationoftheoptimizerwhich f (X |x<γ ≤x<γ )). (37) i n i−1 i leads to even more “critical” rare events in the next iteration. We define the stopping criterion regarding the parameter or We can further rewrite the optimization function regarding the objective value. In practice, we want to start with an θi and π˜i respectively. For π˜i, we have appropriate sampling distribution to get a good solution with 1 (cid:88) max c lnπ˜, (38) laesssamiteprlaintigond.isSteriebuseticotnio.n V-C1 for a discussion of initializing π˜i N n∈Si n i which obtains an analytical form for the optimizer (cid:80) c 1{n∈S } π˜ = n∈Si n i . (39) i (cid:80) c n∈Si n For θ , we have i 1 (cid:88) max c lnexp(θ X −κ(θ )) θi N n∈Si n i n i f (X |γ ≤x<γ ), (40) Fig.3. IterationsofCrossEntropy. i n i−1 i IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 6 which is an exponential change of measure with D only. 2) Adjusting sample size N: The choice of sample size N i We note that we can simplify this optimization function by should not only depend on the total number of rare events rewriting the log term as obtained in each iteration. For each parameter of interest, we need sufficient non-zero c ’s to guarantee the qualification of 1 (cid:88) n max c (lnexp(θ X −κ(θ ))+ theestimation.Wenotethattheparametersestimationdepend θi N n∈Si n i n i onlyontherareeventinthecorrespondingpiece,soweadjust lnfi(Xn|γi−1 ≤x<γi)), (41) sample size N to ensure that each piece with large portion π˜i contains enough rare event samples. which is equivalent to 3) Setting a lower bound for π˜ : When we update π˜ in i i 1 (cid:88) mθaixN n∈Sicn(θiXn−κ(θi)), (42) e(3v9e)n,tisfacmnp=le 0infothrealplience∈, wSie,hmaevaeniπ˜nig=tha0t.tWhehreeniswneohraavree π˜ = 0, the support of the IS distribution will differ from since the latter term does not depend on θ . i i the original distribution. We note that it might cause bias in For a bounded exponential distribution with parameter λ, our simulation analysis. On the other hand, once π˜ hits 0, it the Cross Entropy iteration solves i will be 0 in the following iterations. Therefore, we need to 1 (cid:88) keep π˜ > 0. Setting a low bound for π˜ , for example, 0.01, max c (θ X − i i θi N n∈Si n i n when there is no rare event for piece i, gives an efficient IS distribution while avoiding the problems. e−(λ−θi)γi−1 −e−(λ−θi)γi ln ). (43) 4) Updating parameter θ : The absence of rare event λ−θ i i samples also leads to failures in updating θ . In this case, i For a bounded normal distribution with parameters µ = 0 we use either the value of θ in the last iteration, or we set i andσ,theoptimizationfunctionfortheCrossEntropyiteration it to 0, i.e. reset the distribution as the real distribution. We is note that we can tolerant some inaccurate estimation if π˜ is i max (cid:88) c X θ −((cid:88) c )(σ2θi2+ small, since a small π˜i indicates that this piece might not be θi n∈Si n n i n∈Si n 2 im5p)orCtahnatntgointhgetrraurnecaevtieonntsγ. : The truncations of the Piece- i Φ(γi−θiσ2)−Φ(γi−1−θiσ2) wise Mixture Distribution are fixed throughout the Cross En- ln σ σ ). (44) Φ(γi)−Φ(γi−1) tropymethod.Thus,ifthereisabadselectionoftruncationin σ σ ouroriginaldistributionmodel,theCrossEntropycannotgive C. Discussion on Numerical Implementation an efficient IS distribution. The changing of truncation points WehavepresentedtheoptimizationfunctionsforCrossEn- is hard to implement by optimization, so we use a heuristic tropy iterations, but we cannot reliably apply these equations approach for adjusting the truncation points to emphasize the in practice without considering some of the problematical tail part of the Piecewise IS distribution. numerical details. In this section, we discuss methods to In any iteration, if the number of rare events is not enough overcome these numerical issues. to properly update the parameters, we check π˜ of the current i 1) Initializing Cross Entropy Iterations for Rare Events: sampling distribution. If the π˜ of the tail piece is the largest k Since rare events occur with small probability, using the possible value, we increase the value of the all truncation originaldistributionassamplingdistributiontostarttheCross points except γ with a certain value. Shifting the truncation 0 Entropy iterations it becomes computationally burdensome gives more weight to the tail part. Then by sampling from to sample a single rare event. One possible approach is to the adjusted distribution, we check if the number of events of initialize with guess of sampling distribution. When we have interest is sufficient. We repeat these actions until we obtain some rough knowledge about the optimal IS distribution, enough rare events in the iteration. we can use the knowledge to construct a proper sampling We propose this heuristic approach, since the flexibility of distribution. the Piecewise Mixture Distribution is not fully exploited if Forcaseswherewehavelittleknowledgeabouttheoptimal we cannot change the truncation points. We note that finding IS distribution, we construct adaptive events that gradually a more systematic procedure to locate the knots remains an reduce the rarity. For rare events denoted by ε, we define open question. the sequence of events to be ε ⊃ ε ⊃ ... ⊃ ε ⊃ ε, 1 2 n where ε is not rare for our initializing sampling density. For 1 VI. SIMULATIONANALYSIS each iteration t, we gradually reduce the rare event set ε and t A. Automated Vehicle Model use ε to replace ε in the objective function. Since ε is a t t subsetofε ,theISdistributionforε alsoprovidesmore First,wepresentourPiecewiseMixtureModelsforR−1and t−1 t−1 chances for samples from ε . We use the optimal solution TTC−1 andthencomparetheresultswiththesingleparamet- t in (t−1)th iteration θ∗ as the sampling parameter θ for ric distribution model used in [11]. For both approaches, we t−1 t the next iteration and choose ε to have a relatively larger divide the data of TTC−1 into three segments regarding the t probability to occur under f . Since ε gradually approaches rangeofv.Sincethethreesegmentsaresimilarindistribution, θt t ε as we iterate, eventually we obtain the optimal parameters we only show the results of the segment for v in the range of for ε. 15 to 25 m/s. IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 7 Fig.4. PiecewiseMixtureDistributionfittingforR−1. Fig.6. ComparisonoffittingforR−1. 1) Piecewise mixture models for R−1 and TTC−1: Fig. 4 shows the fitting of R−1 using two bounded exponential dis- tributionsandthreeboundedexponentialdistributions.Adding one more truncation point provides a better fitting to the body part of distribution while having the same fitting of tail. In Fig. 5, we truncated the data into two parts. For the tail part, we use the exponential distribution. For the body part, the mixture of two normal distributions gives a better fit. The Piecewise Mixture Models enable us to use different distributions for the body part and the tail part. 2) Comparison with single parametric distribution models: Fig.6 and Fig.7 compare the new model and the previous Fig.7. ComparisonoffittingforTTC−1 givenvL between15and25m/s. model. We note that Piecewise Mixture Models provide more flexibility in data fitting. lane change model determines whether a crash happens by checking to see if the value of R, the range between two B. Cross Entropy Results vehicles,reaches0.Meanwhile,theTTC alsogoesto0when Here,weusethelanechangemodeltoexemplifytheCross a crash happens. To construct events less rare than a crash, Entropy method. For the three variables R,TTC,v, the dis- we relax the criterion for crash to be either R hits t > 0 R tribution is f(R,TTC,v)=f(v)f(R)f(TTC|v) where f(v) or TTC hits t > 0. By changing these two thresholds, TTC is the empirical distribution. Since we have three conditional t and t as shown in Fig. 8, we construct the adaptive R TTC distributions of TTC regarding the value of v, we find the rare events sequence for the Cross Entropy iterations. We use IS distributions independently for each case. We present the sample size N =1000 for each iteration. results for v from 5 to 15 m/s. Fig. 9 and 10 show the parameters present in each of the Weassumethatwehavelessinformationabouttherelation iterations. We observe that the parameters stabilize gradually. between the distribution of variables and the rare events. Fig. 11 shows how the distribution changes gradually from Our objective is to construct adaptive rare events to help the original distribution to the IS distribution. We note that us approach the IS distribution. We recall that our original the density moves toward the tail part as we iterate. C. Simulation Results In our simulation experiments, we set the convergence criterion as the relative half-width of 100(1−α)% confidence interval drops below β. In this case, we use α=0.2 and β = 0.2 to study the number of samples needed for convergence. OurgoalistocomparetheefficiencyofthePiecewiseMixture Distribution and single exponential distribution models. Fig. 12 shows that both models give a similar estima- tion as the number of experiments grows large, and that the Piecewise Mixture Distribution model converges slightly faster than the single parametric model. The circles show that the relative half-width of the Piecewise Mixture Distribution Fig.5. PiecewiseMixtureDistributionfittingforTTC−1givenvLbetween 15and25m/s. model reaches the target confidence value after 7800 samples, IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 8 Fig.8. CrossEntropyiterationswithsequenceofeventswiththresholdsfor crash.Weleaveiteration1blanktokeepthex-axisconsistentwithFig.9and 10. whereas the single parametric model needs about 13800 sam- Fig.10. CrossEntropyiterationswithsequenceofeventsofTTC−1 forv from5to15m/s. ples.UsingthePiecewiseMixtureDistributionmodelreduced the sample size by 44%. To reduce stochastic uncertainty, we repeat the tests 10 times and calculate the average. It takes 7840 samples on average to obtain a converged estimation using the Piecewise Mixture Distribution model, whereas it takes 12320 samples on average using the single accelerated distribution model to converge. Table I compares the two models with the crude MonteCarlomethod[24].Weestimatethenumberneededfor convergence of crude Monte Carlo by using the fact that the numberofeventsofinterestoccurringisBinomialdistributed. We compute the standard deviation of the crude Monte Carlo estimation Pˆ(x∈ε) by (cid:115) Fig.11. DistributionchangethroughCrossEntropyiterationswithsequence Pˆ(x∈ε)(1−Pˆ(x∈ε)) ofeventsofTTC−1 forv from5to15m/s. std(Pˆ(x∈ε))= , (45) n which allows us to estimate where zα/2 is the (1−α/2) quantile of normal distribution. WecalculatetherequiredsamplesizeN ofcrudeMonteCarlo z2 (1−Pˆ(x∈ε)) in Table I from an estimation Pˆ(x ∈ ε) = 7.4×10−7 with Nˆ = α/2 , (46) β2Pˆ(x∈ε) 80% confidence interval (7.0×10−7,7.8×10−7). Finally, we apply the heuristic approach in Section V-C5 to thedatasegmentwithv from5to15m/s.Werunsimulations with this segment and compare the results with the standard Fig.9. CrossEntropyiterationswithsequenceofeventsofR−1 forvfrom Fig.12. Estimationofcrashprobabilityforonelanechangeusingpiecewise 5to15m/s. andsingleaccelerateddistributions. IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 9 TABLEI TABLEII NUMBEROFSAMPLES(N)NEEDEDTOCONVERGE. COMPARISONOFTHECOMPUTATIONTIMEBETWEENSINGLE PARAMETRICMODELANDPIECEWISEMODEL. Piecewise Single Crude N 7840 12320 5.5×107 Stages Crude Single Piecewise RatiotoPiecewise 1 1.57 7×103 Fitting - 4parametersto 18parametersto estimate estimate 30,000simulations 24,000simulations CrossEntropy - 4parameters 18parameters 5.5×107 12,320 7840 Simulation simulations simulations simulations the same confidence level compared to the single parametric Model. APPENDIXA INVERSECDFOFPIECEWISEMIXTUREDISTRIBUTIONS We can sample from Piecewise Mixture Distribution by the inverse CDF approach. Here, we derive the inverse CDF for Piecewise Mixture Distribution. Fig. 13. Relative half-width of crash probability estimation for one lane changewithleadingvehicle’sspeedinrangeof5to15m/s,comparingsingle, The CDF of Piecewise Mixture Distribution (5) can split piecewiseandheuristicaccelerateddistributions. into  ...  approach for the ewise Mixture Distribution and single para- F(x)= (cid:80)i−1π +π F (x|γ ≤x<γ ) γ ≤x<γ . j=1 j i i i−1 i i−1 i metric distribution models. Fig. 13 shows the convergence of ... confidence half-width. We note that the relative half-width of (47) the heuristic, which is smaller than the standard approach for Therefore the inverse function can be written as the Piecewise Mixture Distribution model, indicates that the  ... latter model’s performance can be further improved. F−1(y)=F−1(y−(cid:80)ij−=11πj|γ ≤x<γ) (cid:80)i−1π ≤y<(cid:80)i π . ...i πi i−1 i j=1 j j=1 j VII. CONCLUSIONS (48) Thispaperproposedanewmodelforacceleratedevaluation where Fi−1 is the inverse conditional CDF of Fi. Below, we of AVs. The Piecewise Mixture Distribution Models provide give two example of inverse conditional CDF. more accurate fitting to the surrounding human-controlled FortheinverseCDFofconditionalexponentialdistribution, vehicle behaviors than the single parametric model used we have in the literature. The proposed model was more efficient F−1(y|F (γ )≤y <F (γ ))= and reduced the evaluation time by almost half than single θ θ 1 θ 2 parametric model. The Cross Entropy procedure described in Fθ−1((Fθ(γ2)−Fθ(γ1))y+Fθ(γ1)), (49) this paper effectively worked in this scenario analysis. We whereF andF−1aretheCDFandinverseCDFofexponential provided practical solutions to deal with the numerical issues distribution. which occurred while calculating the optimal parameters. The For conditional normal distribution, the inverse CDF is heuristic approach exploited the flexibility of the Piecewise Mixture Distribution structure. Testing the proposed model F−1(y|F (γ )≤y <F (γ ))= θ θ 1 θ 2 on a large dataset of cut-in crashes caused by improper lane γ −θσ2 γ −θσ2 changes, the Piecewise Mixture Distribution model reduced σΦ−1((Φ( 2 )−Φ( 1 ))y+ σ σ the simulation cases by about 33% compared with the single γ −θσ2 parametric model under the same convergence requirement. Φ( 1 ))+θσ2. (50) σ Moreover, the proposed model was 7000 times faster than the Crude Monte Carlo method. Table II summarizes the comparison of the computation efforts between the models. We note that using the Piecewise APPENDIXB Mixture Distribution model increases the number of parame- EMALGORITHMFORMIXTUREBOUNDEDNORMAL ters estimated, where the estimation of parameters is almost DISTRIBUTION instant.Inthe CrossEntropystage,the numberofsimulations required for the Piecewise model is not significantly less than Here, we present a numerical MLE algorithm with mixture thesingleparametricmodel,becauseweassumenoknowledge bounded normal distribution. The steps are as follows. about the optimal IS distribution for the Piecewise model. ALGORITHM: Overall,thePiecewisemodelneedsfewersimulationstoreach 1) Initialize {p ,σ }, j =1,...,m. j j IEEETRANSACTIONSONINTELLIGENTTRANSPORTATIONSYSTEMS,VOL.XX,NO.XX,DECEMBER2016 10 2) E step: update [16] R. Y. Rubinstein, “Rare Event Simulation via Cross-entropy and Im- portance Sampling,” in Second International Workshop on Rare Event p f (X |σ ) τj = j j n j . (51) Simulation,1999,pp.1–17. n (cid:80)m p f (X |σ ) [17] J. Aldrich, “RA Fisher and the Making of Maximum j=1 j j n j Likelihood 1912-1922,” Statistical Science, 1997. [Online]. Available: 3) M step: update http://projecteuclid.org/euclid.ss/1030037906 (cid:80)N τj [18] D.CoxandD.Hinkley,TheoreticalStatistics,1979. p = n=1 n (52) [19] S.BoydandL.Vandenberghe,ConvexOptimization,2004. j N [20] A. Dempster, N. Laird, and D. Rubin, “Maximum Likelihood from IncompleteDataviatheEMAlgorithm,”Journaloftheroyalstatistical and society,1977. σ (Φ(γ2)−Φ(γ1)) [21] J. Bucklew, Introduction to Rare Event Simulation. Springer Science σ =argmin−τjln j σj σj . (53) &BusinessMedia,2004. j σj n φ(Xσjn) [[2223]] VD..VPa.pKnirko,esSet,atRis.ticYa.lRLuebairnnsitneginT,haenodryP,.1W99.8G.lynn, The Cross-Entropy 4) Repeat 2 and 3 until L(θ|D) converges. MethodforEstimation. ElsevierB.V.,2013,vol.31. [24] S. Asmussen and P. Glynn, Stochastic Simulation: Algorithms and Analysis. Springer,2007. APPENDIXC [25] P.-T. d. Boer, “A Tutorial on the Cross-Entropy Method,” pp. 19–67, 2005. VANILLACROSSENTROPYMETHOD ALGORITHM:[25] 1) Initialize θ . s 2) Sample {X ,...,X } from f and update 1 N θs θ =argmax 1 (cid:88)N I (X ) f(Xi) lnf (X ). (54) Zprhei-ycuanadnidHatueanPgh.DZ.hisytuuadnenHtuinangIndisusatrisaelcoannddyOepa-r θ N i=1 ε i fθs(Xi) θ i erationsEngineeringattheUniversityofMichigan, AnnArbor.Hisresearchinterestsincludesimulation 3) Update andstochasticoptimization. θ =θ. (55) s 4) Repeat 2 and 3 until θ “converges”. REFERENCES [1] GoogleAutoLLC,“Monthlyreports?GoogleSelf-DrivingCarProject.” [Online].Available:https://www.google.com/selfdrivingcar/reports/ [2] “PreliminaryReport,?HighwayHWY16FH018.” [3] A.Evan,“FatalTeslaSelf-DrivingCarCrashRemindsUsThatRobots Aren’tPerfect,”IEEESpectrum,2016. DingZhaoDingZhaoreceivedthePh.D.degreein [4] FESTA-Consortium,“FESTAHandbookVersion2DeliverableT6.4of 2016 from theUniversity of Michigan, AnnArbor. theFieldopErationalteStsupporTAction,”FESTA,Tech.Rep.,2008. HeiscurrentlyaResearchFellowintheUniversity [5] NHTSA,“TrafficSafetyFacts2014,”DOT,Tech.Rep.,2014.[Online]. of Michigan Transportation Research Institute. His Available:http://www-nrd.nhtsa.dot.gov/Pubs/811620.pdf\nhttp://www- research interest includes evaluation of connected nrd.nhtsa.dot.gov/Pubs/809778.pdf and automated vehicles, vehicle dynamic control, [6] H.PengandD.Leblanc,“EvaluationofthePerformanceandSafetyof driverbehaviorsmodeling,andbigdataanalysis. AutomatedVehicles,”WhitePap.NSFTransp.CPSWork,2012. [7] M.Aust,“EvaluationProcessforActiveSafetyFunctions:Addressing KeyChallengesinFunctional,FormativeEvaluationofAdvancedDriver AssistanceSystems,”2012. [8] Y. Kou, “Development and Evaluation of Integrated Chassis Control Systems,”2010. [9] A.UngorenandH.Peng,“AnAdaptiveLateralPreviewDriverModel,” VehicleSystemDynamics,vol.43,no.4,pp.245–259,42005. [10] W. Ma and H. Peng, “A Worst-case Evaluation Method for Dynamic Systems,”Journalofdynamicsystems,,1999. HenryLamHenryLamreceivedtheB.S.degreein [11] D.Zhao,H.Lam,H.Peng,S.Bao,D.J.LeBlanc,K.Nobukawa,and actuarialsciencefromtheUniversityofHongKong C. S. Pan, “Accelerated Evaluation of Automated Vehicles Safety in in2005,andtheA.M.andPh.D.degreesinstatistics Lane-Change Scenarios Based on Importance Sampling Techniques,” from Harvard University, Cambridge, in 2006 and IEEETransactionsonIntelligentTransportationSystems,2016. 2011. [12] D.Zhao,H.Peng,S.Bao,K.Nobukawa,D.J.LeBlanc,andC.S.Pan, From2011to2014,hewasanAssistantProfessor “Accelerated evaluation of automated vehicles using extracted natural- in the Department of Mathematics and Statistics istic driving data,” in Proceeding for 24th International Symposium of at Boston University. Since 2015, he has been an VehiclesonRoadandTracks,2015. Assistant Professor in the Department of Indus- [13] D.Zhao,X.Huang,H.Peng,H.Lam,andD.J.Leblanc,“Accelerated trial and Operations Engineering at the University EvaluationofAutomatedVehiclesinCar-FollowingManeuvers,”ArXiv, of Michigan, Ann Arbor. His research focuses on p.12,2016. stochastic simulation, risk analysis, and simulation optimization. Dr. Lam’s [14] Z.Huang,D.Zhao,H.Lam,D.J.LeBlanc,andH.Peng,“Evaluation works have been funded by National Science Foundation and National of Automated Vehicles in the Frontal Cut-in Scenario - an Enhanced Security Agency. He has also received an Honorable Mention Prize in the Approach using Piecewise Mixture Models,” 10 2016. [Online]. Institute for Operations Research and Management Sciences (INFORMS) Available:http://arxiv.org/abs/1610.09450 GeorgeNicholsonBestStudentPaperAward,andFinalistinINFORMSJunior [15] D. Bezzina and J. R. Sayer, “Safety Pilot: Model Deployment Test FacultyInterestGroupBestPaperCompetition. Conductor Team Report,” NHTSA, Tech. Rep. June, 2014. [Online]. Available:http://safetypilot.umtri.umich.edu/

See more

The list of books you might like

Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.
DMCA & Copyright: Dear all, most of the website is community built, users are uploading hundred of books everyday, which makes really hard for us to identify copyrighted material, please contact us if you want any material removed.
Copyright @ PDFdrive.to, 2022 | We love our users