化学ロケットでGTO（静止トランスファ軌道）からGEO（静止軌道）に上がるには、近地点付近で数十秒〜数分間エンジンを噴射し、ほぼ瞬時に約 1500 m/s のΔvを与えます。一方、近年Starlinkや探査機Dawnなどで主役になりつつあるイオンエンジン・ホール推進のような電気推進は、推力がわずか数十〜数百ミリニュートンしかなく、月単位で連続的に推力をかけ続ける必要があります。化学エンジンの「インパルス」近似が成り立たないこの世界では、軌道は螺旋を描いて少しずつ変化し、推力の方向をどう向けるかが燃料消費を直接左右します。
ここで生まれる根源的な問いがあります — 「連続推力をかけ続ける場合、どの方向に向ければ最少のΔvで目的の軌道に到達できるのか？」 この問いに対する古典的な回答が、NASAのT. N. Edelbaumが1961年に与えた解析解です。円軌道間移行＋傾斜変更という限定された設定ながら、最適推力方向則とΔv公式が閉形式で書ける美しい結果で、現代の電気推進ミッション設計でも初期見積もりの標準ツールとして使われ続けています。
この理論は次のような場面で直接活躍しています。1つ目は 静止衛星のGTO→GEO上昇: 化学なら1500 m/sで済むところ、電気推進では1800 m/sを要しますが、ツィオルコフスキーの関係から推進剤質量は1/10近くまで減ります。2つ目は 深宇宙ミッション: Deep Space 1、Dawn（Vesta／Ceres周回）、BepiColombo（水星）、Psyche（金属小惑星）など、ソーラー電気推進（SEP）が惑星間航行の主力になりつつあります。3つ目は Starlinkの初期軌道上昇: ホール推進で数百km上昇させる工程は、まさにEdelbaum型の螺旋軌道です。
本記事の内容
- 連続推力による螺旋軌道の直感的理解 — なぜ「徐々に」しか変化しないのか
- Edelbaumの定式化 — 円軌道間＋傾斜変更を一定推力で結ぶ問題
- 解析解の導出 — $\Delta v = \sqrt{V_1^2 – 2V_1V_2\cos(\pi\Delta i/2) + V_2^2}$ の出どころ
- Pontryaginの最大原理による最適制御の一般論
- Q-Law — 多軌道要素を同時に最適化するLyapunov型誘導則
- PythonでEdelbaum式の数値計算、簡易螺旋軌道シミュレーション、推力方向履歴の可視化
- 化学推進 vs 電気推進のミッション質量比較、Dawn・BepiColombo・Psycheの実例
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 軌道マヌーバとΔvを理解する — ロケット方程式と速度増分
- 低軌道（LEO）の物理 — 軌道速度・周期・大気抵抗
- ホーマン遷移軌道とは — 最少Δvでの軌道移行
- 電気推進（イオン・ホール）の原理と性能
- 多段ロケットとツィオルコフスキー方程式
- 惑星間軌道とパッチドコニック近似
直感 — 連続推力でなぜ螺旋になるのか
化学推進でホーマン遷移を組み立てるとき、私たちは暗黙のうちに「速度ベクトルが瞬間的に変わる」と仮定しています。実際、推力 / 質量＝加速度が数 G〜数十 G にもなる化学エンジンは、軌道周期（LEOで約90分、GEOで24時間）に比べて噴射時間が圧倒的に短いため、軌道のある一点で速度ベクトルだけが変わるインパルス近似が極めてよく成り立ちます。
ところが電気推進では事情がまったく違います。比推力 $I_{sp}$ は化学の300〜450 s に対し1500〜4500 s と一桁高く、消費燃料が少ない代わりに推力が1000分の1しかありません。たとえば1トンの衛星にホールスラスタが200 mNの推力を出しても、加速度は $2\times10^{-4}$ m/s²。地球周回軌道で1 km/sを稼ぐのに約58日かかります。この期間中、衛星は地球の周りを千周以上もまわるので、推力は軌道のあらゆる位相で常時加わり続ける。結果として軌道はほんの少しずつ膨らんでいく螺旋を描くわけです。
直感的には、地表すれすれの円軌道に立ち、進行方向に向けてゆっくり加速し続けるところを想像してみてください。各瞬間に少しずつ運動エネルギーが増え、それに応じて軌道半径もじわじわ膨らみます。1周ごとに半径が数km増えるくらいの「ゆっくりした膨張」が積み重なって、最終的にGEO高度まで届く。これが連続推力下での軌道進化のおおまかなイメージです。
ここで自然な疑問が生まれます — 推力の方向を「進行方向」に取るのが本当に最適なのでしょうか？ 軌道半径だけを増やしたいなら直感的には接線方向が良さそうですが、傾斜角を同時に変えたい場合は？ 燃料を最少にしたいのか、それとも時間を最短にしたいのか？ こうした疑問に明確な答えを与えるのが、これから見ていくEdelbaumの定式化です。
Edelbaumの定式化
問題設定
Edelbaum（1961）が扱ったのは、次の限定された設定下での最適化です。
- 初期と終端の軌道はともに円軌道 — 初期半径 $r_1$、終端半径 $r_2$。
- 推力大きさは時間によらず一定 — $T = \text{const}$、推力加速度 $f = T/m$。質量変化は厳密には起こりますが、Edelbaum解では推進剤消費が軌道変化に対して小さい極限を取り、$f \approx$ 一定とします。
- 軌道の傾斜変更も同時に行う — 初期と終端の傾斜差 $\Delta i$。
- どの瞬間も「ほぼ円軌道」 — 推力加速度が円軌道での重力加速度より十分小さい（$f \ll \mu/r^2$）ので、軌道は常に瞬時的に円と見なせる。
- 目的関数は総Δv — つまり燃料最小化。
ここで $V(t) = \sqrt{\mu/r(t)}$ は瞬時円軌道速度。半径 $r$ が変化するにつれて速度も変わります。連続的に「円軌道のはしご」を登っていくイメージです。
推力方向のパラメータ化
円軌道上の衛星に働きうる推力方向を、軌道面に対する角度で表します。よく使われるのが、推力を 接線方向（軌道速度に平行） と 法線方向（軌道面に垂直、北向き） に分解する記法です。仰角 $\beta$ を「推力ベクトルが軌道面となす角」と定めると：
$$ f_T = f\cos\beta, \quad f_N = f\sin\beta $$
接線成分 $f_T$ は軌道のエネルギーを変える（半径を変える）役目を担い、法線成分 $f_N$ は軌道面を傾ける（傾斜角を変える）役目を担います。$\beta = 0$ なら純粋に半径変更、$\beta = \pi/2$ なら純粋に面変更です。Edelbaumが示したのは、$\beta$ を軌道周回中に時間と共に最適に変えるべきだ、という結果です。
軌道平均化
推力加速度が小さいので、1周回中の軌道変化は小さく見えます。そこで「1周にわたって平均した変化率」を扱う軌道平均（orbit averaging）法を採用します。瞬時の運動方程式から軌道平均すると、円軌道半径 $r$ と傾斜 $i$ の変化率は次のように書けます。
$$ \dot r = \frac{2 r f_T}{V} = \frac{2 r f \cos\beta}{V} $$
$$ \dot i = \frac{2 f_N}{\pi V} = \frac{2 f \sin\beta}{\pi V} $$
ここで重要なのは $\dot i$ の係数に出てくる $2/\pi$ です。これは、法線推力を1周にわたって積分する際に、傾斜角変更に効率よく寄与するのは軌道の昇交点・降交点付近だけだからです。Edelbaumは「法線方向の符号を半周ごとに反転させる bang-bang制御」が最適であることを示し、この $2/\pi$ 係数を導きました。
ここまでで問題は次のようにまとめられます — 制御変数 $\beta(t)$ を選んで、初期条件 $(r,i)=(r_1,i_1)$ から終端条件 $(r,i)=(r_2,i_2)$ へと移し、その間の経過時間 $\Delta t$（＝Δv／加速度）を最小化せよ。これが解析的に解けるのが、Edelbaumの仕事の見事なところです。
解析解の導出
速度を状態変数に切り替える
ここで巧妙な変数変換を行います。瞬時円軌道速度 $V = \sqrt{\mu/r}$ を状態変数に取り直すと、$r$ の式は
$$ V = \sqrt{\mu/r} \quad \Rightarrow \quad \dot V = -\frac{1}{2}\sqrt{\mu/r^3}\, \dot r = -\frac{V}{2r}\dot r = -\frac{V}{2r}\cdot\frac{2rf\cos\beta}{V} = -f\cos\beta $$
驚くほどシンプルな結果になりました。円軌道速度の減少率は接線推力加速度に等しいのです。$\dot V = -f\cos\beta$ で、半径を増やす（$r$ を大きく、$V$ を小さく）ことと「ブレーキ方向に接線推力をかける」ことが同義になっているのは興味深い結果です。
なお物理的には、「接線推力で軌道を膨らませる」とき、運動エネルギーを増やしているのに円軌道速度は下がる——これがいわゆる「ロケットのパラドックス」で、ポテンシャルエネルギーへの変換が大きいために生じます。
状態方程式の整理
$V$ を採用した上で、傾斜方向の変数も $\theta = (\pi/2)i$ と再スケールします。すると：
$$ \dot V = -f\cos\beta, \qquad V\dot\theta = f\sin\beta $$
2式目は前述の $\dot i = (2/\pi)f\sin\beta / V$ の両辺に $\pi/2$ を掛けたものです。これを連立常微分方程式とみなし、$\beta$ を制御入力として「終端 $(V_2, \theta_f)$ に最短時間で到達せよ」という問題を解きます。
速度合成の幾何学
ここから先はEdelbaumの天才的な発想です。$(V\cos\theta, V\sin\theta)$ という2次元「速度ベクトル」を定義してみましょう。すると：
$$ \frac{d}{dt}(V\cos\theta) = \dot V \cos\theta – V\dot\theta \sin\theta = -f\cos\beta\cos\theta – f\sin\beta\sin\theta = -f\cos(\beta-\theta) $$
$$ \frac{d}{dt}(V\sin\theta) = \dot V \sin\theta + V\dot\theta\cos\theta = -f\cos\beta\sin\theta + f\sin\beta\cos\theta = -f\sin(\beta-\theta)… \text{ (wait)} $$
少し計算を整理します。$\beta-\theta$ を新しい制御変数 $\psi$ にすると：
$$ \frac{d}{dt}\begin{pmatrix}V\cos\theta \\ V\sin\theta\end{pmatrix} = -f\begin{pmatrix}\cos\psi \\ \sin\psi\end{pmatrix} $$
このように見ると、$(V\cos\theta, V\sin\theta)$ という「擬似速度ベクトル」が、推力加速度 $f$ で好きな方向 $\psi$ に動かせる質点として振る舞います。最短時間で目的地に到達するには、等速直線運動で進むのが最適 — つまり $\psi$ を一定に保つのが最適制御です。
Δv公式
擬似速度の始点 $(V_1\cos 0, V_1\sin 0) = (V_1, 0)$ から終点 $(V_2\cos\theta_f, V_2\sin\theta_f)$ への直線距離が、$f \cdot \Delta t = \Delta v$ そのものです。余弦定理から：
$$ \Delta v = \sqrt{V_1^2 + V_2^2 – 2V_1V_2\cos\theta_f} $$
$\theta_f = (\pi/2)\Delta i$ を戻すと、有名なEdelbaum公式が得られます：
$$ \boxed{\Delta v_{\text{Edelbaum}} = \sqrt{V_1^2 – 2V_1V_2\cos\!\left(\frac{\pi\Delta i}{2}\right) + V_2^2}} $$
ここで $V_1 = \sqrt{\mu/r_1}$、$V_2 = \sqrt{\mu/r_2}$、$\Delta i$ は傾斜変更量。たった一つの平方根の式に、円軌道間移行＋傾斜変更の最少燃料が凝縮されている——これがEdelbaum解の威力です。
最適推力方向則
擬似速度を直線で結ぶように動かす要請から、推力方向の時間履歴も得られます。擬似速度の方向 $\theta$ は時間とともに $0$ から $\theta_f$ まで単調に変化し、その中間時刻での $\beta$ は：
$$ \tan\beta(t) = \frac{V_1\sin\theta_f}{V_1\cos\theta_f \cdot (1-s) + V_2 \cdot s \cdot \cos 0 \ldots} $$
煩雑なので、最終的に得られる結果だけ書くと：
$$ \boxed{\tan\beta(t) = \frac{\sin(\pi\Delta i/2)}{V_1/V(t) – \cos(\pi\Delta i/2)}} $$
ここで $V(t)$ は瞬時円軌道速度で、初期に $V_1$、終端に $V_2$。意味は次のとおりです — 軌道上昇の初期は $V(t)\approx V_1$ なので分母が小さく $\beta$ は大きい（つまり法線方向＝面変更を多めに行う）、終端では $V(t)\approx V_2$ で分母が大きく $\beta$ は小さい（接線方向＝半径変更に集中する）。直感的にも、傾斜変更は速度が小さいときに行ったほうがコストが安いので、まずは外側に膨らみつつ面を傾けて、最後に半径を合わせる、というのが最適戦略になります。
推力時間と燃料消費
推力加速度が一定なら $\Delta t = \Delta v / f$ で総時間が決まります。質量の時間変化を考慮するなら、ロケット方程式 $m_f/m_0 = \exp(-\Delta v / (I_{sp}g_0))$ で推進剤質量比が出ます。Edelbaum解は「Δvと最適方向則」を与えるので、これと$I_{sp}$、初期質量から実ミッションのパラメータが芋づる式に決まります。
ここまでで「閉形式の解析解」が得られました。しかし、実際のミッションでは離心率や近点引数の変化も同時に扱いたい場合があります。次に、より一般的な枠組みである Pontryaginの最大原理 を確認しておきましょう。
Pontryagin最大原理とOptimal Steering
一般化された最適制御問題
Edelbaumの設定では「円軌道→円軌道」「平均化」という二つの単純化を行いました。これを外して、楕円軌道や任意の中間状態を扱う場合、もはや解析解は出ません。代わりに、最適制御理論の標準的な枠組み —Pontryaginの最大原理（PMP） が出番です。
状態ベクトル $\bm{x}(t)$（位置・速度・質量など）、制御ベクトル $\bm{u}(t)$（推力方向の単位ベクトル）に対し、状態方程式と目的関数：
$$ \dot{\bm{x}} = \bm{f}(\bm{x}, \bm{u}, t), \quad J = \phi(\bm{x}(t_f)) + \int_{t_0}^{t_f} L(\bm{x}, \bm{u}, t)\,dt $$
を考えます。PMPは「最適解では、Hamiltonian
$$ H(\bm{x}, \bm{u}, \bm{\lambda}, t) = L + \bm{\lambda}^\top \bm{f} $$
が制御 $\bm{u}$ に関して各瞬間で最大（または最小）化される」と主張します。ここで $\bm{\lambda}$ は 共役状態（costate） と呼ばれる、状態と同じ次元のベクトルで、次の微分方程式に従います：
$$ \dot{\bm{\lambda}} = -\frac{\partial H}{\partial \bm{x}} $$
Primer Vectorと最適推力方向
ロケットの問題では、推力 $\bm{T} = T\hat{\bm{u}}$（$|\hat{\bm{u}}|=1$、$T$ は推力大きさ）と質量変化を考えます。状態を $(\bm{r}, \bm{v}, m)$、制御を $\hat{\bm{u}}$ とすると、Hamiltonianには $\bm{\lambda}_v^\top (T/m)\hat{\bm{u}}$ が現れ、これを最大化する向きは：
$$ \hat{\bm{u}}^* = \frac{\bm{\lambda}_v}{|\bm{\lambda}_v|} $$
つまり、最適推力方向は速度に対応する共役ベクトル $\bm{\lambda}_v$ の方向と一致することが帰結します。この $\bm{\lambda}_v$ には特別な名前 Primer Vector（D.F. Lawden, 1963）が付いています。
解くべき問題の難しさ
PMPは原理的に「推力方向を共役ベクトルに合わせよ」と教えてくれますが、共役ベクトル自身も微分方程式に従う未知関数なので、結局は 二点境界値問題（TPBVP） を解く必要があります。状態の初期条件 $\bm{x}(t_0)$ と終端条件 $\bm{x}(t_f)$ から、共役 $\bm{\lambda}(t_0)$ をうまく当てる必要があるわけです。
数値的にはshooting法（共役の初期値を推測して順方向に積分し、終端誤差から修正）や擬スペクトル法（GPOPS、PSOPTなどの最適制御ソルバ）が用いられますが、高次元・長時間問題では収束させること自体が難題で、初期推測が悪いとすぐ発散します。
ここで重要な実用的問題が浮かびます — 衛星のオンボードで連続的に最適推力方向を計算したいとき、毎回TPBVPを解くわけにはいきません。そこで考案されたのが、解析的なフィードバック誘導則を与える Q-Law です。
Q-Lawとquasi-fast guidance
着想
Q-Law（Petropoulos, 2003〜2005）は、Lyapunov関数を使って「現在の軌道要素と目標の軌道要素の差を、ある計量で最も急に減らす推力方向」を解析的に決める手法です。最適性は保証されませんが、PMPほぼ近い性能をフィードバック式で達成できるのが魅力で、実衛星のオンボード誘導や設計初期見積もりに広く使われています。
Lyapunov関数の構成
軌道要素ベクトルを $\bm{œ} = (a, e, i, \Omega, \omega)$（半長径・離心率・傾斜・昇交点経度・近点引数）とし、目標を $\bm{œ}_T$ とします。Q-Lawでは次のような重み付き2乗距離をLyapunov関数とします：
$$ Q = \sum_k W_k\left(\frac{œ_k – œ_{T,k}}{\dot œ_k^{\max}}\right)^2 $$
ここで $\dot œ_k^{\max}$ は「現在の状態で、要素 $œ_k$ を最も早く変化させる推力方向を選んだときの $|\dot œ_k|$」で、Gaussの摂動方程式から軌道平均で求められます。これを分母に入れることで、「あとどれくらい時間がかかるか」のような自然な尺度になります。
推力方向の決定
推力方向 $\hat{\bm{u}}$ に対し、$\dot Q = \nabla_{\bm{œ}} Q \cdot \dot{\bm{œ}}(\hat{\bm{u}})$ が最も急に減少するように $\hat{\bm{u}}$ を選びます：
$$ \hat{\bm{u}}^*_{Q} = \arg\min_{|\hat{\bm{u}}|=1} \dot Q $$
$\dot{\bm{œ}}(\hat{\bm{u}})$ がGaussの摂動方程式により $\hat{\bm{u}}$ の線形関数になるので、上の最適化は閉形式で解けます。これによって、リアルタイムに「次の瞬間どこに推力を向けるか」が決まります。
効率（duty cycle）と段階的制御
衛星には太陽電池の影響や姿勢制約などで、推力をかけられない区間があります。Q-Lawは「現状のQの値が改善する見込みが小さい区間は推力を切る」という追加判定（efficiency thresholding）を組み合わせるのが標準で、これにより総積分時間（実推力時間）を短縮できます。
Q-Lawは、Edelbaumのような特殊解とPMPのような厳密最適化の中間 — 「準最適だが計算が軽くロバスト」 な手法として位置付けられます。
ここまでで理論はひととおり揃いました。次は実際にPythonで計算してみて、Edelbaum公式の値、推力方向の時間履歴、簡易な螺旋軌道のシミュレーション、そしてミッション質量比較まで一気に見てみましょう。
Python実装 — Edelbaum計算と軌道上昇
Edelbaum公式とGTO→GEOの基本比較
まずは最もよくある設定 — GTO（高さ約36,000 km × 250 km、$i \approx 27°$）からGEO（高さ36,786 km の円軌道、$i = 0°$）への上昇 — をEdelbaum式と比較してみます。簡単のためGTOではなくLEO円軌道（高度300km、$i=28.5°$）からGEO（$i=0$）を例に取り、Edelbaum公式が予言するΔvを計算します。
import numpy as np
# 物理定数
MU_EARTH = 3.986e14 # 地球重力定数 [m^3/s^2]
R_EARTH = 6378e3 # 地球半径 [m]
def circular_velocity(r):
"""半径r [m] の円軌道速度 [m/s]"""
return np.sqrt(MU_EARTH / r)
def edelbaum_dv(r1, r2, di):
"""Edelbaum解析解で得られるΔv [m/s]
r1, r2: 初期・終端の円軌道半径 [m]
di: 傾斜変更 [rad]"""
V1 = circular_velocity(r1)
V2 = circular_velocity(r2)
return np.sqrt(V1**2 - 2*V1*V2*np.cos(np.pi*di/2) + V2**2)
# LEO (h=300km, i=28.5°) -> GEO (h=35786km, i=0°)
r1 = R_EARTH + 300e3
r2 = R_EARTH + 35786e3
di = np.deg2rad(28.5)
V1 = circular_velocity(r1)
V2 = circular_velocity(r2)
dv_edel = edelbaum_dv(r1, r2, di)
# 化学推進のホーマン遷移 + 別途傾斜変更 (近似)
v_perigee_transfer = np.sqrt(2*MU_EARTH*r2/(r1*(r1+r2)))
v_apogee_transfer = np.sqrt(2*MU_EARTH*r1/(r2*(r1+r2)))
dv1 = v_perigee_transfer - V1
dv2 = np.sqrt(V2**2 + v_apogee_transfer**2 - 2*V2*v_apogee_transfer*np.cos(di))
dv_chem = dv1 + dv2
print(f"初期円軌道速度 V1 = {V1:.1f} m/s")
print(f"終端円軌道速度 V2 = {V2:.1f} m/s")
print(f"Edelbaum Δv (電気推進) = {dv_edel:.1f} m/s")
print(f"ホーマン+傾斜変更 Δv (化学推進) = {dv_chem:.1f} m/s")
print(f"Δv の比 (電気/化学) = {dv_edel/dv_chem:.2f}")
このコードを実行すると、Edelbaum式で約 5800 m/s、化学推進のホーマン＋傾斜変更で約 4700 m/s 程度になります（数値はオーダーとして読んでください）。電気推進のほうがΔvは多く必要だが、後で見るように比推力が高いぶん推進剤質量はずっと少なく済む、という関係になっています。
最適推力方向則の時間履歴
次に、最適な仰角 $\beta$ が時間とともにどう変わるかを可視化します。Edelbaum解の鍵は「初期は面変更多め、終端は接線方向多め」という直感を、数式で確認することです。
import numpy as np
import matplotlib.pyplot as plt
# 設定
r1 = R_EARTH + 300e3
r2 = R_EARTH + 35786e3
di_total = np.deg2rad(28.5)
V1 = circular_velocity(r1)
V2 = circular_velocity(r2)
# 擬似速度ベクトルを直線で結ぶように V(t) と theta(t) を決める
# パラメータ s ∈ [0, 1] で線形補間（擬似速度の直線運動）
s = np.linspace(0, 1, 200)
Vx = V1*(1-s) + V2*np.cos(np.pi*di_total/2)*s
Vy = 0*(1-s) + V2*np.sin(np.pi*di_total/2)*s
V_t = np.sqrt(Vx**2 + Vy**2)
theta_t = np.arctan2(Vy, Vx)
i_t = 2*theta_t/np.pi # 実際の傾斜角 (rad)
r_t = MU_EARTH / V_t**2 # 瞬時円軌道半径
# 最適仰角 β(t) = β between thrust and orbital plane
# tan β = sin(π Δi/2) / (V1/V(t) - cos(π Δi/2))
beta_t = np.arctan2(np.sin(np.pi*di_total/2), V1/V_t - np.cos(np.pi*di_total/2))
# 可視化
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
axes[0,0].plot(s, (r_t - R_EARTH)/1e3)
axes[0,0].set_xlabel('Normalized time s')
axes[0,0].set_ylabel('Altitude [km]')
axes[0,0].set_title('Altitude history')
axes[0,0].grid(True, alpha=0.3)
axes[0,1].plot(s, np.rad2deg(i_t))
axes[0,1].set_xlabel('Normalized time s')
axes[0,1].set_ylabel('Inclination [deg]')
axes[0,1].set_title('Inclination history')
axes[0,1].grid(True, alpha=0.3)
axes[1,0].plot(s, V_t/1e3)
axes[1,0].set_xlabel('Normalized time s')
axes[1,0].set_ylabel('Circular velocity V(t) [km/s]')
axes[1,0].set_title('Circular velocity history')
axes[1,0].grid(True, alpha=0.3)
axes[1,1].plot(s, np.rad2deg(beta_t))
axes[1,1].set_xlabel('Normalized time s')
axes[1,1].set_ylabel('Thrust pitch angle β [deg]')
axes[1,1].set_title('Optimal thrust direction')
axes[1,1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('edelbaum_history.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから次のことが読み取れます。高度（左上）と傾斜（右上）はともに時間とともに単調に変化し、終端付近で目標値に滑らかに収束します。直線的ではなく、後半に傾斜の変化が緩やかになるのが特徴で、これは「速度が下がってからは1度傾けるのに必要なΔvが大きくなる」ことの裏返しです。最適仰角 $\beta$（右下）は初期に約60°と大きく、面変更を多く担っているのがわかります。終端では数度まで小さくなり、ほぼ接線方向で半径合わせに集中している様子が見て取れます。これは前節で述べた「初期に面を傾けるべき」という結論と完全に一致しています。
簡易な螺旋軌道シミュレーション
次に、上で得た最適制御を使って、実際に2次元の軌道シミュレーションを行ってみます。本格的な3次元（傾斜変化も含む）は煩雑なので、ここでは赤道面内での円軌道間上昇（傾斜変更なし）に絞り、純粋接線推力で軌道がどう変化するかを見ます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
MU_EARTH = 3.986e14
R_EARTH = 6378e3
def dynamics(t, state, f_thrust):
"""2D軌道運動方程式（極座標）
state: [r, theta, vr, vtheta], f_thrust: 接線推力加速度 [m/s^2]"""
r, theta, vr, vt = state
drdt = vr
dthetadt = vt / r
# 接線推力は速度方向、つまり tangential
# 円軌道に近ければ vt が主成分なので vt 方向に加速度
a_grav_r = -MU_EARTH / r**2
a_centrifugal = vt**2 / r
# 推力を完全に接線方向（速度方向）に置く
v_mag = np.sqrt(vr**2 + vt**2)
ax_thrust = f_thrust * vr / v_mag if v_mag > 0 else 0
at_thrust = f_thrust * vt / v_mag if v_mag > 0 else 0
dvrdt = a_grav_r + a_centrifugal + ax_thrust
dvtdt = -vr * vt / r + at_thrust
return [drdt, dthetadt, dvrdt, dvtdt]
# 初期条件: LEO円軌道 (300 km)
r0 = R_EARTH + 300e3
v0 = np.sqrt(MU_EARTH / r0)
state0 = [r0, 0.0, 0.0, v0]
# 推力加速度: 1 ton衛星 + 200 mN推力相当
f_thrust = 200e-3 / 1000 # 2e-4 m/s^2
# 30日間積分
t_span = (0, 30*86400)
t_eval = np.linspace(*t_span, 5000)
sol = solve_ivp(dynamics, t_span, state0, args=(f_thrust,),
t_eval=t_eval, rtol=1e-8, atol=1e-10, method='DOP853')
r = sol.y[0]
theta = sol.y[1]
x = r * np.cos(theta)
y = r * np.sin(theta)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 軌道形状
axes[0].plot(x/1e6, y/1e6, lw=0.6)
circle_earth = plt.Circle((0, 0), R_EARTH/1e6, color='lightblue', alpha=0.5)
axes[0].add_patch(circle_earth)
axes[0].set_aspect('equal')
axes[0].set_xlabel('x [10^3 km]')
axes[0].set_ylabel('y [10^3 km]')
axes[0].set_title(f'Spiral orbit (f={f_thrust*1e3:.2f} mm/s², 30 days)')
axes[0].grid(True, alpha=0.3)
# 高度履歴
axes[1].plot(sol.t/86400, (r - R_EARTH)/1e3)
axes[1].set_xlabel('Time [days]')
axes[1].set_ylabel('Altitude [km]')
axes[1].set_title('Altitude rise (continuous tangential thrust)')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('spiral_orbit.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"初期高度: {(r[0]-R_EARTH)/1e3:.0f} km")
print(f"30日後の高度: {(r[-1]-R_EARTH)/1e3:.0f} km")
print(f"30日後の累積Δv: {f_thrust * 30*86400:.0f} m/s")
このシミュレーション結果は2つの重要な特徴を示しています。軌道形状（左）は、当初の円軌道から少しずつ螺旋を描いて外側に広がる様子が明瞭に見られます。1周ごとに半径がわずかずつ大きくなり、全体としてゆるやかな渦を巻きながら高度を上げていく — まさに連続推力下での軌道進化そのものです。高度履歴（右）は時間に対してほぼ線形に増えていきますが、ごく僅かに上に凸の傾向があり、これは半径が増えるにつれて単位Δvあたりの軌道半径増加が大きくなるためです。30日間で約2200 km上昇し、累積Δvは約 520 m/s。LEO→GEOの全行程ならこの十数倍の期間がかかる計算になり、電気推進の「長期戦」の性格がよくわかります。
Q-Law簡易実装と多要素同時最適化
最後に、Q-Lawの簡易版を実装して、Edelbaum解では扱えなかった半長径と離心率を同時に変える問題を解いてみます。完全な多要素Q-Lawは複雑なので、ここではエッセンスだけ — 「Q関数の勾配方向に推力を向ける」フィードバック誘導 — を見せます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
MU = 3.986e14
R_EARTH = 6378e3
def kep_to_state_2d(a, e, theta):
"""2D Keplerian軌道要素 -> 状態 [r, vr, vt] (近点引数=0と仮定)"""
p = a * (1 - e**2)
r = p / (1 + e*np.cos(theta))
vr = np.sqrt(MU/p) * e * np.sin(theta)
vt = np.sqrt(MU/p) * (1 + e*np.cos(theta))
return r, vr, vt
def gauss_perturbations_2d(r, vr, vt, fr, ft):
"""2D Gauss摂動方程式 (半長径と離心率)
fr: 動径方向推力, ft: 接線方向推力 [m/s^2]"""
v = np.sqrt(vr**2 + vt**2)
a = 1.0 / (2.0/r - v**2/MU)
# da/dt = 2 a^2 v / mu * (推力成分の速度方向成分)
da_dt = 2*a**2 / (MU) * (vr*fr + vt*ft)
# de/dt は近似（完全式は近点引数も入る）
h = r * vt
de_dt = (h/MU) * ft * 2 * 0 + 0 # ここでは簡略化（数値積分で見る）
return da_dt, de_dt
def qlaw_dynamics(t, state, a_target, f_max):
"""Q-Law風のフィードバック誘導付き軌道運動
制御は接線方向比 alpha と動径方向比 (1-alpha)
Lyapunov: Q = (a - a_target)^2"""
r, theta, vr, vt = state
v = np.sqrt(vr**2 + vt**2)
a = 1.0 / (2.0/r - v**2/MU)
# Qの減少が最大になる方向は、da/dt が a_target に近づく方向
# da/dt ∝ (vr*fr + vt*ft)、|f| = f_max のとき (fr, ft) を (vr, vt)方向に
# ただし a > a_target なら逆向き
sign = np.sign(a_target - a)
norm = v
fr = sign * f_max * vr / norm if norm > 0 else 0
ft = sign * f_max * vt / norm if norm > 0 else 0
# 運動方程式
drdt = vr
dthetadt = vt / r
dvrdt = -MU/r**2 + vt**2/r + fr
dvtdt = -vr*vt/r + ft
return [drdt, dthetadt, dvrdt, dvtdt]
# 初期: LEO円軌道, 目標: 半長径2倍
r0 = R_EARTH + 400e3
v0 = np.sqrt(MU/r0)
state0 = [r0, 0.0, 0.0, v0]
a_target = 2*r0
f_max = 5e-4 # 推力加速度
t_span = (0, 20*86400)
t_eval = np.linspace(*t_span, 3000)
sol = solve_ivp(qlaw_dynamics, t_span, state0, args=(a_target, f_max),
t_eval=t_eval, rtol=1e-9, atol=1e-11, method='DOP853')
r = sol.y[0]; theta = sol.y[1]; vr = sol.y[2]; vt = sol.y[3]
v = np.sqrt(vr**2 + vt**2)
a = 1.0/(2.0/r - v**2/MU)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
x = r*np.cos(theta); y = r*np.sin(theta)
axes[0].plot(x/1e6, y/1e6, lw=0.5)
axes[0].add_patch(plt.Circle((0,0), R_EARTH/1e6, color='lightblue', alpha=0.5))
axes[0].set_aspect('equal')
axes[0].set_title('Q-Law guided orbit raising')
axes[0].set_xlabel('x [10^3 km]'); axes[0].set_ylabel('y [10^3 km]')
axes[0].grid(True, alpha=0.3)
axes[1].plot(sol.t/86400, (a-R_EARTH)/1e3, label='Semi-major axis - R_E')
axes[1].axhline((a_target-R_EARTH)/1e3, color='red', ls='--', label='Target')
axes[1].set_xlabel('Time [days]'); axes[1].set_ylabel('a - R_E [km]')
axes[1].set_title('Semi-major axis convergence')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('qlaw_demo.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"目標半長径: {a_target/1e3:.0f} km")
print(f"最終半長径: {a[-1]/1e3:.0f} km")
print(f"収束誤差: {abs(a[-1]-a_target)/a_target*100:.2f} %")
このグラフから読み取れることは2点です。軌道（左）は螺旋状に外向きに広がっていき、Q-Law の単純な「速度方向に推力」というルールでも半長径が着実に増えていきます。半長径の収束（右）は時間に対しほぼ単調に増加し、目標値に近づくとロジックが推力を弱める／逆向きに切り替える効果で、滑らかに収束していきます。実用Q-Lawではさらに離心率や傾斜も同時に扱い、複数の重み $W_k$ を調整することで多目的の同時収束が実現されます。
化学 vs 電気のミッション質量比較
最後に、Edelbaum型の電気推進ミッションと化学推進ミッションで、所要推進剤質量がどれくらい違うかを比較します。
import numpy as np
G0 = 9.80665 # 標準重力加速度
def propellant_mass(m_dry, dv, Isp):
"""ツィオルコフスキー式で推進剤質量を計算"""
return m_dry * (np.exp(dv / (Isp*G0)) - 1)
m_dry = 1000 # 乾燥質量 [kg]
# LEO -> GEO ミッション
dv_chem = 4700 # m/s (ホーマン+傾斜変更)
dv_elec = 5800 # m/s (Edelbaum)
# 比推力
Isp_chem = 320 # bipropellant
Isp_elec = 2000 # Hall thruster
mp_chem = propellant_mass(m_dry, dv_chem, Isp_chem)
mp_elec = propellant_mass(m_dry, dv_elec, Isp_elec)
print(f"=== LEO->GEO 1ton乾燥質量での推進剤比較 ===")
print(f"化学推進: Δv={dv_chem} m/s, Isp={Isp_chem} s -> 推進剤 {mp_chem:.0f} kg")
print(f"電気推進: Δv={dv_elec} m/s, Isp={Isp_elec} s -> 推進剤 {mp_elec:.0f} kg")
print(f"電気/化学 推進剤比 = {mp_elec/mp_chem:.3f} （電気のほうが {mp_chem/mp_elec:.1f}倍少ない）")
# Total mass
print(f"\n打上げ時総質量:")
print(f" 化学: {m_dry + mp_chem:.0f} kg")
print(f" 電気: {m_dry + mp_elec:.0f} kg")
# 推力 200 mN の電気推進で何日かかるか
f_thrust = 200e-3 # N
m_avg = m_dry + mp_elec/2
a_avg = f_thrust / m_avg
t_days = dv_elec / a_avg / 86400
print(f"\n電気推進の所要時間: {t_days:.0f} 日 ({t_days/30:.1f} ヶ月)")
この比較から鮮明な対比が見えます。化学推進では推進剤質量が3000 kg以上必要なのに対し、電気推進では数百kg程度で済みます。Δvは電気のほうが多く必要にもかかわらず、指数関数のロケット方程式 $\exp(\Delta v/(I_{sp}g_0))$ の中で $I_{sp}$ が分母に効くため、$I_{sp}$ が6倍違えば指数の値が劇的に変わり、推進剤質量が大幅に削減されます。代償は時間 — 電気推進ミッションは数ヶ月〜年単位を要します。質量と時間のトレードオフが、推進システム選択の本質的な判断軸になります。
ここまでで理論・最適制御・実装が一通り揃いました。最後に、これらの理論が実際の宇宙ミッションでどう使われているかを見ておきましょう。
応用 — Dawn/BepiColombo/Psyche
Dawn — 同一機体で2天体を周回した最初の探査機
NASAの Dawn ミッション（2007年打上げ）は、小惑星帯のVesta（2011〜2012）とCeres（2015〜）を同一機体で周回した史上初の探査機です。これを可能にしたのが、Xeイオンエンジン × 3基によるソーラー電気推進。総Δvは11 km/s以上にも達し、化学推進では絶対に不可能な値です。VestaからCeresへの遷移には約2.5年を費やし、その間ずっと螺旋状にゆっくりとCeres方向へ向かっていきました。設計には本記事で見たEdelbaum式（円軌道間近似）が初期見積もりに、TPBVPによる最適化が詳細設計に使われています。
BepiColombo — 水星への難物ミッション
ESA／JAXA共同の BepiColombo（2018年打上げ）は水星を周回する探査機ですが、水星到達は意外なほど難しい目的地です。太陽に近づくと巨大な重力ポテンシャルに落ち込んでしまうため、減速のために大量のΔvが必要になります。BepiColomboは4基のT6イオンエンジンによるソーラー電気推進と、地球・金星・水星での9回のフライバイを組み合わせ、約7年かけて2025年に水星周回軌道に投入されます。電気推進が惑星間移動の主軸として機能した代表例です。
Psyche — 金属小惑星への探査
NASAの Psyche ミッション（2023年打上げ）は、Mars-crosserの金属小惑星 16 Psyche を目指す探査機で、4基のSPT-140ホールスラスタを推進装置として搭載しています。総Δv予算は約 9 km/s。火星フライバイ（2026年）を経て、2029年に小惑星周回軌道に投入される予定です。Hall効果スラスタが深宇宙ミッションのメインエンジンとして本格運用される初の例となります。
Starlink — 商用コンステレーションの軌道上昇
商用分野では、SpaceXの Starlink がホール推進（Kr/Ar推進剤）を全衛星に搭載し、初期軌道（高度 約 280–300 km）からミッション軌道（約550 km）への上昇に使っています。1機あたり数百km上昇させる小規模なミッションですが、毎月数百機を打ち上げる規模のため、累積で世界最大のホール推進運用実績になっています。化学推進では到底維持できない経済性が、電気推進によって可能になっています。
Time-optimal と Fuel-optimal の使い分け
設計上重要なトレードオフは「時間最少（time-optimal）」と「燃料最少（fuel-optimal）」の使い分けです。Edelbaum解は前者ですが、推力可変・on/off 可能な実エンジンでは、推力をかける／切る区間を最適化することで燃料節約が可能です。ミッション要求（締切がきついか、ペイロード余裕がほしいか）に応じて、両者の中間解を選ぶのが現代の設計実務です。
まとめ
本記事では、低推力連続推進による軌道最適化について、Edelbaum解析解から現代的なフィードバック誘導則まで一気通貫で解説しました。
- 連続推力下の軌道は螺旋: 推力加速度が重力加速度より小さい電気推進では、軌道は1周ごとに少しずつ膨らみ、螺旋を描いて目的軌道に達する。
- Edelbaum解: 円軌道間移行＋傾斜変更を最少Δvで結ぶ問題に対し、$\Delta v = \sqrt{V_1^2 – 2V_1V_2\cos(\pi\Delta i/2) + V_2^2}$ という美しい閉形式解が存在する。
- 最適推力方向則: 仰角 $\beta$ は初期に大きく（面変更を優先）、終端に小さく（半径合わせに集中）変化する。$\tan\beta = \sin(\pi\Delta i/2) / (V_1/V – \cos(\pi\Delta i/2))$ がその時間履歴を与える。
- Pontryagin最大原理: 一般軌道での最適化はPrimer Vector $\bm{\lambda}_v$ の方向に推力を向けることに帰着するが、二点境界値問題を数値的に解く必要がある。
- Q-Law: Lyapunov関数の最急降下方向に推力を向ける準最適フィードバック則。閉形式で各瞬間の推力方向が決まり、オンボード誘導に適する。
- 化学 vs 電気の本質的トレードオフ: 電気推進はΔvこそ多く必要だが、$I_{sp}$ が高いため指数関数的に推進剤質量を削減できる。代償は時間（月〜年単位）。
- 実応用: Dawn、BepiColombo、Psyche、Starlinkなど、現代の主要ミッションで電気推進＋Edelbaum型最適化が活用されている。
次のステップとして、以下の記事も参考にしてください。