异构交通流matlab源码

% 源程序
clear all


tic %开始计时，一般与toc连用，表示结束
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 调整参数，影响稳定性
a=1.3;%最大加速度%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%,调整对稳定性有重要影响
b=3;%期望减速度
T=1.6;%安全车头时距%
v00=33;%期望速度
s0=4;%拥堵车距
alpha=0.3;%三个0到1范围内的系数
beta=alpha;
gamma=alpha;
tau1=0.2;  tau2=0;%反应时间
s=5; % 车长
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  调整密度，影响稳定性
road=2000;
N=50;
h=road/N;%平均初始车头间距，对应的密度也不相同
tend=40000;
dt=0.1;

allv=cell(5,1);
alldx=cell(5,1);
allx=cell(5,1);
iN=4;
Be=1;

allv=cell(5,1);
alldx=cell(5,1);
allx=cell(5,1);
iN=4;
%% 随机参数 ，确定网联车的分布
% for
ii=2;
dx=[];v=[];dv=[];x=[];
df=0.2;
Jy=[1 0];%判断是否网联车，1是，0不是
p1=0.3;p2=1-p1; % 网联车的比p1；  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
prob1=[p1 p2];
L=length(Jy);%Jy的长度 1行2列所以为2
cc=zeros(L,1);%2行1列0阵
i=1;
FVm=[]; % 向量，确定每辆车是否是网联车
fprob=prob1;
prob=prob1;
fVm=Jy;   
fnum1=[];
while i<=N
    vm=randsrc(1,1,[fVm; fprob]);   %按照fprob的概率随机生成fVm中指定的数字 也就是0和1  1代表网联车概率为p1
    FVm(i)=vm;
    nn=find(abs(Jy-vm)<10^(-10));%满足条件的所有vm的所在位置
    fnum1(:,i)=nn(1)';%nn的第一个元素赋值给fnum1的第一列
    cc(nn(1))=cc(nn(1))+1;
    c1=find(abs(cc-N*prob')>1/2);
    fVm=Jy(c1);
    fprob=prob(c1)/sum(prob(c1));
    i=i+1;
end
AJc=zeros(1,N);   % 生成1xN的全0阵  向量，确定当前车和前车是否为网联车，以确定网联是否起效。注意：网联起效比例跟理论有一点差距的，理论值属于期望平均值。
%
for i=1:N
    if i==N
        if FVm(i)==1 && FVm(1)==1  %环形链接
            AJc(i)=1;
        else
            AJc(i)=0;
        end
    else
        if FVm(i)==1 && FVm(i+1)==1  %前后链接
            AJc(i)=1;
        else
            AJc(i)=0;
        end
    end
end

AJc=FVm;
%% 稳态车头间距
dx1=h;
aa=0; bb=100;
ve=(aa+bb)/2;
mm1=T*ve;
rf=a* (1-(ve/v00).^4-((s0+mm1)./(dx1-s)).^2);
while abs(rf)>10^(-20)
    if rf>0
        aa=ve;
    else
        bb=ve;
    end
    ve=(aa+bb)/2;
    %     v1=[1:30];
    mm1=T*ve;
    rf=a* (1-(ve/v00).^4-((s0+mm1)./(dx1-s)).^2);
end

%% 初始条件
v1=ones(1,N)*ve; dv1=zeros(1,N)*ve;   %稳态速度代入作为初始速度
dx1=ones(1,N)*h;  
x1=[];
for i=1:N
    x1(i)=sum(dx1(1:i-1));%每一列求和
end
Bi1=ceil(tau1/dt); Bi2=ceil(tau2/dt); %%%%%%向前取整
Be=Bi1+1;
x=[];v=[];dx=[];dv=[];
for i=1:Be
    x(i,:)=x1;
    v(i,1:N)=v1*0;
    dv(i,:)=dv1;
    dx(i,:)=dx1;
end
v(Be,1:2)=v(Be,1:2)+1;

完整代码添加QQ1575304183