基礎電気電子工学実験・課題G1

© 2018,2019 H. Toyoki


概要説明

1本の道路を多数の車が走る場合の数値シミュレーションを行います.直線的な高速道路のようなもっとも単純な道路においてさえ自然渋滞が発生しうることが理解できる数理モデルを数値計算(コンピュータシミュレーション)してみます。

まずは、下記URLの動画を見てください。

http://iopscience.iop.org/1367-2630/10/3/033001/media

最初に等間隔で走行している車列がちょっとしたきっかけで乱れたとき、元に回復せず部分的に混んだ状態(渋滞クラスタ)が生じてしまうことが見て取れます。高速道路での自然渋滞はこのような現象だと考えられています。

このような1本道の車の動きをあらわす数理モデルとして「最適速度モデル」(Optimal Velocity Model, OVモデル)というものがよく知られており、その派生モデルは現在でも研究されています。

この実習では、OVモデルを数値計算して動作を観察してみます。

最初にjupyterノートブック(以前はipython notebookと呼ばれた)の説明と操作練習を少し行い、そのあと最適速度モデルのシミュレーションを行い、レポートを作成します。

レポート作成は、Wordを用いることにしましょう。メモのために、早速、Wordを立ち上げておいてください。

前処理

プログラムは、Webブラウザ内で動くものだけを使用しますが、ブラウザがInternet Exploreだとうまく動作しないことがあります。

念のため、Windowsのメニュー「設定」から、下の図のように「設定>システム>既定のアプリ」でWebブラウザをGoogle Chromeにしてください。

set chorome

pythonによるプログラムとjupyter notebook

プログラム言語pythonは、最近、データ解析や機械学習などに用いられることが多い人気言語で、さまざまな周辺ツールやライブラリがどんどん作られつつあります。ここでは、プログラムと様々な説明(通常の文章)を混在させることができるjupyter noteboookを利用します。

pythonはC等と異なり、インタプリタ型の言語です。スクリプト言語と呼ばれることもあります。Cのように、機械語に翻訳(コンパイル)してから実行するのではなく、実行段階で一行一行プログラムを読み取って実行していきますから速度は遅いですが、プログラムが簡単です。

プログラム上の大きな違いは、変数宣言をせずにいきなり使うことができる点です。変数の型もその場に応じて解釈されます(その解釈も時間がかる要因です)。

jupyter notebookの使い方

jupyter notebookはWindowsメニューから立ち上げます。ノートブックの編集、実行などの操作はWebブラウザ内で行われます。つまり、jupyterはpythonインタプリタの本体(バックエンド)とWebサーバ及び動的なページ(HTML5)を作成するシステムから成り立っています。

このページもjupyterを使って作成しました。

サンプルプログラムのダウンロードとjupyterの立ち上げ

  • 実習用のフォルダを「ドキュメント」フォルダの下に作ってください。
  • 次のリンクは2つの実習用ノートブックファイルを圧縮したものです。クリックしてダウンロードし、ファイルを解凍して実習用フォルダの下に格納してください。

    実習用ファイル

  • Windowsのメニューからjupyterを立ち上げてください。jupyter notebookは、メインメニュー"Anaconda3"のサブメニューにあります。
  • jupyterが立ち上がったらダウンロードしたフォルダに移動し、ipynbファイルを読み込んでみてください。

jupyterの基本的なメニュー 編集、実行など

ここで、練習用ノートブックを使って、jupyterを操作する練習をします。(新規作成はしなくても構いません。)

ノートブックの読み込み、新規作成

  • 既存ノートブックの読み込みは、ファイル一覧から該当ファイルをクリックします。
  • 新規には"New"から"Python2"を選びます。(自分のパソコンで実習する場合はpython3を利用しましょう)
  • 新規ファイルは"Untitled"という名になっているので、そこをクリックして名前を適当に変更します。

メニュー

  • ノートブックは「セル」に分かれています。
  • アイコンメニューは基本的にはセルに対するものです。

Web上にたくさんのチュートリアルがありますので、参考にしてください。

編集上のTips

jupyterにはプルダウンメニュー以外にもたくさんの機能、機能に対するショートカットキーがあります。"Help"→"Keyboard Shortcuts"で一覧が表示されます。 日本語での紹介は、たとえば、http://qiita.com/masafumi_miya/items/6524dbd227705351a00c

  • 編集をもとに戻す(undo)はCtrl-Z
  • セルカット(アイコンメニューのはさみ)をもとに戻すは、ESCキーを押したあとZ
  • 戻したものを撤回(redo)は、Ctrl-Y

レポートについて

作成及び提出方法

jupyter notebookやpythonに関して試したこと、最適速度モデルのシミュレーションをやってみてわかったことを画像付きで報告してください。

作成は、Wordあるいはjupyter notebookで作成してください。Wordで作成した場合は、かならずpdf形式に変換して提出してください。ipynb形式の場合はそのままで結構です。(jupyterはpdf形式に変換する機能も備わっていますが、現在のところ日本語が入っているとうまく処理できません。工夫次第で可能になりますが、教室のパソコンでは無理です。)

提出は、レポートを電子メールに添付してtoyoki@yamanashi.ac.jp宛に送ってください。

  • 必ず件名及び本文を書くこと。(本文は、「xxxのレポートを提出します。 氏名」くらいでよいですが、何もないのは失礼にあたります。)
  • ファイル名の先頭は"t18ee000"のように学籍番号から始めてください。(アルファベットは小文字で)
  • 締め切りは実習日の次の金曜日の昼までとします。

実習

以下に示す最適速度モデルのシミュレーションプログラムを実行し、「基本的なパラメータ」の値によって交通流がどのように変化するか(車の軌跡や空間パターンのスナップショット)を調べ、レポートを作成しなさい。

Tips

レポートへの画像の貼りこみ

文書への画像の貼りこみは、Wordの場合は、「挿入」→「スクリーンショット」を用いるのが簡便です。

Markdownの記法

Markdownの記法はたとえば、http://qiita.com/tbpgr/items/989c6badefff69377da7 を参照してください。

最適速度モデル

実習時間が少ないので、詳細は理解できなくてもよいが、簡単に説明する。

実際の運転者の行動、車の動作は複雑ですが、もっとも本質的に重要なのは次のようなことだと考えられます。

  • 運転者は前方の車との車間距離($\Delta x$)に応じて速度を決めるとする。それを最適速度と呼び$V(\Delta x)$と書く。
  • 最適速度は車間距離の単調増加関数であり、$\Delta x\rightarrow \infty$で一定値(最高速度)に漸近する。 ov_function.png
  • 加速・減速

    • $v(t) < V(\Delta x)$ならば加速
    • $v(t) > V(\Delta x)$ならば減速
    • 上記をまとめて、加速度は$V(\Delta x(t)) - v(t)$に比例すると定義

    車はN台あって、それぞれの速度と位置を$v_n(t)$, $x_n(t)$とし、番号$n$は進行方向に向かってつけることにし、微分方程式を書いてみましょう。まず、位置の変化は速度そのものですから、

    $$ \frac{dx_n}{dt} = v_n $$

とあらわされます。速度の変化は加速度であり、上記の考えを

$$ \frac{dv_n}{dt} = a \left[ V(\Delta x_n) - v_n\right] $$

のように定義します(右辺は[力/質量]、この場合は人間の操作も含めた[車の駆動力/質量]を表します)。最適速度は$V(0)=0$であり、単調増加して一定値に近づく関数として

$$ V(\Delta x_n) = \tanh (\Delta x_n - c) + \tanh(c) $$

を採用します(上の図)。ただし、$\qquad n=0,1,2, .... ,N-1$であり、$n=N-1$の前は$n=1$、つまり周期的とします。

道路も実験のビデオと同様に周回道路とします。

OVモデルのシミュレーション

OVモデルの数値計算本体のクラス定義

次のセルは、OVモデル(2N次元の微分方程式)を数値的に解くプログラムの本体で、下に述べる可視化プログラムで利用します。

この部分はクラス(関数と変数の集合体)として記述してあります。クラスはオブジェクト指向言語C++やJava等に共通する書法(概念)です。

下の方にある可視化プログラムを実行する前に、直下にセルを実行しておいてください。pythonインタプリタに記憶されます。

In [1]:
"""
OVモデルの本体
"""
import numpy as np
from scipy.integrate import odeint

class ov_model:
    def __init__(self, carNum, L, a, dt, initialize=True):
        self.carNum = carNum
        self.dt =dt
        self.L = L
        self.set_params(L,a)
        if initialize:
            self.initialize_x_v()
        
    def set_params(self,L,a):
        self.a = a
        self.L = L
       
    def initialize_x_v(self): # must call this after set_params
        self.x = np.arange(0, self.L, self.L/self.carNum)        # initial position
        self.v = np.random.random(self.carNum)     # initial velocity
        self.x_v = np.concatenate((self.x,self.v))
        self.x_v_trace = np.empty((self.carNum*2))
        self.t_trace = np.zeros(1) # array of time for plot
        self.present_time = 0.0
        self.x_v_trace[:] = self.x_v       # data for plot
        
    def ov_eq(self, x_v, t, carNum, a, L):
        '''
        r.h.s of the ov equation
        '''
        x = x_v[:carNum]
        v = x_v[carNum:]
        x_fw = np.roll(x,-1)
        delx = (x_fw- x + L) % L
        return np.concatenate((v, a*(np.tanh(delx - 2.0) +np.tanh(2.0) -v)),axis=0)

    def update_unit(self, n_steps, store_trace=False):
        '''
        solving the ov equation by "n_steps" iterations
        '''
        t = np.arange(0, n_steps*self.dt, self.dt)
        new_x_v = odeint(self.ov_eq,self.x_v,t, args=(self.carNum,self.a,self.L))
        self.x_v = new_x_v[-1] % self.L
        self.present_time = self.present_time+self.dt*n_steps
        if store_trace:
            self.x_v_trace = np.vstack((self.x_v_trace, self.x_v))
            self.t_trace = np.vstack((self.t_trace,np.array([self.present_time])))
    
    def shift_x_v_t(self, n):
        if n > self.t_trace.size:
            n=-1
        self.t_trace = self.t_trace[n:]
        self.x_v_trace = self.x_v_trace[n:]

軌跡の描画

説明は後回しにして、次のセルをまずは実行してみる。

これは、横軸を時刻、縦軸を道路上の位置として車の軌跡を描くプログラムです。

主要なパラメータをGUI(Graphical User Interface)でセットできるようになっています。

17行目の"%matplotlib auto"の記述を"auto"の代わりに"inline"とするとノートブックの内部に描画されます。

次のことに注意して試してみてください。

  • パラメータをセットして、"Run Interact"ボタンを押すと計算がスタートします。 初期($t=0$)状態は、位置は等間隔、速度は乱数で与えてあります。
  • "Initialze"のチェックをはずして"Run"させると、初期化することなく続きを見ることができます。("%matplotlib auto"で見るときは、一度描画の画面を閉じてから実行してください。)

  • 基本的なパラメータは、加速減速の即応能力(感度と呼ぶ)$a$と車の密度$\\\rho=N/L$です。理論的な解析からこの2つのパラメータによって振る舞いが2つに分かれることが知られています。下の図の"Stable"領域では均一な走行が安定です。"Unstable"領域では均一な状態が不安定になり部分渋滞が発生します。 phase_dig_ov_small.png

  • シミュレーションは、本来は連続な時間、空間を離散的に区切り計算しますから当然計算誤差を含みます。$a$が小さいか密度が高い場合は計算が難しくなります。そのような場合、本来は起こらないはずの「追い越し」がおきます。いったん追い越しがおこると数値計算はめちゃめちゃになり意味をなさなくなります。
In [2]:
%matplotlib inline
from ipywidgets import interact_manual, Checkbox
from IPython.display import clear_output, display, HTML
model = ov_model(50, 100, 1.5, 0.01, True)

def show_trace_ov(carNum, rho, a, initialize, max_t):
    global model
    unit_steps =10
    dt=0.01
    if initialize:
        model = ov_model(carNum, carNum/rho, a, dt, initialize)
    else:
        model.shift_x_v_t(int(max_t/dt))
    for i in range(int(max_t/dt/unit_steps)):
        model.update_unit(unit_steps, True)
    # plot the traces
    import matplotlib.pyplot as plt
    %matplotlib inline
    fig = plt.figure(figsize=(10,5))
    axes = fig.add_subplot(1,1,1)
    axes.set_xlabel("time $t$")
    axes.set_ylabel("position $x$")
    for i in range(carNum):
        axes.plot(model.t_trace, model.x_v_trace[:,i], '.', markersize=3)
    plt.show()
    
interact_manual(show_trace_ov, carNum=50, rho=(0.1,1.5), a=(0.5, 3.0), initialize=Checkbox(value=True),max_t=(100,300))
Out[2]:
<function __main__.show_trace_ov(carNum, rho, a, initialize, max_t)>

車の動きのアニメーション

実験ビデオに対応する車の動きをシミュレーションによって再現してみましょう。

リアルタイムアニメーション(計算しながら変化を描画すること)はWebブラウザ内ではできないので、"%matplotlib auto"で実行します。

アニメーション作成にはmatplotlibに備わっているanimationというライブラリを用いる。ただし、これはインラインGUIと併用できないので、パラメータはプログラムの初めのところに埋め込みます。(4~6行目)

アニメーションは、静止画を連続的に同じ画面に表示することにより動いているように見せかけるものです。36行目のframe数分だけintervalミリ秒間隔で表示します。速すぎるようであればintervalを長くしてみてください。repeatをTrueにするとエンドレスで繰り返し表示されます。

アニメーションは別ウィンドウが開きそこに表示されます。セルの下に空白の四角がでるだけの時は、コード内のコメントにあるように、別のノートブックを開き、そこで試してみてください。)

In [ ]:
"""
OVモデルのアニメーション表示 (円形の道路上を車が移動する)
新しいウィンドウが開かないときは、FileメニューからNew Notebookを開き、
上の「OVモデルの本体」と書かれたセルと、このセルをコピーして実行してみてください。
違うノートブック間のコピー&ペーストは、ブラウザのコピー、ペースト機能を使ってください。
(jupyterのメニューにあるコピー、ペーストは同じノートブックないだけで機能します。)
"""
%matplotlib auto
import matplotlib.pyplot as plt
from matplotlib import animation

carNum=50
rho=0.5
a=1.2
L=carNum/rho
model = ov_model(carNum, L, a, 0.01, initialize=True)

# create the animation frame and execute
radius = L/np.pi/2.0
fig = plt.figure(figsize=(6,6))
axes = fig.add_subplot(1,1,1)
pt, = axes.plot([], [], "o")
frame_size = radius*1.1
axes.set_xlim((-frame_size, frame_size))
axes.set_ylim((-frame_size, frame_size))
axes.set_title('a='+ str(a)+ ', $\\rho=$'+ str(rho))
    
def init_anim():
    pt, = axes.plot([], [], "o")
    model.initialize_x_v()
    return pt,
    
def anim_update(i):
    n_steps = 20
    model.set_params(L,a)
    model.update_unit(n_steps, False)
    arg_fac = 2.0*np.pi/L
    x = radius*np.cos(model.x_v[:carNum]*arg_fac)
    y = radius*np.sin(model.x_v[:carNum]*arg_fac)
    pt.set_xdata(x)
    pt.set_ydata(y)
    return pt,
    
anim=animation.FuncAnimation(fig, anim_update, init_func=init_anim,
                               frames=1000, interval=20, blit=True, repeat=False)
plt.show()

参考

説明図の作成

本文で説明のために使用した図を作成するプログラムを参考のために示します。関数を描画するだけのものです。なれると手軽に使えて重宝します。

安定・不安定領域の図作成プログラム

In [5]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib auto
x = np.arange(0, 4.0, 0.02)
plt.plot(x, 2./(np.cosh(x-2.0)**2))
plt.xlabel("$\Delta x = \\rho^{-1}=L/N$")
plt.ylabel("$a$")
plt.text(0.3, 1.5, "stable", fontsize=12)
plt.text(1.5, 0.5, "unstable", fontsize=12)
plt.show()
Using matplotlib backend: Qt5Agg

最適速度関数の描画

In [3]:
'''
最適速度関数の描画
'''
import numpy as np
import matplotlib.pyplot as plt
%matplotlib auto
x = np.arange(0, 8.0, 0.02)
plt.plot(x, np.tanh(x-2.0)+np.tanh(2.0))
plt.xlabel("$\Delta x$")
plt.ylabel("$V(\Delta x)$")
plt.title("$V(\Delta x) = \\tanh (\Delta x - 2.0) + \\tanh (2.0)$")
plt.show()
Using matplotlib backend: Qt5Agg

車間距離-速度のアニメーション

次のセルは、横軸に前方車との車間距離、縦軸に速度をとって車の状態を描画するプログラムです。

この図は何を物語っているのだろうか。考えてみましょう。現象をつかみやすいように台数を増やし、描画間隔を短くしてあります。14~16行目や52行目のパラメータを変更して試してみてください。(このプログラムも別ウィンドウ上に結果が表示されます。)

In [7]:
class ov_model_ex(ov_model):
    def __init__(self, carNum, L, a, dt, initialize=True):
        super().__init__(carNum, L, a, dt, initialize)
        
    def get_headway(self):
        x = self.x_v[:carNum]
        x_fw = np.roll(x,-1)
        delx = (x_fw- x + self.L) % self.L
        return delx
    
import matplotlib.pyplot as plt
from matplotlib import animation
%matplotlib auto
carNum=200
rho=0.5
a=1.2
L=carNum/rho
model = ov_model_ex(carNum, L, a, 0.01, initialize=True)

# create the animation frame and execute
fig = plt.figure(figsize=(6,6))
axes = fig.add_subplot(1,1,1)
pt, = axes.plot([], [], "o")

axes.set_xlim((0, 4.0))
axes.set_ylim((0.0, 1.0+np.tanh(2.0)))
axes.set_title('a='+ str(a)+ ', $\\rho=$'+ str(rho))
axes.set_xlabel("$\Delta x$")
axes.set_ylabel("$v$")
# plot optimal velocity
x = np.arange(0,4.0,0.01)
y = np.tanh(x-2.0) + np.tanh(2.0)
axes.plot(x,y,"r-")

def init_anim():
    pt, = axes.plot([], [], "o")
    model.initialize_x_v()
    return pt,
    
def anim_update(i):
    n_steps = 20
    model.set_params(L,a)
    model.update_unit(n_steps, False)
    arg_fac = 2.0*np.pi/L
    x = model.get_headway()
    y = model.x_v[carNum:]
    pt.set_xdata(x)
    pt.set_ydata(y)
    return pt,
    
anim=animation.FuncAnimation(fig, anim_update, init_func=init_anim,
                               frames=3000, interval=10, blit=True, repeat=False)
plt.show()       
Using matplotlib backend: Qt5Agg
In [ ]: