09. 月面への鉛直軟着陸

月面へ降りていく着陸機が、地面にぶつからずに速度ゼロで接地するには、エンジンを噴射して減速しなければなりません。 燃料を最も節約するには、早めに少しずつ吹かすのがよいのでしょうか、それとも最後に一気に吹かすのがよいのでしょうか。 この問題には紙と鉛筆で求まる答えがあり、数値解と比べられます。

項目	内容
状態	高度 $h$、速度 $v$（上向きが正）
制御	推力 $T$
目的	推力の時間積分（使う燃料に比例）を最小にする
区間	2区間（各12ノード）、境目の時刻も終了時刻も自由
新しく学ぶこと	積分で表す目的関数（燃料）、境目の時刻を自由にした2区間

状況

着陸機は月面の上空 100 m にいて、毎秒 10 m で降下しています。月の重力加速度は 1.62 m/s² で、何もしなければどんどん速くなります。 エンジンは上向きに最大 4000 N の推力を出せ、推力の大きさは 0 から最大まで自由に変えられます。

ロケットエンジンが使う推進剤の量は、推力 × 噴射時間（力積）に比例します。そこで、推力の時間積分をできるだけ小さくして、 高度 0・速度 0 で接地させることを考えます。着陸機の質量は 1000 kg で、燃料を使っても変わらないとします （降下中に使う燃料が機体の質量に比べて少ない、という近似です）。

この例題では時間を「エンジンを止めて落ちる区間」と「減速する区間」の2つに分けます。 ただし、いつ切り替えるかは決めずに、最適化に選ばせます。

問題設定

運動方程式

$$\begin{aligned} \dot{h} &= v \\[2pt] \dot{v} &= \frac{T}{m} - g \end{aligned}$$

• $h$ は月面からの高度 [m]、$v$ は速度 [m/s] で、上向きを正とします。降下中は $v < 0$ です。
• 推力 $T$ [N] は上向きにはたらき、質量 $m$ で割った $T/m$ が推力による加速度です。
• $g$ は月の重力加速度で、常に下向きにはたらきます。

両方の区間で同じ式を使います。

区間と接続条件

区間1を $0 \le t \le t_1$、区間2を $t_1 \le t \le t_f$ とします。境目の時刻 $t_1$ と終了時刻 $t_f$ はどちらも最適化の変数です。 境目では高度と速度がつながります。

$$h^{(1)}(t_1) = h^{(2)}(t_1), \qquad v^{(1)}(t_1) = v^{(2)}(t_1)$$

推力には接続条件をつけないので、境目で推力が不連続に変わってかまいません。

境界条件

$$\begin{aligned} &h(0) = 100\ \mathrm{m},\quad v(0) = -10\ \mathrm{m/s} \\[2pt] &h(t_f) = 0,\quad v(t_f) = 0 \end{aligned}$$

上空 100 m を毎秒 10 m で降下している状態から始め、高度 0 で速度 0（軟着陸）にします。

制約

両方の区間で、次の範囲を守ります。

$$0 \le T \le 4000\ \mathrm{N}, \qquad 0 \le h \le 100\ \mathrm{m}, \qquad v \le 0$$

推力は上向きのみで、最大推力を超えられません。高度は地面より下に潜らず、出発点より上にも行きません。 速度は常に下向きか 0 で、一度上昇してからまた降りる、という動きは許しません。

目的関数

$$J = \frac{1}{T_{\max}\, t_u} \int_0^{t_f} T\, dt$$

推力の時間積分（力積）を最小にします。前の係数 $1/(T_{\max} t_u)$ は、値を 1 程度にそろえるための定数で、答えには影響しません。 ここで $t_u = \sqrt{h_0/g} \approx 7.86$ s は時間の代表的な大きさです。

これまでの例題の目的関数は「終端での値」（最終時刻、最終高度）でしたが、ここでは時間全体にわたる積分です （目的関数の2つの形）。

パラメータ

記号	意味	値
$h_0$	初期高度	100 m
$v_0$	初期速度（上向き正）	−10 m/s
$g$	月の重力加速度	1.62 m/s²
$m$	着陸機の質量	1000 kg
$T_{\max}$	最大推力	4000 N

内部では、高度を 100 m、速度を $\sqrt{h_0 g} \approx 12.7$ m/s、時間を $t_u$、推力を 4000 N を単位としてスケーリングしています。

初期推定

初期推定には、後で説明する解析解をそのまま与えています。区間1は推力 0 の自由落下、区間2は最大推力での減速で、 区間の時刻も解析解の切り替え時刻・着地時刻にしています。数値解が解析解からずれないかを確かめるための設定です。

結果

高度・速度・推力の時間履歴。青が区間1、橙が区間2で、点は各区間12個のノードでの値、線は区間ごとの多項式です。縦の破線は区間の境目、黒の破線は解析解、赤い点線は最大推力 4000 N です。

降下の様子。左の軌跡は区間1（自由落下）を青、区間2（減速）を橙で描き、炎は噴射を表します。右は速度と推力の時間変化で、黒の破線が解析解、縦の黒線が現在の時刻です。

• 推力は区間1ではずっと 0、区間2ではずっと最大推力 4000 N です。途中の値は使われません。
• エンジンを点火する時刻は 3.632 s で、そのときの高度は 53.0 m、 速度は -15.88 m/s です。
• 着地は 10.305 s で、使った力積は 26,695 N·s です。
• 解析解の点火時刻は 3.632 s、着地時刻は 10.305 s で、数値解と一致しています。

解説

燃料を減らすことは、時間を短くすること

運動方程式 $\dot{v} = T/m - g$ を、出発から着地まで積分してみます。

$$v(t_f) - v(0) = \frac{1}{m}\int_0^{t_f} T\, dt - g\, t_f$$

$v(0) = -10$、$v(t_f) = 0$ を入れて整理すると、使う力積は

$$\int_0^{t_f} T\, dt = m \bigl( |v_0| + g\, t_f \bigr)$$

となります。右辺の $|v_0|$ は「最初の降下速度を打ち消す分」で、どう吹かしても変わりません。 もう一方の $g\, t_f$ は「着地までずっと重力に引かれ続ける分」で、重力損失と呼ばれます。 噴射して空中にとどまっている時間が長いほど、重力に逆らうために燃料を余計に使います。

つまり燃料最小は、着地までの時間 $t_f$ を最短にすることと同じです。 最短時間で、しかも速度 0 で接地するには、ぎりぎりまで加速度が最も下向き（推力 0）で落ち、 最後に加速度が最も上向き（最大推力）で止まるのが一番です。途中の推力を使うと、その分ゆっくり降りることになり、 重力損失が増えます。こうして「自由落下してから最大推力」というバンバン制御が答えになります。

切り替え時刻を手で求める

推力 0 の間の加速度は $-g$、最大推力の間の加速度は $b = T_{\max}/m - g = 2.38\ \mathrm{m/s^2}$（上向き）です。 点火する瞬間の降下速度を $v_s$ とすると、自由落下と減速で進む距離の和が初期高度に等しいことから

$$\frac{v_s^2 - v_0^2}{2g} + \frac{v_s^2}{2b} = h_0$$

が成り立ち、これを解くと

$$\begin{aligned} v_s &= \sqrt{\frac{(v_0^2 + 2 g h_0)\, b}{b + g}} \\[4pt] t_1 &= \frac{v_s - |v_0|}{g} \\[4pt] t_f &= t_1 + \frac{v_s}{b} \end{aligned}$$

となります。数値解はこの値を有効数字の範囲で再現しています。

比べるために、出発直後から一定の推力でじわじわ減速して、ちょうど地面で止まる場合を考えます。 このときの減速度は $v_0^2/(2h_0) = 0.5\ \mathrm{m/s^2}$、着地まで 20 s かかり、力積は $1000 \times (10 + 1.62 \times 20) = 42400$ N·s です。 最適解の約 1.6 倍の燃料を使います。

区間の境目を自由にする意味

推力は点火の瞬間に 0 から 4000 N へ不連続に跳びます。1本の多項式でこれを表すと、例題07のように振動します。 そこで、時間を2つの区間に分け、境目の時刻 $t_1$ を最適化の変数にしました。ソルバーは境目を点火の瞬間に合わせることで、 各区間の推力を一定値として、振動なしに表せます。例題08では境目の時刻を固定しましたが、ここでは境目の位置そのものを最適化で決めています。

試してみよう

• 最大推力 PARAMETERS["max_thrust"] を 3000 N に下げる。減速の加速度が小さくなるので、点火を早める必要があり、着地まで時間がかかって力積が増えることを、上の式とも比べて確かめる。
• 初期速度 PARAMETERS["initial_velocity"] を 0（空中で静止した状態から）にすると、自由落下の区間がどう変わるかを見る。
• 推力の初期推定 thrusts = [0, p["max_thrust"]] を [2000, 2000] にして、解析解から離れた推定から出発しても同じ答えに戻るかを確かめる。

実行する

uv run python examples/09_soft_landing.py

examples/09_soft_landing.py を GitHub で見る