インパルス制御によるハロー軌道への移行とその維持
2009SE066 池山拓弥

2009SE103 加藤圭

指導教員：市川朗

1

はじめに現在，数多くの宇宙探査機が打ち上げられ，それらによ

ここで，変数は表 1 のように定義する．

る様々なミッションが実行されている．例えば，1995 年末
表 1 変数定義には ESA/NASA により，太陽表面・コロナ・太陽風を連続的に観測することを目的とした太陽観測衛星 SOHO(Solor変数

and Heliospheric Observatory) が地球・太陽系のラグランジュ点 (L1 点) 近傍のハロー軌道に打ち上げられた [2]．

D0Me
µ1µ2µ

ラグランジュ点とは二つの天体の引力と宇宙機の遠心力が釣り合う点であり L1 から L5 まで存在し，ラグランジュ点近傍の周期的な軌道をハロー軌道という．特に地球・月系の L2 点近傍のハロー軌道は月に隠れることがないので軌道上の衛星は地球から常に観測可能である [3]．従って，n

衛星を仲介して月の裏側の月面基地と地球との交信などが可能である．今後，ハロー軌道の重要性は増していくと考ρえられる．しかし，この軌道は非常に不安定であるので，制御を行わなければ軌道を維持することが困難である．本

D1

研究では，地球・月系において L2 点に宇宙港が存在するとし，そこから宇宙機を打ち上げることを想定し，ハロー軌道への移行とその軌道を維持することを目的とする．軌
道への移行・維持を行う際には，インパルス制御を用い，

D2

説明

値

地球-月間の距離
地球と月の質量比µ1 = GMe とする定数µ2 = GM とする定数µ = µ1 + µ2 とする定数n = (µ/D03 )1/2 とす

384,748 [km]81.3045M
398,601 [km3 /s2 ]4,887 [km3 /s2 ]403,488 [km3 /s2 ]る角速度

ρ = M/(Me + M ) とする月の質量の割合

ρD0 と す る 地 球 か ら重心までの距離(1 − ρ)D0 とする月から重心までの距離かかる時間 (整定時間) を用いる．

24,674 [km]380,073 [km]万有引力定数

6.67384 ×
10−11[m3 s−2 kg−1 ]D

D = RM − Re-u

円制限三体問題

0.01215G

その制御性能について考察する．また，性能の評価には消費燃料と比例の関係にある総速度変化 (∆V ) と軌道移行に

2.661699 ×
10−6 [rad/s]

状態方程式における入力ベクトル-

地球，月，宇宙機の三体問題を考える．宇宙機の質量が月，地球に比べて十分小さいので，地球と月の運動に影響(1) と (2) よりを与えない．これを制限三体問題と呼ぶ．ニュートンの運D̈ = −

動方程式は

GMe m
GMe M
Me R̈e =D+
reD3
re3
GMe M
GM m
M R̈M = −D+rD3r3(1)G(Me + M )rre
D + Gm( 3 − 3 )3Drreとなり，また m ≪ Me ，M より (3) はD̈ = −(2)となる．

(3)µDD3となる．R を地球・月系の重心からの位置ベクトルとする．mu を推進力とすると，ニュートンの運動方程式より

Ym

宇宙機

mR̈ = −rerR(4)となる．O − {i, j, k} は，原点 O を地球・月系の重心，i地球

Me

GMe m
GM mr + mure −3rer3月XObM

トル，k を鉛直方向の単位ベクトルとする回転座標系である．また，地球・月系の円運動の角速度ベクトルは nk とRm
Reを動径方向の単位ベクトル，j を月の回転方向の単位ベクRMO

図 1 三体問題なる．さらにR = Xi + Y j + Zk
r e = D1 i + R
r = R − D2 i

(5)
(6)
(7)

とおく．(4)，(5)，(6)，(7) を用い，τ = t/(1/n)，X̄ = となる．ここでX/D0 ，Ȳ = Y /D0 ，Z̄ = Z/D0 とし，無次元化をすると，Tx(τ ) = ( x̄ ȳ x̄′ ȳ ′ z̄ z̄ ′ )
001 0 0 0

000 1 0 0 X̄ ′′ − 2Ȳ ′ − X̄
00 2 0 0 
 2σ + 1A=
1−ρρ0−σ + 1 −2 0 0 0 
= − 3 (X̄ + ρ) − 3 (X̄ − 1 + ρ) + ūx(8)000 0 0 1 r̄er̄
000 0 −σ 01−ρρȲ ′′ + 2X̄ ′ − Ȳ = − 3 Ȳ − 3 Ȳ + ūy(9)0 0 0r̄er̄
 0 0 0 1−ρρ
 1 0 0 Z̄ ′′ = − 3 Z̄ − 3 Z̄ + ūz(10)B= 0 1 0 r̄er̄
 0 0 0 0 0 1となる．ここで
1−ρρl−2σx̄−(x̄+l+ρ)−(x̄+l−1+ρ)22[]1/2
 2r̄e3r̄3r̄e = (X̄ + ρ)2 + Ȳ 2 + Z̄ 21−ρρσȳ−ȳ−ȳ[]f(x)=
1/2r̄e3r̄3r̄ = (X̄ − 1 + ρ)2 + Ȳ 2 + Z̄ 21−ρρuxσz̄
−z̄
−z̄
33ūx = 2r̄er̄
n D0reである．
r̄e =D0

インパルス制御4であり，′ は τ による微分を示す．以下では，運動方程式(8)，(9)，(10) を用いる．

3パルス入力において，積分を一定にし時間を無限小，入力を無限大にしたものはインパルス入力と呼ばれる．連続時間のシステムは

ハロー軌道

ẋ(t) = Ax(t) + Bu(t)， x(0) = x0(8)，(9)，(10) は (11) を満たす定常点を持ち，ラグランジュ点として知られている．となる．時刻 [s, s + h] で入力 us /h を与えると，t < s の1−ρρ(X̄ + ρ) + 3 (X̄ − 1 + ρ)r̄e3r̄
ρ1−ρȲ + 3 ȲȲ =3r̄er̄
Z̄ = 0とき

X̄ =x(t) = eAt x0(11)となり，t > s + h のとき

x(t) = eAt x0 + eAt
その解のうち，月の裏側にある点を L2 点と呼び，以下で与えられる [1]．

1h∫ s+he−Ar Bus drsとなる．h → 0 とすれば

x(t) = eAt x0 + eA(t−s) Bus
L2 = (l2 (ρ), 0, 0)， l2 (ρ) = 1.15568となる．s における左右の極限値はここで，L2 点を原点とした座標系 {x̄, ȳ, z̄} を考えると，X̄ ，Ȳ ，Z̄ はそれぞれ x̄ + l2 ，ȳ ，z̄ と表すことができ，その原点で線形化すると以下のようになる．

x(s−) = lim x(t) = eAs x0t↗sx(s+) = lim x(t) = eAs x0 + Bust↘sとなる．ここで，一定時間ごとにインパルス入力を入れる．

x̄′′ − 2ȳ ′ − (2σ + 1)x̄ = ūxȳ ′′ + 2x̄′ + (σ − 1)ȳ = ūyz̄ ′′ + σz̄ = ūz

時間 ks に uk−1 を加え，x(k) = x(ks+) とすると

x1 = x(s+) = eAs x0 + Bu0xk+1 =ここで σ = ρ/|l2 (ρ) − 1 + ρ| + (1 − ρ)/|l2 (ρ) + ρ| =333.19043 である．また，X̄ ，Ȳ ，Z̄ をそれぞれ x̄ + l2 ，ȳ ，z̄で置き換えた場合の (8)，(9)，(10) の状態方程式はẋ(τ ) = Ax(τ ) + Bf (x) + Bu， x(0) = x0

limt↘(k+1)sx(t) = eAs xk + Bukが得られ，状態方程式は漸化式

xk+1 = Ad xk + Bukとなる．ここで Ad = eAs である．

(12)

図 2 より α = 12 で ∆V が最小となる r は 0.2，そのと最適レギュレータ

5きの ∆V の値は 6.9 × 10−1 となり，さらに同じ R で整定システム (12) の最適レギュレータ問題とはJ(u; x0 ) =∞∑時間も 11.4 と最小となった．

( T)
xk Qxk + uTk Ruk

7.2

中継軌道を導入した軌道移行制御

k=0を最小にする問題である．ここで，Q は半正定行列，R は軌道移行制御において中継軌道を導入したシミュレー√正定行列である．(A, B) が可制御，( Q, A) が可観測であるとき，uk は離散時間リッカチ代数方程式 (DARE)ションを行う．本研究では，総速度変化を抑える方法とし

X = Q + ATd XA − ATd XB(R + B T XB)−1 B T XAdr に対する総速度変化のグラフを図 3 に，α = 30 の場合のr に対する総速度変化のグラフを図 4 に示す．の解 X によりて，中継軌道を導入する方法を提案する．また，中継軌道を導入しない場合と比較し検証を行う．α = 12 の場合のuk = −(R + B T XB)−1 B T XAd xk

2
中継軌道3つ中継軌道4つ中継軌道5つ1.8で与えられる．フィードバックゲインを

1.6

1.4
∆ VK = (R + B T XB)−1 B T XAd
とおく．次に，目標軌道を自由運動

1.2

10.8

xf k+1 = Ad xf k ， xf (0) = xf 0

(13)

0.6

0.4-2とし，(12) の軌道を (13) の軌道に追従する制御を考える．

-1.5-1-0.5

00.5

1r

図 3 α = 12 のときの r − ∆V

2 つの軌道の誤差 ek = xk − xf k はek+1 = Ad ek + Buk
0.75を満たす．よって，安定化フィードバック uk = −Kek に中継軌道3つ中継軌道4つ中継軌道5つ0.7より誤差 ek は 0 に収束し，目標の軌道に追従する．

0.65

0.6

ハロー軌道への移行

0.55∆ V60.5

目標とするハロー軌道と制御軌道の方程式は非線
0.45形 で あ る が ，本 研 究 で は ，線 形 フ ィ ー ド バ ッ ク uk =−K (x(ks+) − xf (ks)) を用いる．

シミュレーション結果

70.4

0.35-2-1.5-1-0.500.5

11.52r

重み行列 R と 1 周期ごとのインパルス入力回数 α を変化

図 4 α = 30 のときの r − ∆Vさせ，シミュレーションを行う．ここで重み行列 Q = I6 ，
図 3 より，α = 12 のとき ∆V が最小となる r は，中継R = 10r I3 と置き，r を変化させその軌道特性を調べる． 軌道を 3 つ導入した場合では r = 0.5，4 つ導入した場合

7.1

では r = 0.6，5 つ導入した場合では r = 0.5 となった．ま軌道移行制御

制御軌道の初期値を L2 点におき，ハロー軌道に移行させるための制御を行う．ここで，許容誤差を 10−3 とする．

た，図 4 より，α = 30 のとき ∆V が最小となる r は，中継軌道を 3 つ導入した場合では r = 1.2，4 つ導入した場合では r = 1.4，5 つ導入した場合では r = 1.4 となった．

301.528

定時間を示し，表 3 に α = 30 の場合の ∆V の最小値と1.426

そのときの r，整定時間を示す．また，実際の軌道シミュ

1.3

24

1.2

22

1.1

20

1180.9

160.8

140.7
0.6-3∆V
整定時間-2.5-2-1.5-1-0.50レーションを図 5 に示す．
整定時間

∆V

表 2 に α = 12 の場合の ∆V の最小値とそのときの r，整1.6

表 2 α = 12 の場合の ∆V と整定時間

12

中継 3 つ100.5

中継 4 つr

中継 5 つ図 2 α = 12 のときの軌道移行制御 r − ∆V

中継無し

r
0.5
0.6
0.5
0.2

∆V
5.4×10−1
5.1×10−1
5.9×10−1
6.9×10−1

整定時間

19.1
21.6
31.6
11.4

-5

表 3 α = 30 の場合の ∆V と整定時間

7

整定時間

6

中継 3 つ中継 4 つ中継 5 つ中継無し∆V
3.5×10−1
3.5×10−1
3.8×10−1
4.0×10−1

等間隔インパルス
任意のインパルス間隔

12.3
17.7
21.6
3.55ΔVr
1.2
1.4
1.41.2x 10432L2点中継軌道1中継軌道2中継軌道3中継軌道4目標軌道
制御軌道

0.01z0.005

168101214161820α

図 7 最適な R を用いた軌道維持制御 α − ∆V0

これより，等間隔インパルスでは α が増加したときには-0.005

-0.01
0.1
0.05

0.04
0.02
00-0.05
-0.02
y-0.1

-0.04x

図 5 中継軌道を 4 つ導入した場合の三次元運動∆V が単調減少であり，インパルス時刻を変動させる場合も同様の傾向がみられる．次に α = 12, r = 0.4 で後者の場合のインパルス入力を 1 周期分示したハロー軌道の yz
平面上への射影を図 8 に示した．また，この条件にかかわらず初期値付近に入力をいれない場合 ∆V を抑えられるという結果が得られた．

表 2，表 3 より中継軌道を導入すると，∆V を抑えることができるが，中継軌道を介している為，整定時間は増加-3x 1014することがわかる．また，α = 12 の場合の総速度変化は制御軌道
インパルス入力

12

中継軌道を導入しないときに比べて最大で約 26%，最小で108

約 14% 改善され，α = 30 の場合の総速度変化は中継軌道6を導入しないときに比べて最大で約 13%，最小で約 5% 改z4善された．以上より，入力回数が少ない方が中継軌道を導20入する効果が大きくなる．

-2-47.3

任意のインパルス入力間隔による軌道維持制御-6

制御軌道と目標軌道の方程式を同じ初期値で解き，ハ-8-0.1

-0.05

00.05

0.1

y

ロー軌道を維持するための制御を行う．本節では等間隔にインパルス入力をいれる場合とさらに任意のインパルス時刻の場合にわけて考察を行う．

図 8 α = 12, r = 0.4 の yz 平面

8おわりに中継軌道を導入したシミュレーション結果より，中継軌-5道を導入することによって，総速度変化を抑えつつ目標軌

x 10109

等間隔インパルス
任意のインパルス間隔道への移行が可能であることが確認できた．しかし，中継軌道の設定の仕方や，制御入力の回数を変えた場合など，87条件を変化させたときには，再度シミュレーションを行い，∆V6検証する必要がある．

5軌道維持制御ではインパルス時刻の変動を許し，さらに4初期値付近では自由運動にして制御を行う場合が総速度変32化が抑えられることを示した．本研究で用いたインパルス
1-3-2-101r

図 6 α = 12 のときの軌道維持制御 r − ∆V

入力時刻の探索アルゴリズムはすべての場合を調べてはいないため，アルゴリズムの改善によってはさらに総速度変化が抑えられる可能性がある．

α = 12 と固定した図 6 より等間隔インパルスでは ∆Vが最小となる r は r = 0.4 で，∆V = 3.2×10−5 となる．イ 参考文献
ンパルス入力時刻を任意に変動させる場合では r = −0.9 [1] A.Ichikawa: Recent Developments in Formation Flyのとき ∆V = 1.4 × 10−5 の最小値となっており，総速度変ing, Lecture Notes, 2010.
化は等間隔インパルスの場合に比べて約 56% 改善された．
[2] 歌島昌由：宇宙開発事業団技術報告 ラグランジュ点
上記より α によって ∆V が最小となる R が存在する．
近傍の軌道力学，宇宙開発事業団，1997．
等間隔インパルスの ∆V が最小の R をそれぞれの α に用
[3] 木下宙：天体と軌道の力学，東京大学出版，1998．
いて α − ∆V を図 7 に示す．