08. Goddard問題・区間分割

例題07では、推力がゼロへ急に切り替わるところを1本の多項式で表そうとして、推力のグラフが波打ちました。切り替わる時刻で時間を2つに分け、それぞれを別の多項式で表すと、この振動は消えます。物理の問題は例題07とまったく同じで、変わるのは「時間の分け方」だけです。

項目	内容
状態	高度 $h$ 、速度 $v$ 、質量 $m$
制御	推力 $T$
目的	最終高度 $h(t_f)$ を最大にする
区間	2区間（各25ノード）、境目は $t = 0.075$ に固定、終了時刻は自由
新しく学ぶこと	区間に分ける（ノット）ことで、不連続な制御を表す

状況

状況は例題07と同じです。鉛直に上がるロケットを、決まった量の推進剤で最も高く上げます。空気抵抗のため、最適な推力は「最大推力 → 中間の推力（特異弧）→ ゼロ」と変わります。

例題07の答えでは、推進剤を使い切る時刻はおよそ 0.077 でした。そこで今回は、時間を

区間1： $0 \le t \le 0.075$ （推進剤を燃やす区間）
区間2： $0.075 \le t \le t_f$ （推力ゼロで惰性で上がる区間）

の2つに分け、区間ごとに別々の多項式で状態と推力を表します。区間の境目をノット（knot）と呼びます。推力は区間の境目で不連続に変わってかまいません。量はすべて例題07と同じ無次元量です。

問題設定

運動方程式

両方の区間で、例題07と同じ式を使います。

\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(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$ は排気速度です。各項の意味は例題07を参照してください。

区間と接続条件

区間1の終わりの時刻を $t_1$ とし、次の等式で固定します。

t_1 = 0.075

境目では、高度・速度・質量がすべてつながる（連続である）条件を課します。

\begin{aligned} h^{(1)}(t_1) &= h^{(2)}(t_1) \\[2pt] v^{(1)}(t_1) &= v^{(2)}(t_1) \\[2pt] m^{(1)}(t_1) &= m^{(2)}(t_1) \end{aligned}

上付きの $(1), (2)$ は区間の番号です。ロケットの位置や速度が瞬間的に跳ぶことはないので、状態は連続です。一方、推力 $T$ には接続条件をつけません。区間1の最後と区間2の最初で推力が違ってよいので、「吹かしていた推力が急にゼロになる」ことを表せます。

境界条件

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

例題07と同じです。地表で静止した状態から出発し、推進剤を使い切って上昇が止まった瞬間を終端とします。

制約

\begin{aligned} &0 \le T \le 3.5, \qquad h \ge 1, \qquad v \ge 0 \\[2pt] &0.6 \le m \le 1, \qquad t_1 \ge 0.05, \qquad t_f \ge 0.1 \end{aligned}

推力・高度・速度・質量の範囲は例題07と同じです。区間の時刻にも下限をつけています（ $t_1$ は等式で 0.075 に固定しているので、 $t_1 \ge 0.05$ は実際には効きません）。

目的関数

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
$t_1$	区間の境目の時刻	0.075

初期推定

区間の時刻は $[0, 0.1]$ と $[0.1, 0.3]$ から始めます（ $t_1$ は解く過程で 0.075 に移ります）。

高度は全体で 1 から 1.01 へなめらかに変わる3次式、速度は 0。
質量は区間1で 1 から 0.6 へ直線的に減り、区間2では 0.6 のまま。
推力は区間1で最大推力 3.5、区間2で 0。

「区間1で燃やし、区間2で惰性」という答えの形を、初期推定の段階で区間に割り当てています。

結果

高度・速度・質量・推力の時間履歴。青が区間1、橙が区間2です。点は各区間25個のノードでの値、線は区間ごとの多項式、縦の破線は区間の境目 t = 0.075、赤い点線は燃料切れの質量 0.6 と最大推力 3.5 です。

上昇の様子。左の軌跡は区間1を青、区間2を橙で描いています。右は速度と推力の時間変化で、縦の破線が区間の境目、縦の黒線が現在の時刻です。

最高到達高度は 1.0128 で、例題07の 1.0128 とほぼ同じです。
推力は区間1の終わりまで出ていて、区間2では最初からゼロに保たれます。境目で推力が不連続に落ちています。
区間1の終わりの質量は 0.6000 で、推進剤はほぼ区間1で使い切られています。
最高点に達する時刻は 0.200 です。

解説

振動が消えた理由

擬スペクトル法は、区間全体を1本の多項式で表します。多項式はなめらかな関数なので、推力が急にゼロになる「段差」をうまく表せず、その前後で値が波打ちます（例題07）。

区間を分けると、段差は区間1の最後のノードと区間2の最初のノードの間に置かれます。各区間の中では推力はなめらかに変わるので、それぞれ1本の多項式で無理なく表せます。区間2の推力がすべてのノードでほぼ 0 に保たれているのはこのためです。

例題07と比べると、区間1の中の推力の波打ちも小さくなっています。これは推力が「最大推力 → 中間の推力」と変わる部分を、 25ノードの多項式だけで表せば済むようになったためです。

境目の時刻の決め方

この例題では、境目を例題07の結果を見て $t_1 = 0.075$ に固定しました。例題07の燃焼終了時刻 0.079 より少し早いので、推進剤を少し早めに燃やし切ることになり、最終高度は例題07よりわずかに低くなります（差は小数第5位）。区間1の終わりで推力が再び最大推力近くまで上がっているのも、残りの推進剤を 0.075 までに燃やし切ろうとするためと考えられます。

本来、燃焼終了の時刻は答えを見るまで分かりません。等式 $t_1 = 0.075$ を外せば、境目の時刻も最適化の変数として決められます。ただしそのときは、区間1の終わりで推進剤を使い切るといった条件を加えないと、境目が燃焼終了と一致する保証はありません。

試してみよう

境目の時刻の等式（result.equal(ts0, 0.075)）の 0.075 を 0.06 や 0.09 に変える。0.06 では推進剤を早く燃やし切る必要があり最終高度が下がるか、0.09 では区間1の中に推力の段差が入り再び振動が出るかを確かめる。
各区間のノード数（Problem(..., [25, 25], ...)）を [15, 10] などに減らし、例題07で 50 ノード使ったときと精度を比べる。
等式 result.equal(ts0, 0.075) を消し、境目の時刻を自由にすると、どこに落ち着くか、推力がどう変わるかを見る。

実行する

uv run python examples/08_goddard_1knot.py

examples/08_goddard_1knot.py を GitHub で見る