学习常用模型及算法：2.元胞自动机和交通模型

元胞自动机在数模中的应用十分广泛。

元胞自动机原理

一维元胞自动机。给出任意一个状态，都能知道下一时刻的状态。规则已给出，总共有2^3 = 8种可能。

二维元胞自动机——生命游戏

元胞的状态是有限的——不是有车就是无车，不是死就是活。

可以用简单的规则来模拟复杂的问题。

元胞自动机的构成要素

左图可用来模拟疾病传播情况

最常用的就是正方型网格，其次是六边型网格。三角型网格几乎不出现。

最常用的为左数两种邻居，但具体按照题目要求决定。

边界处理总共有四种类型。其中汽车交通就采取的是吸收型（汽车经过以后就会消失不见）

右侧采取的是周期型，相当于将相同的网格拼在一起。

总和型是合法型的一种特殊情况

第一个模型：森林火灾模型

火灾在森林中蔓延的模型。元胞总共有三种状态：空格、树、火。

其中如果树的任一邻居是火（或者树遭雷劈），它的下一阶段就会变成火。再下一阶段就变成了空格。会以极小概率再次变成树。

系统稳定条件：着火的密度和由空格转化而成的树的密度应该是相等的。

时间尺度分离条件：遭雷劈的概率<<空格转化成树的概率<<树被火烧完的时间尺度

用切片重拼的方法表示上下左右邻居。

这里用RGB分别表示着火、树和空格。

% simulate forest fire with cellular automata
% zhou lvwen: zhou.lv.wen@gmail.com
% August 15 2010 n = 300;
Plight = 5e-6;
Pgrowth = 1e-2;
UL = [n 1:n-1];
DR = [2:n 1];
veg=zeros(n,n);
imh = image(cat(3,veg,veg,veg));
% veg = empty=0 burning=1 green=2
for i=1:3000%nearby fires?sum =            (veg(UL,:)==1) + ...(veg(:,UL)==1)     +      (veg(:,DR)==1) + ...(veg(DR,:)==1);veg = 2*(veg==2) - ...( (veg==2) & (sum>0 | (rand(n,n)<Plight)) ) + ...2*((veg==0) & rand(n,n)<Pgrowth) ;set(imh, 'cdata', cat(3,(veg==1),(veg==2),zeros(n)) )drawnow
end

第二个模型：交通模型

交通模型，这是比较重要的一个模型。

其中最大密度的情况是前车的尾部与后车的头部紧挨，此时每l的长度就有一辆车。

流量的大小为密度x车速。密度为0时，路上没有车，流量自然为0。密度到达上限时，路段完全拥堵，流量也为0。

所以流量与密度是类似抛物线的关系。

如图，以deltaN为中介，可推出流量对x的导数与密度对时间的导数的和为0。

在前面简单交通规则的基础上，加了一些额外的细节。

每个车的初始速度不尽相同；加速时至少比原先速度快一个单位，但要小于规定的最大速度；速度不能快到与前车相撞；减速时至少比原先速度小一个单位，但要在0以上；以及位置更新规则。

临界密度的情况下，每6个格子中有一辆车，车间距为5，每个车刚好可以以最大速度vmax=5行驶而不用担心相撞。

代码分析与实现

ns.m求平均流量；gaplength是用来求两辆车的距离（车间距+1）

12行和19行都是对越界情况的处理。

main.m

% This the main script implements the Nagel Schreckenberg cellular automata
% based traffic model. Flow vs density curves are displayed.
%
% zhou lvwen: zhou.lv.wen@gmail.com
%
density = 0:0.02:1;
roadlength = 100;
vmax = 5;
tmax = 200;
pbrak = 0;
flux = [];
vmean = [];for rho = density[R, J, V] = ns(rho, pbrak, roadlength, tmax, 0, 0);flux = [flux; J];vmean = [vmean; V];
end% ------------------------- density vs. volecity --------------------------
figure
plot(density, vmean,'k.','markersize',15);
hold on
plot(density,min(vmax, 1./density-1),'-r','linewidth',2)
ylim([0,5.55])
legend({'Cellular automata aproach', ...'$v(\rho) = \min\{v_{\max}, 1/\rho-1\}$'}, ...'interpreter','latex')
xlabel('density in vehicles/cell')
ylabel('velocity in cell/time')% --------------------------- density vs. flux ----------------------------
figure
plot(density, flux,'k.','markersize',15);
hold on;
plot(density,min(density*vmax, 1-density),'-r','linewidth',2)
legend({'Cellular automata aproach', ...'$J(\rho) = \min\{\rho\cdot v_{\max}, 1-\rho\}$'}, ...'interpreter','latex')
xlabel('density in vehicles/cell')
ylabel('flux in vehicles/time')

ns.m

function [rho, flux, vmean] = ns(rho, p, L, tmax, animation, spacetime)
%
% NS: This script implements the Nagel Schreckenberg cellular automata based
% traffic model. Car move forward governed by NS algorithm:
%
%   1. Acceleration. If the vehicle can speed up without hitting the speed
%      limit vmax it will add one to its velocity, vn 鈫� vn + 1. Otherwise,
%      the vehicle has constant speed, vn 鈫� vn.
%
%   2. Collision prevention. If the distance between the vehicle and the car
%      ahead of it, dn, is less than or equal to vn , i.e. the nth vehicle
%      will collide if it doesn't slow down, then vn 鈫� dn 鈭� 1.
%
%   3. Random slowing. Vehicles often slow for non-traffic reasons (cell
%      phones, coffee mugs, even laptops) and drivers occasionally make
%      irrational choices. With some probability pbrake, vn 鈫� vn 鈭� 1,
%      presuming vn > 0.
%
%   4. Vehicle movement. The vehicles are deterministically moved by their
%      velocities, xn 鈫� xn + vn .
%
% USAGE: flux = ns(rho, p, L, tmax, isdraw)
%        rho       = density of the traffic
%        p         = probability of random braking
%        L         = length of the load
%        tmax      = number of the iterations
%        animation = if show the animation of the traffic
%        spacetime = if plot the space-time after the simuation ended.
%        flux      = flux of the traffic
%
% zhou lvwen: zhou.lv.wen@gmail.com
%if nargin == 0; rho = 0.25; p = 0.25; L = 100; tmax = 100; animation = true; spacetime = true;
endvmax = 5;                  % maximun speed
% place a distribution with density
ncar = round(L*rho);
rho = ncar/L;x = sort(randsample(1:L, ncar));
v = vmax * ones(1,ncar);   % start everyone initially at vmaxif animation; h = plotcirc(L,x,0.1); endflux = 0;                  % number of cars that pass through the end
vmean = 0;
road = zeros(tmax, L);for t = 1:tmax% accelerationv = min(v+1, vmax);%collision preventiongaps = gaplength(x,L); % determine the space vehicles have to movev = min(v, gaps-1);% random speed dropsvdrops = ( rand(1,ncar)<p );v = max(v-vdrops,0);% update the positionx = x + v;passed = x>L;          % cars passed at time rx(passed) = x(passed) - L;% periodic boundary conditionsif t>tmax/2flux = flux + sum(v/L); %flux = flux + sum(passed); vmean = vmean + mean(v);endroad(t,x) = 1;if animation; h = plotcirc(L,x,0.1,h); end
end
flux = flux/(tmax/2);
vmean = vmean/(tmax/2);if spacetime; figure;imagesc(road);colormap([1,1,1;0,0,0]);axis image; end% -------------------------------------------------------------------------function gaps = gaplength(x,L)
%
% GAPLENGTH: determine the gaps between vehicles
%
ncar = length(x);
gaps=zeros(1, ncar);
if ncar>0gaps = x([2:end 1]) -x;gaps(gaps<=0) = gaps(gaps<=0)+L;
end% -------------------------------------------------------------------------function h = plotcirc(L,x,dt,h)
W = 0.05;  R = 1;
ncar = length(x);theta = [0 : 2*pi/L : 2*pi];
xc = cos(theta);     yc = sin(theta);
xinner = (R-W/2)*xc; yinner = (R-W/2)*yc;
xouter = (R+W/2)*xc; youter = (R+W/2)*yc;xi = [xinner(x);  xinner(x+1); xouter(x+1); xouter(x)];
yi = [yinner(x);  yinner(x+1); youter(x+1); youter(x)];
if nargin == 3color = randperm(ncar);h = fill(xi,yi, color); hold onplot(xinner,yinner, 'k', xouter,youter, 'k','linewidth',1.5)plot([xinner; xouter], [yinner; youter],'k','linewidth',1.5)axis image; axis((R+2*W)*[-1 1 -1 1]); axis off
elsefor i=1:ncar;  set(h(i),'xdata',xi(:,i),'ydata',yi(:,i));  end
end
pause(dt)

横坐标是归一化的密度，纵坐标是车流量。由此可见车流量先随密度单调递增，到达某一值后随密度单调递减，直至为0。