交通流中的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）加速：),1min(max v v v n n +→

规则（1）反映了司机倾向于以尽可能大的速度行驶的特点。

（2）减速：),min(n n n d v v →

规则(2)确保车辆不会与前车发生碰撞。

（3）随机慢化： 以随机概率p 进行慢化，令：)0,

1-min(n n v v → 规则(3)引入随机慢化来体现驾驶员的行为差异，这样既可以反映随机加速行为，又可以反映减速过程中的过度反应行为。这一规则也是堵塞自发产生的至关重要因素。

（4）位置更新：n n n v x v +→ ，车辆按照更新后的速度向前运动。 其中n v ，n x 分别表示第n 辆车位置和速度；l （l ≥1）为车辆长度；11--=+n n n x x d 表示n 车和前车n+1之间空的元胞数；p 表示随机慢化概率；max v 为最大速度。

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.1 NaSch模型运行图

图3.1.2 NaSch模型

3.2流量密度分析

图3.2描述了交通流量与密度的关系，从图中可知，该模型中，当密度为0——0.185时，流量随密度的增加而增加；当密度超过0.185时，流量开始随密度的增加而下降。

图3.2 基于NaSch模型的流量密度图

3.3 NaSch模型时空图分析

图3.3.1和图3.3.2描述了，时间步从11001开始到11500结束，共500个时间步的空间和时间的关系，从图中可以模拟出自发产生的堵塞现象。

图3.3.1 基于NaSch模型的时空图

图3.3.2 基于NaSch模型的时空图

4 模型评价

优点：该程序基本实现了NaSch模型的基本功能，并且最大速度、元胞数量、车辆数量以及运行间隔时间都可以修改，程序很灵活，并且可以清晰的看出每一次运行过程。

缺点：当时间步超过20000步时，内存占用量大。

附件

% 主程序：NaSch_3.m程序代码

% 单车道最大速度3个元胞开口边界条件加速减速随机慢化clf

clear all

%build the GUI

%define the plot button

plotbutton=uicontrol('style','pushbutton',...

'string','Run', ...

'fontsize',12, ...

'position',[100,400,50,20], ...

'callback', 'run=1;');

%define the stop button

erasebutton=uicontrol('style','pushbutton',...

'string','Stop', ...

'fontsize',12, ...

'position',[100,500,50,20], ...

'callback','freeze=1;');

%define the Quit button

quitbutton=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]);

%CA setup

n=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')

axis equal

axis tight

stop=0; %wait for a quit button push

run=0; %wait for a draw

freeze=0; %wait for a freeze（冻结）

while (stop==0 & x<11502)

if(run==1)

%边界条件处理，搜素首末车，控制进出，使用开口条件

a=searchleadcar(cells);

b=searchlastcar(cells);

[cells,v]=border_control(cells,a,b,v,vmax);

i=searchleadcar(cells); %搜索首车位置

for j=1:i

if i-j+1==n

[z,v]=leadcarupdate(z,v);

continue;

else

%======================================加速、减速、随机慢化

if cells(i-j+1)==0; %判断当前位置是否非空

continue;

else v(i-j+1)=min(v(i-j+1)+1,vmax); %加速

%=================================减速

k=searchfrontcar((i-j+1),cells); %搜素前方首个非空元胞位置

if k==0; %确定于前车之间的元胞数

d=n-(i-j+1);

else d=k-(i-j+1)-1;

end

v(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;

end

end

end

cells=z;

memor_cells(x,:)=cells; %记录速度和车辆位置

memor_v(x,:)=v;

x=x+1;

set(imh,'cdata',cells) %更新图像

%update the step number diaplay

pause(0.0001);

stepnumber = 1+str2num(get(number,'string'));

set(number,'string',num2str(stepnumber))

end

if (freeze==1)

run = 0;

freeze = 0;

end

drawnow

end

figure(1)

for l=11001:1:11500

for k=1:1:1000

if memor_cells(l,k)>0

plot(k,l,'k.');

hold on;

end

end

end

xlabel('空间位置')

ylabel('时间（s）')

title('时空图')

for i=1:1:500

density(i)=sum(memor_cells(i,:)>0)/1000;

flow(i)=sum(memor_v(i,:))/1000;

end

figure(2)

plot(density,flow,'k.');

title('流量密度图')

xlabel('density')

ylabel('flow')

%

% /////////////////////////////////////////////////////////////////////// %

%

% 函数：searchlastcar.m程序代码

function [location_lastcar]=searchlastcar(matrix_cells)

%搜索尾车位置

for i=1:length(matrix_cells)

if matrix_cells(i)~=0

location_lastcar=i;

break;

else %如果路上无车，则空元胞数设定为道路长度

location_lastcar=length(matrix_cells);

end

end

% 函数：searchfrontcar.m程序代码

function [location_frontcar]=searchfrontcar(current_location,matrix_cells)

i=length(matrix_cells);

if current_location==i

location_frontcar=0;

else

for j=current_location+1:i

if matrix_cells(j)~=0

location_frontcar=j;

break;

else

location_frontcar=0;

end

end

end

% 函数：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));

for i=1:n

j=z(i);

if j==0

matrix_cells(j)=0;

else

matrix_cells(j)=1;

end

end

matrix_cells_start=matrix_cells;

% 函数：randslow.m程序代码

function [new_v]=randslow(v)

p=0.3; %慢化概率

rand('state',sum(100*clock)*rand(1));%?¨?????ú??×?

p_rand=rand; %产生随机概率

if p_rand<=p

v=max(v-1,0);

end

new_v=v;

% 函数：leadcarrupdate.m程序代码

function [new_matrix_cells,new_=leadcarupdate(matrix_cells,v) %第一辆车更新规则

n=length(matrix_cells);

if v(n)~=0

matrix_cells(n)=0;

v(n)=0;

end

new_matrix_cells=matrix_cells;

new_v=v;

% 函数：searchleadcar.m程序代码

function [location_leadcar]=searchleadcar(matrix_cells)

i=length(matrix_cells);

for j=1:i

if matrix_cells(i-j+1)~=0

location_leadcar=i-j+1;

break;

else

location_leadcar=0;

end

end

% 函数：speadstart.m程序代码

function [v_matixcells]=speedstart(matrix_cells,vmax)

%道路初始状态车辆速度初始化

v_matixcells=zeros(1,length(matrix_cells));

for i=1:length(matrix_cells)

if matrix_cells(i)~=0

v_matixcells(i)=round(vmax*rand(1));

end

end