2013高教社杯全国大学生数学建模竞赛A题第四问求解思路——元胞自动机NaSch交通流仿真模型

学校数模练习，本人编程，根据建模手的模型编程，感觉写的还行，思路清晰，遂发文记录

模型求解：
 

根据我们在模型建立中构建的交通流元胞自动机模型，需要得到车辆排队长度达到上游路口所需要的时间，因此我们给出了以下的算法求解步骤。

算法步骤：

1. 初始化基本参数，包括车辆长度、车道排队长度，设置元胞数组等。
2. 以分钟为外循环，秒钟为内循环进行交通流的仿真，在外循环中初始化车辆流入时间，在内循环中进行车辆规则的更新以及进行程序结束的判断，如果车辆排队长度达到140m则输出仿真系统运行的总时间。
3. 车辆规则的更新需要调用车辆更新规则函数，每次调用函数，首先根据每一车道的车辆占用情况更新车道速度，车道速度分为低、中、高三档速度。
4. 进行车辆行驶行为的更新，先判断车辆是否满足换道条件，如果满足则进行换道，换道过后，这一条道路后续车辆会受到前车换道的影响，前进速度会减缓
5. 假如车辆不满足换道条件，就前进，如果是距离出口较近（小于50m）的车辆，需要减速慢行，将其速度调整为低档速度，前进前先获取前方车辆占用情况，再分为没有车、全是车、有车有空位三种情况处理，如果有车就跟车，没有车就根据速度前进。在前进中还需要判断驶出的特殊情况，由于出口只有1道并且旁边发生车祸，车辆从1道驶出需要3s。
6. 原有车辆更新完之后，进行新车辆流入的更新，根据车辆流入时间进行车辆流入的更新。
7. 有任意道车辆排队长度达到140m则结束仿真并输出仿真系统运行的总时间。

结果分析：

通过MATLAB运行交通流元胞自动机仿真模型得到仿真系统运行总时间12个结果如下：

（可以给个表，也可以给图，需要的话喊我画）

去掉其中的最大最小值6m39s和5m58s并对剩下的10个结果取平均值得到平均仿真系统运行的总时间6m23s，这说明从事故发生开始，经过大约6m23s，车辆排队长度将达到上游路口。

运行代码：

% problem 4
clear
close all;
clc;

%初始化基本参数
global M S interval;
interval = 0; % 计时器变量
car_length = 5;% 全部为小车（包括安全车间距1m，车身4m）
total_unit_length = floor(140/car_length);% 车道单位长度=28 每格5m
flag = zeros(3,total_unit_length);% 定义元胞道路数组，0表示没车，1表示有车
queue_length = zeros(3,1);% 每条道路排队长度，用来限制车辆速度，起始为0
newcar = [7,16,12];% 各个车道每1min更新的车辆数量

%仿真循环
for M = 1:15 % 先以15min为结束条件，实际结束条件为任意车道排队长度达到140/5 = 28
    % 初始化车辆流入时间，一车道为1-60s内中5s有车进入，二车道为1-25s内中11s有车进入，三车道为1-25s内中9s有车进入
    newcar_road1 = randperm(60,newcar(1));
    newcar_road2 = randperm(25,newcar(2));
    newcar_road3 = randperm(25,newcar(3));

    for S = 1:60 % 60s为一个周期
        % 更新规则函数
        [queue_length,flag] = function1(queue_length,flag,total_unit_length,newcar_road1,newcar_road2,newcar_road3);
        % 程序结束判断：当任一车道排队长度达到28，输出仿真运行时间，程序结束
        if any(queue_length == 28)
            sprintf("仿真运行的时间为%d分%d秒",M,S)
            return
        end
    end
end

更新函数：
 

function [queue_length,flag] = function1(queue_length,flag,total_unit_length,newcar_road1,newcar_road2,newcar_road3)
global S interval;
     % 车道速度更新
        advance = [3 3 3]; % 排队长度小于10为3单位长度/s
        for q = 1:3
            if(queue_length(q)>10 && queue_length(q)<=20) % 排队长度大于10小于20为2单位长度/s
                advance(q) = 2;
            elseif(queue_length(q) > 20) % 排队长度大于20为1单位长度/s
                advance(q) = 1;
            end
        end

        % 原有车辆更新
        oldcar_index = find(flag == 1); % 获得原有车辆下标
        % 如果路上没车就不需要更新车辆了，有车才更新
        if(~isempty(oldcar_index))
            l = length(oldcar_index); % 获得总下标数
            for i=1:l
                % 获取车辆下标对应的车道位置
                if(mod(oldcar_index(i),3) == 0)
                    road_index = 3;
                else
                    road_index = mod(oldcar_index(i),3);
                end

                % 更新判断：
                % 换道判断，如果不在1道 并且 右边一格或者右边一格的再前面一格没车 并且 右边车道排队长度比左边小 或者 原车坐标小于等于30触发换道判断，反之只能前进
                % 换道：原有道路排队长度减1，更新道路排队长度加1
                if(road_index~=1)&&(~flag(oldcar_index(i)-1))&&((queue_length(road_index)>=queue_length(road_index-1))||(oldcar_index(i)<=30))&&(oldcar_index(i)<=84)
                    % 如果右边一格没车换道
                    flag(oldcar_index(i)) = 0;
                    flag(oldcar_index(i)-1) = 1;
                    queue_length(road_index) = queue_length(road_index) - 1;
                    queue_length(road_index - 1) = queue_length(road_index - 1) + 1;
                    % 换道导致后面的车缓行：前进速度减1单位长度/s
                    if(advance(road_index-1)>1)
                        advance(road_index-1) = advance(road_index - 1) - 1;
                    end
                else
                    % 3.前进：如果前进方向有车就跟车或是保持不动
                    % 当原有车下标小于30，即离出口50m内，减速慢行,更新速度为1单位长度/s
                    if(oldcar_index(i)<=30)
                        % 由于前进速度为1单位长度/s只用考虑前方一个单位长度的占用情况
                        % 驶出判断：车只能由下标1处驶出
                        % 车辆驶出后，排队长度减1
                        if(oldcar_index(i)==1)
                            interval = interval + 1;
                            if(interval == 3)
                            %驶出
                            flag(1) = 0;
                            queue_length(1) = queue_length(1) - 1;
                            interval = 0;
                            end
                            % 下标为2、3却不能换道只能在原地
                        elseif(oldcar_index(i)==2||oldcar_index(i)==3)
                        else
                            % 未驶出
                            if(~flag(oldcar_index(i)-3))
                                flag(oldcar_index(i)) = 0;
                                flag(oldcar_index(i)-3) = 1;
                            end
                        end
                        % 当原有车下标大于12则服从车道规定速度行驶（不会触发驶出判断）
                    else
                        temp = [];
                        % 将原有车辆坐标最大前进单位长度里的车辆占用情况赋值给temp
                        for j=advance(road_index):-1:1
                            temp = [temp,flag(oldcar_index(i)-3*j)];
                        end

                        % 没车
                        if(all(temp(:) == 0))
                            flag(oldcar_index(i)) = 0;
                            flag(oldcar_index(i)-3*advance(road_index))= 1;
                            % 全是车：不动，不需要操作
                        elseif(all(temp(:) == 1))
                            % 有空出来的单位长度:找到下标最大的占用位置，进而算出新前进单位长度
                        else
                            idx = find(temp == 1);
                            new_advance = length(temp)-max(idx);
                            flag(oldcar_index(i)) = 0;
                            flag(oldcar_index(i)-3*new_advance)= 1;
                        end
                    end
                end
            end
        end

        % 新车辆流入
        % 驶入判断：如果为汽车驶入的时间并且每个车道最后一个位置没有被占用则成功驶入
        if (any(newcar_road1 == S)&&flag(1,total_unit_length)==0)
            flag(1,total_unit_length) = 1;
            queue_length(1) = queue_length(1) + 1;
        end
        if(S<=25)
            if (any(newcar_road2 == S)&&flag(2,total_unit_length)==0)
                flag(2,total_unit_length) = 1;
                queue_length(2) = queue_length(2) + 1;
            end
            if (any(newcar_road3 == S)&&flag(3,total_unit_length)==0)
                flag(3,total_unit_length) = 1;
                queue_length(3) = queue_length(3) + 1;
            end
        end
end