15. 低推力軌道遷移

イオンエンジンのような電気推進は、推進剤の使い方はとても効率的ですが、推力は重力に比べてごくわずかです。そのような宇宙機で、内側の円軌道から外側の円軌道へ移るには、何周も回りながら少しずつ軌道を広げていくことになります。燃料を最も節約するには、いつ、どの向きに吹かせばよいでしょうか。

項目	内容
状態	軌道半径 $r$ 、半径方向速度 $v_r$ 、接線方向速度 $v_t$
制御	推力加速度の正負別の成分 $u_r^+,\ u_r^-,\ u_t^+,\ u_t^-$ （すべて 0 以上）
目的	推力加速度の大きさの時間積分 $\int (\lvert u_r\rvert + \lvert u_t\rvert)\,dt$ を最小にする
区間	1区間、60ノード、終了時刻は自由（上限 55）
新しく学ぶこと	無次元化した軌道の問題、正負に分けた制御による絶対値の目的関数、多数のノードを使う長時間の問題

状況

中心の天体のまわりを半径 1 の円軌道で回っている宇宙機を、半径 4 の円軌道へ移します。推力は各方向に重力（半径 1 での値）の 1 % までしか出せないので、1周では軌道をほとんど広げられません。

量はすべて無次元で表します。長さの単位は最初の軌道半径、速度の単位はその円軌道速度で、重力定数 $\mu$ は 1 になります。このとき最初の軌道を1周する時間は $2\pi \approx 6.28$ 、目標の軌道（半径 4）を1周する時間は $2\pi \times 8 \approx 50.3$ です。遷移にかけてよい時間は最大 55 で、最初の軌道の約 8.8 周分にあたります。

推力は質量あたりの値（推力加速度）で与え、宇宙機の質量の変化は考えません。使う推進剤の量は、推力加速度を時間で積分した速度変化の合計（ $\Delta v$ ）で測ります。

問題設定

座標系

軌道面の中で、位置を中心の天体からの極座標で表します。

$r$ ：中心からの距離（軌道半径）。
$v_r$ ：中心から外向きの速度（半径方向速度）。
$v_t$ ：それに直角で、回っている向きの速度（接線方向速度）。
$u_r,\ u_t$ ：推力加速度の半径方向成分と接線方向成分。

回った角度 $\theta$ は運動方程式の右辺に現れないので、状態に入れていません。図を描くときに $\dot{\theta} = v_t / r$ を積分して求めています。

正負に分けた制御

推力加速度の各成分を、0 以上の2つの量の差で表します。

\begin{aligned} u_r &= u_r^+ - u_r^- \\ u_t &= u_t^+ - u_t^- \\ &u_r^\pm,\ u_t^\pm \ge 0 \end{aligned}

$u_r^+$ は外向き、 $u_r^-$ は内向き、 $u_t^+$ は進む向き、 $u_t^-$ は逆向きの推力です。この分け方の理由は目的関数で説明します。

運動方程式

\begin{aligned} \dot{r} &= v_r \\[2pt] \dot{v}_r &= \frac{v_t^2}{r} - \frac{1}{r^2} + u_r^+ - u_r^- \\[2pt] \dot{v}_t &= -\frac{v_r v_t}{r} + u_t^+ - u_t^- \end{aligned}

$1/r^2$ は重力（ $\mu = 1$ ）、 $v_t^2/r$ は遠心力にあたる項です。 $v_t = 1/\sqrt{r}$ のとき両者がつり合い、円軌道になります。
$-v_r v_t / r$ は、外へ動くと接線方向速度が減ることを表します。推力がなければ角運動量 $r v_t$ が保存されるためです。

境界条件

\begin{aligned} &r(0) = 1,\quad v_r(0) = 0,\quad v_t(0) = 1 \\ &r(t_f) = 4,\quad v_r(t_f) = 0,\quad v_t(t_f) = 0.5 \end{aligned}

始めと終わりはどちらも円軌道です。半径 4 の円軌道速度は $1/\sqrt{4} = 0.5$ で、外側の軌道ほど遅く回ります。軌道の上のどこで終えるか（角度）は自由です。

制約

\begin{aligned} &0 \le u_r^+,\ u_r^-,\ u_t^+,\ u_t^- \le 0.01 \\ &1 \le r \le 4 \\ &1 \le t_f \le 55 \end{aligned}

推力加速度は各成分とも 0.01 までです。大きさではなく成分ごとの上限なので、半径方向と接線方向に同時に吹けば、合計の大きさは最大 $0.01\sqrt{2}$ になります。軌道半径は途中で最初の軌道より内側にも、目標より外側にも出ないようにします（始点と終点では境界条件で値が決まっているので、途中のノードにだけ課します）。

目的関数

J = \int_0^{t_f} \left( u_r^+ + u_r^- + u_t^+ + u_t^- \right) dt

本当に最小化したいのは $\int (\lvert u_r\rvert + \lvert u_t\rvert)\,dt$ です。推力加速度の大きさを積分すると速度変化の合計になり、これが使う推進剤の量に対応します。しかし絶対値 $\lvert u\rvert$ は $u = 0$ で折れ曲がっていて微分できないので、勾配を使う最適化法には扱いにくい形です。

そこで $u = u^+ - u^-$ と分けて、和 $u^+ + u^-$ を最小化します。こちらは各変数について1次式で、なめらかです。同じ $u$ を表す分け方はいくつもありますが（例えば $u^+ = 0.006,\ u^- = 0.001$ でも $u = 0.005$ ）、和を小さくしたいので、最適解では片方が必ず 0 になり、 $u^+ + u^- = \lvert u\rvert$ が成り立ちます。外向きと内向きに同時に吹くのは推進剤のむだなので、最適化が自然に避けてくれるわけです。

なお、この $J$ は半径方向と接線方向の大きさを別々に足しています。2方向に同時に吹くと、推力の本当の大きさ $\sqrt{u_r^2 + u_t^2}$ より多めに数えます。方向ごとに別のスラスタを持つ宇宙機のモデルと考えることもできます。

パラメータ

記号	意味	値
$\mu$	重力定数	1
$u_{\max}$	推力加速度の各成分の上限	0.01
$r_0,\ v_{r0},\ v_{t0}$	最初の円軌道	1, 0, 1
$r_f,\ v_{rf},\ v_{tf}$	目標の円軌道	4, 0, 0.5
$t_{f,\max}$	遷移時間の上限	55

初期推定

終了時刻の初期推定は 50 とします。

軌道半径 $r$ ：1 から 4 へ直線的に増やします。
半径方向速度 $v_r$ ：0 のまま。
接線方向速度 $v_t$ ：各時刻の $r$ での円軌道速度 $1/\sqrt{r}$ 。
推力：進む向きの $u_t^+$ を上限の半分 0.005、ほかは 0。

「各時刻でほぼ円軌道のまま少しずつ広がる」という、実際の解に近い推定です。数周にわたる長い問題では、ノードの間隔が粗いと渦巻きを表せないので、1区間に 60 ノードを使っています。

結果

軌道面での経路。中央の点が中心の天体、灰色の破線は最初の軌道（r = 1）と目標の軌道（r = 4）です。点はノードでの値、線はそれらを通る多項式、橙色の矢印は推力の向き（間引いて表示し、長さは推力の大きさに比例）です。

時間履歴。軌道半径、接線方向速度（実線）と半径方向速度（破線）、半径方向と接線方向の推力加速度（それぞれ正の成分から負の成分を引いた値）。赤い点線は各成分の上限 ±0.01 です。

遷移の様子。右側は軌道半径と接線方向の推力加速度で、縦線が現在の時刻です。

宇宙機は約 3.8 周回りながら軌道を広げ、時刻 55.0 に目標の円軌道に入ります。終了時刻は上限の 55 に達しています。
目的関数（速度変化の合計）は 0.579 です。
接線方向の推力はほとんどの時間で上限 0.01 です。例外は終盤で、時刻 36〜41 ごろに外向きの推力を加え、 43〜48 ごろには推力をすべて止めて慣性で飛んでいます。
最後の数単位時間は、進む向きに加えて内向きにも吹き、半径方向速度を 0 に戻しながら円軌道に入ります。
接線方向速度は、進む向きに押し続けているのに 1 から 0.5 へ減っていきます。

解説

進む向きに押すと遅くなる

推力の大部分は接線方向、つまり進む向きです。進む向きに押すと軌道のエネルギーが増え、軌道が少しずつ外へ広がります。ところが外側の円軌道ほど回る速さは遅い（ $v_t = 1/\sqrt{r}$ ）ので、広がった結果として速さはむしろ落ちていきます。運動方程式でいえば、推力 $u_t$ による加速を、外へ動くことによる $-v_r v_t/r$ の減速が上回っているのです。

時間の上限と燃料の兼ね合い

小さな推力で時間をかけて回るほど、燃料の点では有利になります。理想的に推力が小さい場合、ゆっくり渦巻き状に広がる遷移の速度変化は $1 - 0.5 = 0.5$ （最初と最後の円軌道速度の差）です。一方、大きな推力を2回だけ瞬間的に吹くホーマン遷移なら約 0.449 で済みます。今回の結果 0.579 はどちらよりも大きくなっています。

理由は終了時刻が上限の 55 に張り付いていることです。0.01 の推力でこの時間内に半径 4 まで届くには、ほぼずっと吹き続けるしかなく、推力を止めて節約する余地がほとんどありません。時間の上限が緩ければ、燃料の効率のよい場所でだけ吹き、それ以外は慣性で飛ぶ解が選ばれます。

長い問題とノード数

約 4 周の渦巻きを1本の多項式で表すため、60 ノードを使っています。1周あたり約 15 ノードです。推力の切り替わりの前後では、ノードの間の線が上限 0.01 を越えたり振動したりしています。制約はノードでしか課していないためで、不連続な制御を1本の多項式で表すときの限界です（区間に分ける）。

試してみよう

PARAMETERS["tf_max"] を 55 から 70 に緩める。手元の計算では、終了時刻は上限より手前の約 65 で決まり、目的関数は約 0.50 まで下がって、推力を止める区間が何か所も現れます。
目標の軌道を変えるときは、rf と vtf を一緒に変える必要があります。例えば rf を 2 にするなら vtf は $1/\sqrt{2} \approx 0.707$ です。周回数と目的関数がどう変わるかを見る。
u_max を 0.02 に上げ、周回数が減ることと、目的関数がホーマン遷移の値（約 0.449）に近づくかを確かめる。

実行する

uv run python examples/15_low_thrust_orbit_transfer.py

examples/15_low_thrust_orbit_transfer.py を GitHub で見る