楕円軌道上のフォーメーション
2011SE059 堀田一希 2011SE137 小島竜也
指導教員：市川朗

1

はじめに√(R0 + x)2 + y 2 + z 2 である.
原点 x = y = z = 0 で (3) 式を線形化すると,ここで, R =現在地球を周回する人工衛星の数は 3500 個以上といわれており, その軌道の大きさ, 周期, 軌道と赤道面との傾斜µ)x = ux
R03µÿ + 2θ̇ẋ + θ̈x − (θ˙2 − 3 )y = uyR0µz̈ + 3 z = uz
R0ẍ − 2θ̇ẏ − θ̈y − (θ˙2 + 2

角は様々である. 本研究では, 地球周回楕円軌道上の主衛星とその近くを飛行する従衛星の相対運動について考える.
主衛星の軌道が楕円軌道の場合の線形化した相対運動方程式は,Tschauner-Hempel(TH) 方程式と呼ばれる. この周期解を用いて楕円軌道のフォーメーションの研究を行う.

(4)となる. この (4) 式が Tshauner-Hempel 方程式（TH 方程式）と呼ばれる.

Tx = [ x y ẋ ẏ z ż ]u = [ ux uy uz ]Tとおくと, 状態方程式はẋ = A(t)x + Bu,

図 1 楕円軌道上の主衛星と従衛星

x(0) = x0(5)となる. ここで,

2

楕円軌道
主衛星の軌道はR0 =p, p = a(1 − e2 )1 + ecosθ(1)00 (θ̇)2 + 2µ/R03
A(t) = −θ̈00で与えられるとする. 軌道面内の運動方程式は.

R̈0 − R0 θ̇2 = −µR02R0 θ̈ + 2Ṙ0 θ̇ = 000θ̈(θ̇)2 − µ/R0300

 0(2)

 0 1B= 0
00100−2θ̇00012θ̇00000000−µ/R030 0 0 0 01000
1001で与えられる. また θ̇ = ( pµ3 ) 2 (1 + e cos θ) が成り立つ. 主衛星とともに回転する座標系 {i, j, k} を導入する.

である. 時間 t を無次元化して τ = nt とし,(x̄(τ ), x̄(τ ), x̄(τ )) = (1/a)(x̄(τ /n), x̄(τ /n)x̄(τ /n)) とおiは動径方向, j は飛行方向, kは軌道面外への単位ベクトルr を主衛星に対する従衛星の位置ベクトルとし,r = xi + yj + zkとする.ニュートンの運動方程式よりR̈ + Rµ3 R = 0 が得られ, 次の方程式が得られる.

く. また (x, y, z) をそれぞれ長半径 a によって無次元化する. このとき (2) の方程式は

˙ 2=− 1R̄¨0 − R̄0 (θ̄)
R¯02R̄ θ̄¨ + 2R̄˙ θ̄˙ = 00µµẍ − 2θ̇ẏ − θ̈y − θ˙2 x − 2 = − 3 (x + R0 ) + ux
R0Rµÿ + 2θ̇ẋ − θ̈x − θ˙2 y = − 3 y + uy
Rµz̈ = − 3 z + uz
R(6)0で与えられる. なおこの（6）以降の・は時間 τ に関しての(3)微分を表す.またθ̄ =13(1−e2 ) 2(1 + ecosθ̄)2 が成り立つ.

近地点を固定して離心率 e の値をパラメータとして楕円軌道を表したものを図 2 に示す.

00 
0 0 1 0

て Yamanaka と Ankerson[2] により解かれた.

TH 方程式を解くためには, 軌道方程式 (6) の初期条件は必要である. 無次元化した (1) の楕円軌道は以下のようになる.
R̄0 =1 − e2
1 + ecosθ̄(10)(1 − e2 )esinθ̄θ̄˙esinθ̄
R̄˙0 ==12(1 + ecosθ̄)(1 − e2 ) 2
図 2 楕円軌道（e = 0.1 から 0.9 まで)
よって (6) の初期条件は

˙[R̄0 , θ̄, R̄˙0 , θ̄](τθ)
21−eesinθ̄(1 + ecosθ̄)2=[, θ̄,]1 ,3
1 + ecosθ̄(1 − e2 ) 2 (1 − e2 ) 2

離心率 e = 0 としたとき円になる. このとき図２においてe の値を増やしていくと楕円が大きくなるのは近地点を固定しているためである.
(3),(4) の方程式は.

ここで τθ̄ はθ̄（τθ ）= θ̄ である. 特に˙ 2 x̄ − 1 = − 1 (x̄ + R̄ ) + u¯
¨ − 2θ̄˙ȳ˙ − θ̄¨ȳ − (θ̄)
x̄
0xR̄3R¯02
˙ 2 ȳ = − 1 ȳ + u¯(7)ȳ¨ + 2θ̄˙x̄˙ + θ̄¨x̄ − (θ̄)
yR̄31z̄¨ = − 3 z̄ + u¯zR̄√
1+e](1 − e)3˙[R̄0 (0), θ̄(0), R̄˙0 (0), θ̄(0)]= [1 − e, 0, 0,

˙[R̄0 (τπ/2 ), θ̄(τπ/2 ), R̄˙0 (τπ/2 ), θ̄(τπ/2 )]= [1 − e2 , π/2, e(1 − e2 )− 2 , (1 − e2 )− 2 ]31

˙ 2 + 2 ]x̄ = u¯
¨ − 2θ̄˙ȳ˙ − θ̄¨ȳ − [(θ̄)
x̄
x
R¯03
˙ 2 − 1 ]ȳ = u¯ȳ¨ + 2θ̄˙x̄˙ + θ̄¨x̄ − [(θ̄)
yR¯031z̄¨ + ¯3 z̄ = u¯z
R0˙[R̄0 (τ3π/2 ), θ̄(τ3π/2 ), R̄˙0 (τ3π/2 ), θ̄(τ3π/2 )]= [1 − e2 , 3π/2, −e(1 − e2 )− 2 , (1 − e2 )− 2 ]
31(8)

√˙
[R̄0 (π), θ̄(π), R̄˙0 (π), θ̄(π)]= [1 + e, π, 0,1−e](1 + e)3となる.となる.

ここでは τ に関する微分であることを示し,楕円軌道上の θ = 0 の点を近地点といい,θ = π の点を遠2(ūx (τ ), ūy (τ ), ūz (τ )) = (1/an )(ūx (τ /n), ūy (τ /n),¯ ) = θ(τ /n) とする. 無次元化したūz (τ /n)), θ(τ地点という. 近地点は楕円軌道上で, 地球から一番近いと状態方程式はころを表し, 遠地点は一番遠いところを表す.

x̄˙ = A(τ )x̄ + B ū,
00
 ˙ 2
 (θ̄) + 23R0
A(τ ) = −θ̄¨00x̄(0) = x¯0(9)軌道面外運動は常に周期解となる.(1) の解により軌道面内の運動が周期解であるための必要十分条件は00
θ̄¨100012θ̄˙000

˙ 2− 1(θ̄)
R0300−2θ̄˙
0000000− R130¨ 0 ) + ρ2 ỹ(θ˙ 0 ) = 0 で与えられる
0 (3ρ + e2 − 1)x̃(θ0 ) + esx̃(θ0  .ここでρ = 1 + ecosθ である.この条件はθ0 = 0, π のとき
0  簡単になる. 実際に代入すると,
˙˙0  (2+e)x̃(0)+(1+e)ỹ(0)= 0, (2−e)x̃(π)+(1−e)ỹ(π)=01  となる. これによって, 式 (8) の周期解の初期条件は以下の0ようになる.となる.

TH 方程式は t を θ で置き換え,(x̃(θ), ỹ(θ), z̃(θ)) = (1 + cosθ)(x, y, z) とおくことによっ2+e√x̄(0) (近地点)(1 − e) i − e22−e˙√ȳ(π)=−x̄(π) (遠地点)(1 + e) i − e2˙ȳ(0)=−(11)

3

とおく. 目標軌道の状態方程式はフォーメーション
無次元化した (8) の状態方程式はx̄˙ = Ax̄ + Bh(x̄) + B ū,x¯˙f = Ax¯f + Bh(x¯f )x(0) = x¯0(12)と表せる. ここでは A, B, h(x̄) を次のように与えられる.

 0 0 1 0 0
 0 0 0 1 0
 3 0 0 2 0A=
 0 0 −2 0 0

0 0 0 0 00 0 0 0 −1
 0 0 0  0 1B=
 0
00



001000 0 0 0 10

フォーメーションフライトは従衛星を目標軌道に乗せることを目的とする. 誤差 e=x − x¯f は以下を満たす.

ė = Ae + B(ū + h(x̄) − h(x¯f ))0 0 0 
01ū = −Ke − h(x̄) + h(x¯f )(18)ė = (A − BK)e(19)を用いるととなる.



従って A − BK が安定であれば誤差は 0 に収束する.
フィードバックゲイン K は最適レギュレータのリッカチ
方程式

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

A は, 主衛星の軌道が円軌道であるときのシステム行列である. なお線形化せず非線形項でシミュレーションを行う場合は, ここからさらに (7) と (8) の差を出す.

(7) と (8) の差を f (x̄) とすると,f (x̄) はf (x̄) = (17)

ここで, フィードバック[(θ̇)2 + R23 − 3]x¯1 + θ̄¨x¯2 + [2θ̄˙ − 2]x¯40¨˙2˙1h(x̄) =  θ̄x¯1 + [(θ̄ ) + R03 ]x¯2 + [2 − 2(θ̄)]x¯3 1(1 − R3 )x¯5



(16)

− R2¯3 + R1¯2 − R1¯3 (x̄ + R̄0 )00( R1¯3 − R1¯3 )ȳ0
( R1¯3 − R1¯3 )z̄0(20)を用いて設計する
なお非線形の場合は (17) と (18) がそれぞれė = Ae + B(ū + h(x̄) − h(x¯f ) + f (x̄))

ū = −Ke − h(x̄) + h(x¯f ) − f (x̄)

となる。と表せられる.f (x̄) を (12) に加えると非線形項を含んだ状態方程式は4

シミュレーション
今回非線形の方程式に非線形項を使いフォーメーション設計したものと, 線形化された TH 方程式をフォーメーx̄˙ = Ax̄ + B(h(x̄) + f (x̄)) + B ū,x(0) = x¯0 (13)

ション設計したものを, 離心率 e = 0.2 として比較したところ, 図 3 と図 4 のようにほぼ同じ値を示した, これよりとなる.
主衛星の初期位置が近地点の時、従衛星の初期軌道および目標軌道の初期位置を以下のように置く.

x˜0 = [0.01 0 0

用性があると言えるので, 今回は線形化された TH 方程式のフォーメーション設計でシュミレーションを行う.

−0.01(2 + e)
√0.01 0](1 − e) 1 − e2−0.005(2 + e)
√0.005 0](1 − e) 1 − e2x˜f0 = [0.005 0 0

線形化された TH 方程式でのシミュレーション結果には有(14)また主衛星の初期位置を遠地点としたときは

x˜0 = [−0.01 0 0x˜f0 = [−0.005 0 00.01(2 − e)√
− 0.01 0](1 + e) 1 − e20.005(2 − e)√− 0.005 0](1 + e) 1 − e2

図 3 L1 ノルム比較(15)

図 4 Lz1 ノルム比較

本研究では, フィードバックゲインK の値は最適レギ
ュレータの重み行列をR = I, Q = 10−q I とおき, q を変化させて設計する.

図 5(面内運動の入力の L1 ノルム), 図 6(面外運動の入力の L1 ノルム), 図 7(整定時間 ST ), 図 8(面内運動の入力のL1 ノルムと整定時間 ST ) では離心率 e = 0.2 において,初期位置を近地点とした場合と遠地点とした場合を比較したものである.
図 7 は初期軌道から目標軌道に収束するまでの時間のグラフであり, 整定時間はほぼ同じになった. 図 8 は面内運動の L1 ノルムと整定時間 ST の関係をグラフで表したもので,L1 ノルムは ST と共に減少しやがて ST と共に増加していく.L1 ノルムが最小のとき整定時間は過度に増加し
図 5 面内運動の L1 ノルムていないためここが最適な値と言える.
図 5, 図 6 において, フィードバックゲインは HCW 方程式（主衛星が円軌道の場合の運動方程式）のリッカチ方程式を用いているため e > 0 のとき, ノルムは q と共に減少するが, やがて q と共に増加する. このことから離心率 e に対して, 最小の L1 ノルムを与える q が必ず存在する. 面外運動と面内運動の最小の L1 ノルムを見るとどちらも主衛星の初期位置が近地点より遠地点の方が小さいことが分かる. そのときの q の値（近地点 q = 0.625 遠地点q = 0.125）を入力した制御軌道を図 9 と図 10 に示した.

図 6 面外運動の L1 ノルム

図 9 制御軌道:離心率 e=0.2 図 10 制御軌道:離心率 e=0.2

5終わりに本研究を楕円軌道上のフォーメーションというテーマで進めてきて, 線形化した TH 方程式を使って設計した
図 7 整定時間 ST

フィードバック制御で非線形のシステムを制御することができるということが分かった. さらにその制御での重み qの最適な値を求め地球周回軌道上の主衛星の近くを飛行する従衛星を初期軌道から目標軌道にのせるとき, 制御を開始する初期位置を近地点とするよりも遠地点にした方が燃料効率が良いということが分かった.

参考文献[1] A Ichikawa:Relative motion along an elliptic orbit 講義プリント

図 8 L1 ノルムと整定時間 ST[2] K.Yamanaka and Ankerson,Newstate transi tion
matrix for relative motion an orbitary ellipticalorbit,J.Guidance,Control,Dynamics,vol.25,2002,99.6066.