フィードバックによる楕円軌道上のフォーメーションフライト
M2013SC010 南 佳那
指導教員：市川1

はじめにフォーメーションフライトとは, 複数の宇宙機編隊飛行によりミッションを達成する方法であり, その基礎研究として地球周回軌道上の主衛星とその付近を飛行する従衛星の相対運動が注目されている. フォーメーションフライ
トでは, 従衛星が適切な相対位置関係を維持することが望ましい. 楕円軌道上の主衛星に対する従衛星の相対運動方
程式を原点において線形化したものを Tschauner-Hempel
方程式 (TH 方程式) と呼ぶ. 軌道面内運動はある条件を満たすと周期的となる [1]. 状態空間表現の行列 A は時変である. 本研究では, 楕円軌道上の主衛星とその近傍にある従衛星のフォーメーション形成についての問題を扱い,最少燃料かつ目標軌道への収束が速いフィードバックゲインの設計法を提案する. フィードバックゲインは, 最適レギュレータ理論を用いて設計する. フォーメーション再形成を行なうためのフィードバック制御は, リッカチ微分方程式及び特異リッカチ微分方程式の 2 種類の方程式の周期解を用いて行い, 比較を行う. リッカチ微分方程式
は評価関数の入力への重みを変化させ, 特異リッカチ微分方程式は指数関数のパラメータを変化させる. 2 つの設計法により生成したフィードバックゲインの性能は, 総速度変化 (L1 ノルム) と軌道が収束するまでの時間である整定時間により評価する. また, 主衛星の初期位置の変化に伴う従衛星の消費燃料についても考察を行う.

2

楕円軌道上の相対運動

主衛星の軌道は, R0 = p/(1 + e cos θ) で与えられる. こ
こで, R0 は動径, θ は真近点離角, p = A0 (1 − e2 ) は半直弦, A0 は長半径, e は離心率とする. 軌道周期は T =2π(A30 /µ)1/2 で与えられ, µ は地球の重力定数である. nは, 角速度 θ̇ の平均とする. 次の方程式µR̈0 − R0 θ̇2 = − 2R0(1)R0 θ̈ + 2Ṙ0 θ̇ = 0が 2 体問題における Newton の運動方程式より導かれる[2][3]. 動径 R0 の楕円軌道上主衛星とその近傍の従衛星の相対運動 (図 1) を考えるため, 主衛星の重心を原点とする回転座標系 o − {i, j, k} を考える.

朗

ここで, o は主衛星の質量中心, i は動径方向の単位ベクト
ル, j は飛行方向, k は軌道面外への単位ベクトルである.
主衛星から見た従衛星の位置ベクトル r を r = xi+yj+zk
とおく. 従衛星の地球の質量中心からの位置ベクトルは R = R0 + r であるから Newton の運動方程式より|R| = R, |R0 | = R0 とすると,µµ= − 3 (x + R0 ) + ux
2R0Rµÿ + 2θ̇ẋ + θ̈x − θ̇2 y = − 3 y + uy
Rµz̈ = − 3 z + uz
Rẍ − 2θ̇ẏ − θ̈y − θ̇2 x −(2)が得られる. ここで, ux , uy , uz は宇宙機の制御加速度であり, R = [(R0 + x)2 + y 2 + z 2 ]1/2 である. 原点において線形化すると,µ)x = uxR03
µÿ + 2θ̇ẋ + θ̈x − (θ̇2 − 3 )y = uyR0µz̈ + 3 z = uz
R0ẍ − 2θ̇ẏ − θ̈y − (θ̇2 + 2(3)となり, この方程式は TH 方程式と呼ばれる. 式 (3) の状
態方程式は x = [x y ẋ ẏ z ż]T , u = [ux uy uz ]T とするとẋ = A(t)x + Bu, x(0) = x0である. ここで,
0
0 2µ θ̇ + 2 R30A(t) = 
−θ̈
0

000
θ̈

100012θ̇000θ̇2 − Rµ3−2θ̇000000000− Rµ30


00 
0 

0 
1 000 0 1
B=
 0
 00000100


00 0 0 
0 1で あ る.
時 間 t を τ = t/(1/n) に よ り 置 き 換 え(x̄(τ ), ȳ(τ ), z̄(τ )) = (1/A0 )(x(τ /n), y(τ /n), z(τ /n)) と置
くと式 (1), (2), (3) の無次元化方程式は

図 1 楕円軌道上における主衛星と従衛星

˙ 2=− 1
¨ − R̄ (θ̄)
R̄00R̄02R̄ θ̄¨ + 2R̄˙ θ̄˙ = 000

˙ 2 x̄ − 1 = − 1 (x̄ + R̄ ) + ū¨ − 2θ̄˙ȳ˙ − θ̄¨ȳ − (θ̄)
x̄
0x
R̄02
R̄3
˙ 2 ȳ = − 1 ȳ + ūȳ¨ + 2θ̄˙x̄˙ + θ̄¨x̄ − (θ̄)
yR̄31z̄¨ = − 3 z̄ + ūzR̄˙ 2 + 2 ]x̄ = ū¨ − 2θ̄˙ȳ˙ − θ̄¨ȳ − [(θ̄)
x̄xR̄03
˙ 2 − 1 ]ȳ = ūȳ¨ + 2θ̄˙x̄˙ + θ̄¨x̄ − [(θ̄)
yR̄031z̄¨ + 3 z̄ = ūzR̄03.1

リッカチ微分方程式 (RDE)

Q を状態への重み, R を入力への重みとし, 評価関数を,∫ ∞J(ū; x̄0 ) =(x̄T (τ )Qx̄(τ ) + ūT Rū(τ ))dτ(8)0とすると, これを最小化するゲイン K はリッカチ微分方
程式−Ẋ = AT X + XA + Q − XBR−1 B T X(4)(9)の周期解 X(τ ) を用いてと な る.
こ こ で,(ūx (τ ), ūy (τ ), ūz (τ ))=(1/A0 n2 )(ux (τ /n), uy (τ /n), uz (τ /n)), θ̄(τ ) = θ(τ /n) であり, 式 (4) の状態方程式は

˙ = A(τ )x̄ + B ū, x̄(0) = x̄0x̄(5)となる. この方程式に対するシステム行列は
0010000010 ˙ 2¨˙ (θ̄) + 23θ̄02θ̄0R̄0
A(τ ) = ˙¨˙120−θ̄(θ̄) − R̄3 −2θ̄ 0
0000000000 − R̄13K = R−1 B T X(τ )(10)で与えられる. 制御入力をこのゲイン K を用いて ū =−Kē とすると,ē˙ = (A − BK)ē(11)となり, 誤差を０に収束させることができる.

3.2 特異リッカチ微分方程式 (SRDE)00 
安定化の条件下で, 評価関数
0 
∫ ∞


J(ū;x̄)=exp(γτ )ūT (τ )Rū(τ )dτ
γ00 

01 
を最小化するフィードバックゲイン Kγ は0(12)0で あ る.
TH 方 程 式 は (x̃(θ), ỹ(θ), z̃(θ)) = (1 +e cos θ)(x, y, z) とおくことで Yamanaka と Ankersen により解かれている [4]. 軌道面内の運動が周期的となるための必要十分条件は˙ 0 ) + ρ ỹ(θ0 ) = 0(3ρ + e − 1)x̃(θ0 ) + esx̃(θ2˙2(6)で与えられる. ただし, ρ = 1 + e cos θ であり, s = ρ sin θである. この式 (6) を (x̃, ỹ) = A0 ρ(x̄(τ ), ȳ(τ )) と置き換えると

ρ (1 + ρ)x̄(τ ) − ρ e sin θȳ(τ )(7)3˙ ) + ρ(1 − e2 ) 32 ȳ(τ˙ )=0+e sin θ(1 − e2 ) 2 x̄(τ22を得る. 条件式 (7) は θ = 0 またはπ のとき簡単になる.
このとき式 (7) は各々,˙(2 + e)x̄(0) + (1 − e)(1 − e2 )1/2 ȳ(0)=0(2 − e)x̄(π) + (1 + e)(1 − e )2 1/2 ˙ȳ(π) = 0となる.

3フィードバック設計

システム (5) の周期解を目標軌道 x̄f とし, 初期軌道 x̄
から目標軌道 x̄f へフィードバック制御での移行を考える. 従衛星の制御軌道, 目標軌道の方程式を

˙ = A(τ )x̄ + B ūx̄˙ f = A(τ )x̄fx̄
とおくと, 誤差 ē = x̄ − x̄f は ē˙ = A(τ )ē + B ū, ē(0) = ē0を満たす. 安定化フィードバックゲイン K は最適レギュレータ理論を用いて設計する. ゲイン K の設計には, 2 種類の評価関数を最小とする方法を用いる.

Kγ = R−1 B T Xγ(13)で与えられる. この Xγ は特異リッカチ微分方程式−Ẋγ = AT Xγ + Xγ A + γXγ − Xγ BR−1 B T Xγ (14)を満たす周期解 Xγ (τ ) である [5].

4シミュレーション

状態方程式 (5) において, リッカチ微分方程式の入力への重み R と指数関数のパラメータ γ を変化させ, 2 つの制御器の性能を入力の L1 ノルム, 整定時間により比較する. 主衛星の初期位置 x̄l0 , 従衛星の初期軌道の初期位置
x̄0 は,e sin θ̄ (1 + e cos θ̄)2 T
1 − e2] (15)θ̄
13
1 + e cos θ̄ (1 − e2 ) 2(1 − e2 ) 2ce sin θ̄ + ċρx̄0 = 0.01[cos θ̄ − (1 + 1/ρ) sin θ̄3(1 − e2 ) 2
ρ(e − 2c) − (1 + 1/ρ)se sin θ̄ 1 + ecos θ̄3ρ(1 − e2 ) 2x̄l0 = [−sin θ̄
1(1 − e)(1 − e2 ) 2]T(16)とする. c は c = ρ cos θ であり, 従衛星の目標軌道の初期位置 x̄f 0 は, x̄f 0 = 1/2x̄0 とする. 主衛星が近地点, 遠地点にきた場合, 上記の θ̄ に θ̄ = 0, π を代入したものを用いる. θ = π/2, 3π/2 の場合は θ̄(τθ ) = θ̄ となる時間
τπ/2 , τ3π/2 を求め, そのときの θ̄ を用いて初期位置を決定した. 整定条件とする目標軌道との軌道誤差 ē は 10−4 以内とする. 異なるフィードバックゲインの設計法を比較する際の主衛星の軌道の離心率は, e = 0.3 とした. 総速度変化は面内と面外で分け, L1 ノルム, L1z ノルムとし,∫ ST
L1 = 0 (ū2x + ū2y )1/2 dτ∫ ST
L1z = 0 |ūz |dτとする. ただし, ST は整定時間である.

4.1

制御開始位置での比較

4.2

主衛星の初期位置が近地点 (θ = 0), 遠地点 (θ = π),
θ = π/2, θ = 3π/2 である場合の従衛星の L1 ノルム, L1zノルムの比較を, 2 つの設計法について行う.

4.2.1

RDE を用いた際の比較

図 3 は L1, L1z ノルムを各 θ の場合について比較したグラフであり, 図 4 は各 θ の場合について, 総ノルムを比較したものである.

RDE と SRDE の比較

RDE を用いて設計するフィードバックゲインは, 式 (9)を満たす周期解 X(τ ) を用いて式 (10) で与えられる. このときの評価関数は式 (8) である. 面内・面外の状態への重み Q を I4 , I2 に固定し, R を Rin = 10r I2 , Rout = 10rとして r を増加させていくことで入力の重みを変化させ,L1 ノルム, L1z ノルムの変化を調べる. SRDE を用いて設計するフィードバックゲインは, 式 (14) を満たす周期解 Xγ (τ ) を用いて式 (13) で与えられる. この時の評価関数は式 (12) である. 面内・面外の入力への重み R をRin = I2 , Rout = 1 に固定し, γ を γ = 10−ρ として ρ を増加させていき, L1 ノルム, L1z ノルムの変化を調べる.
主衛星及び従衛星の初期位置は式 (15), (16) に θ̄ = 0 を代入した値 (近地点) とする.

(a)L1 ノルム(b)L1z ノルム

図 3 整定時間-L1, L1z ノルム

図 4 総ノルムの比較 (RDE)(a) 面内運動(b) 面外運動

図 2 整定時間-L1, L1z ノルム (主衛星：近地点)
図 2 は, 上記の内容で RDE と SRDE の整定時間と L1 ノルム, L1z ノルムのグラフの比較を行っている. 総ノルムが 0.0150(86.4m/s) 以内の範囲で整定時間が最小となる
r, γ をとり, 整定時間の比較を行う. このとき r = 0.6250,
γ = 1.0000 である. 各整定時間は RDE が 11.55, SRDEが 6.48 である. 整定時間について有次元化すると, 主衛星が長半径 A0 = 12, 000km の楕円軌道を飛行するとしたとき, RDE は 6.7 時間, SRDE は 3.7 時間であり, 差は3.0 時間で SRDE の方が整定時間が短い. 以上から, 総ノルムが等しいとき, SRDE を用いてフィードバックゲインを設計した方が収束の速い制御器を設計可能であると言える. 遠地点 (θ̄ = π) においても同様の条件で比較を行う. このとき r = 0.8125, γ = 10−0.1250 である. 各整定時間は RDE が 15.06, SRDE が 9.47 である. 整定時間について有次元化すると, 主衛星が長半径 A0 = 12, 000km の楕円軌道を飛行するとしたとき, RDE は 8.7 時間, SRDEは 5.5 時間であり, 差は 3.2 時間で SRDE の方が整定時間が短い. よって, 総ノルムが等しいとき, 近地点同様遠地点でも SRDE を用いてフィードバックゲインを設計した方が収束の速い制御器を設計可能であると言える.

図 4 より, r が −1 から −0.1875 以内においては θ = 3π/2の場合が総ノルムが最も小さい. しかし, −0.1875 < r ≤1.3125 の場合には, θ = π/2 の場合が最も総ノルムが小さい. また, r > 1.3125 より大きくなると近地点の場合に総ノルムが小さくなり, θ = π のとき最小ノルムをとる.
よって, 用いる r の値により, 主衛星が θ = 0, π/2, π, 3π/2であるとき消費燃料が少ないと言える.

(a) 整定時間 (面内)(b) 整定時間 (面外)

図 5 整定時間
整定時間については, r の増加に伴い入力 ū が小さくなる
ため, 整定時間が大きくなることが分かる (図 5). 各初期位置の場合で, 面内・面外ともに θ の値による整定時間の差はほぼない. また, r が小さいときは面内の整定時間と面外の整定時間の差は小さいが, r の増加に伴い差が広がり, 面外の方が長くなる. ここで, 整定時間が 10 時間以内の範囲で総ノルムが最大となる r をとり, 総ノルムの比較を行うと, 表 1 となる. このとき主衛星は A0 = 12, 000kmの楕円軌道を飛行するとした. 表 1 より, 整定時間が 10時間以内と等しいときには θ = π/2 から制御を開始した方が良いと言える.
表 1 各主衛星の位置での総ノルム
主衛星の位置r
総ノルム 有次元 (m/s)θ=0θ = π/2θ=πθ = 3π/2
4.2.2

0.93751.06250.93750.93750.01380.01300.01440.0139

80.0
75.0
83.0
80.1

主衛星の初期位置に依存せずほぼ同じ値をとる. また, 面外の整定時間よりも面内の整定時間が長く, ρ の増加に伴
い面内と面外の整定時間の差が大きくなる. ここで, 整定時間が 10 時間以内の範囲で総ノルムが最大となる ρ をとり, 総ノルムの比較を行うと, 表 2 となる. このとき主衛星は A0 = 12, 000km の楕円軌道を飛行するとした. 表 2より, 整定時間が 10 時間以内と等しいときには近地点から制御を開始した方が良いと言える.
表 2 各主衛星の位置での総ノルム
主衛星の位置ρ
総ノルム 有次元 (m/s)θ=0θ = π/2θ=πθ = 3π/2

SRDE を用いた際の比較

図 6 は L1, L1z ノルムを各 θ の場合について比較したグラフであり, 図 7 は各 θ の場合について, 総ノルムを比較したものである.

(a)L1 ノルム(b)L1z ノルム

図 6 r-L1, L1z ノルム

50.4375
0.37500.43750.43750.0123
0.01280.01260.0124

70.9
73.8
72.6
71.5

おわりに

本研究では, 楕円軌道上の漸近的フォーメーション形成問題を考察した. 主衛星の軌道が円軌道である場合 HillClohessy-Wiltshire 方程式と同様に, 楕円軌道である場合
TH 方程式も特異リッカチ微分方程式で設計したフィードバックの方が, 等しい総ノルムを用いる場合, 整定時間が短いため, リッカチ微分方程式で設計したフィードバックよりも整定時間の面では性能がよいことが分かった. また, 主衛星が θ = 0, π/2, π, 3π/2 から制御を開始した場合の L1 ノルム, L1z ノルムを比較した結果, r, γ の値により最小の総ノルムをとる場合は変化すると分かった. また, 例えば整定時間が 10 時間以内としたときには, 特異リッカチ微分方程式で設計したフィードバックゲインを用い, かつ主衛星が θ = 0 にきた場合に従衛星の制御を開始する場合に総ノルムが最小となることが分かった. また, リッカチ微分方程式・特異リッカチ微分方程式において総ノルムを最小にする場合は θ = π, 0 の場合であった.

参考文献
図 7 総ノルムの比較 (SRDE)
図 7 より, ρ の値により総ノルムが最小となるときの θ の値は変化すると分かる. −0.6875 ≤ ρ < 0.1875 においては θ = 3π/2 の場合の総ノルムが最も小さく, それ以外の範囲では近地点の場合が最も小さい.

[1] Bando, M., and Ichikawa, A. : Graphical Generationof Periodic Orbits of Tschauner-Hemple Equations,
Journal of Guidance, Control, and Dynamics, Vol.
35, No. 3, 2012, pp. 1002-1007.
[2] Wie, B. : Space Vehicle Dynamics and Control,
AIAA, Reston, VA, 1998, pp. 282-285.
[3] Alfriend, K.T., Vadali, S.R., Gurfil, P., How, J.P.,
and Breger, L.S. : Spacecraft Formation Flying, Dynamics, Control and Navigation, Elsevier, Amsterdam, 2010, pp. 103-121.

(a) 整定時間 (面内)(b) 整定時間 (面外)

図 8 整定時間
整定時間については, ρ の増加に伴い整定時間が大きくなることが分かる (図 8). 各位置の場合で, 面外については[4] Yamanaka, K., and Ankersen, F. : New state transition matrix for relative motion on an arbitrary elliptical orbit, Journal of Guidance, Control, and Dynamics, Vol. 25, No. 1, 2002, pp. 60-66.
[5] Zhou, B., Duan, G., and Lin, Z. : A Parametric Lyapunov Equation Approach to the Design of Low GainFeedback, IEEE Transaction on Automatic Control,Vol. 52, No. 6, 2008, pp. 1548-1554.