円制限三体問題におけるハロー軌道への移行とその維持
2010SE046 廣村達哉 2010SE257 宇佐美元啓
指導教員：市川朗
Y

はじめに

1

re

月と地球, もしくは, 地球と太陽といった二天体と宇宙機による三体問題では, 天体の質量に対して宇宙機の質量よる既存の重力場内を運動していると考えることができ

rRは無視出来るほど微小である. したがって宇宙機が天体の運動に与える影響も無視できるため，宇宙機は二天体に

mMe

obRmReる. この重力場の力学的平衡点はラグランジュ点と呼ばれL1 から L5 まで存在することが知られており, これらの点の周りをまわる周期軌道をハロー軌道と呼ぶ. ハロー軌
道を利用した世界初の人工衛星は 1978 年に打ち上げられた探査衛星 International Sun-Earth Explorer-3(ISEE-3)である. ISEE-3 は太陽-地球系の L1 点近傍で観測活動を行った [1]. ハロー軌道上の宇宙機には観測だけでなく第MXRM

O - { I,J,K }図 1 三体問題

球と月の重心間距離であり G は万有引力定数である.

二天体の裏側との通信や深宇宙探査衛星の中継基地としてµ = G(Me + M ) = 403, 488[km3 /s2 ]の役割も期待されており, ISEE-3 以後様々な宇宙機が打µ1 = GMe = 398, 601[km3 /s2 ]ち上げられている. しかし, ハロー軌道の方程式は解析的µ2 = GM = 4, 887[km3 /s2 ]D0 = 384, 748[km]D1 = ρD0 = 4, 674[km]D2 = (1 − ρ)D0 = 380, 073[km]に周期解を求められないため，数値的に求めた初期値による周期解を用いて目標軌道を設計し, 実際の運動のずれを補正しなければならない. 本研究では，地球-月系における円制限三体問題の L2 点ハロー軌道を用いる. この軌道は地球から見て月の裏側に位置する L2 点を中心とする最大径約 70,000[km] の軌道であり, 軌道上の宇宙機には, 月の裏側の観測や L2 点に配置した宇宙港と地球との通信などの役割が期待されている. L2 点に宇宙港のような起点を仮定しハロー軌道への移行制御を行う. 収束後は維持制御

n = (µ/D03 )1/2 = 2.661699 × 10−6 [rad/s]
地球の重心と回転座標系の原点 Ob との距離を D1 , Obと月の重心との距離を D2 , 宇宙機の制御加速度を u =[ux uy uz ]T とすると, 宇宙機の運動方程式は座標を用いることで,を行う. これは精度の高い近似周期解を用いても自由運動
では 1 周期目で 35[km] 程度, 2 周期目で 66,000[km] 程度のずれが生じるためである. 制御は最適レギュレータ理論に基づくフィードバック制御により実現する. 制御性能の評価は軌道収束の早さ (整定時間) と燃料消費に比例する速度変化の総和 (L1 ノルム) を用いる. また, 維持制御中の宇宙機のスラスタがトラブルにより停止した場合を仮定し, 目標軌道との誤差が規定内で飛行可能な停止日数の限界とトラブルによって増加したノルムを示すことで, スラスタ故障が宇宙機に与える影響を考察する.

2軌道方程式

2.1

三体問題

三体問題において, 地球と月の共通重心を中心として二天体が円運動するとき, 円制限三体問題という. ここで, 図

1 の通り, Me を地球の質量, M を月の質量, m を宇宙機の質量とする. 質量比は Me : M = 81.3045 : 1 とし, 二天体の総質量に対する月の質量比率は ρ = M/(Me + M ) =0.01215 とする. また, O を慣性座標系 O − {I, J , K} の原点とする. ここで, モデル化に用いた値を示す. D0 は地GMe(X + D1 )re3
GM− 3 (X − D2 ) + uxr
GMGMeŸ + 2nẊ − n2 Y = − 3 Y − 3 Y + uyrerGMeGMZ̈ = − 3 Z − 3 Z + uzrer

Ẍ − 2nẎ − n2 X = −(1)となる. ここで re = [(X + D1 )2 + Y 2 + Z 2 ]1/2 , およびr = [(X − D2 )2 + Y 2 + Z 2 ]1/2 とした. 次に, τ = t/(1/n),
X̄ = X/D0 , Ȳ = Y /D0 , Z̄ = Z/D0 とおき, (1) を無次元化すると,1−ρ(X̄ + ρ)r̄3ρe− 3 (X̄ − 1 + ρ) + ūxr̄
1−ρρ000Ȳ + 2X̄ − Ȳ = − 3 Ȳ − 3 Ȳ + ūyr¯er̄
1−ρρ00Z̄ = − 3 Z̄ − 3 Z̄ + ūzr̄er̄

X̄ 00 − 2Ȳ 0 − X̄ = −(2)0となる.

は, τ による一階微分を表す. ここで r¯e =[(X̄ + ρ) + Ȳ 2 + Z̄ 2 ]1/2 , r̄ = [(X̄ − 1 + ρ)2 + Ȳ 2 + Z̄ 2 ]1/2 ,ūx = ux /n2 D0 , ūy = uy /n2 D0 , ūz = uz /n2 D0 , r̄e =re /D0 , r̄ = r/D0 , である. (2) は, 次の等式を満たすラグ
ランジュ点と呼ばれる点をもつことが知られている [1].
23.2

誤差方程式

目標軌道の初期値と制御軌道の初期値をそれぞ れ, xf0[x0 y0 x00y00= [xf0 yf0 x0f0 yf0 0 zf0 zf0 0 ]T , x0 =z0 z00 ]T とすると, 目標軌道と制御軌道の方程式はそれぞれ

ρ1−ρ
X̄ =(X̄ + ρ) + 3 (X̄ − 1 + ρ)r¯e 3r̄1−ρρȲ =Ȳ + 3 Ȳ
3r¯er̄
Z̄ = 0x0f = Axf + Bh(xf ), xf (0) = xf0x0 = Ax + Bh(x) + Bu, x(0) = x0となる. ここで軌道誤差 e = x − xf は,

円制限三体問題において, ラグランジュ点は 5 つ存在する.

()e0 = Ae + B h(x) − h(xf ) +Bu

そのうち, 地球から見て月の裏側に位置するラグランジュを満たす. 最適レギュレータの状態の重み行列を Q, 入力

点 L2 =(l2 (ρ), 0, 0) において (2) を線形化すると,の重み行列を R とし, リッカチ方程式000X̄ − 2Ȳ − (2σ + 1)X̄ = ūx00

AT X + XA + Q − XBR−1 B T X = 00

Ȳ + 2X̄ + (σ − 1)Ȳ = ūy(3)00

Z̄ + σ Z̄ = ūzの解 X を求める. 求めた X より, フィードバックゲイン

K = R−1 B T X が得られる. 制御入力をとなる. ここで σ = ρ/(l2 (ρ)−1+ρ)3 +(1−ρ)/(l2 (ρ)+ρ)3u = −Ke − h(x) + h(xf )である.
とおくと3制御系設計

3.1

状態方程式の導出となり, 軌道誤差を 0 に収束させることができる.

(X̄, Ȳ , Z̄) = (x + l2 (ρ), y, z) とし, 座標の原点を L2 とする. 線形化方程式 (3) を x = [x y x0 y 0 z z 0 ]T , u =[ux uy uz ]T を用いて, 状態方程式で表すと
x0 = Ax + Bu

(4)



0000
 2σ + 10
A=01−σ0000
 0 0 0  0 0 0  1 0 0 
B= 0 1 0 
0 0 00 0 1100−2000 01 02 00 00 00 −σ4.1

評価方法

ラグランジュ点 L2 から, ハロー軌道へ移行させる制御,ハロー軌道に維持させる制御を行う. 移行制御, 維持制御
力の重み行列を変化させながらシミュレーションを行うこ0 0 0 
10とにより, L1 ノルムの変化を調べる. L1 ノルムは加速度の絶対積分で表される. L1 ノルムは, 燃料消費と比例関係にあるため, これを最小化する q, r が最適な qopt , ropt である. 移行制御においては, 収束までの整定時間 ST も考慮した. 本研究では微分方程式の数値解をルンゲクッタ法を用いて解き, そのときの刻み幅を 10−4 とした.

4.2

移行制御

ラグランジュ点 L2 から, ハロー軌道へ移行させる制御を行う. 移行に用いた状態の重み行列 Qpre , 入力の重み行列 Rpre は, Qpre = diag[10 10 qpre qpre 10 qpre ],x0 = Ax + Bh(x) + Buとなる. ここでl2 − 2σx − 1−ρ
r̄e3 (x + l2 + ρ)− r̄ρ3 (x + l2 − 1 + ρ)ρσy − 1−ρr̄ 3 y − r̄ 3 yeρσz − 1−ρr̄e3 z − r̄ 3 zである.

シミュレーションと性能評価

0 である. また, 非線形方程式 (2) はh(x) = 4それぞれにおける最適レギュレータの状態の重み行列, 入となる. ここで行列 A, B は

e0 = (A − BK)eRpre = diag[10rpre 10rpre 10rpre ] である. Qpre の位置成分を 10 として, 速度成分 qpre を変化させ, L1 ノルムを最小にする qpreopt を求めた. qpre の変化と各 qpre における rpre の変化, そのときの L1 ノルムの推移の 3D シミュレーションを図 2 に示す. ここで, ST = 2.1191 であり, 収束するまでの整定時間である. 最適な値は, qpreopt =2.000, rpreopt = −0.5000 である.
∫ ST
L1 =
(ux 2 + uy 2 + uz 2 )1/2 dτ(5)0表 1 K による L1 ノルムの比較
0.65

0.7

0.6

0.65
0.6

維持制御

0.55

0.55

norm

K = KI
0.24771.6679 × 10−4-移行制御

0.5

0.50.45

K = Kpre0.2284(92.2%)
1.5495 × 10−4(92.9%)

K = Khalo
0.2371(95.7%)
1.5446 × 10−4(92.6%)

0.45

0.4
0.4

0.35

0.35

0.3
0.25

0.30.2210-1r-2-302486100.25

q pre

維持可能条件5

本章では, 目標軌道との軌道誤差 ±10−4（384.7480[km]）以内で維持制御中の宇宙機のスラスタがトラブルにより停止した場合を仮定し, 軌道誤差 ±10−3 (3847.48[km]) 以内で飛行可能な日数の限界を示す. また, スラスタの復旧に限界停止日数までかかった場合と, 正常な維持制御が行わ

図 2 qpre , rpre の変化に対する L1 ノルムの推移れた場合とのノルム変化を示す. スラスタは 3 軸方向に出力可能であり, それぞれが独立であるとする. また, 故4.3

障した軸以外はフィードバック制御が行われているものと維持制御する.

ハ ロ ー 軌 道 に 維 持 さ せ る 制 御 を 行 う.

維持に用いた 状 態 の 重 み 行 列 Qhalo , 入 力 の 重 み 行 列 Rhalo は,Qpre = diag[10 10 qhalo qhalo 10 qhalo ], Rhalo =diag[10rhalo 10rhalo 10rhalo ] である. 移行制御同様, Qhalo
は, 位置成分を 10 として, 速度成分 qhalo を変化させ, L1ノルムを最小にする qhaloopt を求めた. qhalo の変化と各qhalo における rhalo の変化, そのときの L1 ノルムの推移の 3D シミュレーションを図 3 に示す. 最適な値は,qhaloopt = 0.800, rhaloopt = −0.3750 である.
∫ TL1 =(ux 2 + uy 2 + uz 2 )1/2 dτ(6)5.1

限界停止日数と増加ノルム

ハロー軌道上の宇宙機は, 周期的に速度が変化するため,スラスタ停止時のタイミングを考慮する必要がある. そこで, 図 4 のように, ハロー軌道の 1 周期 (約 14.8488 日) を15 分割し, 1 日ずつ停止位置をずらしながら停止位置からの限界停止日数を求めた.

x 10z
x 10-4

day : 550x 10

維持軌道
L2点
故障位置

-310

day : 100-5-445

day : 00.1
0.05
04.5

ynorm

4

3.5

-0.05
-0.1

-0.04

-0.02

00.02

0.04

x3

図 4 故障位置

32.5

221.5
2

10-1

r-2-30246810q halo

図 3 qhalo , rhalo の変化に対する L1 ノルムの推移

4.4性能評価

状態に対する重み行列 Q に単位行列を用いたときのフィードバックゲイン KI Qpre を用いたときの Kpre ,Qhalo を用いたときの Khalo のそれぞれの場合の 移行制御の L1 ノルム と, 維持制御の L1 ノルムの 1 周期平均を表1 に示す. また, KI を用いたときの L1 ノルムを 100 として, Kpre と Khalo の L1 ノルムの比率を括弧内に示す.

シミュレーション結果を表 2 に示す.
表 2 故障日, 故障軸別の限界停止日数
0 日目
1 日目
2 日目
3 日目
4 日目
5 日目
6 日目
7 日目
8 日目
9 日目
10 日目
11 日目
12 日目
13 日目
14 日目

X軸
7981011101112121211101098

XY 軸
87912121213141312121110109XZ 軸
7981011101112121211101098XYZ 軸
87912121213141312121110109Y軸
79486312076
7476
8686858483828180

YZ 軸
794863120767487
86868584
83
82
81
80

Z軸
452451464463462
461
460
459
458
457
456
455
454
453452

次に, 故障したスラスタが限界停止日に復旧し, その後正常な維持制御となるまでに必要となったノルムの増加分を表 3 に示す.

表 3 故障日, 故障軸別の増加ノルム
0 日目
1 日目
2 日目
3 日目
4 日目
5 日目
6 日目
7 日目
8 日目
9 日目
10 日目
11 日目
12 日目
13 日目
14 日目

X軸
0.006286
0.009190
0.009742
0.009733
0.010132
0.008085
0.007132
0.005482
0.005169
0.006621
0.005530
0.004648
0.007508
0.007157
0.007000

XY 軸
0.007101
0.006286
0.010432
0.009569
0.007745
0.005829
0.004554
0.004910
0.004333
0.003907
0.005972
0.004826
0.004198
0.008328
0.008095

XZ 軸
0.006292
0.009163
0.009719
0.009641
0.010173
0.008072
0.007150
0.005466
0.005151
0.006612
0.005522
0.0046400.007516
0.007165
0.007009

XYZ 軸
0.007124
0.006285
0.010499
0.009730
0.007823
0.005862
0.004530
0.004886
0.004319
0.003896
0.005974
0.004829
0.004200
0.008357
0.008124

Y軸
0.022948
0.021565
0.020419
0.026866
0.020376
0.019571
0.020731
0.024050
0.023278
0.022859
0.023141
0.023301
0.023282
0.023232
0.023238

YZ 軸
0.021873
0.020757
0.019542
0.026705
0.019617
0.018812
0.029017
0.022942
0.022184
0.021785
0.022055
0.022208
0.022190
0.022142
0.022147

Z軸
0.010285
0.009565
0.010423
0.010422
0.010346
0.010285
0.010274
0.010289
0.010310
0.010322
0.010326
0.010327
0.010327
0.010328
0.010328

図 9 Y 軸故障

図 10 YZ 軸故障

図 9, 図 10 より, Y 軸故障と YZ 軸故障の波形に大きな違いが見られないことから Z 軸入力の影響は最も小さいことがわかる. また, Y 軸故障関連群では Y 軸故障の影響を

表 2, 表 3 より, X 軸 XY 軸 XZ 軸 XYZ 軸故障, Y 軸

YZ 軸故障の各故障日における限界停止日数および増加ノルムが非常に近い値であることが分かる. また, 停止日数の観点からは X 軸 Y 軸 Z 軸の順に故障の影響が大きく,ノルム増加の観点からは Y 軸 X 軸 Z 軸の順に故障の影響が大きいといえる.
5.2

X 軸入力が補う形になっているため, 限界停止日数は X 軸故障関連群より長くなるものの, 増加ノルムは最も大きくなることがわかる. 最後に, Z 軸故障の入力波形を図 11 に示す.

軸別の故障影響評価

故障軸ごとの限界日数が最小の時の入力波形を示し, 故障の影響評価を行う. 前節の結果から X 軸 XY 軸 XZ 軸XYZ 軸故障を X 軸故障関連群, Y 軸 YZ 軸故障を Y 軸故障関連群とし, 比較を行った. また, Z 軸故障の波形から Z

図 11 Z 軸故障軸スラスタの有用性を検討する. X 軸故障関連群の入力波形を図 5 から図 8 に示す.

図 11 より, Z 軸スラスタ故障の影響は軽微であるといえる. よって, X 軸スラスタと Y 軸スラスタのみで維持制御を行っても X 軸故障関連群の数 10 倍の維持日数となる.
ノルムの増加も Y 軸故障関連群の 1/2 程度である. また,Z 軸スラスタが復旧した時の入力の最大値は X 軸故障関連群や Y 軸故障関連群の入力の最大値より小さい. よって,Z 軸入力は毎周期入力する必要はなく, スラスタも出力の小さいものでよい.
図 5 X 軸故障

図 6 XY 軸故障

6

おわりに地球-月系の円制限三体問題における L2 点近傍のハロー軌道上を運動する宇宙機の制御系設計を行った. 最適レギュレータの状態の重み行列, 入力の重み行列を変化させながらシミュレーションを行うことにより, ハロー軌道への移行制御と維持制御にかかる L1 ノルムを最小化するフィードバックゲイン Kpre , Khalo を示した. さらに, 維
持制御中のスラスタ故障を仮定し, 故障軸スラスタの限界停止日数および故障の影響によるノルムの増加分を求め
図 7 XZ 軸故障

図 8 XYZ 軸故障た. その結果から, スラスタの故障影響評価を故障軸別に行った.

図 5 から図 8 より, X 軸故障関連群では Y 軸スラスタやZ 軸スラスタからの入力があっても, XYZ 軸の全てが故障した場合とほぼ同じ入力波形となることがわかる. これは,X 軸故障の影響力が大きく, YZ 軸のスラスタで補完することができないためであると考えられる. したがって, X軸故障関連群の限界停止日数および増加ノルムが近い値になると考えられる. 次に, Y 軸故障関連群の入力波形を図9, 図 10 に示す.

参考文献
[1] M. Utashima., ”Orbital Mechanics Near Lagrange’sPoints”,1997[2] A. Ichikawa., M. Bando.,”PERIODIC ORBITSAND FORMATION FLYING NEAR THE LIBRATION POINTS”, 2012.