07. Goddard問題・単一区間

鉛直に打ち上げるロケットを、決まった量の推進剤で最も高く上げるには、推力をどう出せばよいでしょうか。 ロケットの父 R. H. Goddard が提起した、最適制御の代表的な問題です。

項目	内容
状態	高度 $h$、速度 $v$、質量 $m$
制御	推力 $T$
目的	最終高度 $h(t_f)$ を最大にする
区間	1区間、50ノード、終了時刻は自由
新しく学ぶこと	質量が減る運動方程式、最大化問題の書き方、不連続な制御を1本の多項式で表したときの限界

状況

ロケットは推進剤を燃やして推力を得ます。推進剤を使うほど機体は軽くなり、同じ推力でもよく加速するようになります。 一方、速く飛ぶと空気抵抗が速度の2乗で増え、推進剤を抵抗に打ち勝つために浪費します。 ただし空気は高度が上がると急に薄くなるので、高いところでは抵抗はあまり問題になりません。

「最初から全開で吹かす」と低空で速度が出すぎて抵抗で損をし、「ゆっくり吹かす」と重力に逆らう時間が長くなって損をします。 最適な推力はこの間のどこかにあり、それは時刻とともに変わります。

計算を簡単にするため、量はすべて地球半径・地表の重力・初期質量を 1 とする無次元量で表します。

問題設定

運動方程式

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

右辺に現れる抵抗 $D$ と重力 $g$ は次のとおりです。

$$D(h, v) = D_c\, v^2 \exp\!\left(-H_c\,\frac{h - h_0}{h_0}\right), \qquad g(h) = g_0 \left(\frac{h_0}{h}\right)^{2}$$

• 抵抗 $D$ は速度の2乗に比例し、空気の密度にあたる指数関数の項によって、高度が上がると急に小さくなります。
• 重力 $g$ は地球中心からの距離の2乗に反比例します。
• 質量は推力に比例して減ります。$c$ は排気速度で、大きいほど同じ推力で推進剤の消費が少なくなります。

境界条件

$$h(0) = 1,\quad v(0) = 0,\quad m(0) = 1, \qquad v(t_f) = 0,\quad m(t_f) = 0.6$$

地表で静止した状態から出発します。終端では推進剤を使い切って（質量 0.6）、上昇が止まった（速度 0）瞬間とします。 その瞬間が最高到達点です。

制約

$$0 \le T \le 3.5, \qquad h \ge 1, \qquad v \ge 0, \qquad 0.6 \le m \le 1$$

推力の上限は初期重量の 3.5 倍です。

目的関数

$$J = -h(t_f)$$

最終高度を最大にしたいので、マイナスをつけて最小化します。

パラメータ

記号	意味	値
$h_0,\ g_0,\ m_0$	初期高度・地表の重力・初期質量	1
$D_c$	抵抗の係数	310
$H_c$	空気が薄くなる速さ	500
$c$	排気速度	0.5
$T_{\max}$	最大推力	3.5
$m_f$	推進剤を使い切ったときの質量	0.6

初期推定

高度は 1 から 1.01 へ、質量は 1 から 0.6 へなめらかに変わる3次式、速度は 0、推力は最大から 0 へ下がる3次式で与えます。 終了時刻の初期推定は 0.3 です。

結果

高度・速度・質量・推力の時間履歴。点は50個のノードでの値、線はそれらを通る多項式です。

上昇の様子。炎の長さは推力の大きさを表します。推力がゼロになった後は2倍速で再生しています。

• 最高到達高度は 1.0128 で、地球半径の約 1.3 % の高さまで上がります。
• 推進剤は時刻 0.079 ごろに尽き、その後は惰性で上昇して時刻 0.199 に最高点に達します。
• 推力は「最大推力 → 中間の推力 → ゼロ」の3段階に変わります。

解説

中間の推力（特異弧）

最初は重力に逆らって早く速度を得るため、最大推力で吹かします。速度が上がって抵抗が大きくなると、 推力をあえて落とし、抵抗と重力の損失がつり合う速度を保つように調整します。この推力が最大でもゼロでもない区間を 特異弧（singular arc）といいます。高度が上がって空気が薄くなると、残りの推進剤を使い切ってから惰性で上昇します。

推力の上下限だけで答えが決まる「バンバン制御」（例題02）と違い、特異弧では推力の値そのものを 物理のつり合いから決める必要があり、間接法では扱いが難しいことで知られます。 直接法では特別な扱いなしに求まります。

推力の振動

推力のグラフを見ると、中間の推力の区間で値が波打ち、燃焼終了の前後でなめらかに（実際より遅く）ゼロへ落ちています。 本当の最適解では推力はゼロへ不連続に切り替わります。50ノードの1本の多項式で不連続な関数を表そうとしたため、 多項式が振動しているのです。擬スペクトル法が不連続を苦手とする典型例です。

例題08 では、燃焼終了の時刻で区間を分けることで、この振動を解消します。

試してみよう

• ノード数（Problem([0.0, 0.3], [50], ...) の 50）を 30 や 80 に変え、推力の振動と最高到達高度がどう変わるか比べる。
• 抵抗の係数 PARAMETERS["Dc"] を小さくすると、特異弧が短くなり、最大推力で吹かす時間が長くなることを確かめる。
• 最大推力 T_max を下げると、推力の上限で吹かす時間が長くなる。

実行する

uv run python examples/07_goddard_0knot.py

examples/07_goddard_0knot.py を GitHub で見る