模型的建立基本来自于:​​http://www.doc88.com/p-2078634086043.html​

花了一天半的时间用新学会的matlab实现了一下。

 

───────────────────────────────────────────────────────────────────────────────────────

2018-1-31更新:

居然真的有人会看我的博文诶,那我不能就这么不负责任的直接甩代码跑路了,稍微评价一下我自己的代码吧……

 

1、首先是优点,就是能用(囧),调调参数,几个车道,多少元胞,多少辆车,自动驾驶车手动驾驶车比例,基本都能改。

2、缺点一:各种规则十分不清晰明确,正确性有待商榷……

  特别是变道规则,又臭又长,我自己写完之后都不想看第二遍,而且我怀疑跟论文里写的规则有出入(虽然我说实话,论文里的变道规则也说的不是很清晰)

3、缺点二:速度贼慢,慢的跟爬一样。因为这个代码是为了数学建模敲的,所以也多考虑优化速度什么的,直接就简单粗暴地写了。

  这导致了,我在出结果的时候,车道数不敢设大于3,元胞数不敢设大于200,循环时常不敢超过300,;

  即使这样,初始放入车辆的时候,元胞占用率大于0.3的时候,程序跑的就肉眼可见的慢,如果想跑个几十上百组的数据求平均的话,估计可以打局王者荣耀再来看看结果……

4、稳定性贼差,我算结果的时候,几百组几百组求平均(元胞占用率0.0几的情况下,跑的还是挺快的),结果都稳定不下来……导致我们敏感性分析基本等于瞎扯。

 

感觉自己是真鸡儿的菜,写了这么久的代码,敲个稍微大点的程序,还是敲得跟坨shit一样基于元胞自动机NaSch模型的多车道手动-自动混合驾驶仿真模型的Matlab实现_javascript

───────────────────────────────────────────────────────────────────────────────────────

2018-5-7更新:

emmm发现居然少贴了一个函数的代码,补上了calc_Vave()函数;

然后关于这个模型怎么用……

最简单的用法是:

更改global变量的值(车辆数、车道数、自动驾驶车辆占比等等等等),更改模拟时长T,然后运行NaSch()这个函数,返回的ans答案是整个时长T内所有车辆的平均速度

───────────────────────────────────────────────────────────────────────────────────────

 

NaSch模型主函数:



function [ret_Vave] = NaSch()
rng('shuffle');

global CellsNum; CellsNum = 200; %The number of cells in each lane
global LanesNum; LanesNum = 3; %The number of Lanes
global D; D = 0.02; %The cell occupancy rate
global L; L = 5; %The length of one cell
global CarsNum; CarsNum = round( LanesNum .* CellsNum .* D ); %The number of cars
global SelfRatio; SelfRatio = 0.5; %Self-driving vehicle rate
global Vmin; Vmin = 1; %Minimum speed limit
global Vmax; Vmax = 6; %Maximum speed limit

%Meet the lane changing conditions...
global ChangeRightP; ChangeRightP = 0.5; %The probability of changing into the right lane
global ChangeLeftP; ChangeLeftP = 0.8; %The probability of changing into the left lane

[z,cars] = RoadInit();
cells = z;
cars = CarInit(cars);

T = 1000; %The duration of simulation
Vave = zeros(1,T); %Average velocity per unit time

global loopCars; loopCars = [];
%Start the simulation
for t = 1:1:T

[cells,cars] = EntryControl(cells,cars);

%Calculate each distance
[d,d_safe,d_l_for,d_l_bac,d_r_for,d_r_bac] = calcDist(cells,cars);

%Get the next state
for i = 1:1:CarsNum
if(cars(i).M==1) %Leftmost lane
[cells,cars] = L_NextState(i,cells,cars,d,d_safe,d_l_for,d_l_bac,d_r_for,d_r_bac);
elseif(cars(i).M==LanesNum) %Rightmost lane
[cells,cars] = R_NextState(i,cells,cars,d,d_safe,d_l_for,d_l_bac,d_r_for,d_r_bac);
else %Middle lane
[cells,cars] = M_NextState(i,cells,cars,d,d_safe,d_l_for,d_l_bac,d_r_for,d_r_bac);
end
end

cars = RandSlow(cars);
[cells,cars] = Update(cells,cars);
cells = cells(1:LanesNum,1:CellsNum); %Control the size of Cells Matrix
Vave(t) = calc_Vave(cars); %Calculate current average velocity
end

ret_Vave = mean(Vave); %Calculate average velocity
end


 

初始化函数:



function [cell_mat,car_mat] = RoadInit()
global LanesNum;
global CellsNum;
global CarsNum;

mat = zeros(LanesNum,CellsNum);

Car(1,CarsNum) = struct('ID',[],'V',[],'M',[],'N',[],'Self',[]);
%Vehicle Information: ID, Speed, Position, Whether self driving or not

for i = 1:1:CarsNum
m=round( 1 + ( LanesNum - 1 ) * rand(1) );
n=round( 1 + ( CellsNum - 1 ) * rand(1) );
while( mat(m,n) ~= 0 )
m = round( 1 + ( LanesNum - 1 ) * rand(1) );
n = round( 1 + ( CellsNum - 1 ) * rand(1) );
end
mat(m,n) = i;

Car(i).ID = i;
Car(i).M = m;
Car(i).N = n;
end

cell_mat = mat;
car_mat = Car;
end



function car_mat = CarInit(Car)
global CarsNum;
global Vmin;
global Vmax;
global SelfRatio;

for i = 1:1:CarsNum
Car(i).V = round( Vmin + ( Vmax - Vmin ) * rand(1) );
if(rand(1)<=SelfRatio) Car(i).Self = 1;
else Car(i).Self = 0;
end
end

car_mat = Car;
end


 

入口控制:



function [new_cells,new_cars] = EntryControl(cells,cars)
global loopCars;
global LanesNum;
global Vmin;

while(StackSize(loopCars) > 0)

if( all( cells(:,1) ~= 0 ) == 1 )
break
end

for i = 1:1:LanesNum
if(cells(i,1) == 0)
id = topStack(loopCars); popStack(loopCars);
cells(i,1) = id;
cars(id).V = Vmin;
cars(id).M = i;
cars(id).N = 1;
break
end
end

end

new_cells = cells;
new_cars = cars;
end


 

计算各种车距:



function [Ret_d,Ret_d_safe,Ret_d_l_for,Ret_d_l_bac,Ret_d_r_for,Ret_d_r_bac] = calcDist(cells,cars)
global LanesNum;
global CellsNum;
global CarsNum;

d = linspace(-1,-1,CarsNum);
d_safe = linspace(-1,-1,CarsNum);
d_l_for = linspace(-1,-1,CarsNum);
d_l_bac = linspace(-1,-1,CarsNum);
d_r_for = linspace(-1,-1,CarsNum);
d_r_bac = linspace(-1,-1,CarsNum);

for i = 1:1:CarsNum

m = cars(i).M; %The m_th lane
n = cars(i).N; %The n_th cell

%Calculate the distance to front car
d(i) = inf;
for f = n+1:1:CellsNum
if(cells(m,f) ~= 0)
d(i) = f - n;
break
end
end

%If the left lane exist
if( m-1 >= 1 )
%Forward
d_l_for(i) = inf;
for p = n:1:CellsNum
if( cells( m-1 , p ) ~= 0 )
d_l_for(i) = p - n;
break
end
end
%Backward
for p = n:-1:1
d_l_bac(i) = inf;
if( cells( m-1 , p ) ~= 0 )
d_l_bac(i) = n - p;
break
end
end
end

%If the right lane exist
if( m+1 <= LanesNum )
%Forward
d_r_for(i) = inf;
for p = n:1:CellsNum
if( cells( m+1 , p ) ~= 0 )
d_r_for(i) = p - n;
break
end
end
%Backward
d_r_bac(i) = inf;
for p = n:-1:1
if( cells( m+1 , p ) ~= 0 )
d_r_bac(i) = n - p;
break
end
end
end

%The d_safe of two kinds of cars is different
if(cars(i).Self == 1)
d_safe(i) = cars(i).V;
else
d_safe(i) = 2 .* cars(i).V;
end

end

Ret_d = d;
Ret_d_safe = d_safe;
Ret_d_l_for = d_l_for;
Ret_d_l_bac = d_l_bac;
Ret_d_r_for = d_r_for;
Ret_d_r_bac = d_r_bac;
end


 

加减速规则实现:



function new_cars = Accelerate(id,cars,d,d_Safe)
global Vmax;
global Pa;
Pa = 0.8; %Acceleration probability

if( d(id) > d_Safe(id) && cars(id).V < Vmax )
if(rand(1) <= Pa)
cars(id).V = cars(id).V + 1;
end
end

new_cars = cars;
end



function new_cars = Decelerate(id,cars,d,d_safe)
global Vmax;
global Vmin;

if( d(id) < d_safe(id) && d(id) >= Vmax - Vmin )
cars(id).V = max( cars(id).V - 1 , Vmin );
elseif( d(id) < Vmax - Vmin )
cars(id).V = Vmin;
end

new_cars = cars;
end


 

关于变道的状态更新:



function [new_cells,new_cars] = L_NextState(i,cells,cars,d,d_safe,d_l_for,d_l_bac,d_r_for,d_r_bac)
global Vmax;
global LanesNum;
global loopCars;

if( d_r_for(i) > d_safe(i) ) %与右前方车辆距离大于安全距离
ok1 = 1;
else
ok1 = 0;
end

if( d_r_bac(i) > 1 + Vmax - min(cars(i).V+1,Vmax) ) %与右后方车辆安全
ok2 = 1;
else
ok2 = 0;
end

if( cars(i).M+2 <= LanesNum && cells(cars(i).M+2,cars(i).N) ~= 0 ) %右边相隔一个车道平行位置有车
id = cells(cars(i).M+2,cars(i).N); %获取该车编号
if( min( d_l_for(id) , d_l_bac(id) ) < d_safe(id) ) %判断该车不能向左变道
ok3 = 1;
else
ok3 = 0;
end
else
ok3 = 1;
end

ChangeLaneOK = ok1 && ok2 && ok3; %得到是否有变道条件

if( d(i) < d_safe(i) ) %与前车距离小于安全距离,能变道则必须变道
if(ChangeLaneOK)

cells(cars(i).M,cars(i).N) = 0;
cars(i).M = cars(i).M + 1;
if( cells( cars(i).M , cars(i).N ) ~= 0 ) %撞车控制
loopCars = pushStack(loopCars,cars(i).ID);
else
cells(cars(i).M,cars(i).N) = cars(i).ID; %变道
end

cars = Accelerate(i,cars,d,d_safe); %加速
else
cars = Decelerate(i,cars,d,d_safe); %减速
end
else %与前车距离大于等于安全距离
ChangeLaneDesire = getCLD(i,cars,'R'); %计算变道欲望
if(ChangeLaneOK && ChangeLaneDesire)

cells(cars(i).M,cars(i).N) = 0;
cars(i).M = cars(i).M + 1;
if( cells( cars(i).M , cars(i).N ) ~= 0 ) %撞车控制
loopCars = pushStack(loopCars,cars(i).ID);
else
cells(cars(i).M,cars(i).N) = cars(i).ID; %变道
end

cars = Accelerate(i,cars,d,d_safe); %加速
else
cars = Accelerate(i,cars,d,d_safe); %加速
cars = Decelerate(i,cars,d,d_safe); %减速
end
end

new_cells = cells;
new_cars = cars;
end



function [new_cells,new_cars] = M_NextState(i,cells,cars,d,d_safe,d_l_for,d_l_bac,d_r_for,d_r_bac)
global Vmax;
global LanesNum;
global loopCars;
%计算中间车道上车辆的下一时刻状态


if( d_l_for(i) > d_safe(i) ) %与左前方车辆距离大于安全距离
ok1 = 1;
else
ok1 = 0;
end

if( d_l_bac(i) > 1 + Vmax - min(cars(i).V+1,Vmax) ) %与左后方车辆安全
ok2 = 1;
else
ok2 = 0;
end

if( cars(i).M-2 >= 1 && cells(cars(i).M-2,cars(i).N) ~= 0 ) %左边相隔一个车道平行位置有车
id = cells(cars(i).M-2,cars(i).N); %获取该车编号
if( min( d_r_for(id) , d_r_bac(id) ) < d_safe(id) ) %判断该车不能向右变道
ok3 = 1;
else
ok3 = 0;
end
else
ok3 = 1;
end

ChangeLeftLaneOK = ok1 && ok2 && ok3; %得到是否有向左变道条件


if( d_r_for(i) > d_safe(i) ) %与右前方车辆距离大于安全距离
ok1 = 1;
else
ok1 = 0;
end

if( d_r_bac(i) ~= inf && d_r_bac(i) ~= -1 && cells( cars(i).M + 1 , cars(i).N - d_r_bac(i) ) > 0 ) %右后方有车
id = cells( cars(i).M + 1 , cars(i).N - d_r_bac(i) ); %得到该车的编号
if( d_r_bac(i) > d_safe(id) ) %与右后方车辆安全
ok2 = 1;
else
ok2 = 0;
end
else
ok2 = 1;
end

if( cars(i).M+2 <= LanesNum && cells(cars(i).M+2,cars(i).N) ~= 0 ) %右边相隔一个车道平行位置有车
id = cells(cars(i).M+2,cars(i).N); %获取该车编号
if( min( d_l_for(id) , d_l_bac(id) ) < d_safe(id) ) %判断该车不能向左变道
ok3 = 1;
else
ok3 = 0;
end
else
ok3 = 1;
end

ChangeRightLaneOK = ok1 && ok2 && ok3; %得到是否有向右变道条件


if( d(i) < d_safe(i) ) %与前车距离小于安全距离

if(ChangeRightLaneOK) %能够向右变道则向右变道

cells(cars(i).M,cars(i).N) = 0;
cars(i).M = cars(i).M + 1;
if( cells( cars(i).M , cars(i).N ) ~= 0 ) %撞车控制
loopCars = pushStack(loopCars,cars(i).ID);
else
cells(cars(i).M,cars(i).N) = cars(i).ID; %变道
end

cars = Accelerate(i,cars,d,d_safe); %加速

else %不能向右变道则尝试向左变道

if(ChangeLeftLaneOK == 1) %能向左变道则向左变道

cells(cars(i).M,cars(i).N) = 0;
cars(i).M = cars(i).M - 1;
if( cells( cars(i).M , cars(i).N ) ~= 0 ) %撞车控制
loopCars = pushStack(loopCars,cars(i).ID);
else
cells(cars(i).M,cars(i).N) = cars(i).ID; %变道
end

else %两边都不允许变道
cars = Decelerate(i,cars,d,d_safe); %减速
end

end

else %与前车距离大于等于安全距离

ChangeLaneDesire = getCLD(i,cars,'R'); %计算变道欲望

if(ChangeRightLaneOK && ChangeLaneDesire) %能够向右变道,且有变道欲望

cells(cars(i).M,cars(i).N) = 0;
cars(i).M = cars(i).M + 1;
if( cells( cars(i).M , cars(i).N ) ~= 0 ) %撞车控制
loopCars = pushStack(loopCars,cars(i).ID);
else
cells(cars(i).M,cars(i).N) = cars(i).ID; %变道
end

cars = Accelerate(i,cars,d,d_safe); %加速
else
cars = Accelerate(i,cars,d,d_safe); %加速
cars = Decelerate(i,cars,d,d_safe); %减速
end

end

new_cells = cells;
new_cars = cars;
end



function [new_cells,new_cars] = R_NextState(i,cells,cars,d,d_safe,d_l_for,d_l_bac,d_r_for,d_r_bac)
global Vmax;
global loopCars;
%计算最右侧车道上车辆的下一时刻状态

if( d_l_for(i) > d_safe(i) ) %与左前方车辆距离大于安全距离
ok1 = 1;
else
ok1 = 0;
end

if( d_l_bac(i) > 1 + Vmax - min(cars(i).V+1,Vmax) ) %与左后方车辆安全
ok2 = 1;
else
ok2 = 0;
end

if( cars(i).M-2 >= 1 && cells(cars(i).M-2,cars(i).N) ~= 0 ) %左边相隔一个车道平行位置有车
id = cells(cars(i).M-2,cars(i).N); %获取该车编号
if( min( d_r_for(id) , d_r_bac(id) ) < d_safe(id) ) %判断该车不能向右变道
ok3 = 1;
else
ok3 = 0;
end
else
ok3 = 1;
end

ChangeLaneOK = ok1 && ok2 && ok3; %得到是否有变道条件

if( d(i) < d_safe(i) ) %与前车距离小于安全距离
ChangeLaneDesire = getCLD(i,cars,'L'); %计算变道欲望
if(ChangeLaneOK && ChangeLaneDesire)
cells(cars(i).M,cars(i).N) = 0;
cars(i).M = cars(i).M - 1;
if( cells( cars(i).M , cars(i).N ) ~= 0 ) %撞车控制
loopCars = pushStack(loopCars,cars(i).ID);
else
cells(cars(i).M,cars(i).N) = cars(i).ID; %变道
end
else
cars = Decelerate(i,cars,d,d_safe); %减速
end
else %与前车距离大于等于安全距离
cars = Accelerate(i,cars,d,d_safe); %加速
cars = Decelerate(i,cars,d,d_safe); %减速
end

new_cells = cells;
new_cars = cars;
end


 

计算变道欲望:



function ret = getCLD(id,cars,Type)
%Calculate the desire to change the lane
global ChangeRightP;
global ChangeLeftP;

if(Type == 'R') %To right
if(cars(id).Self == 1) %Self driving
ret = 1;
else %Manual driving
if(rand(1) <= ChangeRightP)
ret = 1;
else
ret = 0;
end
end;
else %To left
if(cars(id).Self == 1) %Self driving
ret = 1;
else %Manual driving
if(rand(1) <= ChangeLeftP)
ret = 1;
else
ret = 0;
end
end;
end

end


 

随机慢化:



function new_cars = RandSlow(cars)
global CarsNum;
global Vmin;

p=0.1; %手动车慢化概率
q=0.03; %自动驾驶车慢化概率

for i = 1:1:CarsNum
if(cars(i).Self==1)
if(rand(1) <= q)
cars(i).V = max( cars(i).V - 1 , Vmin );
end
else
if(rand(1) <= p)
cars(i).V = max( cars(i).V - 1 , Vmin );
end
end
end

new_cars = cars;
end


 

状态更新:



function [new_cells,new_cars] = Update(cells,cars)
global CarsNum;
global CellsNum;
global loopCars;

for i = 1:1:CarsNum
cells( cars(i).M , cars(i).N ) = 0;
if( cars(i).N + cars(i).V > CellsNum ) %Drive out of the road
loopCars = pushStack(loopCars,cars(i).ID);
else
cells( cars(i).M , cars(i).N ) = cars(i).ID;
end;
end

new_cells = cells;
new_cars = cars;
end


 

计算某一时刻的平均车速:



function Ret = calc_Vave(cars)
global CarsNum;

sum = 0;
for i = 1:1:CarsNum
sum = sum+cars(i).V;
end

Ret = sum./CarsNum;
end


 

栈模拟:



function newStack = pushStack(Stack,val)
newStack = [val,Stack];
end



function newStack = popStack(Stack)
newStack = Stack(2:end);
end



function ret = StackSize(Stack)
ret = size(Stack,2);
end



function Val = topStack(Stack)
Val = Stack(1);
end