NS元胞自动机模型--python实现

实现：

]# -*- coding: utf-8 -*-
''' NS模型
    场景：
        周期型边界
        道路长度:cell=1000个元胞
        车辆初始分布为均匀分布
        初始速度:v0=vmax=5
        随机慢化概率:p=0.1
        仿真时步为2000时步，从500时步开始采样
        -1表示元胞，其他值表示车辆
    要求：绘制车辆加速度的分布图(密度：0.05, 0.2, 0.4, 0.6)
'''

#import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import os

#设置车辆初始均匀分布,车辆密度ρ=0.05, 0.2, 0.4, 0.6
def InitDist(v0=5, ρ=0.05, cell=1000):
    arra = np.zeros(1000, dtype=int)#创建1000个元胞
    arra = np.array([-1] * cell)
    v_count = int(ρ*cell)           #车辆数
    temp = np.linspace(0, 999, num=v_count).astype(int)#生成均匀分布数
    for i in temp:
        arra[i] = v0
    return v_count,arra #返回车辆初始分布与车辆数


#NS建模
def NS(cell=1000, n=InitDist()[0],v0=5, vmax=5, p=0.1,times=2000, start=500):
    cell_arr = np.zeros((times,cell), dtype=int)#NS模型演化过程表
    df = pd.DataFrame(cell_arr).apply(lambda x: x.iloc[:] - 1)
    df.loc[0,:] = InitDist()[1]                 #设置车辆初始状态

    #根据车辆数n生成相应的任意随机慢化概率,支持文件导入
    p_rand = pd.DataFrame()
    if(os.path.exists('p_array.csv')==False):
        p_rand = pd.DataFrame(np.random.randint(0, 10, size=(times, n)))
        p_rand = p_rand.apply(lambda x: x.iloc[:] / 10)
        p_rand.to_csv("p_array.csv",index=True,header=True)
    else:
        p_rand = pd.DataFrame(pd.read_csv('p_array.csv',index_col=0))

    #NS演化规则
    i = 0
    for i in range(times - 1):
        index = []                      #保存车辆所在元胞的下标
        d_temp = list(df.iloc[i,:])     #每次取出当前路段元胞的车辆分布情况
        p_temp = list(p_rand.iloc[i,:]) #每次取出当前路段车辆的换道概率
        for j in range(cell):
            if d_temp[j] >= 0:
                index.append(j)         #将车辆的当前位置保存

        #对每一辆车进行处理
        for k in range(n):
            v = d_temp[index[k]]        #获取当前速度

            #加速
            if(v < vmax):
                v = min(v+1, vmax)
            #减速
            if index[(k+1)%n] > index[k]:
                gap = index[(k + 1) % n] - index[k] - 1   #两车相隔gap个元胞
            else:
                gap = cell - index[k] + index[(k+1)%n] - 1#周期型边界，需处理边界问题
            if (gap < vmax):
                v = min(v,gap)
            #随机慢化
            if(p > p_temp[k]):
                v = max(v-1,0)
            #print(v)
            #运动
            s = index[k] + v
            s = s % cell         #边界处理
            df.loc[i+1,s] = v    #更新下一轮车辆状态变化值

    df.to_csv("result.csv")
    return df

if __name__ == "__main__":
    df = NS()
    print(df)