元胞自动机在数模中的应用十分广泛。
元胞自动机原理
一维元胞自动机。给出任意一个状态，都能知道下一时刻的状态。规则已给出，总共有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;
end
vmax = 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 vmax
if animation; h = plotcirc(L,x,0.1); end
flux = 0; % number of cars that pass through the end
vmean = 0;
road = zeros(tmax, L);
for t = 1:tmax
% acceleration
v = min(v+1, vmax);
%collision prevention
gaps = gaplength(x,L); % determine the space vehicles have to move
v = min(v, gaps-1);
% random speed drops
vdrops = ( rand(1,ncar)<p );
v = max(v-vdrops,0);
% update the position
x = x + v;
passed = x>L; % cars passed at time r
x(passed) = x(passed) - L;% periodic boundary conditions
if t>tmax/2
flux = flux + sum(v/L); %flux = flux + sum(passed);
vmean = vmean + mean(v);
end
road(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>0
gaps = 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 == 3
color = randperm(ncar);
h = fill(xi,yi, color); hold on
plot(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
else
for i=1:ncar; set(h(i),'xdata',xi(:,i),'ydata',yi(:,i)); end
end
pause(dt)
横坐标是归一化的密度，纵坐标是车流量。由此可见车流量先随密度单调递增，到达某一值后随密度单调递减，直至为0。
6348
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
