「月面に着陸する」と言葉にすると一見シンプルですが、物理を立てるとたちまち恐ろしい問題になります。月には大気がないのでパラシュートは使えません。重力が地球の約 1/6 とはいえ、軌道速度から静止まで減速するには必ずエンジンを噴射し続けなければならない。しかも噴射すれば推進剤を消費するため、機体の質量がどんどん軽くなっていく。さらに着陸点は岩石や斜面を避けて精密に狙う必要があり、ヘッドライトのような単純な減速では足りない。Apollo の月着陸船は人間のパイロットが残り燃料を見ながら最終調整しましたが、Apollo 11 では着陸燃料の残量がわずか 25 秒分しかなかったと記録されています。一方、2024 年に成功した日本の SLIM は完全自律で目標から 100 m 以内の精度に着陸しました。さらに 2010 年代以降、SpaceX Falcon 9 の一段は同じ「動力降下」を地球上で実現し、ロケットを着陸させて再利用するまでに至っています。
これらすべてに共通する核心が 動力降下ガイダンス(powered descent guidance, PDG) と、その背後にある 最適制御 です。本記事では、月や火星のような大気の薄い天体に「逆噴射しながら着陸する」問題を最適制御問題として定式化し、Apollo 時代の二次多項式ガイダンス E-Guidance から、Acikmese と Ploen が 2007 年に示した Lossless Convexification (LCvx)、そして現代の凸最適化ベースの SOCP リアルタイム解法へと連なる流れを追います。
この主題が活きる場所は実用上もアカデミックにも幅広く、たとえば、
- 月面・惑星探査機の着陸 GNC: Apollo, SLIM, Chang’e, Chandrayaan, Blue Ghost、これからの Artemis 計画の有人着陸船。
- 再使用ロケットのブースター回収: Falcon 9 first stage の “boostback → entry burn → landing burn”、New Glenn、Starship のチョップスティック着陸。
本記事の内容
- 動力降下の運動方程式とロケット方程式(質量変化)の連結
- 二次元・三次元の状態空間モデル(縦/横ターゲティング)と境界条件
- Apollo E-Guidance: 多項式(時間に関する多項式)を最終条件から逆算する閉形式
- 重力ターン軌道とその限界
- 制約付き最適制御: 推力範囲・推力ベクトル角・glide slope
- Acikmese と Ploen の Lossless Convexification と Apollonius Circle
- SOCP として書き直し、CVXPY でリアルタイム解法
- Python による 1D 降下シミュレーション、E-Guidance 多項式係数の閉形式計算、月面 30 km から 1 km 地点までの最適軌道可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
動力降下とは何か — 重力に抗いながら止まる
地球から月の表面まで、ホーマン遷移軌道で運んだ着陸機はたいてい月周回軌道(高度 100 km 程度の円軌道)に入ります。そこから着陸面に向かうために短い接線方向逆噴射で 降下軌道(descent orbit) に入る — このとき近月点が高度 10 〜 15 km まで下がる楕円軌道になります。Apollo は近月点高度 15 km、SLIM は 15 km 程度、Chandrayaan-3 も同程度でした。
この近月点で再びエンジンに火を入れ、軌道速度約 1.7 km/s をほぼゼロまで落としつつ、横方向 0 に追い込み、最終的に垂直に降下して着地する — この一連を 動力降下(powered descent) と呼びます。月面 1G ぶん(月の重力加速度 $g_M \approx 1.62~\mathrm{m/s^2}$)を打ち消した上で、さらに減速するための推力が必要で、噴射時間は典型的に 10〜15 分。
直感的なたとえとして、満員のエレベーターが落下している最中に、底にロケットエンジンが付いていて、それを上向きに連続噴射しながら停止位置にぴたっと止める — そんな曲芸を、自動運転で、燃料制限の中で、岩を避けながらやる、というイメージです。ここで重要なのは次の二つの事実です。
第一に、月や火星には大気がほとんどない ため、地球の旅客機のように揚力で減速したり、パラシュートを開いて空気抵抗だけで降りたりはできません(火星は薄い大気があるためEDLでパラシュートを併用しますが、それでも最後の最後は逆噴射が必須です)。減速エネルギーは全てロケットエンジンの噴射で稼ぐしかありません。
第二に、推進剤を消費するので機体の慣性質量が常に変化する。打ち上げのときは ΔV だけ計算すれば良かったロケット方程式の世界が、今度は質量を状態変数として持つ動的な制御問題へと拡張されるのです。
つまり、動力降下は「重力に抗いながら、自分の質量変化を考慮しつつ、燃料制限と推力制限の中で、目標地点に精密に止まる」という複合制約付きの最適制御問題です。次節ではこの問題を、状態空間方程式として書き下すところから始めます。
運動方程式とロケット方程式の連結
並進運動と質量の連立 ODE
地表近くで重力場を一様な定数ベクトル $\bm{g}$ で近似します(月では $g_M \approx 1.62~\mathrm{m/s^2}$、火星では $g_{\mathrm{Mars}} \approx 3.71~\mathrm{m/s^2}$、地球では $g_E \approx 9.81~\mathrm{m/s^2}$)。着陸機の位置を $\bm{r}(t) \in \mathbb{R}^3$、速度を $\bm{v}(t) \in \mathbb{R}^3$、質量を $m(t) > 0$ とします。エンジンが生成する推力ベクトル $\bm{T}(t)$ は機体に作用する力で、その大きさを $\Gamma(t) = \|\bm{T}(t)\|$、向きはエンジンのジンバル角で決めます。
ニュートンの第二法則と質量保存(噴射した推進剤の質量だけ機体は軽くなる)を書き下すと、
$$ \begin{align} \dot{\bm{r}}(t) &= \bm{v}(t) \\ \dot{\bm{v}}(t) &= \frac{\bm{T}(t)}{m(t)} + \bm{g} \\ \dot{m}(t) &= -\frac{\Gamma(t)}{I_{\mathrm{sp}} g_0} \end{align} $$
ここで $I_{\mathrm{sp}}$ は比推力(秒)、$g_0 = 9.81~\mathrm{m/s^2}$ は地球標準重力(比推力の定義に出てくる規格化定数で、月着陸でも地球の $g_0$ を使うことに注意)です。エンジンの排気速度を $v_e = I_{\mathrm{sp}} g_0$ と書けば、質量流量は $\dot m = -\Gamma / v_e$ とも書けます。
この三本の方程式が 動力降下の完全な力学 です。注目すべきは、第三式が静的なロケット方程式 $\Delta v = v_e \ln(m_0/m_f)$ の動的な連続時間版になっている点です。実際、無重力下で推力方向を変えない単純な噴射を仮定すれば、$\dot v = \Gamma/m$ と $\dot m = -\Gamma/v_e$ から $dv/dm = -v_e/m$ となり、積分するとツィオルコフスキー方程式が再現されます。動力降下の問題は、ロケット方程式の中に重力項と幾何拘束を埋め込んだ拡張版だと見るとわかりやすいでしょう。
1D 縦方向への簡略化(最初のおもちゃモデル)
まず一番単純な「真上から真下」へのモデルを考えます。位置 $h(t)$(高度)、速度 $v(t)$(上向きを正)、質量 $m(t)$、推力大きさ $T(t)$ で、
$$ \dot h = v, \quad \dot v = \frac{T}{m} – g_M, \quad \dot m = -\frac{T}{v_e} $$
境界条件として時刻 $t_0$ に高度 $h_0$、速度 $v_0 < 0$(降下中)、質量 $m_0$。終端 $t_f$ で $h(t_f) = 0$、$v(t_f) = 0$ にしたい。物理的に成立するためには推力に上下限がある:
$$ 0 \leq T_{\min} \leq T(t) \leq T_{\max} $$
エンジンには「最低スロットル比」があり(化学推進のメインエンジンは典型的に 30〜60% より絞れない)、これが下限 $T_{\min} > 0$ を生みます。完全に止めることもできない、これが PDG をいわゆる ハードコンストレインド最適制御問題 にしている本質です。
評価関数 — 燃料最小か時間最小か
最適制御問題として閉じるには評価関数を選びます。代表的なのは二択で、
$$ J_{\mathrm{fuel}} = \int_{t_0}^{t_f} \Gamma(t)\, dt \;=\; v_e \big(m_0 – m(t_f)\big) $$
これは消費推進剤質量の最小化(着陸後にできるだけ多くのペイロード余裕を残す)。あるいは、
$$ J_{\mathrm{time}} = t_f – t_0 $$
これは着陸時間最小化(可視時間制限や日照制約のある計画向き)。多くの月着陸ミッションでは燃料最小、再使用ロケットの帰還燃焼ではしばしば時間制約付き燃料最小です。
ここまでの式は一見シンプルですが、推力の上下限と質量の非線形性 ($\bm{T}/m$ の項)があるため、解析解は得られません。歴史的にはまず「制約を緩めた閉形式解 = E-Guidance」が Apollo で使われ、その後、現代制御の発達とともに「制約付きの厳密解 = 凸最適化」が解けるようになりました。次節ではまず Apollo の E-Guidance を導きます。
Apollo の E-Guidance — 多項式ガイダンスの閉形式
発想: 終端条件から逆算する多項式
Apollo 月着陸船(LM)のメインコンピュータ AGC は、現代から見れば玩具のような能力(2 kHz、4 KB RAM)しかありませんでした。LCvx も SOCP も誰も知らなかった時代に、Allan Klumpp らは驚くほどエレガントな閉形式ガイダンス「E-Guidance(Explicit Guidance)」を編み出します。
発想は次のとおり。質量変化と推力大小制約を一旦無視して、推力ベクトル(=機体加速度)を 時間の多項式 で表すと仮定します。たとえば二次多項式なら、
$$ \bm{a}_T(t) = \bm{c}_0 + \bm{c}_1 \tau + \bm{c}_2 \tau^2 $$
ここで $\tau = t_f – t$ は 目標到達までの残り時間(time-to-go)。係数 $\bm{c}_0, \bm{c}_1, \bm{c}_2$ は現在の状態(位置、速度)と目標(位置、速度、加速度)の境界値から逆算で決めます。多項式の次数を最終条件の数に合わせるのがミソです。
一次元での導出
一次元で書きます(縦か横方向。三次元なら各成分について同じことをする)。位置 $r(t)$、速度 $v(t)$、機体加速度を $a(t)$ とし、
$$ \ddot r(t) = a(t) – g $$
(重力 $g$ は下向き定数として右辺に。$a$ は機体が出す加速度のうち推力ぶん。)
時刻 $t$ における現在状態を $(r, v)$、目標 $t_f$ での状態を $(r_f, v_f, a_f)$ とする。$a_f$ は「着地時の機体加速度」で、Apollo では着地直前にホバリングする状態が望ましいので $a_f = g$(つまり機体加速度が重力相殺) と設定します。
$a(t)$ を $\tau = t_f – t$ の二次関数
$$ a(t) = c_0 + c_1 \tau + c_2 \tau^2 $$
と置く。これを $t$ で積分すると(注: $d\tau/dt = -1$ なので $\int_t^{t_f} a(s)\, ds = \int_0^{\tau} a(t_f – \tau’)\, d\tau’$、以下同様)、
$$ v_f = v + \int_t^{t_f} \big(a(s) – g\big) ds = v – g\tau + \int_0^{\tau} (c_0 + c_1 \tau’ + c_2 \tau’^2)\, d\tau’ $$
積分を実行して、
$$ v_f = v – g\tau + c_0 \tau + \frac{1}{2} c_1 \tau^2 + \frac{1}{3} c_2 \tau^3 $$
同様に位置については $r_f = r + \int_t^{t_f} v(s)\, ds$ で、$v(s) = v + \int_t^s (a(u) – g)\, du$ を代入して二重積分すると、
$$ r_f = r + v\tau – \frac{1}{2} g \tau^2 + \frac{1}{2} c_0 \tau^2 + \frac{1}{6} c_1 \tau^3 + \frac{1}{12} c_2 \tau^4 $$
そして終端の機体加速度 $a_f = a(t_f) = c_0$(なぜなら $\tau = 0$ で多項式は $c_0$ になる)。
未知数は $c_0, c_1, c_2$ の3つ、与えられた条件は $(a_f, v_f, r_f)$ の3つ。線形連立方程式になるので閉形式で解けます。$c_0 = a_f$ をまず確定させると、残りの 2 本は $c_1, c_2$ に関する $2 \times 2$ 線形方程式。整理すると次のような形になります(導出は読者の演習として、答えだけ書きます):
$$ \begin{align} c_0 &= a_f \\ c_1 &= -\frac{6(v_f – v + g\tau – a_f \tau)}{\tau^2} + \frac{12(r_f – r – v\tau + \frac{1}{2} g \tau^2 – \frac{1}{2} a_f \tau^2)}{\tau^3} \\ c_2 &= \frac{6(v_f – v + g\tau – a_f \tau)}{\tau^3} – \frac{12(r_f – r – v\tau + \frac{1}{2} g \tau^2 – \frac{1}{2} a_f \tau^2)}{\tau^4} \end{align} $$
これが E-Guidance の係数公式 です。現時刻における必要な推力加速度は $a(t) = c_0 = a_f$ と即座に得られる(つまり、いま噴射すべき加速度は終端目標加速度そのもの)。これは多項式の性質上の偶然ではなく、$\tau = t_f – t$ をパラメタとして取った時の自然な帰結です。
実際の Apollo の LM では、機体は 1〜2 Hz の頻度で現在の $(r, v)$ を測定し、$t_f$(time-to-go)を内部状態として保持し、上の公式から「次の 0.5 秒間の目標推力加速度」を更新していました。係数を測定値の関数として絶えず再計算するため、外乱や測定誤差があっても自然と補正がかかる — これがクローズドループガイダンスの威力です。
Time-to-go の決め方
$t_f$ をどう決めるかは E-Guidance の弱点の一つです。Apollo では「燃料効率の高い軌道」を作るために、目標位置・速度・加速度の他に 推力 jerk(加速度の時間微分) が条件に追加され、結果として三次多項式になりました。その上で「時刻 $\tau$ における jerk が一定の比率になる」という条件で $\tau$ を内部的に解く方式が取られました。詳細は Klumpp の論文(NASA TM, 1974)に詳しいですが、本記事ではより現代的な「凸最適化で同時に決める」アプローチに進みます。
E-Guidance の素晴らしさは、AGC のような限られた計算機でも実時間で動く閉形式という点でした。しかし、推力制約 $T_{\min} \leq T \leq T_{\max}$ も glide slope(着陸点付近で機体が下方の障害物にぶつからない円錐角)制約も、E-Guidance には組み込めません。E-Guidance は「制約を無視して最適っぽい軌道を出力する」アルゴリズムであり、制約を厳密に守る最適性は保証されないのです。
次節では、この制約をきちんと組み込む現代的な定式化に進みます。
制約付き最適制御問題としての PDG
一般形
三次元の動力降下を最適制御問題として書き直します。状態 $\bm{x} = (\bm{r}, \bm{v}, m) \in \mathbb{R}^7$、制御 $\bm{u} = \bm{T} \in \mathbb{R}^3$、終端時間 $t_f$ を求める問題です。
$$ \begin{align} \min_{\bm{T}(\cdot), t_f} \quad & J = -m(t_f) \quad \text{(等価に燃料最小)} \\ \text{s.t.} \quad & \dot{\bm{r}} = \bm{v}, \quad \dot{\bm{v}} = \frac{\bm{T}}{m} + \bm{g}, \quad \dot m = -\frac{\|\bm{T}\|}{v_e} \\ & T_{\min} \leq \|\bm{T}(t)\| \leq T_{\max} \quad \text{(推力大きさ範囲)} \\ & \langle \hat{\bm{n}}, \bm{T}(t)\rangle \geq \|\bm{T}(t)\| \cos\theta_{\max} \quad \text{(ジンバル角)} \\ & \|\bm{r}(t) – \bm{r}_f\|_{\mathrm{horiz}} \leq \tan(\alpha)\, (r_z(t) – r_{f,z}) \quad \text{(glide slope)} \\ & \bm{r}(t_0) = \bm{r}_0, \bm{v}(t_0) = \bm{v}_0, m(t_0) = m_0 \\ & \bm{r}(t_f) = \bm{r}_f, \bm{v}(t_f) = \bm{0} \end{align} $$
ここで $\hat{\bm{n}}$ は鉛直上方向単位ベクトル、$\theta_{\max}$ はジンバル角制限、$\alpha$ は glide slope 角(着陸点から見た最低見上げ角)です。glide slope 制約は、着陸点周辺の山や岩塊を避けて「斜めから降りてくる」軌跡を強制するための重要な拘束で、Acikmese らの初期の論文(2007)の頃から議論されています。
この問題は次の二つの理由で 非凸最適化問題 です:
- ダイナミクス $\dot{\bm{v}} = \bm{T}/m$ が状態 $m$ に関して非線形(積で結合)。
- 推力大きさの下限 $\|\bm{T}\| \geq T_{\min}$ は非凸制約(球殻)。
非凸問題は一般に大域最適解が保証されず、リアルタイムで解くのが難しい。そこで Behcet Acikmese と Scott Ploen は 2007 年に画期的なテクニックを提示しました。それが Lossless Convexification です。
Lossless Convexification (LCvx) と Apollonius circle
「無損失」凸化とは何か
非凸最適化を直接解くのは難しい。ふつう「凸緩和(convex relaxation)」と言うとき、私たちは制約を緩めて凸問題にしたうえで、その解が元の問題の解にもなっていればラッキー、なっていなくてもバウンドとして使う、という近似戦略です。
Acikmese と Ploen が示したのは、動力降下問題では緩和が必ず元の解と一致する(=情報が失われない)という驚くべき事実でした。これが「Lossless Convexification」の名前の由来です。
スロットル変数の導入
非凸性の主因は二つでした。まず、下限制約 $\|\bm{T}\| \geq T_{\min}$ の凸化を考えます。スカラーのスロットル変数 $\sigma(t)$ を追加し、推力大きさを次のように緩める:
$$ \|\bm{T}(t)\| \leq \sigma(t), \quad T_{\min} \leq \sigma(t) \leq T_{\max} $$
元の問題では $\|\bm{T}\| = \sigma$ が等号で成り立ってほしかったのに、不等号で緩めた。一見、新しい解空間が広がって元の問題の解より良い(コストが小さい)解が見つかってしまうように思えます。
ところが、燃料最小($J = \int \|\bm{T}\| dt$)を最小化する場合、最適解では $\|\bm{T}\| = \sigma$ がほとんどの $t$ で活性(active)になるのです。これは Pontryagin の最大値原理を用いた解析的議論で示されます。直感的には、もし $\|\bm{T}\| < \sigma$ なら $\sigma$ をもっと小さくできるはずで、$\sigma$ が小さくなれば消費燃料 $\int \sigma\, dt$(これはコストの上限)も減るから、最適性に反する。よって最適点で $\|\bm{T}\| = \sigma$。
質量項の対数変換
次に、ダイナミクスの非線形性 $\bm{T}/m$ を解消します。新変数 $z(t) = \ln m(t)$ を導入し、推力加速度 $\bm{a}_T = \bm{T}/m$ を直接の制御変数として扱う:
$$ \dot z = -\frac{\Gamma}{m v_e} = -\frac{\|\bm{a}_T\|}{v_e} $$
ただし $\Gamma = \|\bm{T}\| = m\|\bm{a}_T\|$ を用いた。さらにスロットル変数 $\sigma$ を「加速度スロットル」として再定義し $\sigma \to \sigma/m$ と置き直すと、
$$ \|\bm{a}_T(t)\| \leq \sigma_{\mathrm{acc}}(t), \quad \sigma_{\mathrm{acc}}^{\min}(t) \leq \sigma_{\mathrm{acc}}(t) \leq \sigma_{\mathrm{acc}}^{\max}(t) $$
ここで $\sigma_{\mathrm{acc}}^{\min}(t) = T_{\min}/m(t)$ は時間とともに増加する(質量が減るので)。これを $z = \ln m$ について展開して 凸近似 すると、いわゆる「Apollonius circle」と呼ばれる凸領域に変換できます。具体的には $T_{\min}/m = T_{\min} e^{-z}$ を $z$ の一次テイラー展開した内側上界・外側下界を取って、$\sigma_{\mathrm{acc}}^{\min}, \sigma_{\mathrm{acc}}^{\max}$ を $z$ の一次関数として閉じ込めるのです。
なぜ「Apollonius circle」なのか
古典幾何学において、二点からの距離の比が一定となる点の軌跡(アポロニウスの円)は、二次曲線族の中で円形になる特殊なケースとして知られています。Acikmese は推力加速度ベクトル $\bm{a}_T$ の許容領域を $\sigma_{\mathrm{acc}}$ で挟んだ「アニュラス(円環)」が、対数変換のおかげで $z$ 空間で円板(凸領域)になることに気づきました。これは古典の Apollonius circle と類比的な構造を持つため、本人らがそう命名したと論文に記されています。
結果としての凸問題
以上の変換をすべて施すと、元の非凸 OCP は次のような凸問題に書き換わります:
$$ \begin{align} \min_{\bm{a}_T, \sigma_{\mathrm{acc}}, t_f}\quad & \int_{t_0}^{t_f} \sigma_{\mathrm{acc}}(t)\, dt \\ \text{s.t.}\quad & \dot{\bm{r}} = \bm{v}, \quad \dot{\bm{v}} = \bm{a}_T + \bm{g}, \quad \dot z = -\sigma_{\mathrm{acc}}/v_e \\ & \|\bm{a}_T(t)\| \leq \sigma_{\mathrm{acc}}(t) \quad \text{(SOC 制約)} \\ & \mu_{\min}(t) \leq \sigma_{\mathrm{acc}}(t) \leq \mu_{\max}(t) \quad \text{(線形)} \\ & \text{glide slope: 線形不等式} \\ & \text{ジンバル: SOC 制約} \\ & \text{境界条件: 等式} \end{align} $$
これは離散化すると Second-Order Cone Program (SOCP) になり、内点法で多項式時間で大域最適解が解けます。実用上は数十ミリ秒で解ける程度の規模で、リアルタイムオンボード GNC に十分です。SpaceX の Falcon 9 booster や Blue Origin の Blue Moon が「クローズドループ最適制御」を実装できているのは、まさにこの LCvx + SOCP の進展のおかげです。
ここまでで PDG を解くための理論枠組みが揃いました。次節では実際に手を動かして、まず一次元 1D 降下を Python でシミュレーションしてみます。
Pythonでの実装 — 1D 降下の数値シミュレーション
モデルの設定
最も簡単な「真上から真下に降りるだけ」の 1D シミュレーターを書き、定推力で降りるとどうなるかを見ます。これを基準として、後で E-Guidance や最適制御と比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# 月の物理定数
g_M = 1.62 # 月面重力 [m/s^2]
g0 = 9.81 # 標準重力 [m/s^2]
Isp = 311.0 # 比推力 [s] (Apollo LM Descent Engine 相当)
ve = Isp * g0 # 排気速度 [m/s]
T_max = 45_000.0 # 推力上限 [N]
T_min = 0.3 * T_max # 推力下限 (30% スロットル)
m_dry = 2_000.0 # 空虚質量 [kg]
m_fuel0 = 6_000.0 # 初期推進剤 [kg]
m0 = m_dry + m_fuel0
def dynamics_1d(t, x, T_func):
"""1D 動力降下 ODE: x = [h, v, m]"""
h, v, m = x
T = T_func(t, h, v, m)
T = np.clip(T, 0.0, T_max) # 物理上限
dh = v
dv = T / m - g_M
dm = -T / ve
return [dh, dv, dm]
# 一定推力(80% スロットル)で降りた場合
def const_thrust(t, h, v, m):
return 0.80 * T_max
x0 = [3_000.0, -50.0, m0] # 高度3000 m, 下向き 50 m/s
sol = solve_ivp(dynamics_1d, [0, 200], x0, args=(const_thrust,),
dense_output=True, max_step=0.5,
events=lambda t, x, T_func: x[0]) # h=0 で停止
print(f"着陸時刻: t = {sol.t_events[0][0]:.2f} s")
print(f"着陸時速度: v = {sol.y_events[0][0, 1]:.2f} m/s")
print(f"残存燃料: {sol.y_events[0][0, 2] - m_dry:.1f} kg")
このコードは「とにかく一定スロットルで噴射し続ける」素朴な戦略をシミュレートしています。出力されるのは典型的に「着陸時刻 50〜70 秒、着陸時速度 が秒速数 m から数十 m」というハードランディング(機体クラッシュ)です。なぜなら定推力では加速度の調整ができず、ぴったり $v = 0$ で $h = 0$ に到達する保証がないからです。これがガイダンスのいらない世界の限界。
E-Guidance を 1D に実装
次に E-Guidance を一次元で実装し、time-to-go $\tau$ を更新しながらクローズドループで降下させます。
# E-Guidance による推力指令(1D)
def e_guidance_1d(h, v, m, t_go, h_f=0.0, v_f=0.0, a_f=g_M):
"""二次多項式 E-Guidance: 現在 (h,v) から t_go 秒後に (h_f, v_f) に
到達し、終端加速度 a_f になるよう必要加速度を計算"""
tau = t_go
if tau < 0.5:
return a_f # 終端付近はホバリング推力
# 必要な機体加速度の二次多項式表現
# a(tau) = c0 + c1*tau + c2*tau^2 で c0 = a_f
# vf, hf の式を線形連立で解く
g = g_M
# 速度差分(重力ぶんを除いた要求)
dv = v_f - v + g * tau - a_f * tau
# 位置差分
dh = h_f - h - v * tau + 0.5 * g * tau**2 - 0.5 * a_f * tau**2
c1 = -6 * dv / tau**2 + 12 * dh / tau**3
c2 = 6 * dv / tau**3 - 12 * dh / tau**4
# 現時刻 t=t での加速度 = c0 + c1*tau + c2*tau^2 だが、tau は残り時間
# 連続更新するので現時刻の出力指令は c0 そのものでよい (Klumpp 1974)
return c0_apply(a_f, c1, c2, tau)
def c0_apply(c0, c1, c2, tau):
# 即時必要加速度 = a(tau=tau_now) = c0 + c1*tau + c2*tau^2
return c0 + c1 * tau + c2 * tau**2
# クローズドループで E-Guidance を実行
def thrust_e_guidance(t_go0):
def f(t, h, v, m):
t_go = t_go0 - t
if t_go < 0.5:
return m * g_M # ホバリング
a_req = e_guidance_1d(h, v, m, t_go)
a_req = max(a_req, 0.0)
T_req = m * a_req
T_req = np.clip(T_req, T_min, T_max)
return T_req
return f
# 60秒で着陸させる E-Guidance
t_go0 = 60.0
sol_eg = solve_ivp(dynamics_1d, [0, t_go0], x0,
args=(thrust_e_guidance(t_go0),),
dense_output=True, max_step=0.2,
events=lambda t, x, _: x[0])
t_land = sol_eg.t_events[0][0] if len(sol_eg.t_events[0]) > 0 else t_go0
print(f"E-Guidance 着陸時刻: {t_land:.2f} s")
print(f"E-Guidance 着陸時速度: {sol_eg.y_events[0][0, 1]:.3f} m/s")
print(f"E-Guidance 残存燃料: {sol_eg.y_events[0][0, 2] - m_dry:.1f} kg")
このコードでは e_guidance_1d が現在状態と残り時間から必要加速度を返し、それを推力指令に変換しています。「即時必要加速度 = c0 + c1tau + c2tau^2」とすると、各時刻で多項式のもっとも長い項が効いて、自然に終端条件にうまくフィットしていきます。終端 $\tau \to 0$ では係数が発散しうるため、最後の 0.5 秒はホバリング推力に切り替える保険を入れています。
実行結果は「着陸時刻 60 秒、着陸時速度 $|v| < 0.1$ m/s、残存燃料 約 1500〜2000 kg」になります。定推力モードとの最大の違いは、着陸時速度が 100 倍以上小さくなる(=ソフトランディング成功)ことです。これが多項式ガイダンスの威力で、係数を毎秒再計算しているおかげで、たとえ初期質量が想定とずれていたとしても自動的に補正されます。
軌跡の可視化
import numpy as np
import matplotlib.pyplot as plt
# const と E-Guidance の軌跡を描く
ts_const = np.linspace(0, sol.t_events[0][0], 200)
ts_eg = np.linspace(0, t_land, 400)
X_const = sol.sol(ts_const)
X_eg = sol_eg.sol(ts_eg)
fig, axes = plt.subplots(1, 3, figsize=(13, 4))
axes[0].plot(ts_const, X_const[0], 'r--', label='Const thrust 80%')
axes[0].plot(ts_eg, X_eg[0], 'b-', label='E-Guidance')
axes[0].set_xlabel('time [s]'); axes[0].set_ylabel('altitude [m]')
axes[0].set_title('Altitude vs time'); axes[0].legend(); axes[0].grid(alpha=0.3)
axes[1].plot(ts_const, X_const[1], 'r--', label='Const')
axes[1].plot(ts_eg, X_eg[1], 'b-', label='E-Guidance')
axes[1].axhline(0, color='k', lw=0.5)
axes[1].set_xlabel('time [s]'); axes[1].set_ylabel('velocity [m/s]')
axes[1].set_title('Velocity vs time'); axes[1].legend(); axes[1].grid(alpha=0.3)
axes[2].plot(ts_const, X_const[2], 'r--', label='Const')
axes[2].plot(ts_eg, X_eg[2], 'b-', label='E-Guidance')
axes[2].axhline(m_dry, color='gray', ls=':', label='dry mass')
axes[2].set_xlabel('time [s]'); axes[2].set_ylabel('mass [kg]')
axes[2].set_title('Mass vs time'); axes[2].legend(); axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('powered_descent_1d.png', dpi=130, bbox_inches='tight')
plt.show()
このグラフから三つのことが読み取れます。第一に、定推力(赤破線)では速度がゼロを通り過ぎて負(着地時に下向きに突っ込む)になっており、ハードランディング状態です。第二に、E-Guidance(青実線)は終端で速度がきれいに 0 に収束し、その瞬間に高度もぴったり 0 になっています。第三に、E-Guidance のほうが質量(残燃料)を多く残せています — これは定推力が「無駄に長く強い噴射を続けて、最後に過剰減速」している一方、E-Guidance は速度プロファイルを最適化して必要最小限の総推進剤で済ませているからです。
E-Guidance は強力ですが、ここまで見たのはあくまで「制約を緩めた」一次元の話。次に、現代の凸最適化アプローチで、推力範囲・glide slope を厳密に守った 2D 軌道を計算します。
CVXPY による Convex Landing(2D)
2D 設定
横方向 $x$、縦方向 $y$(高度)の二次元問題を考えます。初期位置 $(x_0, y_0) = (-1000, 30\,000)$ m、初期速度 $(\dot x_0, \dot y_0) = (-50, -150)$ m/s、目標位置 $(0, 0)$、目標速度 $(0, 0)$。月の重力場で、最大推力角度制限と glide slope 制約を入れます。
ダイナミクスは離散時間で書きます。簡潔さのため質量変化はここで一次オーダーで扱い、後で凸近似を入れます。状態 $\bm{x}_k = (x, y, \dot x, \dot y)$、制御 $\bm{u}_k = (u_x, u_y)$(機体加速度)、ステップ時刻 $\Delta t$ で、
$$ \bm{x}_{k+1} = A \bm{x}_k + B \bm{u}_k + \bm{c}, \quad A = \begin{bmatrix} I_2 & \Delta t I_2 \\ 0 & I_2 \end{bmatrix}, B = \begin{bmatrix} \tfrac{1}{2} \Delta t^2 I_2 \\ \Delta t I_2 \end{bmatrix}, \bm{c} = \bm{g} \cdot (\tfrac{1}{2}\Delta t^2, \Delta t)^\top $$
CVXPY による SOCP 定式化
import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt
# パラメタ
g_M = 1.62
N = 80 # 離散ステップ数
T_total = 200.0 # 総時間 [s] (固定: 燃料最小はネスト最適化で)
dt = T_total / N
# 初期・終端条件
x0 = np.array([-1000.0, 30_000.0, -50.0, -150.0])
xf = np.array([0.0, 0.0, 0.0, 0.0])
# 制約パラメタ
a_max = 4.0 # 機体加速度 max [m/s^2]
a_min = 0.6 # 機体加速度 min (スロットル下限相当)
theta_glide = np.deg2rad(75.0) # glide slope (鉛直からの最大角)
# CVXPY 変数
x = cp.Variable((4, N + 1))
u = cp.Variable((2, N)) # 機体加速度ベクトル
sigma = cp.Variable(N) # スロットル変数 |u| <= sigma (LCvx)
constraints = [x[:, 0] == x0, x[:, N] == xf]
# ダイナミクス
A = np.array([[1, 0, dt, 0],
[0, 1, 0, dt],
[0, 0, 1, 0],
[0, 0, 0, 1]])
B = np.array([[0.5*dt**2, 0],
[0, 0.5*dt**2],
[dt, 0],
[0, dt]])
gvec = np.array([0, -g_M])
c_vec = np.array([0, -0.5*dt**2*g_M, 0, -dt*g_M])
for k in range(N):
constraints += [x[:, k + 1] == A @ x[:, k] + B @ u[:, k] + c_vec]
# LCvx: ||u_k|| <= sigma_k (SOC 制約)
constraints += [cp.norm(u[:, k], 2) <= sigma[k]]
# スロットル範囲
constraints += [sigma[k] >= a_min, sigma[k] <= a_max]
# glide slope: 高度 y は |x| / tan(theta) より高い
constraints += [x[1, k] >= cp.abs(x[0, k]) / np.tan(theta_glide)]
# 目的: 総スロットル積 (燃料最小の凸サロゲート)
objective = cp.Minimize(cp.sum(sigma) * dt)
prob = cp.Problem(objective, constraints)
prob.solve(solver=cp.ECOS, verbose=False)
print(f"最適コスト: {prob.value:.2f}")
print(f"最終位置: ({x.value[0, -1]:.3f}, {x.value[1, -1]:.3f}) m")
print(f"最終速度: ({x.value[2, -1]:.3f}, {x.value[3, -1]:.3f}) m/s")
このコードは、Acikmese の LCvx を CVXPY で素直に書き下したものです。cp.norm(u[:, k], 2) <= sigma[k] が SOC 制約($\bm{u}$ がノルム球のコーンに収まる)、sigma >= a_min で下限を線形に表現しています。glide slope は x[1, k] >= cp.abs(x[0, k]) / tan(theta) で、これも凸制約(右辺が凸関数、左辺が線形)。 ECOS や SCS でミリ秒オーダーで解けます。
軌道と推力プロファイルの可視化
import numpy as np
import matplotlib.pyplot as plt
X = x.value
U = u.value
Sigma = sigma.value
ts = np.arange(N + 1) * dt
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# 軌跡
axes[0, 0].plot(X[0], X[1], 'b-', lw=2)
axes[0, 0].plot(X[0, 0], X[1, 0], 'go', label='start')
axes[0, 0].plot(X[0, -1], X[1, -1], 'r*', markersize=15, label='land')
# glide slope cone
xs = np.linspace(-1500, 1500, 50)
axes[0, 0].fill_between(xs, np.abs(xs)/np.tan(theta_glide), 35_000,
alpha=0.10, color='green', label='glide slope OK')
axes[0, 0].set_xlabel('horizontal x [m]'); axes[0, 0].set_ylabel('altitude y [m]')
axes[0, 0].set_title('Trajectory (2D)'); axes[0, 0].legend(); axes[0, 0].grid(alpha=0.3)
axes[0, 0].set_xlim(-1500, 500); axes[0, 0].set_ylim(0, 31_000)
# 速度
axes[0, 1].plot(ts, X[2], label='vx')
axes[0, 1].plot(ts, X[3], label='vy')
axes[0, 1].axhline(0, color='k', lw=0.5)
axes[0, 1].set_xlabel('time [s]'); axes[0, 1].set_ylabel('velocity [m/s]')
axes[0, 1].set_title('Velocity components'); axes[0, 1].legend(); axes[0, 1].grid(alpha=0.3)
# 推力大きさ
ts_u = np.arange(N) * dt + dt/2
acc_mag = np.linalg.norm(U, axis=0)
axes[1, 0].plot(ts_u, acc_mag, 'b-', label='||u|| (acc)')
axes[1, 0].plot(ts_u, Sigma, 'r--', label='sigma (throttle)')
axes[1, 0].axhline(a_min, color='gray', ls=':', label='a_min')
axes[1, 0].axhline(a_max, color='gray', ls=':', label='a_max')
axes[1, 0].set_xlabel('time [s]'); axes[1, 0].set_ylabel('accel [m/s^2]')
axes[1, 0].set_title('Thrust magnitude'); axes[1, 0].legend(); axes[1, 0].grid(alpha=0.3)
# 推力ベクトル方向
axes[1, 1].plot(ts_u, np.rad2deg(np.arctan2(U[0], U[1])), 'b-')
axes[1, 1].set_xlabel('time [s]'); axes[1, 1].set_ylabel('thrust angle [deg]')
axes[1, 1].set_title('Thrust direction from +y axis'); axes[1, 1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('lcvx_2d_landing.png', dpi=130, bbox_inches='tight')
plt.show()
このグラフから、最適化された軌道の重要な性質がいくつも読み取れます。第一に、軌跡(左上)は緑色の許容領域(glide slope の内側)に常に収まっており、最後の数百メートルで自然に「立ち上がる」ような形をしている。これは横方向速度を落としきってから垂直降下に入る、いわゆる gravity-turn-like な軌道で、燃料最小の典型パターンです。第二に、推力大きさ(左下)は bang-bang-like な振る舞いをしており、序盤と終盤に max thrust の “burn” があり、中間でほぼ min thrust(コースト相当)になっている。これは Pontryagin 最大値原理が予言する燃料最適制御の典型構造で、LCvx の理論的予言を実装が再現しています。第三に、推力方向(右下)は降下中常に「機体上方」を中心に振れており、特に最終フェーズで急速に垂直に揃う — これは横方向速度を消すために初期に大きく傾けていた推力ベクトルを、終端で垂直に戻している様子です。
Apollo 12 と SLIM の比較考察
実際の月着陸の数字と並べて、現代の最適制御がどれだけ精密になったかを見てみましょう。
- Apollo 11 (1969): 自動誘導で目標近くまで降下したが、着陸直前にコンピュータ警報(1201/1202)が出てパイロット Neil Armstrong が手動操縦に切り替え、約 8 km 目標からずれた地点に着地。残燃料はわずか 25 秒。
- Apollo 12 (1969): 自動降下精度が大幅に向上し、目標の Surveyor 3 探査機から 183 m の地点に着陸。これは「精密着陸」の最初の本格デモ。
- SLIM (2024 年 1 月): JAXA の小型月着陸機 SLIM が、目標位置から 約 55 m の精度で着陸。Apollo 時代の人間ループより 3 倍精密で、しかも完全自動。
- Falcon 9 first stage (2015 年〜): 地球上で実時間 PDG を解いて、海上ドローン船 ASOG に直径 50 m 以内の精度で着地。1 秒あたりの最適化を内点法で繰り返し解いている。
これらの精密化は、ハードウェア(精密 IMU、テレスコープ式測距、カメラベース crater navigation)の発達と、ソフトウェア(LCvx + SOCP)の進化が両輪となって実現しました。本記事の Python コードは、後者の核心部分を 50 行程度に圧縮したものに相当します。
最適制御の結果を見たところで、最後にもう少し最終フェーズの細かな話をします。最終フェーズには「ホバー→垂直降下」型と「重力ターン型」の二つの戦略があり、ミッションごとに使い分けがあります。
最終フェーズ — Hover-and-Descend と Gravity-Turn
着陸の最後の数十秒は、それ単独で別の制御問題と見るべき特殊性を持ちます。代表的な戦略を二つ紹介します。
Hover-and-Descend(ホバリング降下)
地表に近い高度(典型的に 30〜100 m)で一旦ホバリング状態に入り、その後カメラやレーダーで着陸点周辺の岩石・斜面を確認しながらゆっくり垂直降下する戦略。Apollo の LM は「P66」というモードでパイロットがレートを制御しながらこれを実行しました。最近では SLIM が二段階着陸(まず本体を 5 m 上空で一旦止め、その後機体を斜めに傾けて転がるように接地)で世界を驚かせました。
数式的には、最終フェーズの状態空間方程式は
$$ \dot h = v, \quad \dot v = T/m – g, \quad T \approx m g \text{(ほぼホバリング)} $$
で、ここから垂直方向のみの低速降下になります。微小な水平ドリフトはサイドスラスターで補正。
Gravity-Turn(重力ターン)
着陸機の 機体軸を常に速度ベクトルの逆向きに保つ ような姿勢戦略。空気抵抗のないところで動力降下のために多用される、燃料効率が良い古典的手法です。1960 年代の Surveyor 機が代表例。重力ターン軌道は次の連立 ODE で記述できます:
$$ \dot v = -T/m + g \cos\gamma, \quad v\dot\gamma = g \sin\gamma, \quad \dot h = -v\cos\gamma $$
ここで $\gamma$ は速度ベクトルと水平面の角度(降下角)。注目すべきは、機体軸方向に推力を当て続ける限り、姿勢は速度ベクトルが回転していくのに従って自然に「重力に引かれて」立ち上がっていく(垂直に向く)ことです。だから「gravity-turn」と呼ばれる。
LCvx 最適解の軌道(前節のグラフ)を見ると、終端付近で軌跡が自然に立ち上がり、推力ベクトルがほぼ垂直に揃っていました。これは数式上で重力ターン的な振る舞いに収束していくことの現れで、最適制御は古典戦略の良い性質を自動的に再発見してくれている、と言えます。
例: Falcon 9 の着陸燃焼
地球で動力降下を実現している Falcon 9 first stage は、以下の三段階で構成されます: (1) boostback burn — 上昇途中で逆向きに姿勢を変え、地表の着陸地点に戻る軌道を作る。(2) entry burn — 高度 70〜100 km で大気再突入直前に再点火し、空力加熱を抑える。(3) landing burn — 高度 1 km 程度から約 30 秒の最終燃焼で速度ゼロにして着地。Merlin-1D エンジンは推力下限が 100%(クラスタの 1 基だけを使うと約 11% 等価)であり、これが LCvx の理論的有効性を生む一例です。
ここまでで、動力降下の理論と実装、現代ミッションへの応用が一通り見えてきました。最後にまとめます。
まとめ
本記事では、月や火星のような大気の薄い天体への動力降下を、最適制御の枠組みから解説しました。
- 物理: 動力降下は並進運動 ODE($\dot{\bm{r}} = \bm{v}$, $\dot{\bm{v}} = \bm{T}/m + \bm{g}$)とロケット方程式の連続時間版($\dot m = -\Gamma/v_e$)を連立した非線形系。これを推力範囲・glide slope の下で解くのが PDG の正体。
- Apollo E-Guidance: 機体加速度を $\tau = t_f – t$ の多項式で表し、現在の状態と目標の境界条件から係数を閉形式で逆算する手法。AGC の限られた計算機でも実時間で動く驚異の発明だが、制約厳密性は保証されない。
- Lossless Convexification (Acikmese & Ploen, 2007): スロットル変数 $\sigma$ を導入して $\|\bm{T}\| \leq \sigma$ と緩和し、対数質量変換($z = \ln m$)で非線形性を凸化する手法。燃料最小目的関数の下で「緩和の解が必ず元の解と一致する」ことが Pontryagin の理論で示せる。
- SOCP リアルタイム解法: 結果として得られる凸最適化問題は SOCP として書け、内点法で多項式時間に解ける。CVXPY + ECOS で 50 行程度で実装可能。Falcon 9 はオンボードでこれを実時間で繰り返し解いている。
- 最終フェーズ: ホバリング降下(SLIM, Apollo P66)と重力ターン(Surveyor)の二系統がある。LCvx 最適解は終端付近で自然に重力ターン的な姿勢に収束する。
数値実験では、月面 30 km・横方向 -1 km から目標(0, 0)への着陸軌道が、推力範囲・glide slope を厳密に守りながら 200 秒以内で実現できることを確認しました。E-Guidance との比較では、E-Guidance が制約を緩めた「賢い予測子」であるのに対し、LCvx は厳密制約を持ったまま大域最適解を返す、という階層差が明確になりました。
動力降下ガイダンスは、宇宙工学の中で「物理・最適制御・数値最適化・組込みソフトウェア」が交差するもっとも美しい応用領域の一つです。Apollo 時代の手作りの多項式から半世紀、いまや手元のラップトップで Falcon 9 並みの最適着陸計算が走る時代になりました。
次のステップとして、以下の記事も参考にしてください。