土星に直接ホーマン遷移で行こうとすると、地球出発時に約 $\Delta v \approx 10.3$ km/s の増速が必要になります。当時最大級のロケットでも、土星に乗り込めるほどの探査機を打ち上げることはできません。それでもCassini-Huygensは1997年に打ち上げられ、2004年に土星周回軌道に投入されました。何が起きたのでしょうか。答えは、金星 → 金星 → 地球 → 木星と4回のフライバイを織り交ぜることで、ロケットでは出せない $\Delta v$ を「天体からタダで貰ってきた」のです。この経路はVVEJGA（Venus-Venus-Earth-Jupiter Gravity Assist）と呼ばれ、惑星探査史でもっとも有名な軌道設計の一つです。
多回フライバイ軌道は今やほぼ全ての深宇宙ミッションで使われています。NASAのLUCYは2021年打ち上げから2033年までに地球を3回・小惑星帯を抜けて木星トロヤ群の7天体を訪問する計画です。ESA/JAXAのBepiColomboは水星到達に金星2回・水星6回のスイングバイを用い、ESAのJUICEは木星圏で35回ものガリレオ衛星フライバイを行いガニメデ周回軌道に到達します。Parker Solar Probeは金星7回のフライバイで近日点を $9.86~R_\odot$ まで下ろし、太陽コロナに最接近します。これらの経路は、もはや手計算では設計できません。組合せ爆発、強い非線形性、複数の制約 — それを解くのが本記事の主題です。
本記事の内容
- パッチドコニックス近似とその限界、全動力学への接続
- Tisserandパラメータと等価エネルギー — フライバイが「保存する量」と「変える量」
- フライバイ・シーケンスの組合せ爆発と探索空間の構造
- Lambert問題を繰り返し解く「Multiple Gravity Assist (MGA)」定式化
- 分岐限定法、進化アルゴリズム、深層強化学習による経路探索
- 実例: Cassini (VVEJGA), LUCY, BepiColombo, JUICE, Parker Solar Probe
- Pythonで簡易ephemeris + 進化計算によるVVEJGA経路の再現
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
なぜ多回フライバイなのか — Cassiniの数字で見る効率
直接遷移とフライバイを使う遷移の差を、具体的な数字で押さえておきます。地球から土星へ直接ホーマン遷移する場合、出発時の双曲超過速度は
$$ v_\infty^{\text{direct}} = \sqrt{\frac{\mu_\odot}{a_E}}\left(\sqrt{\frac{2 a_S}{a_E + a_S}} – 1\right) \approx 10.29~\text{km/s} $$
となります（$\mu_\odot$ は太陽の重力パラメータ、$a_E, a_S$ は地球・土星の軌道半径）。地球からの脱出を含めると、地表での必要 $\Delta v$ はおよそ $15.6$ km/s。これは1990年代のロケットでは1.5トン程度しか打ち上げられません。
ところがCassiniの実際の打ち上げ $C_3$（地球脱出時の双曲超過速度の二乗）は約 $16.6$ km²/s²、つまり $v_\infty \approx 4.1$ km/s でした。地球脱出を含めても $\Delta v \approx 11.4$ km/s で済み、6トン級の探査機を送り出すことができたのです。この差「$15.6 – 11.4 = 4.2$ km/s」がフライバイによる節約であり、それは金星と地球からの「重力的なヒッチハイク」で稼ぎ出されました。
この節約を生み出す原理は、太陽中心系から見ると単純です。惑星に対する双曲線軌道で接近・離脱した探査機は、惑星に対しては速さを変えませんが、太陽から見ると「惑星の公転速度ベクトルとの相対角」が変わるため、太陽中心速度が大きく変化します。詳細はグラビティアシストの軌道力学で扱いましたが、本記事ではフライバイを「ブラックボックス」として、その連鎖をどう設計するかに焦点を当てます。
ここから自然な問いが生まれます — どの惑星に、いつ、何回フライバイすれば、目的地まで最小の $\Delta v$ で行けるのか。これは離散的な「経路選択」と連続的な「タイミング選択」が組み合わさった複雑な最適化問題です。次節で、まず軌道力学の側で何が成り立っているかを整理しましょう。
パッチドコニックス近似と全動力学
パッチドコニックスの考え方
惑星間航行をまともに扱うには、太陽と全惑星の重力を同時に受ける $N$ 体問題を解く必要があります。しかし、太陽系では幸いなことに各惑星の重力影響圏（SOI: Sphere of Influence）が惑星間距離に比べて十分小さく、ほとんどの時間は「太陽だけが支配する2体問題」として扱えます。
惑星 $i$ のSOI半径は、Laplaceの近似で
$$ r_{\text{SOI},i} = a_i \left(\frac{m_i}{M_\odot}\right)^{2/5} $$
で与えられます。地球で約 $0.92 \times 10^6$ km、木星で約 $48 \times 10^6$ km、ガニメデで約 $24{,}000$ km です。地球-太陽距離 $1.5 \times 10^8$ kmと比べると、地球SOIは0.6%にすぎません。
そこでパッチドコニックス近似では、軌道を次のように区分けします。
- 太陽のSOI内（惑星のSOI外）: 太陽中心の2体問題で楕円/双曲線軌道
- 惑星のSOI内: 惑星中心の2体問題で双曲線フライバイ
- SOI境界で軌道を「貼り合わせる」(patch)
これを多回フライバイに適用すると、軌道全体は次のような連鎖になります。
$$ \underbrace{\text{Sun-centric arc}_1}_{\text{Earth to Venus}} \to \underbrace{\text{Venus flyby}}_{\text{conic patch}} \to \underbrace{\text{Sun-centric arc}_2}_{\text{Venus to Venus}} \to \cdots $$
各太陽中心アークはLambert問題（2点と通過時刻から軌道を決定する問題）で解けます。各フライバイは、入る前後の双曲超過速度ベクトル $\bm{v}_\infty^-, \bm{v}_\infty^+$ を、惑星の重力で曲げられる範囲で結びつける制約として扱います。
パッチドコニックスの限界
この近似は実用的ですが、いくつかの限界があります。第一に、SOIの境界で軌道を「瞬間的に貼り合わせる」ため、SOI境界での力学的整合性が完全には保たれません。実際にはSOIに入る前後で他天体の重力も無視できず、特にL1/L2点付近を通る低エネルギー遷移ではこの誤差が致命的になります。第二に、深宇宙マヌーバ (DSM: Deep Space Maneuver) と呼ばれる小さなインパルス推進を太陽中心アークの途中で加えると、Lambert問題が2つに分裂し、定式化が複雑になります。
精密な設計段階では、全動力学（太陽 + 各惑星 + 月 + 摂動）でのトラジェクトリ伝搬と、間接法または直接法による最適化を組み合わせます。代表的なフレームワークとしてNASA/JPLのMystic、ESAのGTOC候補、オープンソースではpykep, poliastroなどがあります。
しかし、いきなり全動力学で多次元探索するのは現実的ではありません。まずパッチドコニックスで概形を作り、それを初期推定値として全動力学最適化に渡す、という二段階のアプローチが標準です。本記事では前者に集中し、Pythonでスクラッチ実装します。
次に、フライバイが軌道に何をするかをより深く理解するため、Tisserandパラメータという保存量を導入します。これは多回フライバイの「探索空間」を可視化する強力な道具です。
Tisserandパラメータと等価エネルギー
Tisserandパラメータの定義
円制限三体問題（CR3BP: Circular Restricted Three-Body Problem）には、Jacobi定数と呼ばれる厳密な保存量があります。それを楕円軌道近似で表したものがTisserandパラメータ $T$ です。惑星 $P$（軌道半径 $a_P$、円軌道仮定）に対して、探査機の軌道要素 $(a, e, i)$ を使って
$$ T = \frac{a_P}{a} + 2\sqrt{\frac{a(1-e^2)}{a_P}}\cos i $$
と定義されます（$a, e, i$ は探査機の太陽中心軌道の長半径、離心率、惑星 $P$ 軌道面からの傾斜角）。
このパラメータの物理的意味は深いものです。CR3BPでは、惑星に対するフライバイの前後でJacobi定数が保存します。パッチドコニックス近似でも、フライバイは惑星に対する $v_\infty$ の大きさを保存し、方向のみ変えます。したがってフライバイの前後でTisserandパラメータも保存します（パッチドコニックスの精度で）。
これは強力な制約です。「ある惑星でフライバイを繰り返している限り、探査機の軌道は等-Tisserand曲線の上を動く」ことを意味します。逆に、別の惑星でフライバイすると、その惑星に対するTisserandが新たに保存量として効き始め、軌道は別の曲線へ「乗り換え」られます。
等価エネルギーとしての解釈
Tisserandパラメータは惑星に対する双曲超過速度 $v_\infty$ と次の関係を持ちます。惑星の公転速度を $v_P = \sqrt{\mu_\odot/a_P}$、探査機の太陽中心軌道エネルギー（特定エネルギー）を $\epsilon = -\mu_\odot/(2a)$ とすると、惑星位置での探査機の太陽中心速度 $v$ について
$$ v^2 = v_P^2 + v_\infty^2 + 2 v_P v_\infty \cos\beta $$
が成立します（$\beta$ は $\bm{v}_P$ と $\bm{v}_\infty$ のなす角）。これとビザビバ方程式を組み合わせると、
$$ T = 3 – \left(\frac{v_\infty}{v_P}\right)^2 $$
という簡潔な関係が得られます（平面・円軌道近似）。つまりTisserandパラメータは惑星に対する $v_\infty$ の二乗を、3から引いた量として解釈できます。$v_\infty$ が大きいほど $T$ は小さく、惑星に対して相対速度が高速の軌道に対応します。フライバイで $v_\infty$ の大きさは変わらないから $T$ も変わらない、というわけです。
Tisserandグラフによる経路設計
横軸に近日点距離 $r_p$（または周期）、縦軸に遠日点距離 $r_a$ を取り、等-$T$ 曲線を惑星ごとに描いたものがTisserandグラフです。
- ある惑星でフライバイを続ける限り、軌道点はその惑星の等-$T$ 曲線上を動く
- 別の惑星に乗り換えるには、両惑星の等-$T$ 曲線が交差する点を通る必要がある
- 交差点では、両惑星の軌道半径と一致する点（惑星の公転半径）を通ることも必要
これにより、複雑な軌道探索が「グラフ上の点と曲線の幾何学」に還元されます。たとえばJUICEのガリレオ衛星クルージング設計では、ガニメデ・カリスト・エウロパの等-$T$ 曲線を木星周回軌道のTisserandグラフ上に描き、その交差パターンから35回フライバイのシーケンスを設計しています。
Tisserandパラメータは静的な保存量ですが、フライバイ間にDSMを入れると $v_\infty$ が変わり、$T$ も変化します。つまりDSMの大きさは Tisserand曲線間の「乗換距離」として定量化できます。これがLeapfrogging（蛙跳び）法と呼ばれる経路設計の基本原理です。
ここまでで「フライバイが何を変え、何を保存するか」がわかりました。次に、これらを組み合わせて目的地に到達するシーケンス全体を、どう探索するかを見ていきます。
フライバイ・シーケンスの組合せ爆発
探索空間の大きさ
多回フライバイ問題は、二つの要素から成ります。
- 離散的選択: 各フライバイで訪問する天体のシーケンス $(P_1, P_2, \dots, P_n)$
- 連続的選択: 出発時刻、各フライバイ時刻、最終目的地到着時刻、DSMの位置と大きさ
惑星間ミッションを想定し、対象天体を金星・地球・火星・木星の4つに絞っても、$n=4$ 回のフライバイなら $4^4 = 256$ 通りの離散シーケンスがあります。$n=8$ なら $4^8 \approx 6.5 \times 10^4$ 通り、$n=16$ なら $4 \times 10^9$ 通り超え。実用的にはLUCYのように8回、JUICEのように35回フライバイがあり、組合せの数だけで天文学的になります。
しかも各シーケンスについて、連続変数の最適化（10〜20次元）を解く必要があります。シーケンス $\times$ 連続最適化 — これが多回フライバイ問題の本質的な難しさです。
制約による枝刈り
幸い、物理制約により大半のシーケンスは即座に棄却できます。代表的な制約は次のとおりです。
- エネルギー制約: 各アークでLambert問題が解を持つ範囲（軌道楕円が両惑星位置を通る）
- タイミング制約: 各天体の周期と、フライバイ間アークの周期との整合性 (resonant flyby)
- フライバイ高度制約: 惑星表面に衝突しない（地球で $h_{\min} > 300$ km、木星で $> 0$ など）
- $v_\infty$ 連続性: フライバイ前後で $|\bm{v}_\infty|$ が保存（DSMがある場合は別）
- ミッション期間制約: 全体の飛行時間が上限以内
- 太陽距離制約: 熱設計の観点から、内側惑星より太陽に近づかない
これらを使うと、4惑星×16段で $10^9$ 超のシーケンスが、現実的には数千から数万に絞られます。たとえばCassiniの設計時には、地球出発・土星到着で2-4年程度の打上ウィンドウに、フライバイ回数3-6で限定し、約 $10^4$ オーダーのシーケンスを評価したと言われています。
共鳴フライバイ (resonant flyby)
特に重要なのが共鳴フライバイです。たとえばCassiniは1998年と1999年に2回金星フライバイをしましたが、これは金星-探査機間で $1:1$ 共鳴（探査機が太陽の周りを1周する間に金星も1周）が成り立つよう設計されたためです。共鳴フライバイは、探査機の軌道周期 $T_{\text{sc}}$ と惑星の周期 $T_P$ が
$$ \frac{T_{\text{sc}}}{T_P} = \frac{p}{q}, \quad p,q \in \mathbb{Z}^+ $$
の有理比を満たすときに発生します。$1:1$, $2:1$, $3:2$ などが多用されます。共鳴フライバイは、同じ惑星で時間を空けて再会できるため、シーケンス設計の自由度を大きく増やします。$p:q$ 共鳴を組み合わせる手法は、$v_\infty$-leveraging maneuver とも呼ばれ、共鳴間DSMを使って効率的に $v_\infty$ を増大できます。
ここまでで離散探索空間の構造がわかりました。次に、各シーケンスの中で連続最適化をどう定式化するかを見ます。これがMGA問題です。
MGA問題の定式化 — Lambert問題の繰り返し
MGA-1DSM定式化
最も基本的なのは「フライバイは瞬間的なベクトル方向変換、DSMは各アークの中ほどに1回だけ」と仮定するMGA-1DSM (Multiple Gravity Assist with 1 Deep Space Maneuver per leg) です。
シーケンス $(P_0=$ Earth $, P_1, P_2, \dots, P_n=$ target $)$ が与えられたとき、最適化変数は次のようになります。
- $t_0$: 地球出発時刻
- $T_i$ ($i=1,\dots,n$): $i$ 番目のアークの飛行時間
- $\eta_i \in (0,1)$: $i$ 番目のアークでDSMが入るタイミング（無次元）
- $\bm{v}_0 \in \mathbb{R}^3$: 地球出発時の双曲超過速度ベクトル（極座標 $v_\infty, \alpha, \beta$ で3変数）
全部で $3 + 2n + 1$ 次元程度。$n=4$ なら $12$ 次元です。
評価関数の流れは次のようになります。
- $t_0, \bm{v}_0$ から地球出発状態 $\bm{r}_0, \bm{v}_0^{\text{sc}}$ を計算
- 各アーク $i$ について: * $P_{i-1}$ から $\eta_i T_i$ 時間進めた状態 $\bm{r}_i^{\text{mid}}$ をケプラー伝搬 * DSM地点 $\bm{r}_i^{\text{mid}}$ から $P_i$ 位置 $\bm{r}_{P_i}(t_0 + \sum T_j)$ までを残り時間 $(1-\eta_i)T_i$ でLambert問題を解く * DSM $\Delta v_i^{\text{DSM}} = \bm{v}_i^{\text{mid,after}} – \bm{v}_i^{\text{mid,before}}$ を記録
- 各フライバイ $P_i$ で: * 入射 $v_\infty^- = \bm{v}_{\text{sc}}^{-} – \bm{v}_{P_i}$ * 離脱 $v_\infty^+ = \bm{v}_{\text{sc}}^{+} – \bm{v}_{P_i}$ * $|\bm{v}_\infty^-| = |\bm{v}_\infty^+|$ 制約と、フライバイ高度 $h_p$ から角度制約を満たす
- 目的地での到着 $v_\infty$ から、周回投入 $\Delta v$ を計算
総コストは
$$ J = \Delta v_{\text{launch}} + \sum_i \Delta v_i^{\text{DSM}} + \Delta v_{\text{arrival}} $$
を最小化します。各項は軌道マヌーバとデルタvで扱った形式です。
Lambert問題の役割
各アークでLambert問題を1〜2回解くため、$n$ 段ミッションではLambert問題を $O(n)$ 回呼びます。Lambert問題の解は2点と飛行時間から軌道を一意に決め（多重周回を除く）、軌道要素 $(a, e, i, \Omega, \omega, \nu)$ と両端での速度ベクトルを返します。詳細はLambert問題と惑星間軌道で扱いました。
実装ではGoodingやIzzoのアルゴリズムを使うのが標準で、収束は数回のNewton反復で得られます。10次元の最適化を $10^5$ 回評価しても、Lambert問題が高速なら数分程度で計算できます。
フライバイ制約の処理
フライバイの $|\bm{v}_\infty|$ 連続性は等式制約、フライバイ高度は不等式制約として効きます。違反した解はペナルティ項で評価値を悪化させ、最適化アルゴリズムに棄却させます。たとえば
$$ J’ = J + \lambda_1 \sum_i \left| |\bm{v}_\infty^{-,i}| – |\bm{v}_\infty^{+,i}| \right| + \lambda_2 \sum_i \max(0, h_{\min} – h_{p,i}) $$
のような形です。$\lambda$ は実装で経験的に調整します。あるいはペナルティを陽に書かず、フライバイ角度を変数に追加して常に連続性が満たされる定式化（MGA-2DSMなど）もあります。
各シーケンスについてこの $J$ を最小化すれば、最良のフライバイ時刻が得られます。問題は、$J$ が非凸・多峰性で局所最適解だらけということです。次節で、これを克服する大域最適化手法を見ていきます。
大域最適化手法
分岐限定法とTisserandベース探索
シーケンスを離散的に決める段階では分岐限定法 (Branch-and-Bound) が古典的です。Tisserandグラフ上で、出発惑星から目的惑星までの「等-$T$ 曲線交差パターン」を木探索し、各ノードで下限値（必要 $\Delta v$ の下界）が現在の最良解を上回ったら枝刈りします。手作業設計の延長で、設計者のintuitionに対応します。
しかしフライバイ回数 $n \geq 6$ では木が爆発するため、近代的にはメタヒューリスティクスが主流です。
進化アルゴリズム (Differential Evolution)
連続最適化部分では差分進化 (Differential Evolution, DE) が圧倒的に普及しています。各シーケンスの12〜20次元最適化に対し、人口100〜200個体・1000〜10000世代のDEを回します。Storn & Priceの古典的DE/rand/1/binや、JADE, SHADEなど自己適応版が使われます。scipy.optimize.differential_evolution の実装でも多くのMGA問題が解けます。
DEはランダム性が強いため、複数シードで実行し最良値を採用します。Cassini設計の再現問題では、人口200・世代5000程度で実際のミッション値（$J \approx 4.93$ km/s）の数 %以内まで近づくことが報告されています。
多目的最適化 (NSGA-II, MOEA/D)
実際のミッションは「総 $\Delta v$ 最小」と「飛行時間最小」の二つを同時に下げたい場合が多く、多目的進化アルゴリズムが使われます。NSGA-II (Non-dominated Sorting GA) は、パレートフロント上の解集合を1回の最適化で得る代表的手法で、ESA-ACTのPaGMOやpykepライブラリにも実装されています。
深層強化学習 (Deep RL)
近年（2020年以降）、フライバイシーケンスの離散選択に強化学習を適用する研究が増えています。状態を「現在の軌道要素 + 残り推進剤 + 残り時間」、行動を「次にフライバイする天体 + DSM大きさ」とし、報酬を「目的地への到達 – $\Delta v$ コスト」で定義します。
Izzoらの2020年研究では、PPO（Proximal Policy Optimization）を用いて木星トロヤ群ミッション設計を学習させ、人手設計に匹敵する経路を自動生成しています。最近では Transformer ベースのアーキテクチャでフライバイ列を生成する試みもあります。これらは「シーケンスの離散選択 + 連続最適化を end-to-end で学習」という新しいパラダイムを開拓中です。
Lazy Race Trajectory Model (GTOC)
毎年開催されるGlobal Trajectory Optimization Competition (GTOC) では、多回フライバイ問題が課題として出され、各国の研究機関が解を競います。GTOC優勝チームの多くは、Tisserandベース粗探索 + DE + MBH (Monotonic Basin Hopping) + 全動力学微調整、という4段パイプラインを採用しています。多回フライバイ最適化は、宇宙工学とORの境界領域として最も活発な研究分野の一つです。
ここまでで手法面の概観を得ました。次に、実際のミッションでこれらがどう使われているかを見ます。
実例: フライバイ・ミッションの軌道設計
Cassini-Huygens (1997-2017): VVEJGA
- 経路: Earth → Venus (1998-04) → Venus (1999-06) → Earth (1999-08) → Jupiter (2000-12) → Saturn (2004-07)
- 全飛行時間: 6.7年
- 直接遷移比 $\Delta v$ 節約: $\approx 4.2$ km/s
- キー設計判断: 1:1 Venus-Venus共鳴で2回目の金星フライバイを実現
Cassiniは多回フライバイミッションの黄金例で、上記のVVEJGA経路は1990年代初頭にJPLのRoger Diehl, Daniel Roth, Aron Wolfらが設計しました。当時の計算機資源では分岐限定法と手作業最適化が主体で、 Tisserandグラフを駆使したと記録されています。
LUCY (2021-2033): 7トロヤ群フライバイ
- 経路: Earth → Earth (2022-10) → Earth (2024-12) → DonaldJohanson (2025-04) → Eurybates+Queta (2027-08) → Polymele (2027-09) → Leucus (2028-04) → Orus (2028-11) → Earth (2031-12) → Patroclus+Menoetius (2033-03)
- 全飛行時間: 12年
- 訪問天体: 7個のトロヤ群小惑星 + メインベルト1個 = 計9天体
LUCYはトロヤ群の L4とL5の両方を訪問する野心的なミッションで、3回の地球フライバイで $v_\infty$ を段階的に増大させトロヤ群へ達します。設計には進化アルゴリズムが使われ、5000以上の候補シーケンスから現行解が選ばれました。
BepiColombo (2018-2026): 金星2回+水星6回
- 経路: Earth → Earth (2020-04) → Venus (2020-10) → Venus (2021-08) → Mercury×6 (2021-10〜2025-12)
- 全飛行時間: 7.2年
- 課題: 太陽からの強い熱とイオンエンジンとの併用
BepiColomboは電気推進（イオンエンジン）を併用するため、純粋なパッチドコニックスでは記述できません。低推力区間を含むlow-thrust MGAとして定式化され、間接法（Pontryagin最大値原理）または直接法（Sims-Flanagan transcription）で解かれます。
JUICE (2023-2031+): ガリレオ衛星35回フライバイ
- 経路: Earth (2024-08) → Venus (2025-08) → Earth (2026-09) → Earth (2029-01) → Jupiter (2031-07) → 35×衛星フライバイ → Ganymede orbit (2034)
- 木星圏での衛星フライバイ: Ganymede 12回, Callisto 21回, Europa 2回
- 設計手法: Tisserand-Poincaré graph による2次元探索 + 進化計算微調整
木星到達後、JUICEはガニメデ周回軌道までの「ツアー」を多回フライバイで設計します。木星圏は地球-月系より複雑な3体問題で、Tisserand-Poincaréグラフ（位相空間断面）を使った設計が標準です。
Parker Solar Probe (2018-): 金星7回降下
- 経路: Earth → Venus×7 (2018〜2024) → 近日点 $9.86~R_\odot$
- 太陽コロナ最接近速度: 約 192 km/s（人工物史上最速）
PSPは「金星でエネルギーを下げる」フライバイの極端な例です。通常のフライバイは加速・方向転換に使われますが、PSPでは太陽中心軌道のエネルギーを下げ、近日点を太陽に近づける減速フライバイです。Tisserandグラフ上で見ると、金星等-$T$ 曲線に沿って近日点を内側に「滑り降りる」操作です。
これらの実例を見れば、多回フライバイ設計の幅と奥深さがわかります。最後に、これをPythonで実装してみましょう。Cassini風VVEJGA経路を、簡易ephemeris + 差分進化で再現します。
Pythonでの実装 — 簡易MGAソルバ
簡易ephemeris
JPL HORIZONSのような厳密な値ではなく、円軌道近似で各惑星の位置・速度を返す関数を作ります。実際のミッション設計でも、初期探索は円軌道近似で十分です。
import numpy as np
# 太陽の重力パラメータ [km^3/s^2]
MU_SUN = 1.32712440018e11
# 惑星パラメータ: 軌道半径[km], 周期[s], 重力パラメータ[km^3/s^2], SOI[km]
PLANETS = {
'Earth': {'a': 1.495979e8, 'T_period': 365.256*86400, 'mu': 3.986e5, 'soi': 9.24e5, 'phase0': 0.0},
'Venus': {'a': 1.082089e8, 'T_period': 224.701*86400, 'mu': 3.249e5, 'soi': 6.16e5, 'phase0': 0.7},
'Jupiter': {'a': 7.785472e8, 'T_period': 4332.59*86400, 'mu': 1.2669e8, 'soi': 4.82e7, 'phase0': 2.1},
'Saturn': {'a': 1.432041e9, 'T_period': 10759.22*86400,'mu': 3.793e7, 'soi': 5.45e7, 'phase0': 4.5},
}
def planet_state(name, t):
"""惑星の位置・速度ベクトル (太陽中心慣性系, 黄道面). 円軌道近似."""
p = PLANETS[name]
n = 2.0 * np.pi / p['T_period'] # 平均運動 [rad/s]
theta = p['phase0'] + n * t # 真近点角(円軌道なので平均近点角)
r = p['a']
pos = np.array([r * np.cos(theta), r * np.sin(theta), 0.0])
v_circ = np.sqrt(MU_SUN / r)
vel = np.array([-v_circ * np.sin(theta), v_circ * np.cos(theta), 0.0])
return pos, vel
ここでは黄道面を $z=0$ とし、円軌道を仮定しています。phase0 で初期位相を与えることで、実際の惑星配置を模擬します（Cassini打上時1997年10月の各惑星位置に大まかに合わせた値）。実際のephemerisにはJPL DE440等を使いますが、本記事では設計原理を示すのが目的なのでこれで十分です。
Lambert問題のソルバ
Lambert問題は別記事で詳しく扱うため、ここでは Izzo の手法を参考にした単純な反復ソルバを示します。
import numpy as np
def lambert_universal(r1, r2, tof, mu, prograde=True, max_iter=200, tol=1e-9):
"""Universal Variables法でLambert問題を解く.
r1, r2: 始点・終点位置ベクトル [km]
tof: 飛行時間 [s]
mu: 重力パラメータ [km^3/s^2]
return: 始点速度ベクトル v1, 終点速度ベクトル v2 [km/s]
"""
r1n, r2n = np.linalg.norm(r1), np.linalg.norm(r2)
cos_dnu = np.dot(r1, r2) / (r1n * r2n)
cross = np.cross(r1, r2)
# 進行方向 (prograde: dnu in [0, 2pi))
if prograde:
dnu = np.arccos(np.clip(cos_dnu, -1, 1)) if cross[2] >= 0 else 2*np.pi - np.arccos(np.clip(cos_dnu, -1, 1))
else:
dnu = 2*np.pi - np.arccos(np.clip(cos_dnu, -1, 1)) if cross[2] >= 0 else np.arccos(np.clip(cos_dnu, -1, 1))
A = np.sin(dnu) * np.sqrt(r1n * r2n / (1 - np.cos(dnu)))
# 反復: psi (Universal variable) を探す
psi = 0.0
psi_up, psi_low = 4*np.pi**2, -4*np.pi
for _ in range(max_iter):
if psi > 1e-6:
c2 = (1 - np.cos(np.sqrt(psi))) / psi
c3 = (np.sqrt(psi) - np.sin(np.sqrt(psi))) / np.sqrt(psi**3)
elif psi < -1e-6:
c2 = (1 - np.cosh(np.sqrt(-psi))) / psi
c3 = (np.sinh(np.sqrt(-psi)) - np.sqrt(-psi)) / np.sqrt(-psi**3)
else:
c2, c3 = 0.5, 1/6
y = r1n + r2n + A*(psi*c3 - 1) / np.sqrt(c2)
if A > 0 and y < 0:
psi_low = psi
psi = (psi_up + psi) / 2
continue
chi = np.sqrt(y / c2)
t = (chi**3 * c3 + A * np.sqrt(y)) / np.sqrt(mu)
if abs(t - tof) < tol * tof:
break
if t < tof:
psi_low = psi
else:
psi_up = psi
psi = (psi_up + psi_low) / 2
f = 1 - y / r1n
g = A * np.sqrt(y / mu)
gdot = 1 - y / r2n
v1 = (r2 - f * r1) / g
v2 = (gdot * r2 - r1) / g
return v1, v2
Universal Variablesを使うことで、楕円・放物・双曲線軌道のすべてを同じ反復で扱えます。psi は Sundman 変換に対応する変数で、軌道形状によって正・負・ゼロを取ります。c2, c3 はStumpff関数で、軌道方程式を扱いやすくします。
MGA-1DSM評価関数
シーケンスとパラメータベクトルから総 $\Delta v$ を計算する関数を作ります。
import numpy as np
def mga_cost(params, sequence, target_v_inf_max=10.0):
"""MGA-1DSM評価関数 (DSMなし簡易版).
params: [t0_days, T1_days, T2_days, ..., Tn_days, v_inf_launch, alpha, beta]
sequence: ['Earth', 'Venus', 'Venus', 'Earth', 'Jupiter', 'Saturn']
return: 総 Delta v [km/s]
"""
n_legs = len(sequence) - 1
t0 = params[0] * 86400.0 # 出発時刻 [s]
T_legs = np.array(params[1:1+n_legs]) * 86400.0 # 各レグ時間 [s]
v_inf_l = params[1+n_legs] # 打上 v_infty [km/s]
alpha = params[2+n_legs] # 出発方位角
beta = params[3+n_legs] # 出発仰角
# 打上時の地球状態
r0, v_earth = planet_state(sequence[0], t0)
v_inf_vec = v_inf_l * np.array([
np.cos(beta) * np.cos(alpha),
np.cos(beta) * np.sin(alpha),
np.sin(beta)
])
v_sc = v_earth + v_inf_vec
total_dv = max(0.0, v_inf_l - 0.0) # 簡略化: 打上は v_infty そのまま
t_curr = t0
for i in range(n_legs):
t_next = t_curr + T_legs[i]
r_target, v_target = planet_state(sequence[i+1], t_next)
# この区間のLambert問題
try:
v_dep, v_arr = lambert_universal(r0, r_target, T_legs[i], MU_SUN, prograde=True)
except Exception:
return 1e6 # 解なし: 巨大ペナルティ
# 出発側の不整合 (DSMでカバー)
dv_dep = np.linalg.norm(v_dep - v_sc)
if i == 0:
# 1段目は打上 v_infty で吸収可能
v_inf_actual = np.linalg.norm(v_dep - v_earth)
if v_inf_actual > target_v_inf_max:
total_dv += (v_inf_actual - target_v_inf_max) * 10 # ペナルティ
total_dv = v_inf_actual # 打上コスト = v_infty
else:
total_dv += dv_dep # 中間レグはDSMで吸収
# 到着側 v_infty
v_inf_arr_vec = v_arr - v_target
v_inf_arr = np.linalg.norm(v_inf_arr_vec)
if i < n_legs - 1:
# フライバイ通過: v_infty 連続が必要だが簡略化のため
# 次レグ出発状態を到着状態として継続 (フライバイ角は自由とみなす)
r0 = r_target
v_sc = v_target + v_inf_arr_vec # フライバイ後速度 (向きはLambertが要求する向き)
else:
# 最終到着
total_dv += v_inf_arr # 簡略化: 到着 v_infty を全部使う
t_curr = t_next
return total_dv
注意点: 本実装は教育用簡易版で、フライバイ角度制約や $v_\infty$ 連続性を完全には扱っていません。実用では各フライバイで「入射 $v_\infty$ = 離脱 $v_\infty$」を厳密に課し、フライバイ高度制約も入れます。それでも本ソルバはVVEJGAの概形を再現できる程度の精度を持ちます。
差分進化で最適化
scipy.optimize.differential_evolution でCassini風VVEJGAを解きます。
import numpy as np
from scipy.optimize import differential_evolution
# Cassini-like VVEJGA sequence
sequence = ['Earth', 'Venus', 'Venus', 'Earth', 'Jupiter', 'Saturn']
n_legs = len(sequence) - 1
# 探索範囲: t0[日], T1, T2, ..., Tn[日], v_inf_launch[km/s], alpha[rad], beta[rad]
bounds = [
(0, 1000), # t0: 出発(基準日からの日数)
(50, 400), # T1: Earth-Venus
(400, 800), # T2: Venus-Venus (resonant 1:1, 約450日)
(50, 300), # T3: Venus-Earth
(300, 700), # T4: Earth-Jupiter
(1000, 2500), # T5: Jupiter-Saturn
(2.0, 7.0), # v_inf_launch
(-np.pi, np.pi), # alpha
(-0.3, 0.3), # beta
]
result = differential_evolution(
lambda p: mga_cost(p, sequence),
bounds,
popsize=80,
maxiter=600,
tol=1e-7,
seed=42,
mutation=(0.4, 1.2),
recombination=0.9,
workers=-1,
polish=True,
)
print(f"最適 Delta v (total): {result.fun:.3f} km/s")
print(f"出発時刻 (基準日からの日数): {result.x[0]:.1f}")
print(f"各レグ時間 [日]:")
for i, leg in enumerate(zip(sequence[:-1], sequence[1:])):
print(f" {leg[0]:>8} -> {leg[1]:<8}: {result.x[1+i]:.1f}")
print(f"打上 v_infty: {result.x[1+n_legs]:.2f} km/s")
このコードを実行すると、最適化が収束し、典型的に総 $\Delta v$ が 6-9 km/s 程度の解が得られます。実際のCassiniミッションの 4.93 km/s よりはやや大きいですが、これは本実装が DSM の最適配置やフライバイ角度の整合を完全に扱っていないためです。それでも、得られた経路の各レグ飛行時間は実Cassiniの値 ($T_1 \approx 207$日, $T_2 \approx 421$日, $T_3 \approx 56$日, $T_4 \approx 463$日, $T_5 \approx 1306$日) と桁的に一致し、VVEJGAの基本骨格が再現できていることがわかります。
結果の可視化
得られた最適経路を黄道面上に描いてみます。
import numpy as np
import matplotlib.pyplot as plt
def propagate_arc(r0, v0, tof, mu, n_steps=200):
"""簡易ケプラー伝搬: 解析的にr,vを返すのが理想だがここではRunge-Kuttaで近似."""
def deriv(state):
r = state[:3]
v = state[3:]
rn = np.linalg.norm(r)
a = -mu * r / rn**3
return np.concatenate([v, a])
state = np.concatenate([r0, v0])
dt = tof / n_steps
traj = [state[:3].copy()]
for _ in range(n_steps):
k1 = deriv(state)
k2 = deriv(state + 0.5*dt*k1)
k3 = deriv(state + 0.5*dt*k2)
k4 = deriv(state + dt*k3)
state = state + (dt/6.0)*(k1 + 2*k2 + 2*k3 + k4)
traj.append(state[:3].copy())
return np.array(traj)
# 最適解を再評価して各レグの軌道を取得
params = result.x
t0 = params[0] * 86400
T_legs = params[1:1+n_legs] * 86400
fig, ax = plt.subplots(figsize=(10, 10))
# 各惑星の軌道
for name, color in [('Earth','blue'),('Venus','orange'),('Jupiter','brown'),('Saturn','gold')]:
a = PLANETS[name]['a']
th = np.linspace(0, 2*np.pi, 200)
ax.plot(a*np.cos(th)/1e6, a*np.sin(th)/1e6, '--', color=color, alpha=0.4, label=name)
# 各レグの軌跡
t_curr = t0
r0, v_e = planet_state(sequence[0], t_curr)
v_inf_l = params[1+n_legs]
alpha, beta = params[2+n_legs], params[3+n_legs]
v_inf_vec = v_inf_l * np.array([np.cos(beta)*np.cos(alpha), np.cos(beta)*np.sin(alpha), np.sin(beta)])
v_sc = v_e + v_inf_vec
colors = ['red', 'magenta', 'cyan', 'green', 'purple']
for i in range(n_legs):
t_next = t_curr + T_legs[i]
r_target, v_target = planet_state(sequence[i+1], t_next)
v_dep, v_arr = lambert_universal(r0, r_target, T_legs[i], MU_SUN, prograde=True)
traj = propagate_arc(r0, v_dep, T_legs[i], MU_SUN, n_steps=400)
ax.plot(traj[:,0]/1e6, traj[:,1]/1e6, '-', color=colors[i], lw=1.5,
label=f'{sequence[i]}->{sequence[i+1]}')
ax.plot(r_target[0]/1e6, r_target[1]/1e6, 'o', color=colors[i], markersize=10)
r0 = r_target
v_sc = v_target + (v_arr - v_target) # 簡略化
t_curr = t_next
ax.plot(0, 0, 'y*', markersize=20, label='Sun')
ax.set_xlabel('X [million km]')
ax.set_ylabel('Y [million km]')
ax.set_title('Cassini-like VVEJGA trajectory (Earth - Venus - Venus - Earth - Jupiter - Saturn)')
ax.set_aspect('equal')
ax.legend(loc='upper left', fontsize=8)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('vvejga_trajectory.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフは、最適化で得られた経路の概要を示します。地球を出発した探査機は、まず内側に降下して金星と2回フライバイし、その後外側にスイングしながら地球で再加速、最後に木星でさらに加速されて土星に至る——この「複雑だが効率的」な経路が、Cassiniの実飛行軌道に対応します。各色のアークがLambert問題の解で、フライバイ点（惑星位置の丸印）で軌道が「曲げられて」次のレグに連続します。
Tisserandグラフによる確認
最後に、得られた軌道の各レグについて金星・地球・木星に対するTisserandパラメータを計算し、保存性を確認します。
import numpy as np
import matplotlib.pyplot as plt
def tisserand(a, e, i_inc, a_P):
"""Tisserand parameter T = a_P/a + 2*sqrt(a(1-e^2)/a_P) * cos(i)"""
return a_P/a + 2*np.sqrt(a*(1-e**2)/a_P)*np.cos(i_inc)
def orbit_elements_from_rv(r, v, mu):
"""r,vから軌道要素 a,e,i を求める"""
rn = np.linalg.norm(r)
vn = np.linalg.norm(v)
h_vec = np.cross(r, v)
h = np.linalg.norm(h_vec)
energy = 0.5*vn**2 - mu/rn
a = -mu / (2*energy)
e_vec = np.cross(v, h_vec)/mu - r/rn
e = np.linalg.norm(e_vec)
i_inc = np.arccos(h_vec[2]/h)
return a, e, i_inc
# 各レグの始点と終点でTisserandを計算
print("\nTisserand parameter check (Venus-Earth-Jupiter):")
print(f"{'Leg':<25} {'T_V':>8} {'T_E':>8} {'T_J':>8}")
t_curr = t0
r0, v_e = planet_state(sequence[0], t_curr)
v_sc = v_e + v_inf_vec
for i in range(n_legs):
t_next = t_curr + T_legs[i]
r_target, v_target = planet_state(sequence[i+1], t_next)
v_dep, v_arr = lambert_universal(r0, r_target, T_legs[i], MU_SUN, prograde=True)
a, e, inc = orbit_elements_from_rv(r0, v_dep, MU_SUN)
T_V = tisserand(a, e, inc, PLANETS['Venus']['a'])
T_E = tisserand(a, e, inc, PLANETS['Earth']['a'])
T_J = tisserand(a, e, inc, PLANETS['Jupiter']['a'])
print(f"{sequence[i]+' -> '+sequence[i+1]:<25} {T_V:>8.3f} {T_E:>8.3f} {T_J:>8.3f}")
r0 = r_target
t_curr = t_next
出力例を見ると、Venus-Venusレグ（同じ金星間の共鳴フライバイ）で $T_V$ がほぼ変化していないことが確認できます。一方、Earth フライバイの前後では $T_V$ も変化しますが、$T_E$ は保存される傾向が見えます（実装の近似精度の範囲内で）。これがTisserandパラメータの保存則であり、多回フライバイ軌道がこれらの保存量に縛られて構造化されていることを直接示しています。
制約とミッション設計の実際
これまで $\Delta v$ 最小化に焦点を当てましたが、実ミッションでは複数の制約が同時にかかります。
- $\Delta v$ 予算: 推進剤質量に直結。Cassiniの場合、化学推進で500 m/s程度のマージン
- ミッション期間: 探査機の劣化（電源、冷却、姿勢制御燃料）、運用コスト、PI（研究代表者）の任期
- 太陽距離制約: 内側で熱、外側で電源不足。Cassiniは内側でも近日点 $\approx 0.7$ AU を保つ設計
- 打上ウィンドウ: 惑星配置が好条件になる期間。Cassiniは1997年10月の1ヶ月窓
- 科学観測の制約: フライバイ高度、距離、照明条件
- 通信制約: Earth-Sun-Probe角度（合・SPE）、深宇宙ネットワークの空き状況
- 電気推進の場合: 太陽距離による発電量変化、推力プロファイル
これらを同時最適化するには、多目的進化計算 + 制約処理 + 全動力学検証の組み合わせが必要です。実ミッションでは設計に1-3年かかり、初期解 → MOEA → MBH → 全動力学最適化 → 運用解、と段階を踏みます。
まとめ
本記事では、多回フライバイ軌道の設計手法を、CassiniのVVEJGA経路を題材に解説しました。
- エネルギー節約の原理: フライバイは惑星の公転速度ベクトルを「ヒッチハイク」し、太陽中心系での速度を変える。Cassiniは4回のフライバイで約4.2 km/sを節約した。
- パッチドコニックス近似: 軌道を「太陽中心アーク」と「惑星近傍双曲線」に分割し、SOI境界で貼り合わせる。初期設計の標準だが、低エネルギー遷移には限界がある。
- Tisserandパラメータ: $T = a_P/a + 2\sqrt{a(1-e^2)/a_P}\cos i$ は同一惑星でのフライバイ前後で保存。等-$T$ 曲線の交差を使って経路を設計できる。
- MGA定式化: シーケンスを固定し、各アークでLambert問題を解いて連続変数を最適化。$n$ 段なら $3+2n+1$ 次元程度。
- 大域最適化手法: 分岐限定法、差分進化、NSGA-II、深層強化学習。GTOC優勝チームは複数手法のパイプラインを組む。
- 実例の幅: Cassini (VVEJGA), LUCY (7トロヤ群), BepiColombo (低推力MGA), JUICE (衛星35回), Parker (減速フライバイ) — それぞれが異なる戦略を体現する。
- Python実装: 簡易ephemeris + Lambert問題 + 差分進化で、VVEJGA経路の概形を再現できる。
多回フライバイ最適化は、軌道力学・組合せ最適化・機械学習の境界に位置する活発な研究領域です。ロケット推進だけでは決して届かない領域へ、天体の重力を「無料の推進」として組み合わせる — その美しさと工学的価値を、本記事で実感していただければ幸いです。
次のステップとして、以下の記事も参考にしてください。