低推力スラスタを用いた静止衛星の軌道制御

三菱電機 北村憲司

第13回宇宙科学奨励賞 記念公演
2021.04.27
1

今回の受賞の研究概要
研究の概要
静止軌道遷移●最小エネルギー軌道遷移の解析
平均化法や正準変換を用いた同一平面内のエネル
ギー最小遷移軌道の解析●最短時間軌道遷移のためのフィードバック制御則
アクチュエータ制約や姿勢制約を考慮した最短時間
軌道遷移のためのフィードバック制御則の提案

静止軌道保持GEO●限定された本数のスラスタによる高精度な軌道保持方式，解析的な平均軌道要素の活用．
GEO

2

静止衛星
1rev/day

36000km

1rev/day

Geostationary orbit(GEO)

地表から高度約36000km上空の衛星．軌道周期が地球の自転周期(約24h)と等しいため，地上の観測者からみると上空で静止しているように見える．
実用上とても価値の高い軌道• 気象衛星（ひまわり8/9号など)• 通信衛星（スカパーJSATなど）• 測位衛星（準天頂衛星など）• 観測衛星（GISATなど)[1] http://www.mitsubishielectric.co.jp/society/space/satellite/observation/himawari8-9.html[2] https://qzss.go.jp/overview/index.html

ひまわり8/9号[1]

みちびき

3

静止衛星の運用
静止衛星の運用は大きく2種類に大別される：
1週間～半年10～15年

軌道遷移フェーズ

軌道保持フェーズ
（ミッションフェーズ）

打ち上げ

GEO到着

燃料枯渇
廃棄

ロケット投入軌道
軌道保持ボックス（0.05°程度）

北極
North

East
赤道面
グリニッジ子午線方向

GEO

±𝜙𝑚𝑎𝑥

West
±𝜆𝑚𝑎𝑥

4

軌道遷移のアニメーション

5軌道遷移方式の比較
化学推進による遷移

電気推進による遷移

GEO

軌道GTO

Δ𝒗

(Geostationarytransfer orbit)

連続制御

インパルス制御
遷移日数

0.5日（実際は1～2週間)

数ヶ月～半年

速度増分

1500 m/s

2300 m/s

有効排気速度

2-4 km/s

15-30 km/s

推薬消費量10.15-0.30

運用負荷小大6低推力軌道遷移の課題
計算負荷の小さい準最適な軌道制御則の構築
今後の衛星の運用負荷の低減と自律化をにらみ，オンボードでも実行可能な，計算負荷の小さく，準最適な軌道制御則の構築を目指す．
オンボードでの制御量算出
緊急時のみ
制御コマンド送信

制御コマンド最低限のテレメトリ

詳細なテレメトリ
データ

地上局

••
•
•

多くの地上運用者
高い運用負荷・コスト
長期にわたる制御計画立案
豊富な計算リソース

地上局

••
•
•

最低限のモニタリング
運用負荷・コストの削減
衛星の自律化
少ない計算リソース(オンボード)

7

低推力遷移軌道の最適化手法の課題
数値的手法

間接法や直接法などの手法を使って最適制御問題を数値的に解く．
計算コストが大きく，拘束条件が増えると局所解に陥りやすい
間接法：Sackett [‘75], Krier [‘15]Direct Collocation: Kluever [‘98], Graham [‘15]微分動的計画法：Aziz [‘18]

計算負荷が高いため，地上計算が不可欠．衛星搭載機でのオンボードの計算には不適．
衛星の自律化を目指し，オンボードで計算可能な遷移手法について検討する．
解析的手法

平均化法や最適制御理論を用いて解析的に軌道を求める．
任意の境界条件，評価関数について解を求めることは困難
離心率をゼロに固定：Edelbaum [‘62], Kechichian [‘97, ‘03]
推力方向を軌道接線方向に固定：Bonnard [‘06]より実務的に価値のある楕円軌道(GTO)から円軌道(GEO)への2D-エネルギー最小遷移を解析的に取扱う半解析的手法問題の構造を理解した上で，あるルールに従って推力方向を決定する．
Lyapunovフィードバック則: Petropoulos [‘03]軌道要素の誤差に応じて推力方向を決定: Ruggierro [‘11]

簡易かつ準最適な3D-最短時間遷移問題のためのFB則を提案
姿勢制約，アクチュエータ制約等も考慮(3次元)8並進のダイナミクス(惑星方程式)
Modified Equinoctial 軌道要素 𝑝, 𝑒𝑥 , 𝑒𝑦 , ℎ, 𝑘, 𝐿 [1]で衛星の状態を表現
𝑑𝑝 2𝑝 𝑝=𝑎𝑑𝑡𝑤 𝜇 𝑑𝑠

𝑹ሷ = −𝜇𝑑𝑒𝑥 1 𝑝=[+𝑤 sin 𝐿 𝑎𝑑𝑟 +𝑑𝑡𝑤 𝜇𝑤 + 1 cos 𝐿 + 𝑒𝑥 𝑎𝑑𝑠 −𝑒𝑦 ℎ sin 𝐿 − 𝑘 cos 𝐿 𝑎𝑑𝑛 ],𝑑𝑒𝑦 1 𝑝=[−𝑤 cos 𝐿 𝑎𝑑𝑟 +𝑑𝑡𝑤 𝜇𝑤 + 1 sin 𝐿 + 𝑒𝑦 𝑎𝑑𝑠 +𝑒𝑥 ℎ sin 𝐿 − 𝑘 cos 𝐿 𝑎𝑑𝑛 ],

𝑹+ 𝒂𝑑
𝑅3

面内運動𝑑𝐿𝑤= 𝜇𝑝
𝑑𝑡𝑝2+1 𝑝ℎ sin 𝐿 − 𝑘 cos 𝐿 𝑎𝑑𝑛 ,𝑤 𝜇𝑑ℎ 𝑠 2 𝑝=cos 𝐿 𝑎𝑑𝑛 ,𝑑𝑡 2𝑤 𝜇𝑁

面外運動𝑑𝑘𝑠2 𝑝=sin 𝐿 𝑎𝑑𝑛 ,𝑑𝑡 2𝑤 𝜇𝑅[1] Walker, M. J. H., Ireland, B., and Owens, J. “A Set Modified Equinoctial Orbit Elements.” Celestial Mechanics, Vol. 36, No. 4, 1985, pp. 409–419

𝑆𝒂𝑑9最適制御問題のハミルトニアン
Cost Function𝑡𝑓

1𝐽 = න 𝒂𝑇 𝒂𝑑𝑡
2 𝑡0Dynamics𝒙ሶ = 𝐹 𝒙, 𝐿 𝒂
𝐿ሶ = 𝑔0 𝐿, 𝒙

Boundary Conditions
𝑡0 , 𝑡𝑓 : Fix
𝒙 𝑡0 = 𝒙0 = 𝑝0

𝑒𝑥0

𝑒𝑦0 𝑇 (Fix)𝑝𝑓 = 𝑎𝑓 = 𝑎𝐺𝐸𝑂 , 𝑒𝑥𝑓 = 𝑒𝑦𝑓 = 0 (Fix)𝒂 = 𝑎𝑟𝑎𝜃 𝑇

𝒙 = 𝑝, 𝑒𝑥 , 𝑒𝑦𝑇

𝐿 𝑡0 = 𝐿0 Fix ,GTO
GEO𝐿 𝑡𝑓 = Fix

 最適制御問題のハミルトニアン:

1
𝐻 = 𝝀𝑥 𝐹 𝒙, 𝐿 𝒂 + 𝜆𝐿 𝑔0 𝒙, 𝐿 − 𝒂𝑇 𝒂2
 最適性原理:

𝜕𝐻=𝟎
𝜕𝒂𝒂 = 𝑭 𝐿, 𝒙 𝑇 𝝀𝑥

∴ 最適な遷移軌道に沿って，ハミルトニアンは以下を満たす:
1𝐻 𝑡, 𝒙, 𝐿, 𝝀𝑥 , 𝜆𝐿 = 𝝀𝑇𝑥 𝐹𝐹 𝑇 𝝀𝑥 + 𝜆𝐿 𝑔0 .
210

ハミルトニアンの時間平均問題を解析的に取り扱うために，ハミルトニアンを軌道1周回で平均化する[1]：ഥ 𝑡, 𝒙, 𝐿, 𝝀𝑥 , 𝜆𝐿𝐻

1 𝑡+𝑇= න 𝐻 𝑡, 𝒙, 𝐿, 𝝀𝑥 , 𝜆𝐿 𝑑𝑡
𝑇 𝑡
1 2𝜋𝑑𝑡= න 𝐻 𝑡, 𝒙, 𝐿, 𝝀𝑥 , 𝜆𝐿𝑑𝐿
𝑇 0𝑑𝐿
1 2𝜋 𝝀𝑇𝑥 𝐹 𝒙, 𝐿 𝐹 𝒙, 𝐿 𝑇 𝝀𝑥
= න+ 𝜆𝐿 𝑑𝐿
𝑇 02𝑔0 𝒙, 𝐿ഥ 𝐻は以下のように陽に書き下すことができる．

ഥ 𝑡, 𝒙, 𝐿, 𝝀𝑥 , 𝜆𝐿 =𝐻𝑛ത

5𝑝ҧ 22 −3 + 5𝜁 −2
ҧ 𝜆ҧ2𝑝 + 𝑒𝑦ҧ 2 + 5𝜁 2ҧ 𝜆ҧ2𝑒𝑥
ൣ4𝑝
ҧ
𝜇 3 4𝜁 5ҧ

[1] B. Bonnard, and J. B. Caillau, Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust, Elsevier, Science Direct,
Vol. 24, 2007, p. 395-411.

11平均化されたハミルトニアンの正準変換
ഥ=𝐻𝑛ത

5𝑝ҧ 2𝜇ҧ
3 4𝜁 5ҧ 𝜆2ҧ𝑝 + 𝑒𝑦ҧ 2 + 5𝜁 2ҧ 𝜆2ҧ𝑒𝑥ൣ4𝑝ҧ 2 −3 + 5𝜁 −21𝜇 3𝑝ҧ = 21 − 𝑒2 ,𝑛ഥ=𝐻𝑒𝑥ҧ = 𝑒 cos Θ ,𝑒𝑦ҧ = 𝑒 sin Θ215−(𝜇𝑛) 3 18𝑛2 𝜆2𝑛 + 5 1 − 𝑒 2 𝜆2𝑒 + 2 − 4 𝜆2Θ + 𝑛𝜆𝐿
4𝑒𝑟=ഥ = 𝜅 𝜆2𝑟 +𝐻

12(𝜏𝑛)3 ,𝜙=2/5 sin−1 𝑒1 2𝑟3𝜅𝜁 𝜙 2𝜆+𝜆+𝜆Θ
𝐿2𝑟2 𝜙𝑟2 2𝜏𝜉 = 𝑟 cos 𝜙 ,ഥ = 𝜅 𝜆2𝜉 + 𝜆2𝜂 +𝐻

𝜂 = 𝑟 sin 𝜙32𝜅𝜁 𝜉, 𝜂 2𝜆𝐿 𝜉 2 + 𝜂 2 2 + 2𝜆4𝜏𝜉 + 𝜂2 Θ12平均化された正準方程式
ഥ=𝜅
𝐻

𝜆2𝑟 +1 2𝑟3𝜅𝜁 𝜙 2𝜆+𝜆+𝜆Θ𝐿2𝑟2 𝜙𝑟2 2𝜏ഥ𝑑𝑟 𝜕𝐻== 2𝜅𝜆𝑟𝑑𝑡 𝜕𝜆𝑟

ഥ 2𝜅𝜆2𝜙𝑑𝜆𝑟
𝜕𝐻3 2=−= 3 −𝑟 𝜆𝐿𝑑𝑡𝜕𝑟𝑟2 2ഥ2𝜅𝜆𝜙𝑑𝜙
𝜕𝐻== 2𝑑𝑡 𝜕𝜆𝜙𝑟

ഥ𝑑𝜆𝜙
𝜕𝐻=−=0𝑑𝑡
𝜕𝜙𝜆𝜙 = const.

ഥ𝑑𝐿 𝜕 𝐻
1==𝑟3𝑑𝑡 𝜕𝜆𝐿 2 2𝜏ഥ𝑑𝜆𝐿
𝜕𝐻=−=0𝑑𝑡
𝜕𝐿𝜆𝐿 = const.

ഥ𝑑Θ 𝜕𝐻2𝜅𝜁 𝜙==𝜆Θ𝑑𝑡 𝜕𝜆Θ𝑟2ഥ𝑑𝜆Θ𝜕𝐻=−=0𝑑𝑡𝜕Θ𝜆Θ = const.

∴3 2𝜅𝜆𝐿2𝜏𝑟ሷ = 2𝜅𝜆ሶ 𝑟 = 𝑟𝜙ሶ 2 − 𝑐𝑟 2𝑐≡𝑟 2 𝜙ሶ = ℎ = constℎ ≡ 2𝜅𝜆𝜙距離の2乗に比例する中心力場 𝑓 𝑟 = −𝑐𝑟 2 での運動と等価になる
13境界値問題の解釈
距離の2乗に比例する中心力場において，始点と終点を，指定された時間と周回数で結ぶように，2点間の軌道と中心力の重力定数𝑐を決定せよ．
（距離の2乗に比例する中心力場におけるLambert問題）

中心力が距離だけに依存するので，運動は求積法によって可解𝜂GTO

斥力 𝑐 < 0→周回数小(𝑐 = 0)

𝑓 = −𝑐𝑟 2引力(𝑐 > 0)→周回数大𝑟(𝑡)

𝜙(𝑡)

𝜉GEO
14

軌道に沿ったポテンシャルエネルギー

GEO

𝜇𝑈 𝑅 =−𝑅GTOGEO

𝑐 3𝑢 𝑟 = 𝑟 (𝑐 < 0)
3𝑌𝑋GTO𝜂𝜉
15

周回数を変化させたときの軌道𝜂GTO斥力(𝑐 < 0)𝑓 = −𝑐𝑟2(𝑐 = 0)引力(𝑐 > 0)

GEO

GTO

60 day, 52 Rev
𝐽 = 0.687 [m2 /s 3 ]GTO

60 day, 87 Rev
𝐽 = 0.300 [m2 /s 3 ]𝜉GTO

60 day, 113 Rev
𝐽 = 0.416[m2 /s 3 ]16解析解と数値解の比較
GTO→GEOのエネルギー最小遷移(60day, 58Rev)について，解析解と数値最適解(疑似スペクトル法)の結果を比較．両者は良く一致．

解析解

数値解 (疑似スペクトル法)

17解析解の活用例(1)エネルギー最小問題の解析解を他の問題の初期解として活用することができる．
たとえば，以下のようなホモトピー法的アプローチにより，エネルギー最小解からΔV最小解へ移行することが可能．
𝑡𝑓𝐽=න

𝒂 2−𝛾 𝑑𝑡

𝛾: 0(エネルギー最小) → 1(ΔV最小)𝑡0

エネルギー最小ΔV最小[mm/s 2 ]コースティング𝑎 の上限は0.5 mm/s2

18解析解の活用例(2)最適な制御加速度は，以下の状態フィードバックの形で表現することができる．

 最適性原理:

𝜕𝐻=𝟎
𝜕𝒂𝒂 = 𝑭 𝐿, 𝒙 𝑇 𝝀𝑥
衛星の状態(軌道要素)

オンボードで随伴変数𝝀𝑥 を定期的に更新することで，力学モデルに誤差があったとしても，目標軌道(静止軌道)へ漸近させることが可能．

GTO

𝝀𝑥 を更新しない場合の軌道𝝀𝑥 の更新

GEO
GEO

Time22

おわりに• 実務的に意義のある軌道設定(GTOからGEOへの遷移)について，制御理論・力学系理論(正準変換，最適制御，平均化法)を応用し，エネルギーコストに関して最適な解を導出した．
• 正準変数のとり方を工夫し，エネルギー最小問題を解析的に算出することを可能とした．
→ 搭載計算ソフトウェアでの状態フィードバック制御の可能性を示した．
低推力かつ多周回の軌道遷移は理論も実運用も複雑となりますが，非常に面白い問題構造をもち，今後も活発な研究がなされる分野と考えます．
電気推進によってペイロード輸送効率が高まることにより，商用衛星においては通信の大容量化が進み，ユーザエクスペリエンスの向上が期待されます．また探査ミッションでもミッションペイロードの拡充によって探査の質と量が高まり，新しい発見へとつながることを期待します．
研究を進めるにあたって，(元)大阪大学の山田克彦教授をはじめ，三菱電機先端技術総合研究所，鎌倉製作所の皆様には貴重な助言を数多く頂きました．
この場を借りてお礼申し上げます．
20