ラグランジュ点近傍における小ハロー軌道の生成
2010SE072 伊藤凌平 2010SE105 小出和直

2010SE122 宮沢龍太郎

指導教員：市川朗

Yre

はじめに

1mr
RobMe

人工衛星は通信や測位などに利用されており, 我々の生Rm

活に不可欠である. 近年では惑星の観測技術においての大Re

MX

RM

きな進展をみせている. 人工衛星は与えられたミッション

Oによってさまざまな軌道を採用しており, 用途に応じて適切な軌道を設計することが観測精度や燃費の向上につながる. ここで月面観測について考える. 月の自転周期は地球

図 1 三体問題の周囲を回る公転周期と同期しているため, 月は地球に対して常に同じ面を向けている. したがって, 月の裏側は地球から観測することが出来ない. この問題に対し, 本研究

慣性座標系における原点 O から地球, 月, 宇宙機の位置

では月の裏側を宇宙機によって効率的に観測することを可ベクトルをそれぞれ Re , RM , Rm とする. ここで, Me ,能とする小ハロー軌道をラグランジュ点近傍に生成する.

M , m はそれぞれ地球, 月, 宇宙機の質量であり, G は万ラグランジュ点とは, 質量差のある二つの天体が共通重心

有引力定数で月と地球の重心を原点とする回転座標系の原の周りをそれぞれ円軌道を描いて回っているとき, この二点を ob − {ib , j b , kb } とする. また, r e と r はそれぞれ地天体に比べて質量が無視できるほど小さな宇宙機をある速球, 月から宇宙機への位置ベクトルであり, D=RM -Re ,度を与えてこの軌道面内に置くと, 最初の二天体との相対

D=|D|, re =|r e |, r=|r| である. このとき地球と月の運動位置を変えずに回り続けることができる位置のことであ方程式は,り, 五つの点が存在する. ラグランジュ点では二天体が作

GmGM
D + 3 re ,3DreGMeGmR̈M = − 3 D + 3 r,Drる重力場が遠心力と釣り合っているため宇宙機は二天体にR̈e =対して相対的に不動のままでいることができる. 本研究においては, 地球-月系のラグランジュ点のうち, 地球から見て月の裏側に位置する L2 点 近傍で小ハロー軌道を生成する. ハロー軌道とは, ラグランジュ点の周りを周回する周期軌道である. ここでは, 二天体の重力と宇宙機の向心加となり, 相対ベクトル D は,速度が複雑に関係した三体問題において実現されるものでD̈ = −あり, 宇宙機の運動方程式を L2 点近傍で線形化することにより小ハロー軌道を生成する. また, 宇宙機を L2 点から小ハロー軌道へ移行し, 維持するために, 最適レギュレータを用いて状態フィードバックを設計する. この状態フィードバック制御の性能は制御に必要な総速度変化（ノルム）

G(Me + M )rre
D + Gm( 3 − 3 ).
D3rreを満たす. 制限三体問題において, m≪Me , M であり, 地球と月の二体問題の方程式は,により評価する.

軌道方程式の導出

2D̈ = −µD,D3となる. ここで, µ=G(Me + M ) である. ob から宇宙機への位置ベクトルを R とする. このとき, 宇宙機の運動方
程式

2.1

宇宙機の運動方程式

R̈ = −GMe
GMr e − 3 r + u,
3rer

(1)

慣性座標系の原点を O-{I,J ,K} とし, 地球と月の制限三体問題を考える.

が得られる. ここで, u は制御加速度である. 地球と月の共通重心のまわりの円運動であると仮定し, パラメータのūz = uz /n2 D0 , r̄e = re /D0 , r̄ = r/D0 とし, (5) 式を無次元化する. ′ は τ に関しての微分を意味する. (5) 式より,値をD0 = 384, 748 [km]Me = 81.3045M
3ρ
1−ρ
(X̄ + ρ) − 3 (X̄ − 1 + ρ) + u¯x ,r¯e 3r̄
1−ρρȲ ′′ + 2X̄ ′ − Ȳ = −Ȳ − 3 Ȳ + ūy ,(6)3r¯er̄
1−ρρZ̄ ′′ = −Z̄ − 3 Z̄ + ūz ,r¯e 3r̄

X̄ ′′ − 2Ȳ ′ − X̄ = −2µ1 = GMe = 398, 601 [km /s ]µ2 = GM = 4887 [km3 /s2 ]n = (µ/D03 )1/2 = 2.661699 × 10−6 [rad/s]
ρ = M/(Me + M ) = 0.01215
D1 = ρD0 = 4674 [km]が得られる. ここで,r¯e ,r̄ を

D2 = (1 − ρ)D0 = 380, 073 [km]r¯e = [(X̄ + ρ)2 + Ȳ 2 + Z̄ 2 ]1/2 ,とする. D0 は, 地球の中心から月の中心までの距離であr̄ = [(X̄ − 1 + ρ)2 + Ȳ 2 + Z̄ 2 ]1/2 ,り, n は円運動の角速度である. ここで, ib は地球から見た月の方向であり, 角速度ベクトルは nkb である.

とおく.

R = Xib + Y j b + Zkb ,r e = D1 ib + R,
r = R − D2 ib ,2.3(2)方程式の線形化

以下では, 小ハロー軌道の方程式を導出する.

Yとすると, (1) 式より,{}GMGMeR̈ = − 3 (X + D1 ) − 3 (X − D2 ) ibrr{ e}GMeGM+ − 3 Y − 3 Y jbrer{}
GM
GMe+ − 3 Z − 3 Z kb + urer

1(3)-1が得られる. ここで, D1 と D2 はそれぞれ重心と地球, 月M Lob

1X-1の間の距離である. ここで回転座標系の ib ,j b ,kb の一階微分は, nj b , −nib , 0 となり, (2) 式の R の二階微分は,

図 2 ラグランジュ点

R̈ =(Ẍ − 2nẎ − n2 X)ib(4)+ (Ÿ + 2nẊ − n2 Y )j b + Z̈kb ,(6) 式はラグランジュ点と呼ばれる静止点 Li を持ち,1−ρρ(X̄ + ρ) + 3 (X̄ − 1 + ρ),r¯e 3r̄
1−ρρȲ =Ȳ + 3 Ȳ ,3r¯er̄Z̄ = 0,となる. (3) 式と (4) 式の係数比較をすると,X̄ =GMe(X + D1 )re3
GM− 3 (X − D2 ) + ux ,rGMe
GM2Ÿ + 2nẊ − n Y = − 3 Y − 3 Y + uy ,rer
GMeGMZ̈ = − 3 Z − 3 Z + uz ,rerẌ − 2nẎ − n2 X = −(5)が得られる. ここでは, re , r を,

re = [(X + D1 )2 + Y 2 + Z 2 ]1/2 ,とおく.

(7)で 与 え ら れ る. 月 の 裏 側 の ラ グ ラ ン ジ ュ 点 を L2 =(l2 (ρ), 0, 0) とおく. (7) 式は, ρ の値を変えることによって, 一般の三体システムの定常点も求められ, 地球-月系では ρ=0.01215 であり, (l1 , l2 , l3 )=(0.83692,1.15568,1.00506) となる. (6) 式を L2 点周りで 1 次の項までテイラー展開し, 線形化した方程式を以下に示す.

X¯′′ − 2Y¯′ − (2σ + 1)X̄ = ūx ,Y¯′′ + 2X̄ ′ + (σ − 1)Ȳ = ūy ,r = [(X − D2 )2 + Y 2 + Z 2 ]1/2 ,2.2Meoe(8)Z¯′′ + σ Z̄ = ūz ,方程式の無次元化

これを用いて, 小ハロー軌道を生成する. ここで, σ =計算を簡単にする為に, τ = t/(1/n), X̄ = X/D0 ,22Ȳ = Y /D0 , Z̄ = Z/D0 , ūx = ux /n D0 , ūy = uy /n D0 ,

ρ/(l2 (ρ) − 1 + ρ)3 + (1 − ρ)/(l2 (ρ) + ρ)3 = 3.19043 である.

3

周波数制御
面内運動, 面外運動はそれぞれ周期解をもつが, (8) 式の3 次元運動の周期解は存在しない.
本研究では常に地球の裏側を安定的に観測するために, 解の周波数を面内運動の ωxy , 面外運動 のωz , または任意の周波数 ωxyz , ω = 1, 1/2, 2, 3 に合わせ, 三次元周期解を生成する制御を行う. ここで ωxyz =ωxy +ωz2= 1.8244 とする. (8) 式の状態方程式は,ẋ = Ax + Bu, x(0) = x0 ,図 3 リサージュ軌道である. ここで,[x = X̄ūyu = [ūxX̄ ′Ȳ

0
0 2σ + 1A=0
00Y¯′]TZ̄

Z̄ ′01200000000−σā=0.0091 のとき, 最大の振幅は 3500[km] となり, この軌道上の宇宙機は地球から観測することができる. なお,(Ȳ , Z̄) 運動は, 地球からの概観である. 周期解を求めるため, まず面外運動の周波数を面内運動の周波数に合わせる制御をおこなう. フィードバック制御により, 特殊解が (8)式の第２式を満たすためには,,Tūz ] ,0001−σ00100−20000 0 ,0 
10ūy = f Ȳ ,となる. ここで,f=

≡ diag[Ain , Aout ],
0 0 0
 0 0 0 
 1 0 0 
B =  0 1 0 ,
 0 0 0 0 0 1である. 面外運動の周波数を面内運動の周波数に合わせるには,ūz = −(ω 2 − ωz2 )Z̄,とおく. 以上より,≡ diag[Bin , Bout ],X̄(τ ) = −(ā/k) sin ωτ,とする. 面内運動, 面外運動のシステム行列の特性方程Ȳ (τ ) = −ā cos ωτ,式は,Z̄(τ ) = ā sin ωτ,|λI − Ain | = λ4 − (σ−2)λ2 − (2σ + 1)(σ − 1) = 0,となる.

|λI − Aout | = λ2 + σ = 0,と な る.

面内運動の周波数は解の虚数根である√

ωxy =1.8627 であり, 面外運動は周波数 σ=ωz =1.7862の正弦波である. 周波数制御をするにあたり, 本研究では7 つの周波数に合わせる制御を行った. ここでは, 任意の周波数 ω での周波数制御を例にあげる. 初めに, (8) 式の解を,(9)Z̄(τ ) = ā sin ωz τ,の形で求めると,k=ω 2 + 2σ + 12ωとなる. このとき特殊解は (8) 式の自由運動を満たし, 図 3のようなリサージュ軌道を形成する.

4シミュレーション

4.1

最適レギュレータ(8) 式の状態方程式より,ẋ = Ax + Bu, x(0) = x0 ,である. そのときの評価関数は以下に表わす.

∫ ∞X̄(τ ) = −(ā/k) sin ωxy τ,Ȳ (τ ) = −ā cos ωxy τ,

2ω+ σ − 1 + ω2 ,kJ(u; x0 ) =[xT Qx + uT (t)Ru(t)]dt,0

Q は状態にかかる重みであり, 6 × 6 の半正定対称行列である. また, R は入力にかかる重みであり, 3 × 3 の正定対称行列である. 本研究では Q を固定して考えているため,評価関数を最小にする u(t) を求める. そのために以下の代数リッカチ方程式より適切な X を求める.
AT X + XA + Q − XBR−1 B T X = 0,フィードバックゲイン K は上の式より,

K = R−1 B T X,となる [1].

4.2

誤差方程式

制御軌動を ẋ = Ax+Bu, 目標軌道を x˙f = Axf +Bufとするとき, 制御軌道と目標軌道の誤差は,

e = x − xf

図 5 維持制御ė = Ax + Bu − Axf − Bufωxy での維持制御における R とフィードバックゲイン
K を以下に示す.

となり,= Ae + B(u − uf ),R = 1014.8750 × I3 ,となる. ここで, フィードバックを u − uf = −Ke とす
ると,T
12.6919
7.99900.0000−2.3739 −1.4961
0.0000 
2.33950.0000 
 3.7120
K =  2.3395
1.47440.0000  0.0000
0.00003 × 10−15 0.0000
0.0000 1.475 × 10−7ė = (A − BK)e,となる.

4.3

維持制御

以下のシミュレーションでの目標軌道（維持軌道）と制御軌道の初期値はそれぞれ,x0 = [0ā āω/kxf 0 = [0000

R = 1013.5000 × I3 ,T00

ωxyz , ωz , ω = 1/2, 1, 2, 3 での維持制御における R とフィードバックゲインを以下に示す.

āω]
0T0]

 10.4407

であり, Q = 10 × I6 とする. 図 4 は, 周波数制御をおこ

−1.9528
 3.0536
K= 1.9245

0.00000.0000

なった目標軌道であり 4 周期目までは収束しているが, 4周期目以降で発散していくことが分かった.

6.5802−1.23071.92451.21290.00000.0000T0.0000
0.0000
0.00000.0000
4.93756.4391

それぞれの周波数におけるノルムのグラフと値を以下に示す.

図 4 小ハロー軌道 (4 周目以降発散)目標軌道を維持させるため, 4 周期までに目標軌道の値を初期値にもどす制御を行う必要がある. このときの軌道のシミュレーション結果を図 5 に示す.

図 6 維持制御にかかる周波数毎のノルムの比較

表 1 それぞれの周波数における最小のノルム

ω
xyz
2xyz
11/2
3

ノルム

5.4585 × 10−3
4.2261 × 10−3
2.2103 × 10−2
1.8269 × 10−2
9.0267 × 10−2
2.0996 × 10−1
6.5837 × 10−2
図 9 ノルムと整定時間の比較

図 6 より, 維持制御において r を変化させてもノルムの変化は極めて微小であることがわかる. また, 維持制御にかかるノルムにおいて ωz が最小であることが示された.

4.4移行制御

次に, 制御対象を L2 点から目標軌道へ移行させる制御を行う. ここでは 7 つの周波数に対して r を変えたときのノルムと整定時間の比較を行う. ここでは, 誤差が 10−5 以内のとき収束とした. また, 今回は整定時間の上限を小ハロー軌道 3 周期分の ST = 12.4350 とした.

図 7, 図 8 では, ωxy は r = 0.6250, ωxyz , ωz , ω = 2, 3 はr = 0.1250, ω = 1, 1/2 は r = 0.2500 で最小のノルムとなり, それ以降では発散している. また図 9 より, ノルムと整定時間の両方が小さくなる r が存在することがわかる. このときのそれぞれの R とフィードバックゲイン K を以下に示す.

R = 100.6250 × I3 ,
 12.736

−1.8847
 4.6579
K= 1.18560.00000.0000T0.0000
0.0000
0.0000
0.0000
0.35221.75385.3073−0.2885
1.18562.26520.00000.0000R = 100.1250 × I3 ,
 13.855

4.4795
−1.9953 0.6932
 5.6133 0.7747
K= 0.7737 3.3740

0.0000 0.0000
0.0000 0.0000

0.0000
0.0000
0.00000.0000
1.01413.0866T

R = 100.2500 × I3 ,図 7 ノルムの比較

 13.541

4.7067
−1.9599 0.3747
 5.3291 0.8820
K= 0.8820 3.0204

0.0000 0.0000
0.0000 0.0000

0.0000
0.0000
0.00000.0000
0.78482.6820T

それぞれの周波数におけるノルムを以下に示す.
表 2 それぞれの周波数における最小のノルム

図 8 整定時間の比較

ωxyz2xyz
11/2
3

ノルム

1.770 × 10−2
1.750 × 10−2
1.810 × 10−2
1.760 × 10−2
1.790 × 10−2
2.000 × 10−22.160 × 10−2

表より, 移行制御において最小のノルムの周波数は ω =z であった. 以上より, 移行制御と維持制御において最小のノルムの周波数は ωz となることがわかる.

4.5ミッションにおける燃費評価

小ハロー軌道を生成するにあたり, 周波数が整数の場合は月の公転周期に小ハロー軌道の周期が比例するため観測時に都合がよい. そこで周波数が整数の場合で最もノルムの低い ω = 2 で維持制御を行うが好ましいと考える. ここ
図 12 移行制御 ωz

では, 移行制御において維持制御と同様に周波数 ω = 2 を用いる場合と, 移行制御において最もノルムの低い ωz を用いる場合の全時間にかかるノルムを比較し, どのようにこのときの全時間にかかる入力を図 13 に示す.

差が出るかを検証していく. まず, 周波数を変えずに小ハロー軌道を生成するシミュレーションを示す. 初めに, L2点から目標軌道に移行させるシミュレーションを図 10 に示す.

図 13 全時間での入力 (ωz → ω = 2)

シミュレーションより, ωz から ω = 2 に変化させることで, 大きな入力が見られた. 移行から維持にかけた全時間に
図 10 移行制御 ω=2かかるノルムは 0.1474 となった. 周波数を変化させる場
合と周波数を変化させない場合を比較すると, 周波数を変化させる場合の方が周波数を変化させない場合に比べて小このときの全体にかかる入力を図 11 に示す.

ハロー軌道 1.5 周分のノルムが節約できることがわかる．

5

おわりに
本研究では月の裏側に小ハロー軌道を生成し, 宇宙機をその軌道に乗せて観測を行うために, 周波数制御と収束周期内に加える制御を用いて目標軌道を維持させ, L2 点から人工衛星を目標軌道へ移行させる制御を行った. また 7 つの周波数に合わせた制御を行い, 制御にかかるノルムを比較し, より良い周波数を検討した. 周波数制御を行う中で目標軌道を維持させるのに最もノルムの低い設定をし, 移行制御から維持制御にかけて周波数を変化させることでノルムを抑えることを可能にした. その結果, ミッションを
図 11 全時間での入力 (ω = 2 → ω = 2)行うにあたり, 観測しやすくかつ燃料消費も抑えられる制御方法を見つけることができた.

参考文献
移行から維持にかけた全時間にかかるノルムは 0.1832となった. 次に, 移行制御から維持制御にかけて周波数を

ωz から ω = 2 に変えて, 小ハロー軌道を生成するシミュ
レーションを以下に示す. 図 12 は L2 点から目標軌道までを ωz で移行したシミュレーションである.

[1] 森泰親, 現代制御理論, 森北出版,pp94-96(2013)[2] A. Ichikawa., M. Bando.,”Formation Flying near theLibration Points in the Elliptic Restricted ThreeBody Problem”.