楕円軌道上の相対運動方程式とその周期解
2010SE103 小林史弥 2010SE104 小林弘明
指導教員：市川朗

1

ルは R = R0 i + r である. ニュートンの運動方程式

はじめにR̈ +フォーメーションフライトとは複数の宇宙機による編隊飛行のことである. その基礎研究として, 地球周回軌道上µR=0R3より次の方程式が得られる.の衛星とその近くを飛行する従衛星の相対運動が注目されている. フォーメーションフライトでは従衛星は適切な相対位置関係を維持することが望ましい. 主衛星の軌道が円軌道の場合, その近傍の従衛星の相対運動方程式を原点周りで線形化した Hill-Clohessy-Wiltshire(HCW) 方程式が用いられる. 主衛星の軌道が楕円軌道の場合の線形化した相対運動方程式は Tschauner-Hempel(TH) 方程式と呼ばれる. 本研究では TH 方程式の周期解を用いてフォーメーションフライトの研究を行う.

2µµ= − 3 (x + R0 ) + ux
R02Rµ(3)ÿ + 2θ̇ẋ + θ̈x + θ˙2 x = − 3 + uy
Rµz̈ = − 3 z + uz
R√2
ここで,R= (R0 + x) + y 2 + z 2 である. 原点 x=y=z=0で (3) を線形化すると,
ẍ − 2θ̇ẏ − θ̈y − θ˙2 x −µ)x = ux
R0 3µÿ + 2θ̇ẋ + θ̈x − (θ̇2 −)x = uy
R0 3µz = uzz̈ +R0 3ẍ − 2θ̇ẏ − θ̈y − (θ̇2 + 2

楕円軌道上の相対運動方程式

xr

従衛星

主衛星となる. (4) は Tshauner-Hempel 方程式とよばれる.

yR(4)R0θx = [x y ẋ ẏ z ż]T = [x1 x2 x3 x4 x5 x6 ]T

地球中心u = [ux uy uz ]Tとおくと, 状態方程式はẋ = A(t)x + Bu,x(0) = x0

ここで,
00
 23 θ̇ + 2µ/R0
A(t) = −θ̈
00

図 1 楕円軌道上の主衛星と従衛星
主衛星の軌道はR0 =p,
1 + e cos θ

p = a(1 − e2 )(1)

100−2θ̇00

 00 0 0 0 
0
1000
100

 0
 1B=
 0
00で与えられる楕円軌道とする. ここで

R0 :動径, µ:地球の重力定数, θ:真近点離角, a:長半径,√e:離心率, p:半直弦, 平均運動を n= µ/a3 とする.
楕円軌道の方程式は以下の 2 式で与えられる.
µR̈0 − R0 θ˙2 = − 2RR0 θ̈ + 2Ṙ0 θ̇ = 000θ̈θ̇2 − µ/R0 300012θ̇00000000−R0 3 /µ


00 

0 0 1 0である. 時間 t を無次元化し τ =nt とする. また (x,y,z) をそれぞれ長半径 a により無次元化し(2)主衛星とともに回転する座標系 {i, j, k} を導入する.

i は動径方向, j は飛行方向, k は軌道面外への単位ベクトルである.r を主衛星に対する従衛星の位置ベクトルとし,
r=xi+yj+zk とする. 従衛星の地球中心からの位置ベクト(x̄(τ ), ȳ(τ ), z̄(τ )) = (1/a)(x(τ /n), y(τ /n), z(τ /n))とおく. ここで, 式 (2), 式 (4) を無次元化すると以下のようになる.

R̄0′′ − R̄0 (θ̄′ )2 = −R̄0 θ̄

′′+ 2R̄′0 θ̄′=01R̄02(5)[]
2′x̄′′ − 2θ̄y − θ̄′′ ȳ − (θ̄′ )2 +2 x̄ = ūxR̄0[]1ȳ ′′ + 2θ̄′ x̄′ + (θ̄′ )2 x̄ − (θ̄′ )2 +2 ȳ = ūyR̄01z̄ ′′ = − z̄ + ūzR̄0 0 0 1B=
 0

00(6)

0 0 0 
0 

01000
100(12) [][](θ′ )2 + R23 − 3 x̄1 + θ̄′′ x̄2 + 2θ¯′ − 2 x̄4


[0]
[]
h(x̄) =  θ̄′′ x̄1 + (θ¯′ )2 + R13 x̄2 + 2 − 2(θ¯′ ) x̄3  (13)
( 0 )
1 − R13 x̄5

ここで ′ は τ に関する微分を表す. さらに,[ūx (τ ) ūy (τ ) ūz (τ )] = (1/an2 )[ux (τ /n), uy (τ /n), uz (τ /n)]0θ̄(τ ) = θ(τ /n) とおく (4) を無次元化した状態方程式は
x̄′ = A(τ )x̄ + B ū,00
 ′ 230A(τ ) =  (θ̄ ) −+θ̄2/R′′
00x̄(0) = x̄000θ̄′′θ̄′2 − 1/R0 300

100−2θ̄′000
12θ̄′000(7)

00000−1/R0 3

主衛星の初期位置が近地点のとき従衛星の初期軌道および目標軌道の初期値をそれぞれ以下のようにおく.

[]−0.01(2 + e)
√f0 = 0.01 0 0x0.01 0(1 − e) 1 − e2[]−0.005(2 + e)
√gx=0.005
000.0050f0(1 − e) 1 − e2



00 
0 
 (8)
0 
1 0(14)

主衛星の初期位置が遠地点のときとなる.

][−0.01(2 − e)√f0 = 0.01 0 0x
0.01 0(1 + e) 1 − e2[]−0.005(2 − e)√gx
0.005 0f 0 = 0.005 0 0(1 + e) 1 − e2θ̄′′ は (11) 式の第２式を用いて,
θ̄′′ = −2R̄0′ θ̄R̄0と表せる.

(15)

目標軌道の方程式は以下のように表せる.

TH 方程式は t を θ で置き換え,(ex(θ),ey (θ),ez (θ))=(1+ecos θ)(x,y,z) と お く こ と に よ りYamanaka と Ankerson [1] により解かれた. 軌道面外運動は正弦関数で与えられ周期解となる. 軌道面内の運動がフォーメーションフライトは従衛星を目標軌道に乗せるこ

周期解になるための必要十分条件はとを目的にする. 誤差 e=x − x̄f は以下を満たす.

x̄′f = A(τ )x̄f = Ax̄f + Bh(x̄f )e′ = Ae + B(ū + h(x̄) − h(x̄f ))(3ρ + e2 − 1)ex(θ0 ) + ese
x′ (θ0 ) + ρ2 ye′ (θ0 ) = 0で与えられる. ここで ρ = 1 + e cos θ である. この条件は

ū = −Ke − h(x̄) + h(x̄f )(18)e′ = (A − BK)e(19)′(2+e)ex(0)+(1+e)ey (0) = 0, (2−e)ex(π)+(1−e)ey (π) = 0を用いればとなる. これより, 式 (5) の周期解の初期条件は以下のようになる.

このとき,A − BK が安定であれば誤差は 0 に収束する.

2+e√x̄(0) 近地点(1 − e) 1 − e22−e√ȳ ′ (0) = −x̄(0) 遠地点(1 + e) 1 − e2

フィードバックゲイン K は最適レギュレータのリッカチ

ȳ ′ (0) = −3(17)ここで, フィードバック

θ0 =0,π のとき簡単になる. 実際に代入すると,′(16)

方程式

(9)
A′ X + XA + Q − XBR−1 B ′ X = 0を用いて設計する.

フォーメーション

4

無次元化した方程式 (5) の状態方程式は

x̄′ = Ax̄ + Bh(x̄) + B ū, x(0) = x¯0(10)と表せる. ここで A,B,h(x̄) は次のようになる.

 0 0
 0 0
 3 0A=
 0 0

0 00 0(20)100−20001200000000−10 0 0 0 10

シミュレーション結果
フィードバックゲイン K は最適レギュレータの重み行列を Q ＝ I, R = 10r I とおき, r を変化させて設計する.
はじめに主衛星の初期位置を近地点とする. 図 2 は面内運動の入力の L1 ノルム 図 3 は面外運動の入力の L1 ノルムである.

(11)

ここで離心率 e をパラメータとして表示している. e=0 は主衛星の軌道が円軌道の場合に対応し, L1 ノルムは r と共に単調に減少する. フィードバックゲインはHCW 方程式のリッカチ方程式を用いているため e>0 の場合, L1 ノルムは始めは r と共に減少するが, やがて r とr の依存性は, 主衛星の初期値が近地点の場合と同様で共に増加する. したがって各離心率 e に対して, 最小の L1ある. また整定時間もほぼ同じである. L1 ノルムの大きさノルムを与える r が必ず存在する. 図 4 は目標軌道に収は, 遠地点からの制御を行ったほうが小さくなる.

束するまでの時間のグラフである. 整定時間は r とともに単調に増加し, 離心率に依存しない. これらのグラフを元0.035
e=0.3に L1 ノルム, 整定時間の仕様を満たすフィードバックが0.03設計できる.
0.045
e=0.3
0.04

L1 ノルム

0.025

0.02
e=0.2
0.015

0.035
0.010.03

L1 ノルムe=0.10.0050.025
e=00.02
0-2

e=0.2

-1.5

-1-0.500.5

11.5

2r0.0150.01

図 5 遠地点:面内運動の L1 ノルム

e=0.10.005
e=00-2-1.5-1-0.500.5

11.52r0.04

図 2 近地点:面内運動の L1 ノルム

e=0.3
0.035

0.03

e=0.3
0.035

Lz1ノルム0.040.025
e=0.2

0.02

0.030.015

Lz1ノルムe=0.10.025
e=0.2

0.01

e=00.020.005-2-1.5-1-0.500.5

11.5

2r0.015e=0.1

図 6 遠地点:面外運動の L1 ノルム

0.01

e=00.005-2-1.5-1-0.5

00.5

11.52

また, 実際の軌道シミュレーションを求める.

r

離心率 e=0.3 の場合に図 (6) から最小の L1 ノルムを与図 3 近地点:面外運動の L1 ノルム

える r=0.125 を入力した制御軌道を図 (7) に示した.

制御軌道
初期軌道
目標軌道

90-3x 1080

1070605STz5004030-5200.01-0.04

10-0.0200-2-1.5-1-0.5

00.5

11.50

20.02r

図 4 整定時間 ST
次に主衛星の初期位置を遠地点とした場合の L1 ノルムを図 5, 図 6 に与える.

x-0.010.04y

図 7 制御軌道:離心率 e=0.3このため面内運動の L1 ノルムも急速に増加する.

簡易フィードバックによる制御

5

図 11, 図 12 は主衛星の初期位置が遠地点の時の L1 ノル

ここではフィードバック (18) を線形フィードバックū=−Ke に置き換えた時のシミュレーション結果を紹介する. 簡易フィードバックは離心率が大きくなると収束判定である目標軌道と制御軌道との差が 10−5 以下の条件を満たさなくなるので, 制御可能な r の値を用いる. 図 8, 図 9は, 主衛星の初期位置が近地点の場合の L1 ノルムと L1 ノルムであり,ムである. 整定時間は近地点の場合とほぼ同じである.L1 ノルムは遠地点から制御を行ったほうが小さくなる. フィー
ドバック (18) と比較すると L1 ノルムはやや大きくなる
が, 線形フィードバックでも良好な制御を行うことができる.
0.14

0.12

e=0.3

0.12
0.1

L1 ノルム0.08

L1 ノルム

0.1

e=0.3

0.08

0.06

e=0.2

0.04
0.06

0.020.04
e=0.1
e=0.2
e=00-2-1.5-1-0.5

00.5

11.5

20.02
r
e=0.1
0-2

e=0-1.5-1-0.500.5

11.5

図 11 遠地点:面内運動の L1 ノルム

2r

図 8 近地点:面内運動の L1 ノルム-3x 1012110.014100.013

Lz1ノルム90.01280.011Lz1ノルム7
e=00.016
e=0.1

e=0.30.009

e=0.2
5
e=0.3
e=0.2

0.008
4-2-1.5-1-0.5

00.5

11.5

2

e=0.1
r0.007e=00.006-2-1.5-1-0.500.5

11.5

図 12 遠地点:面外運動の L1 ノルム

2r

図 9 近地点:面内運動の L1 ノルム

6終わりに主衛星の初期位置を遠地点とした場合 L1 ノルムの大きさが小さくなる. 簡易フィードバックは離心率 e=0.1 程度

250

ならばフルフィードバックとあまり差はなかったが, それ200以上の離心率 e では性能が落ちるため, 適切な入力をしな

e=0.3

くてはならない.

ST150
e=0.2

7100

参考文献

e=0.1

50
e=00-2-1.5-1-0.500.5

11.52r

図 10 整定時間 ST
図 10 は整定時間のグラフである.r が 0.5 を超えると目標軌道への収束が遅くなり, やがて収束しなくなる.

1. K.Yamanaka and Ankerson,New state transition matrix for relative motion an orbitaryelliptical orbit,J.Guidance,Control,Dynamics,
vol.25,2002,pp.60-66.
2. Bando,M. and Ichikawa, Akira,”Graphical Generation of Periodic Orbits of the Tshauner-Hempel
Equations Journal of Guidance,Control,andDynamics,Vol.35,No.3,2012,pp-1002-1007