質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.44%
Pythonista

Pythonistaは、iOS上でPythonプログラミングができる開発アプリです。さらに、Pythonの関数・変数などを自動で補完する便利なコードエディタや、PythonスクリプトをiOS上で多様な形で機能させる各種機能も内包しています。

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

MacOS(OSX)

MacOSとは、Appleの開発していたGUI(グラフィカルユーザーインターフェース)を採用したオペレーションシステム(OS)です。Macintoshと共に、市場に出てGUIの普及に大きく貢献しました。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

1回答

1596閲覧

渋滞学 python シミュレーション

araharu

総合スコア1

Pythonista

Pythonistaは、iOS上でPythonプログラミングができる開発アプリです。さらに、Pythonの関数・変数などを自動で補完する便利なコードエディタや、PythonスクリプトをiOS上で多様な形で機能させる各種機能も内包しています。

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

MacOS(OSX)

MacOSとは、Appleの開発していたGUI(グラフィカルユーザーインターフェース)を採用したオペレーションシステム(OS)です。Macintoshと共に、市場に出てGUIの普及に大きく貢献しました。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2023/01/31 06:26

編集2023/01/31 08:18

前提

puthonで交通流の数理モデルからシミュレーション動画を作っています。
下記サイトのものを参考に、下記サイトの真ん中にある動画のようなものを作りたいです。
車の渋滞を表現するモデルで、渋滞が解消されないようにパラメーター(数字が具体的に決められていないもの。Δの部分など)をいじればできるらしいのですが、渋滞がずっと続く形を作れません。
どこをどう直せばいいかなどのアドバイスが欲しいです。

実現したいこと

サイト
https://qiita.com/jabberwocky0139/items/e2526fc5ee3b0dbf144b 
こちらのサイトのコードをまんま使ったものだと、動画は33秒ほどのものでした。

サイトの式と違いですが、以下の数式を使っています。
イメージ説明
Nは台数です

発生している問題・エラーメッセージ

該当のソースコード

python

1import numpy as np 2 3_list = [-1, 2, 0, -3, 5] 4 5def zerobound(x): 6 return max(x, 0) 7 8 9_list = list(map(zerobound, _list)) 10 11 12# 各定数 13V1, V2, C1, C2, Lc = 6.75, 7.91, 0.13, 1.57, 5.0 14λ = 0.5 15 16 17# -> dydx 18def fvdm(v, s_diff): 19 # 各定数 20 V1, V2, C1, C2, Lc = 6.75, 7.91, 0.13, 1.57, 5.0 21 λ = 0.8 22 """ 23 微分方程式の定義(odeint) 24 f = [x, v] 25 """ 26 #インデックススライス 27 # L = [0,1,2,3,4,5] 28 # print(L[2]) 29 # print(L[:2])2まで(終了ポイント) 30 # print(L[1:])1以降(スタートポイント) 31 # print(L[::2])2個ごと(等間隔) 32 33 def V(s): 34 # V1, V2, C1, C2, Lc = 6.75, 7.91, 0.13, 1.57, 5.0 35 return V1 + V2 * np.tanh(C1 * (s - Lc) - C2) 36 37 diff_v = np.append(v, v[0]) 38 # print(v[-1]) 39 diff_v = np.diff(diff_v) 40 # print('diff_v', diff_v) 41 42 """ 43 三項演算子:ifを1行で使える。 44 boolean(ブール)の頭文字は大文字(python) 45 46 result = True if 1 > 0 else Folse 47 1>0ならTrue、他ならFolse 48 49 L = [i for i in range(6)] 50 0~5までのリスト 51 52 """ 53 54 # s_diff = np.array([i if i > 0 else L + i for i in s_diff]) 55 56 57 zero_bound_v = V(s_diff) - v 58 # print(zero_bound_v) 59 # print('boud v', zero_bound_v) 60 # print(V(s_diff)) 61 # print(v) 62 # zero_bound_v = np.array(list(map(lambda x: max(x, 0), zero_bound_v))) 63 64 65 v_dot = K * zero_bound_v + λ*diff_v 66 return v_dot 67 68def fvdm_t(v, s_diff, t): 69 dvdt = fvdm(v, s_diff) 70 71 72 73 74def diff(s): 75 s_diff = np.append(s, s[0]) 76 s_diff = np.diff(s_diff) 77 s_diff %= L 78 s_diff = np.array([i if i > 0 else L + i for i in s_diff]) 79 return s_diff 80 81# t s v 82def runge_kutta(v, s, h): 83 # if h == k and h is None: 84 # x, y = x + h/2, y + k /2 85 86 # k1 87 # dydx = func(x, y) 88 s_diff = diff(s) 89 dydx = fvdm(v, s_diff) 90 # y2 = dydx 91 k1 = dydx * h 92 dydx = fvdm(v + (k1/2), s_diff + (k1 / 2)) 93 k2 = dydx * h 94 k3 = fvdm(v + k2/2, s_diff + k2 /2) * h 95 k4 = fvdm(v + k3/2, s_diff + k3 /2) * h 96 97 y2 = (k1 + 2* k2 + 2* k3 + k4) 98 y2 = y2 / 6 99 y2 = v + y2 100 # for j, i in enumerate(y2): 101 # if i > 15: 102 # print(j, s_diff[j], i) 103 # print((k1 + 2* k2 + 2* k3 + k4) / 6) 104 105 106 return y2 107 108 109def fvdm_solver(x, y, s, h=0.001): 110 y2 = np.array([y]) 111 s2 = np.array([s]) 112 x2 = np.array([0]) 113 114 for i, _ in enumerate(np.arange(0, x, h)): # 0.1 115 dy = runge_kutta(y2[i], s2[i], h) 116 # print(dy) 117 y = y + dy 118 119 120 s_diff = diff(s) 121 # print(s_diff) 122 # print(y) 123 for j, items in enumerate(zip(y, s_diff)): 124 if items[0] > items[1]: 125 # print(items[0], items[1]) 126 if items[0] - items[1] < 0.0: 127 print(items) 128 raise Exception 129 130 131 # if items[i] < items[i] + s_diff 132 if items[1] >= 0: 133 y[j] = items[1] 134 else: 135 y[j] = 0 136 137 138 if items[0] > 30: 139 y[j] = 30 140 if (items[0] < 0.0): 141 y[j] = 0 142 143 # if sum(y) > L: 144 # y = (y / sum(y)) * L 145 # print(y) 146 # print(sum(y)) 147 148 # print(y) 149 150 # print(y) 151 for x, j in enumerate(y): 152 if j < 0: 153 print(x,j) 154 raise Exception 155 156 s = s + y * h 157 s = np.array([i if i < L else i - L for i in s]) 158 159 if all(map(lambda x: x >= 0, diff(s))) == False: 160 raise Exception 161 162 x += h 163 y2 = np.vstack([y2, y]) 164 s2 = np.vstack([s2, s]) 165 x2 = np.append(x2, x) 166 return (x2, y2, s2) 167 168#台数、、、距離 169N, K, c = 20, 0.85, 2 170L = N * 1 171 172# 初期位置 173y = np.arange(N) * L / N 174 175v = [0] + [10] * (N-1) 176# v = np.array([np.random.randint(0, 10) for i in range(N)]) 177 178t = 100 179 180h = 0.005 181 182x2, y2, s2= fvdm_solver(t, v, y, h) 183 184# for i, y in enumerate(s2): 185# print(i, y[0], y[-1]) 186# for i, y in enumerate(y2): 187# # print(max(y)) 188# print(i, y[0]) 189 190# for i in range(100): 191# print(V1 + V2 * np.tanh(C1 * (i - Lc) - C2)) 192 193import matplotlib.animation as animation 194import matplotlib.pyplot as plt 195import seaborn as sbn 196 197r = L / (2 * np.pi) 198fig = plt.figure(figsize=(6, 6)) 199sbn.set_style('ticks') 200ims = [] 201 202for index, x in enumerate(s2): 203 theta = x / r 204 im = plt.plot(r * np.sin(theta), r * np.cos(theta), '.k', markersize=24) 205 ims.append(im) 206 207Writer = animation.writers['ffmpeg'] 208# writer = Writer(fps=15, bitrate=1800) 209writer = Writer(fps=60) 210ani = animation.ArtistAnimation(fig, ims, interval=10) 211 212plt.legend() 213ani.save("output.mp4", writer=writer) 214plt.show()

試したこと

hやtの値を変えて試したりしましたが上手くいかず。。。

補足情報(FW/ツールのバージョンなど)

ここにより詳細な情報を記載してください。

自分の環境は、macbook pro(2019年購入)、コードはVSCで書いてローカル環境で作っています。
時間微分は4次のルンゲクッタ法で求めるというもの(数式理解のためscipy、RK45などのパッケージは使用禁止)。
ソースコードは.py ファイルもしくは.ipynb
動画はmp4形式にしたいです。

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

can110

2023/01/31 06:42

> どこをどう直せばいいかなどのアドバイスが欲しいです。 とりあえずは「コードを理解すればよい」となるかと思うのですが そのコードはご自身で書かれたものでしょうか。あるいはどの程度理解されているでしょうか。 たとえば、コード上の各定数というものが何を表しているかは理解されているでしょうか。
araharu

2023/01/31 06:45 編集

canさま
araharu

2023/01/31 06:46

理解していると思います。勘違いがある可能性などはあるかもしれないですが、自分ではもう何時間も考えてわからないので他人の目からどこがおかしいかを聞きたいというものです
can110

2023/01/31 07:09

ではまずは自身がどこまで理解できているかを他人に伝えるべきかと。 それには他人に説明するのが一番かと思います。 (説明できない=理解できていない) たとえば、コード上の各定数というものが何を表しているか、その数値の単位についてなどを 質問本文に記載すると、他人が理解しやすくなります。
araharu

2023/01/31 07:11

では、説明するので、どこがわからないかを教えていただいていいでしょうか
can110

2023/01/31 07:16

> では、説明するので、どこがわからないかを教えていただいていいでしょうか それは他人に確認するまでもなくやったほうがよいかと思います。 (自分自身の理解度を客観視できるように)
araharu

2023/01/31 07:18

量が多いため、割愛しています。 全部の説明できますよ、という意味ですし、とりあえず勝手に理解できていないというように決めつけているのをやめた方がいいのだと思いますよ^^ 具体的な内容に入れない、説教したいだけならやめていただけると。。
can110

2023/01/31 07:30

「決めつけていない」ので尋ねているだけです。 「全部説明できますよ」といわれても、実際に説明されないと他人には分かりません(よね?) なお具体的なアドバイス(渋滞が発生しやすい条件)は回答に記載していますし 本当にこの問題を解決したいなら、それにそった行動をとることをお勧めします。
araharu

2023/01/31 07:32

>>「全部説明できますよ」といわれても、実際に説明されないと他人には分かりません(よね?) ですので、量が多く割愛しているので具体的にわからない部分を提示していただけると、、、という話です^^
can110

2023/01/31 07:39

とりあえず重要そうなのは「パラメーターをいじればできる」というところだと思うので 「いじれるパラメーター」ってどれで、ぞれぞれどのような意味と単位ですか? といったところでしょうか。
araharu

2023/01/31 07:43

確かにそうですね、すいません
araharu

2023/01/31 07:50

修正しました
can110

2023/01/31 08:02

> 「いじれるパラメーター」ってどれで、ぞれぞれどのような意味と単位ですか? 念のため確認なのですが、修正された内容によって 上記についての説明ができているとお考えでしょうか?
araharu

2023/01/31 08:04

いじれるパラメータがどれか。それの意味は書いてると思いますが。 単位をどう回答して欲しいかわかりませんが、距離とか書いている時点でわかると思います
can110

2023/01/31 08:13

分かりません。 単位によって結果はまったく変わるので数値シミュレーションではそれを明確にするのは重要です。 あと「K=~」の行について、なぜか説明がまったくないですが。気づいてない感じでしょうか。
araharu

2023/01/31 08:16

単位って何を聞いているんですか?メートルとか答えて欲しいのですか??? いじれるパラメータの意味を聞かれているのに、何を答えるのでしょうか。 具体的な数字が決まっているところはいじれませんが。
araharu

2023/01/31 08:16

>>「いじれるパラメーター」ってどれで、ぞれぞれどのような意味と単位ですか?
can110

2023/01/31 08:47

> 単位って何を聞いているんですか?メートルとか答えて欲しいのですか??? はい。「mm」「cm」「m」あたりかとは思いますが、あるいは単位なし(無次元数)ですか? > 具体的な数字が決まっているところはいじれませんが。 提示画像には「2.パラメータなど」と書かれていて どれがパラメータでどれがそうでないかが明確には分からないと思います (もしかしたらわざとそう書いてるのかもしれませんが)。 つまり「K=~」の行の値は「いじれるパラメータ」ではないので説明もしていないということですね。
araharu

2023/01/31 08:51

>>つまり「K=~」の行の値は「いじれるパラメータ」ではないので説明もしていないということですね。 そうですね、 ちなみに画像は以前に自分のためにまとめたものなので、今回「パラメータ」についての話になったのは想定していなかったため紛らわしかったかもですね。 単位はなしです サイト参照して貰えばわかりやすいかもですね
can110

2023/01/31 08:58

了解しました。 なお私は「s,v,x,t」はシミュレーションによって動的に変わる「状態量」であって これら自体を直接いじりようがないモノと解釈しています。 この解釈からが異なるようなので、この質問についてはこれ以上、質問・修正依頼しません。 これまでのやり取りから少しでもご自身の理解度および問題が具体的に第三者に伝わるように修正されていれば幸いです。(なっていなかったらごめんなさい)
araharu

2023/01/31 09:02

>>なお私は「s,v,x,t」はシミュレーションによって動的に変わる「状態量」であって これら自体を直接いじりようがないモノと解釈しています。 そうですよね、なので初期値的なことを言われているのだと思っています。 それについては違うでしょうか
melian

2023/01/31 09:36

処理時間がかなり長いので、 t = 100 h = 0.005 を、 t = 10 h = 0.1 に変更して、 theta = x / r を theta = x / np.pi に変更して実行すると、時計で言うと0時〜1時付近に渋滞らしき?状態が発生していることが分かるかと思います。
araharu

2023/01/31 09:41

この続き?が見たいイメージです。伝わりづらかったらすいません。。。
melian

2023/01/31 09:47

アニメーションを見ていると、影響しているのは以下の様な気もします。(他にもあるかもしれませんが) s = np.array([i if i < L else i - L for i in s])
srsnsts

2023/02/02 05:31

こんにちわ。 すみません、ちょっと気になることがあるので質問してもいいですか? 2つある数式のうちの、最初の1つめの式でラムダが付いている項に、 △がありますけど、これ、ラプラス演算子ですか?
guest

回答1

0

渋滞が解消されないようにパラメーターをいじればできるらしいのですが、渋滞がずっと続く形を作れません。

  • 近くの状況、問題しか見ておらず、そこで発生した問題だけを手っ取り早く雑に解決しようとする。(前1台の状況しか見ない。急加減速する。)
  • 全体の流れや状況を理解せずに、とにかく目的に早く行きたがる。(適度な車間距離をとらない。)

といった行動をとるとダメだったと思います。
あえてそのような行動をとるようにパラメータを調整すればよいでしょう。

ダメというのは渋滞が発生しやすいということです。念のため。
参考:渋滞学

投稿2023/01/31 07:10

編集2023/01/31 13:21
can110

総合スコア38317

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

まだベストアンサーが選ばれていません

会員登録して回答してみよう

アカウントをお持ちの方は

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.44%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問