コンテンツにスキップ

05. 東京・大阪の最速降下

東京から大阪まで地下にトンネルを掘り、摩擦のない乗り物を重力だけで走らせます。下りで加速し、上りで減速して、大阪にちょうど止まって着きます。 最も早く着くトンネルはどんな形で、何分かかるでしょうか。

項目内容
状態東京からの水平距離 xx、深さ yy、速さ vv
制御トンネルの向き θ\theta(真下から測った角度)
目的到着時刻 tft_f を最小にする
区間1区間、30ノード、終了時刻は自由
新しく学ぶこと実際の単位(m、s)で問題を書くこと、桁の大きく違う変数を set_unit_* で同じくらいの大きさにそろえること(スケーリング)

例題03 の最速降下を、東京・大阪間の実際の距離 600 km に当てはめた思考実験です。 例題03と違い、終点は出発点と同じ高さの地表です。重力で下って加速し、同じだけ上って減速するので、 エネルギー保存により大阪にはちょうど速さ 0 で着きます。

計算を簡単にするため、次のように大きく単純化しています。

  • 地面は平らで、重力は深さによらず一定の 9.8 m/s² とします(地球が丸いことや、地下深くで重力が変わることは無視します)。
  • 摩擦も空気抵抗もありません。

数値の大きさにも注意が必要です。距離は数十万 m、時間は数百 s、角度は 1 rad 前後と、変数によって桁がまったく違います。 この例では、それぞれの変数の**代表的な大きさ(単位)**をライブラリに教えて、計算しやすくする方法を学びます。

深さ yy は下向きを正、角度 θ\theta は真下から測ります(θ=90∘\theta = 90^\circ が水平、それより大きいと上り)。

x˙=vsin⁡θy˙=vcos⁡θv˙=gcos⁡θ\begin{aligned} \dot{x} &= v \sin\theta \\[2pt] \dot{y} &= v \cos\theta \\[2pt] \dot{v} &= g \cos\theta \end{aligned}

例題03と同じ式です。速さ vv を向きに応じて水平・鉛直に分けたものが位置の変化で、重力のうちトンネルに沿った成分 gcos⁡θg\cos\theta が速さを変えます。 上り(cos⁡θ<0\cos\theta < 0)では減速します。エネルギー保存から、速さは深さだけで決まります(v=2gyv = \sqrt{2 g y})。

x(0)=0,y(0)=0,v(0)=0x(tf)=l,y(tf)=0\begin{aligned} x(0) &= 0, & y(0) &= 0, & v(0) &= 0 \\[2pt] x(t_f) &= l, & y(t_f) &= 0 \end{aligned}

東京(原点)の地表で静止した状態から出発し、l=600l = 600 km 離れた大阪の地表に着きます。

0≤ x≤ly≥00≤ θ≤π\begin{aligned} 0 \le\ & x \le l \\[2pt] & y \ge 0 \\[2pt] 0 \le\ & \theta \le \pi \end{aligned}

どれも途中のすべてのノードで守る経路制約です。東京より手前や大阪より先へは行かず、地表より上にも出ず、逆戻りする向きも取りません。 最適解ではどれも効いていませんが、最適化の途中で変な道筋を探しに行かないよう範囲を絞っています。

J=tfJ = t_f
記号意味値
gg重力加速度 PARAMETERS["g"]9.8 m/s²
ll東京・大阪間の距離 PARAMETERS["l"]600 000 m
hh初期推定で使う最大深さ PARAMETERS["h"]300 000 m

最適化ソルバーは、すべての未知数がだいたい同じ大きさ(1 前後)であるときに最もよく働きます。 600 000 m の距離と 1 rad の角度を同じ感覚で動かそうとすると、ある変数には変化が小さすぎ、別の変数には大きすぎて、うまく収束しません。 そこで各変数の単位を次のように指定します。

変数単位最適解での大きさ(内部の値)
水平距離 xx300 000 m0 〜 2
深さ yy100 000 m0 〜 約 1.9
速さ vv3000 m/s(= 300 000 m ÷ 100 s)0 〜 約 0.65
向き θ\theta1 rad0 〜 約 3.1
時刻 tt100 s0 〜 約 6.2

ライブラリは内部で各変数をこの単位で割った値を最適化します。運動方程式・境界条件・結果はすべて元の単位(m、s)のまま書けて、 単位の指定は答えを変えず、解きやすさだけを変えます(単位をそろえる)。 速さの単位を「距離の単位 ÷ 時刻の単位」にしているのは、x˙=vsin⁡θ\dot{x} = v\sin\theta の両辺の大きさをそろえるためです。

時刻は 0 から 600 s までとします。

  • 水平距離:0 から 600 km まで時刻に比例して増える直線
  • 深さ:前半の15ノードで 0 から 300 km まで直線的に深くなり、後半の15ノードで 300 km から 0 まで戻る「V字」
  • 向き:0 から π\pi(真下から真上)まで時刻に比例して増える
  • 速さ:初期推定を与えず、0 から始める
トンネルの形。横軸は東京からの距離、縦軸は深さで、下ほど深くなります。緑の線が地表、左右の印が東京と大阪です。点は30個のノードでの値、青の線はそれらを通る多項式、黒の破線はサイクロイド(解析解)です。
トンネルの形。横軸は東京からの距離、縦軸は深さで、下ほど深くなります。緑の線が地表、左右の印が東京と大阪です。点は30個のノードでの値、青の線はそれらを通る多項式、黒の破線はサイクロイド(解析解)です。
深さ・速さ・向き(制御)の時間履歴。点はノードでの値、青の線は多項式、黒の破線はサイクロイドの解析解です。深さは km、速さは km/s、向きは度で表示しています。
深さ・速さ・向き(制御)の時間履歴。点はノードでの値、青の線は多項式、黒の破線はサイクロイドの解析解です。深さは km、速さは km/s、向きは度で表示しています。
乗り物がトンネルを進む様子。左はトンネルの形と現在の位置、右は深さと速さの時間履歴で、縦線が現在の時刻です。
  • 所要時間は 632.7 s(約10分20秒)です。解析解は 620.2 s です。
  • トンネルの最も深い点は中間の 300 km 地点で、深さ 186.9 km(解析解 191.0 km)です。
  • 最も深い点で速さは最大の 1.91 km/s(解析解 1.93 km/s)に達します。
  • 向きは真下(0°)から真上(180°)まで、時間にほぼ比例して変わります。東京ではまっすぐ落ち始め、大阪ではまっすぐ上って地表に出ます。

例題03で見たとおり、最速降下曲線はサイクロイド(円が転がるときの円周上の1点の軌跡)です。

x=a(ϕ−sin⁡ϕ)y=a(1−cos⁡ϕ)\begin{aligned} x &= a(\phi - \sin\phi) \\[2pt] y &= a(1 - \cos\phi) \end{aligned}

例題03では最も深い点(ϕ=π\phi = \pi)で終わりましたが、今回は出発点と同じ高さで終わるので、 ϕ=0\phi = 0 から 2π2\pi まで、ひと山ぶんを通ります。そこで x=2πa=lx = 2\pi a = l から円の半径は

a=l2π≈95.5 kma = \frac{l}{2\pi} \approx 95.5\ \text{km}

と決まります。最大深さ・所要時間・最大速さは

ymax⁡=2a=lπ≈191 kmtf=2πag=2πlgvmax⁡=2g ymax⁡\begin{aligned} y_{\max} &= 2a = \frac{l}{\pi} \approx 191\ \text{km} \\[4pt] t_f &= 2\pi\sqrt{\frac{a}{g}} = \sqrt{\frac{2\pi l}{g}} \\[4pt] v_{\max} &= \sqrt{2 g\, y_{\max}} \end{aligned}

です。l=600l = 600 km、g=9.8g = 9.8 m/s² では tf≈620t_f \approx 620 s です。

所要時間が l\sqrt{l} に比例することに注目してください。距離を4倍にしても時間は2倍にしかなりません。 長いトンネルほど深く掘ることになり、そのぶん速く走れるからです。

30ノードの数値解は、所要時間で約 0.04 s(0.01 % 以下)、最大深さで約 1 km の差で解析解と一致しています。 600 km の問題を、単位を指定しただけで例題03と同じように解けています。

単位を指定しないと、ソルバーから見た未知数は 0 から 600 000 までの距離と、0 から 3 程度の角度が混ざったものになります。 距離を 1 m 動かしても角度を 1 rad 動かしても「同じ1歩」とみなされるので、探索がうまく進みません。 単位を指定すると、どの未知数も 1 前後の値になり、例題03と同じ感覚で解けます。

  • PARAMETERS["l"] を 400000(東京・名古屋くらい)にし、所要時間が 2π⋅400000/9.8≈506\sqrt{2\pi \cdot 400000 / 9.8} \approx 506 s、最大深さが約 127 km になることを確かめる。
  • make_problem の unit_x・unit_y・unit_t をすべて 1 にしてスケーリングを無効にする。答えは変わらないはずなのに、ソルバーが失敗(Inequality constraints incompatible)して止まることを確かめる。
  • PARAMETERS["h"](初期推定の最大深さ)を 150000 や 500000 に変え、どちらも所要時間が約 620.2 s(差は 0.03 s 程度)の同じ答えにたどり着くことを確かめる。
ターミナルウィンドウ
uv run python examples/05_brachistochrone_tokyo_osaka.py

examples/05_brachistochrone_tokyo_osaka.py を GitHub で見る