元胞自动机NaSch模型及其MATLAB代码

元胞自动机NaSch模型及其MATLAB代码作业要求根据前面的介绍，对NaSch模型编程并进行数值模拟：模型参数取值：Lroad=1000，p=0.3，Vmax=5。边界条件：周期性边界。数据统计：扔掉前50000个时间步，对后50000个时间步进行统计,需给出的结果。基本图（流量-密度关系）：需整个密度范围内的。时空图（横坐标为空间，纵坐标为时间，密度和文献中时空图保持一致,画500个时间步即可）。指出NaSch模型的创新之处，找出NaSch模型的不足，并给出自己的改进思路。流量计算方法：密度＝车辆数/路长；流量flux=density×V_ave。在道路的某处设置虚拟探测计算统计时间T内通过的车辆数N；流量flux=N/T。在计算过程中可都使用无量纲的变量。1、NaSch模型的介绍作为对184号规则的推广，Nagel和Schreckberg在1992年提出了一个模拟车辆交通的元胞自动机模型，即NaSch模型（也有人称它为NaSch模型）。时间、空间和车辆速度都被整数离散化。道路被划分为等距离的离散的格子，即元胞。每个元胞或者是空的，或者被一辆车所占据。车辆的速度可以在（0～Vmax）之间取值。2、NaSch模型运行规则在时刻t到时刻t+1的过程中按照下面的规则进行更新：（1）加速：vnmin(vn1,vmax)规则（1）反映了司机倾向于以尽可能大的速度行驶的特点。（2）减速：vnmin(vn,dn)规则(2)确保车辆不会与前车发生碰撞。（3）随机慢化：以随机概率p进行慢化，令：vnmin(vn-1,0)规则(3)引入随机慢化来体现驾驶员的行为差异，这样既可以反映随机加速行为，又可以反映减速过程中的过度反应行为。这一规则也是堵塞自发产生的至关重要因素。（4）位置更新：，车辆按照更新后的速度向前运动。其中vn，xn分别表示第n辆车位置和速度；l（l≥1）为车辆长度；p表示随机慢化概率；表示n车和前车n+1之间空的元胞数；vmax为最大速度。3、NaSch模型实例根据题目要求，模型参数取值：L=1000，p=0.3，Vmax=5，用matlab软件进行编程，扔掉前11000个时间步，统计了之后500个时间步数据，得到如下基本图和时空图。3.1程序简介初始化：在路段上，随机分配200个车辆，且随机速度为1-5之间。图3.1.1是程序的运行图，图3.1.2中，白色表示有车，黑色是元胞。图3.1.1NaSch模型运行图图3.1.2NaSch模型3.2流量密度分析图3.2描述了交通流量与密度的关系，从图中可知，该模型中，当密度为0——0.185时，流量随密度的增加而增加；当密度超过0.185时，流量开始随密度的增加而下降。图3.2基于NaSch模型的流量密度图3.3NaSch模型时空图分析图3.3.1和图3.3.2描述了，时间步从11001开始到11500结束，共500个时间步的空间和时间的关系，从图中可以模拟出自发产生的堵塞现象。图3.3.1基于NaSch模型的时空图图3.3.2基于NaSch模型的时空图4模型评价优点：该程序基本实现了NaSch模型的基本功能，并且最大速度、元胞数量、车辆数量以及运行间隔时间都可以修改，程序很灵活，并且可以清晰的看出每一次运行过程。缺点：当时间步超过20000步时，内存占用量大。附件%主程序：NaSch_3.m程序代码%单车道最大速度3个元胞开口边界条件加速减速随机慢化clfclearall%buildtheGUI%definetheplotbuttonplotbutton=uicontrol('style','pushbutton',...'string','Run',...'fontsize',12,...'position',[100,400,50,20],...'callback','run=1;');%definethestopbuttonerasebutton=uicontrol('style','pushbutton',...'string','Stop',...'fontsize',12,...'position',[100,500,50,20],...'callback','freeze=1;');%definetheQuitbuttonquitbutton=uicontrol('style','pushbutton',...'string','Quit',...'fontsize',12,...'position',[100,600,50,20],...'callback','stop=1;close;');number=uicontrol('style','text',...'string','1',...'fontsize',12,...'position',[20,400,50,20]);%CAsetupn=1000;%数据初始化z=zeros(1,n);%元胞个数z=roadstart(z,200);%道路状态初始化，路段上随机分布200辆cells=z;vmax=5;%最大速度v=speedstart(cells,vmax);%速度初始化x=1;%记录速度和车辆位置memor_cells=zeros(3600,n);memor_v=zeros(3600,n);imh=imshow(cells);%初始化图像白色有车，黑色空元胞set(imh,'erasemode','none')axisequalaxistightstop=0;%waitforaquitbuttonpushrun=0;%waitforadrawfreeze=0;%waitforafreeze（冻结）while(stop==0&x11502)if(run==1)%边界条件处理，搜素首末车，控制进出，使用开口条件a=searchleadcar(cells);b=searchlastcar(cells);[cells,v]=border_control(cells,a,b,v,vmax);i=searchleadcar(cells);%搜索首车位置forj=1:iifi-j+1==n[z,v]=leadcarupdate(z,v);continue;else%======================================加速、减速、随机慢化ifcells(i-j+1)==0;%判断当前位置是否非空continue;elsev(i-j+1)=min(v(i-j+1)+1,vmax);%加速%=================================减速k=searchfrontcar((i-j+1),cells);%搜素前方首个非空元胞位置ifk==0;%确定于前车之间的元胞数d=n-(i-j+1);elsed=k-(i-j+1)-1;endv(i-j+1)=min(v(i-j+1),d);%==============================%减速%随机慢化v(i-j+1)=randslow(v(i-j+1));new_v=v(i-j+1);%======================================加速、减速、随机慢化%更新车辆位置z(i-j+1)=0;z(i-j+1+new_v)=1;%更新速度v(i-j+1)=0;v(i-j+1+new_v)=new_v;endendendcells=z;memor_cells(x,:)=cells;%记录速度和车辆位置memor_v(x,:)=v;x=x+1;set(imh,'cdata',cells)%更新图像%updatethestepnumberdisplaypause(0.0001);stepnumber=1+str2num(get(number,'string'));set(number,'string',num2str(stepnumber))endif(freeze==1)run=0;freeze=0;enddrawnowendfigure(1)forl=11001:1:11500fork=1:1:1000ifmemor_cells(l,k)0plot(k,l,'k.');holdon;endendendxlabel('空间位置')ylabel('时间（s）')title('时空图')fori=1:1:500density(i)=sum(memor_cells(i,:)0)/1000;flow(i)=sum(memor_v(i,:))/1000;endfigure(2)plot(density,flow,'k.');title('流量密度图')xlabel('density')ylabel('flow')%%///////////////////////////////////////////////////////////////////////%%%函数：searchlastcar.m程序代码function[location_lastcar]=searchlastcar(matrix_cells)%搜索尾车位置fori=1:length(matrix_cells)ifmatrix_cells(i)~=0location_lastcar=i;break;else%如果路上无车，则空元胞数设定为道路长度location_lastcar=length(matrix_cells);endend%函数：searchfrontcar.m程序代码function[location_frontcar]=searchfrontcar(current_location,matrix_cells)i=length(matrix_cells);ifcurrent_location==ilocation_frontcar=0;elseforj=current_location+1:iifmatrix_cells(j)~=0location_frontcar=j;break;elselocation_frontcar=0;endendend%函数：roadstart.m程序代码function[matrix_cells_start]=roadstart(matrix_cells,n)%道路上的车辆初始化状态，元胞矩阵随机为0或1，matrix_cells初始矩阵，n初始车辆数k=length(matrix_cells);z=round(k*rand(1,n));fori=1:nj=z(i);ifj==0matrix_cells(j)=0;elsematrix_cells(j)=1;endendmatrix_cells_start=matrix_cells;%函数：randslow.m程序代码function[new_v]=randslow(v)p=0.3;%慢化概率rand('state',sum(100*clock)*rand(1));%?¨?????ú??×?p_rand=rand;%产生随机概率ifp_rand=pv=max(v-1,0);endnew_v=v;%函数：leadcarrupdate.m程序代码function[new_matrix_cells,new_=leadcarupdate(matrix_cells,v)%第一辆车更新规则n=length(matrix_cells);ifv(n)~=0matrix_cells(n)=0;v(n)=0;endnew_matrix_cells=matrix_cells;new_v=v;%函数：searchleadcar.m程序代码function[location_leadcar]=searchleadcar(matrix_cells)i=length(matrix_cells);forj=1:iifmatrix_cells(i-j+1)~=0location_leadcar=i-j+1;break;elselocation_leadcar=0;e