廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

シンプルな LWR モデルの近似解

モデル

渋滞の数理モデルの一つに LWR モデルというのがあり、以下の偏微分方程式で記述される。

\displaystyle \frac{\partial}{\partial t}\rho(x,t) + \frac{\partial}{\partial x} q(x,t)=0

x: 位置
t: 時間
\rho(x,t): 自動車の密度
q(x,t): 自動車の通過量

このモデルでは一次元の道路を通行する自動車の密度を考えている。

f:id:abrahamcow:20170923055005p:plain

自動車の密度 \rho(x,t) は、単位長さあたりの自動車の台数(台数/道のり)で、自動車の通過量 q(x,t) は、単位時間あたりの自動車の台数(台数/時間)なので、以下の等式が成り立つ。

\displaystyle \frac{\partial}{\partial t}\rho(x,t) = - \frac{\partial}{\partial x} q(x,t)

自動車の通過量(出ていく量)が増えるほど、自動車の密度は薄くなるので、右辺の符号はマイナスにしている。
以上が LWR モデルの大雑把な導出である。

上記の式を以下のように書き換える。

\displaystyle \frac{\partial}{\partial t}\rho(x,t) = - \frac{\partial q(\rho)}{\partial \rho} \frac{\partial \rho(x,t)}{\partial x}

\partial q(\rho)/\partial \rho が一定(=q_0)とすると、

\displaystyle \frac{\partial}{\partial t}\rho(x,t) = - q_0 \frac{\partial \rho(x,t)}{\partial x}

となる。

これを一階差分で近似すると、

\displaystyle \rho^{m+1}_{n} - \rho^{m}_{n} = - q_0 \frac{\Delta t}{\Delta x} (\rho^{m}_{n} - \rho^{m}_{n-1})

となる。ここで時間については上付きの添え字で、位置については下付きの添え字で表している。

\displaystyle - q_0 \frac{\Delta t}{\Delta x} = r

と置くと、

\displaystyle \rho^{m+1}_{n} = (1+r)\rho^{m}_{n} - r \rho^{m}_{n-1}

という式が導け、ある時点の位置ごとの密度から、次の時点の位置ごとの密度が求められる。

R によるシミュレーション

初期条件として、\rho^{0}_{0}=1\rho^{0}_{1}=1、あとは 0 とした場合の解を見る。

パラメータは q_0=0.1 とした.

f:id:abrahamcow:20170923055517p:plain

これは横軸に位置、縦軸に密度を取り、時間を色で表したプロット。

時間の経過とともに密度がだんだん薄まっていき、釣り鐘型になることがわかる。

f:id:abrahamcow:20170923101419p:plain

これは横軸に時間、縦軸に密度を取り、位置を色で表したプロット。

f:id:abrahamcow:20170923055538p:plain

これは横軸に位置、縦軸に時間を取り、密度を色で表したプロット。

library("tidyverse")
diffeq <- function(dx, dt, q0, time_len, u0){ 
  n <- length(u0)
  m <- time_len
  r <- -q0*dt/dx
  
  A <- diag(1+r,n)
  A[cbind(2:n,1:(n-1))] <- -r

  b <- matrix(NA,m,n)
  b0 <- u0
  b[1,] <- u0
  for(i in 2:m){
    b1 <- A %*% b0
    b0 <- b1
    b[i,] <- b0
  }
  colnames(b) <- seq(0,by=dx,length.out=n)
  cbind(time=seq(0,by=dt,length.out=time_len),b)
}

u0=numeric(15)
u0[1:2]<-rep(1,2)
out <-diffeq(dx=1,dt=1,q0=1/10,time_len=50,u0=u0)
out_df <-as_data_frame(out) %>% 
  gather(x,rho,-time) %>% 
  mutate(x=as.numeric(x))

ggplot(out_df,aes(x=x,y=rho,colour=time,group=time))+
  geom_line()

ggplot(out_df,aes(x=time,y=rho,colour=x,group=x))+
  geom_line()

ggplot(out_df,aes(x=x,y=time,fill=rho))+
  geom_tile()+
  scale_y_reverse()

参考文献(読みたい論文)

Polson, N., & Sokolov, V. (2015). Bayesian analysis of traffic flow on interstate I-55: The LWR model. The Annals of Applied Statistics, 9(4), 1864-1888. https://arxiv.org/pdf/1409.6034.pdf

abrahamcow.hatenablog.com