ハロー軌道の特性を利用した宇宙機の軌道維持システムの低燃費化
M2014SC020 宇佐美元啓
指導教員：大石泰章

1

はじめに

本研究では月，地球と，宇宙機による三体問題を扱う．
月や地球の質量に対して宇宙機の質量は無視できるほど微小であり，宇宙機は月と地球により生成される重力場内を運動していると考えることができる．さらに，地球と月の運動を，それらの共通重心を中心とした円運動であると仮定し，いわゆる円制限三体問題を考える．円制限三体問題における重力場の力学的平衡点はラグランジュ点と呼ばれ，L1 から L5 まで存在することが知られており，これらの点の周りをまわる周期軌道をハロー軌道と呼ぶ．
本研究で用いるハロー軌道は，地球から見て月の裏側に位置する L2 点を平衡点とする最大径約 70,000[km] の軌道であり，このハロー軌道に沿って宇宙機が運動するよ
図 1 地球，月と宇宙機による三体問題
うな制御を考える．軌道上の宇宙機には，月の裏側の観測や L2 点に配置した宇宙港と地球との通信などの役割が期待されている [1]．しかし，ハロー軌道は不安定周期軌 本研究では，L2 = (l2 , 0, 0) の周りのハロー軌道を考える道であり，制御なしではこの軌道に沿った運動を維持で ため，原点を L2 に移している．さらに，式 (1) は月と地きない．
球の重心間距離 D0 = 3.847480 × 105 [km] を 1，月の公
本稿では，最適レギュレータ理論に基づくフィードバッ 転周期 n = 2.661699 × 10−6 [rad/s] を 2π として無次元化クゲインの設計法を比較する．1 つ目の設計法は，L2 点で されている．
「′ 」は無次元化時間 τ による微分である．
線形化したシステム行列を用いたゲインの設計法であり，これについては既に [2] で考察している．2 つ目は，軌道1−ρ上で線形化した時変システム行列を用いた周期ゲインのx′′ − 2y ′ − x = l2 − 3 (x + l2 + ρ)r
設計法，3 つ目はフロケの理論におけるモノドロミー行列
ρ e− 3 (x + l2 − 1 + ρ) + ux ,により，軌道の発散特性を解析し，発散度合いを反映さrせた周期ゲインの設計法である．本稿では連続入力の場
ρ1−ρ
(1)y ′′ + 2x′ − y = − 3 y − 3 y + uy ,合について上記 3 つの設計法を検討し，さらに，宇宙機
rerの軌道制御法としてより一般的な，パルス入力を用いた
ρ1−ρz ′′ = − 3 z − 3 z + uz .
場合への応用も行う．制御性能の評価には燃料消費に比rer
例する速度変化の総和を用い，最適レギュレータの入力重みを変化させることで，燃料消費の最も少ないフィードバックゲインを求める．

32

軌道方程式

ハロー軌道上の時変システムと特性解析

式 (1) の周期解を解析的に得ることはできないが，適切
図 1 に地球，月，宇宙機からなる系を示す．円制限三 な初期値を与えることで数値的に周期軌道を得ることがで体問題では，地球と月はそれらの共通重心を中心として， きる．この周期軌道をハロー軌道と呼ぶ．ただし，ハロー円運動をしている．ここで，座標の原点を，二天体の共 軌道は不安定周期軌道であり，軌道上の位置によって発散
通重心 Ob とする．慣性座標を Ob -(X, Y, Z) とし，特に 度合いが変化する．本章では，軌道上で線形化した時変シ
X-Y 平面は二天体の運動面であり，Z 軸は運動面に垂直 ステムを求め，フロケの理論を用いて，ハロー軌道上の各であるようにする．また，月と地球の運動にあわせて回 時刻における発散度合いを調べる．本研究におけるハロー−20 0 1.770374×転する座標を Ob -(x̂, ŷ, ẑ) とし，その原点は同様に Ob と 軌道の初期値は xf0 = [−3.563945×10−1−3T− 6.900000 × 100] であり，周期は T = 3.4150
し，図 1 のように x̂ 軸，ŷ 軸をとる．x̂ 軸は，地球から 10見た月の方向であり，地球から観測できるのは， ŷ-ẑ 平 である．
面である．回転座標の角速度は，月の公転角速度であり， 3.1 変分方程式と時変システム行列n[rad/s] とする．ここで，Me は地球の質量，M は月のはじめに，ハロー軌道近傍の微小なずれの変化を表す質量，m は宇宙機の質量である．また，二天体の総質量変分方程式を，軌道上で線形化することで，ハロー軌に対する月の質量比率を ρ とする．
系の力学的平衡点のうちの一つである L2 点の座標を 道の時変システム行列を求める．ハロー軌道を，xf =′ ′′ T
(l2 , 0, 0)，L2 点を基準とした宇宙機の相対座標を (x, y, z)， [xf yf xf yf zf zf ] とかくと，式 (1) は，T制御加速度を u = [ux uy uz ] とすることで，円制限三d体問題の運動方程式は，式 (1) のように表される [4]．た
xf = f (xf )(2)dτ2
だし，式 (1) の re , r はそれぞれ，re = [(x + l2 + ρ) +y 2 + z 2 ]1/2 , r = [(x + l2 − 1 + ρ)2 + y 2 + z 2 ]1/2 である．

と表すことができる．ここで，f (xf ) は，
x′fyf′

モノドロミー行列の固有値は，周期系の安定性解析において有用である [5]．モノドロミー行列の固有値のうち，複素平面上の単位円内に存在するものは安定であり，単位円外に存在する固有値は不安定である．また，モノドロミー行列の固有値は，初期時刻 τ0 に依存せず一定であり，固有ベクトルのみが τ0 により変化する [6]．ハ
ロー軌道における，初期時刻を τ としたモノドロミー行列の固有値を，λ1 (τ ), λ2 (τ ), ..., λ6 (τ )，それらに対応する固有ベクトルを，v1 (τ ), v2 (τ ), ..., v6 (τ ) と定める．ただし|λ1 (τ )| ≧ |λ2 (τ )| ≧ |λ3 (τ )| ≧ |λ4 (τ )| ≧ |λ5 (τ )| ≧ |λ6 (τ )|とする．図 2 は各時刻 τ における固有値の絶対値を示す．
図 2 より，ハロー軌道では，1 つの安定固有値，1 つの不

 2y ′ + xf + l2 − 1−ρ(x+l+ρ)
f2re3 fρ−(x+l−1+ρ)+uf2xr3
f (xf ) = 
 −2x′ + y − 1−ρ y − ρ y + u fy f
re3 fr3 f
zf′1−ρρ− r3 zf − r3 zf + uzeである．軌道 xf からの微小なずれを ∆x とすると，ハロー軌道の近傍点は，x = xf + ∆x と表すことができる．
よって，近傍点のふるまいは，d(xf + ∆x) = f (xf + ∆x)dτ(3)となり，式 (3) から式 (2) を引くことにより，次の式が得られる：

d(∆x) = f (xf + ∆x) − f (xf ).
dτ(4)式 (4) は微少誤差のふるまいであり，ハロー軌道上の点xf の周りのテイラー展開により，線形近似を行うと，

d∂f(∆x) ≃∆x = A(τ )∆xdτ
∂x x=xf (τ )となる．ここで，∆x の係数行列 A(τ ) は時変システム行列であり，次のように表される：

001 0 0 0
 000 1 0 0  ∂f3
 3 ∂f30 2 ∂f0  ∂xf ∂yf∂zf
A(τ ) =  ∂f4 ∂f4.
∂f ∂xf ∂yf −2 0 ∂z4f 0 
 000 0 0 1 ∂f6∂f6
60 0 ∂f0∂xf∂yf∂zf3.2

図 2 モノドロミー行列の固有値推移
安定固有値および，4 つの中心固有値があることが見て取れる．前段落で述べた通り，時間 τ によらず，固有値の大きさはほぼ一定である．最終時刻付近で固有値の変化が見られるのは，シミュレーションで用いたハロー軌
道が完全な周期軌道ではないことによるものであり，初期値の精度がよいほど，変化は小さいと考えられる．
唯一の不安定固有値 λ1 (τ ) は実数なので，対応する固有ベクトル v1 (τ ) は，実ベクトルに選べる．次に，各基本ベクトル ϵi (ϵi は i 番目の成分のみが 1 で，他は零であるような 6 次元ベクトル) が v1 (τ ) 方向の成分をどの程度持つかを調べる．ϵi を v1 (τ ), v2 (τ ), ..., v6 (τ ) の 1 次結合

ϵi = α1i (τ )v1 (τ ) + α2i (τ )v2 (τ ) + ... + α6i (τ )v6 (τ )で表す．α1i (τ ) は，必ず実数になる．また，α1i (τ ) は基本

フロケの理論を用いた軌道特性解析

前節で求めた A(τ ) を用いて，次の周期システムを考える:

x′ (τ ) = A(τ )x(τ ).

(5)

式 (5) の解は，時刻 τ0 における状態 x(τ0 ) を使って，以下のように書ける．

x(τ ) = ΦA (τ, τ0 )x(τ0 ).
ここで，ΦA (τ, τ0 ) は状態遷移行列であり，Φ′A (τ, τ0 ) = A(τ )ΦA (τ, τ0 ), ΦA (τ0 , τ0 ) = Iを満たす．特に，1 周期後の状態遷移行列は τ0 におけるモノドロミー行列と呼ばれる:ΨA (τ0 ) = ΦA (τ0 + T, τ0 )．

図 3 軌道特性の変化ベクトル ϵi の v1 (τ ) 方向の成分であり，ϵi の発散度合いを表すと考えられる．図 3 に軌道 1 周期分の，|α11 |, ..., |α16 |の計算結果を示す．この場合，|α11 | が大きいので，ϵ1 =[1 0 0 0 0 0]T の方向が最も発散しやすく，特に初期値付近と最終値付近の発散度合いが大きい，といったことが理解できる．

4制御系設計4.1

状態空間表現

原点 L2 =(0, 0, 0) において式 (1) を線形化し，状態変数を x = [x y x′ y ′ z z ′ ]T ，制御入力を u = [ux uy uz ]Tとすることで，以下の状態空間表現を得る．

x′ = Ax + Bu(6)ここで行列 A，B は，
00
00 2σ + 10
A=01−σ
0000100−2000 01 02 00 00 00 −σ00 00 0 ,B =  1
 00 
 01 00000100

00 0 0 
0 1である．ただし σ = ρ/(l2 − 1 + ρ) + (1 − ρ)/(l2 + ρ)である．また，非線形方程式 (1) は，33x′ = Ax + Bh(x) + Buと表せる．ここで，h(x) は，非線形項であり，次のようになる [7]：l2 − 2σx − 1−ρ
re3 (x + l2 + ρ)− rρ3 (x + l2 − 1 + ρ)h(x) = .
1−ρρσy − r3 y − r3 yeρσz − 1−ρr3 z − r3 ze

4.2

誤差方程式の安定化

本研究では，初期値 x0 = [x0 y0 x′0 y0′ z0 z0′ ]T から始まる軌道 x を，入力 u を適切に与えることによって，目標軌道 xf に漸近させることを考える．ここで，xf は 3 章′′ T
で用いた初期値 xf0 = [xf0 yf0 x′f0 yf0zf0 zf0] から始まるハロー軌道である．ここで，軌道誤差を e = x − xf とかくと，()e′ = Ae + B h(x) − h(xf ) +Buを満たす．フィードバックゲイン K を用いて，制御入力をu = −Ke − h(x) + h(xf )(7)と定め，K を適切に選ぶことで，軌道誤差 e を零に収束させることができる．

4.3連続入力の場合のゲイン設計法

本節では，文献 [3] で考察した連続入力の場合の制御法
を 3 つ紹介する．
はじめに，原点で線形化した状態方程式 (6) を用いたフィードバックゲインの設計法を考える．設計法は，最適レギュレータ理論に基づく．評価関数を，次のように定める:∫ ∞J(u; x) =(uT Ru + xT Qx)dτ.
(8)の安定化解 P を用いて，K = R−1 B T P となる．ただし，
式 (9) における行列 A は，原点 L2 において線形化したシステム行列であり，軌道上での動特性やその時間変化が考慮されていない．
上記の問題を解決するために，2 番目の設計法として，軌道上で線形化した時変システム行列 A(τ ) を用いた方法を考える．ハロー軌道の周期は T であるから，A(τ +T ) =A(τ ) である．よって，評価関数を最小化する時変フィード
バックゲイン K(τ ) も，周期 T の周期関数となる．K(τ )は，微分リッカチ方程式（Riccati Diﬀerential Equation,RDE）−P ′ (τ ) =AT (τ )P (τ ) + P (τ )A(τ )+ Q − P (τ )BR−1 B T P (τ )の周期解 P (τ ) を用いることで K(τ ) = R−1 B T P (τ ) として得られる．周期解の数値計算法は文献 [3] で示されている．
3 番目の設計法として，3 章で求めた発散度合いを，レギュレータの状態の重みとして利用し，発散度合いに対応した周期ゲインを設計する方法を考える．3 章で求めた各基本ベクトル ϵi の発散度合い α1i (τ ) を用いて，発散度合い行列，

Γ(τ ) = diag[|α11 |, |α12 |, |α13 |, |α14 |, |α15 |, |α16 |] ∈ R6×6を定義する．この Γ(τ ) により，評価関数 (8) の Q を置き換える．これにより，軌道位置によって変化する発散度合いを反映した時変ゲインの設計ができ，発散度合いの大きい（小さい）位置に宇宙機がある場合は，同じ誤差であっても，より大きな（小さな）制御入力を入れることができるようになる．評価関数を最小化する時変フィー
ドバックゲイン K(τ ) の導出方法は，2 番目の方法と同様であり，文献 [3] で示されている．

4.4

パルス制御への応用

本節では，連続入力の場合に紹介した 3 つの制御方法を，宇宙機の制御法としてより一般的な，パルス制御に拡張する．本研究では，パルス入力の幅を w に固定し，1周期に k 回入力することを考える．ハロー軌道の周期はT であるから，パルスの周期は T /k であり w < T /k をみたすように w, k を決めればよい．
ゲインを設計するためには，連続システムを用いて最適化したゲイン K をデューティー比で割る．すなわち，連続システムの場合で紹介した 3 つの方法でゲインを求め，パルス制御におけるゲイン Kd を以下のように定めた：

Kd = K ×T.
wk

(10)

周期ゲイン K(τ ) の場合は，パルス開始時刻 τ0 のゲインを用いて K = K(τ0 ) とし，式 (10) に代入して得られた
Kd を時刻 τ0 + w まで用いる．その後，時刻 τ0 + T /k までは入力を零とする．また，パルス制御を行う場合は，非線形項 −h(x) + h(xf ) を無視する．

0

ここで Q は状態にかかる重み，R は入力にかかる重みである．評価関数を最小化するフィードバックゲインは，代数リッカチ方程式（Algebraic Riccati Equation, ARE）

0 = AT P + P A + Q − P BR−1 B T P(9)5制御性能比較

本章では，宇宙機をハロー軌道上に維持した場合のシミュレーション結果と各制御器を用いた場合の性能評価の結果を示す．

5.1

シミュレーション条件

シミュレーションにおける数値計算の刻み幅は 10−3 とした．また，ハロー軌道の初期値は 3 章で用いた xf0 と同じである．ただし，初期値が与える軌道は完全な周期軌道ではないため，１周期後の値に微小なずれが発生し，最終値を２周期目の初期値として使用できない．そこで，シミュレーションにおいては，ハロー軌道の周期 T ごとに初期値を xf0 に戻し，１周期分の軌道を繰り返し目標軌道として使用する．また，ハロー軌道の初期値と宇宙機の初期位置は同じ値としたため，１周期目は自由運動となり，入力は 0 である．そこで，2 周期目以降の 100 周
図 5 r の変化に伴う ∆V の推移期にかかった燃料消費の，１周期平均を結果として示す．
パルス入力では，パルス幅を w = 50 (1.8596 × 104 [sec]),した．連続入力の場合と同様に，3 種の制御器において，1 周期の入力回数を k = 30 とした．
それぞれ ∆V が最小となるときの r が，その制御器における最適な r である．図 5 より，時変システム行列 A(τ )5.2 燃料消費による評価法を用いることで燃料消費が抑えられることが確認できる．
フィードバックゲインの性能評価を行い，燃料消費のしかし，パルス入力の場合，時変重み Γ(τ ) を用いること少ないシステムを設計する．連続制御の評価は，加速度による低燃費化の効果は確認できない．の絶対積分である L1 ノルムを用いる．L1 ノルムは，次また，図 4，図 5 より連続入力よりもパルス入力を用いのように表わされる:
た場合の方が，燃料消費量が大きくなってしまうことが∫ T見て取れる．これは，パルス入力間の自由運動中に，誤(ux 2 + uy 2 + uz 2 )1/2 dτ.
JL1 =差が拡大してしまうことにより，より多くの制御入力が0必要になることや，入力において非線形項の考慮をして一方，パルス制御の場合は，加速度の総和である ∆V を いないためであると考えられる．
用いる．∆V は，次のように表わされる:

6J∆V =k∑(uxj 2 + uyj 2 + uzj 2 )1/2 w.

j=1

いずれも燃料消費と比例関係にあり，これを小さくするゲインを求めることで，燃料消費の少ないシステムを設計する．シミュレーションでは，評価関数 (8) の入力重み行列を R = 10r × I3 とおき，r を変化させることで L1ノルムおよび ∆V の変化を調べた．

5.3

シミュレーション結果

本節では，各制御器の燃料消費の比較を行う．ただし，値は全て有次元量に戻している．連続入力を用いて，ハロー軌道に宇宙機を 1 周期維持した場合，L1 ノルムは rの変化に伴い，図 4 のように推移した．3 種の制御器におわりに

本稿では，連続入力とパルス入力を用いてハロー軌道の維持制御を行った．連続入力の場合は，軌道の発散特性を利用することで最も燃料消費が抑えられた．一方，パルス入力では時変システム行列 A(τ ) を用いることで燃料消費を抑えられたものの，軌道特性を時変重み Γ(τ ) として用いる方法は，必ずしも効果的ではなかった．

参考文献[1] M. Utashima, Orbital Mechanics Near Lagrange’sPoints, NASDA-TMR-960033, (1997)[2] T. Ikeyama, M. Usami, A. Ichikawa, Halo orbit maintenance in earth-moon circular-restricted
three-body problem, 24th Workshop on Astrodynamics and Flight Mechanics, (2014)[3] 宇佐美元啓， 坂井祐介， 大石泰章， 地球-月系円制
限 3 体問題におけるハロー軌道の維持制御，第 58 回
自動制御連合講演会予稿集，(2015)[4] B. Wie, Space Vehicle Dynamics and Control,AIAA, Reston, (1998)[5] C. Simo, G. Gomez, J. Llibre, R. Martinez, J. Rodriguez, On the optimal station keeping controlof halo orbits, Acta Astronautica, Vol.15, No.6/7,pp.391–397, (1987)

図 4 r の変化に伴う L1 ノルムの推移[6] S. Bittanti, P. Colaneri, Periodic Systems Filteringand Control, Springer, London, (2009)

おいて，それぞれ L1 ノルムが最小となるときの r が，その制御器における最適な r である．図 4 より，時変シス
テム行列 A(τ ) と時変重み Γ(τ ) を用いることで，最も燃料消費を抑えられることが示された．
パルス入力を用いて，ハロー軌道に宇宙機を 1 周期維持した場合，∆V は r の変化に伴い，図 5 のように推移[7] M. Bando, A. Ichikawa, Formation flying along
halo orbit of circular-restricted three-body problem, Journal of Guidance, Control, and Dynamics,
Vol. 38, No. 1, pp.123–129, (2015)