Nagatani modelの交差点渋滞問題への適用

モデルの定義

# 流入口の流入量 (モデル設定として与える)
import math
def injection(t, tm, J0=0.2, Jm=.5, sigma2 = 8.0): # 流入量を返す
    # tm ピークの時刻、 J0 平常時、Jm ピーク、sigma2 混雑時間帯の幅の2乗*2
    if math.fabs(t-tm) > math.sqrt(sigma2)*3:
        return J0
    return J0 + (Jm - J0)/math.exp((t - tm)*(t-tm)/sigma2)

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# 速度の計算 (引数：密度、最高速度)
def velocity(rho, v_c):
    rho_c = 1.0/(v_c + 1.0)
    if rho  < rho_c:
        return v_c
    else:
        return (1.0 - rho)/rho
    
# ノード間の流れ （引数 密度、v_max, v_b, 速度計算関数）
#  返り値： ノードペア間の流れ：密度（ノード）の次元 x 密度（ノード）の次元
def current(rho, v_std, v_fast, f, velocity):
    '''
     v_fast, v_std 道路の最高速度
     v_out1, v_out2 出口側の速度（一定値として与える）
     f 交差点での分岐車数割合（飯田通りとアルプス通りへの流入割合）
    '''
    J = np.zeros((len(rho), len(rho)))
    
    f = 0.5      # 1to4, 2to5, 3to6への割合    
    J[0,3] = rho[0] * velocity(rho[3], v_fast)
    J[1,4] = rho[1] * velocity(rho[4], v_std)
    J[2,5] = rho[2] * velocity(rho[5], v_std)
    J[3,6] = f * rho[3] * velocity(rho[6], v_std)
    J[3,7] = (1. - f) * rho[3] * velocity(rho[7], v_fast)
    J[4,6] = f * rho[4] * velocity(rho[6], v_std)
    J[4,7] = (1. - f) * rho[4] * velocity(rho[7], v_fast)
    J[5,6] = f * rho[5] * velocity(rho[6], v_std)
    J[5,7] = (1. - f) * rho[5] * velocity(rho[7], v_fast)
    
    J[6,8] = rho[6] * velocity(rho[8], v_std)
    J[8,10] = rho[8] * velocity(rho[10], v_std)
    J[10,12] = rho[10] * velocity(rho[12], v_std)
    
    J[7,9] = rho[7] * velocity(rho[9], v_fast)
    J[9,11] = rho[9] * velocity(rho[11], v_fast)
    J[11,13] = rho[11] * velocity(rho[13], v_fast)
    
    return J

# 微分方程式の右辺 x: 密度
def current_eq(rho_in, t, v_std, v_fast, v_out1, v_out2, f,
               peak_time0, peak_time1, peak_time2,
               velocity):
    # 微分方程式の定義
    J = current(rho_in, v_std, v_fast, f, velocity) # velocityは速度-密度関数
    rho = np.zeros(len(rho_in))
    # 流入口
    rho[0] = injection(t, tm=peak_time0) - J[0,3]# 第２引数はピークの時刻
    rho[1] = injection(t, tm=peak_time1) - J[1,4]
    rho[2] = injection(t, tm=peak_time2) - J[2,5]
    # print([t,rho[0],rho[1],rho[2]])
    # 流入部
    rho[3] = J[0,3] - J[3,6] - J[3,7]
    rho[4] = J[1,4] - J[4,6]- J[4,7]
    rho[5] = J[2,5] - J[5,6]- J[5,7]
    # 交差点部
    rho[6] = J[3,6] + J[4,6] + J[5,6] - J[6,8]
    rho[7] = J[3,7] + J[4,7] + J[5,7] - J[7,9]
    # 飯田通り、アルプス通り
    rho[8] = J[6,8] - J[8,10]
    rho[9] = J[7,9] - J[9,11]
    rho[10] = J[8,10] - J[10,12]
    rho[11] = J[9,11] - J[11,13]
    # 出口
    rho[12] = J[10,12] - v_out1 # 定常値
    rho[13] = J[11,13] - v_out2
    
    return rho

具体的な計算及び描画関数と実行

# 時間変化の数値計算（引数：密度, v_max, v_b）と描画
def nagatani_crossSection(init_rho, v_std=1.0, v_fast=1.5, f=0.2, v_out1=0.7, v_out2=0.7,
                         peak_time0=10.0, peak_time1=12.0, peak_time2=14.0):
    # 道路パラメータ
    #f = 0.5 # 飯田通りへの流入割合 (1-fはアルプス)
    #v_std  その他の通り
    #v_fast アルプス通り
    
    dt= 0.002
    n_steps =15000
    t = np.arange(0, n_steps*dt, dt) # 時刻のリスト
    # 密度の初期値
    x = np.full(14, init_rho)
    # 微分方程式の数値計算
    y = odeint(current_eq,x,t, args=(v_std, v_fast, v_out1, v_out2, f,
                                     peak_time0, peak_time1, peak_time2,
                                     velocity)) 
    
    # 描画
    plt.figure(figsize=(15, 12))
    plt.subplots_adjust(wspace=0.2, hspace=0.3) # h = height, w = width
    
    # (0) 流入量（初期パラメータ）のプロット
    ax = plt.subplot(3,2,1)
    in_flow = np.zeros((len(t),3))
    in_flow[:,0] = list(map(lambda x: injection(x, peak_time0), t[:]))
    in_flow[:,1] = list(map(lambda x: injection(x, peak_time1), t[:]))
    in_flow[:,2] = list(map(lambda x: injection(x, peak_time2), t[:]))
    ax.plot(t, in_flow, label=["to 0", "to 1", "to 2"])
    ax.set_xlabel("time")
    ax.set_ylabel("injection") 
    
    # (1) 密度 (density)
    # 上流
    ax = plt.subplot(3,2,3)
    ax.plot(t, y[:, np.array([3,4,5])])
    ax.legend(["アルプス北",
                "富士見",
                "飯田西"])
    ax.set_xlabel("time")
    ax.set_ylabel("Density (initial mean density = $\\rho$")
    ax.set_title("Density")
    
    # 下流
    ax = plt.subplot(3,2,5)
    ax.plot(t, y[:, np.array([6,7])])
    ax.legend([ "飯田東",
                "アルプス南"])
    ax.set_xlabel("time")
    ax.set_ylabel("Density (initial mean density = $\\rho$")
    ax.set_title("Density")
    
    # (2) 流量(current) = (密度) * (速度)
    # 飯田方面
    ax2 = plt.subplot(3,2,4)
    J = np.array(list(map(lambda x: current(x, v_std, v_fast, f, velocity ), y[:])))
    # 個別に書く
    ax2.plot(t, J[:, 3,6], label="アルプス-飯田")
    ax2.plot(t, J[:, 4,6], label="富士見-飯田")
    ax2.plot(t, J[:, 5,6], label="飯田-飯田")
    ax2.legend()
    ax2.set_xlabel("time")
    ax2.set_ylabel("J")
    ax2.set_title("Current（飯田通り向き）")
    # 下流
    ax2 = plt.subplot(3,2,6)
    J = np.array(list(map(lambda x: current(x, v_std, v_fast, f, velocity ), y[:])))
    # 個別に書く
    ax2.plot(t, J[:, 3,7], label="アルプス-アルプス")
    ax2.plot(t, J[:, 4,7], label="富士見-アルプス")
    ax2.plot(t, J[:, 5,7], label="飯田-アルプス")
    ax2.legend()
    ax2.set_xlabel("time")
    ax2.set_ylabel("J")
    ax2.set_title("Current (アルプス通り向き)")
    

nagatani_crossSection(0.2)