出力レギュレーションによるアクティブフォーメーション
2009SE032 後久慎治

2009SE205 西尾壌

指導教員：市川1

はじめに人工衛星は現在の生活に欠かせない通信, 測位, 地球観測などで運用されているとともに, 将来の宇宙開発や技術開発を促進するために惑星探査衛星, 技術試験衛星など様々な衛星が開発, 運用されている. それぞれの人工衛星にはそれぞれの軌道があり, 楕円軌道, 双曲線軌道などがある. もしこれらの軌道上で衛星に故障が起きた場合, 故障した衛星を修復するために, この衛星を目標として別の宇宙機を接近させなければならないこともある. 本研究では, 円軌道上の主衛星に対して, 出力レギュレーション理論を用い, 従衛星を周期の短い目標相対軌道に乗せるフィードバックを設計する. また, 角速度を大きくしていったときの変化を考察する. ここでの評価関数は燃料消費を表す総速度変化とする.

2

朗が 得 ら れ る. こ の 方 程 式 は Hill-ClohessyW iltshire(HCW ) 方 程 式 と よ ば れ る [1],[2],[3]. 推 力を u = 0, 初期値を [x0 y0 ẋ0 ẏ0 ]T とする軌道面内運動の解は

2ẏ0
3x0 + 2ẏ0ẋ0−cos nt +sin nt,nnn
2ẋ0
2ẋ0
6x0 + 4ẏ0y(t) = y0 −+cos nt +sin nt
nnn− (6nx0 + 3ẏ0 )t,x(t) = 4x0 +ẋ(t) = ẋ0 cos nt + (3nx0 + 2ẏ0 ) sin nt,ẏ(t) = (6nx0 + 4ẏ0 ) cos nt − 2ẋ0 sin nt − (6nx0 + 3ẏ0 )となる. また, 初期値を [z0 ż0 ]T とする軌道面外運動の解はż0sin nt,nż(t) = −nz0 sin nt + ż0 cos ntz(t) = z0 cos nt +相対運動の方程式

半径 R0 の円軌道上の主衛星とその近傍の従衛星の相対運動を考えるため, 主衛星の重心を原点とする図 1 の回転座標系 o − {i, j, k} を考える.
このとき相対位置ベクトルを r = xi + yj + zk とする

x(t) = 2c + a cos(nt + α),y(t) = d − 3nct − 2a sin(nt + α),

rxoRz(t) = b cos(nt + β)
Leaderなる. ここで[R0Center of EarthCircular Orbit

図 1 円軌道上の主衛星と, ニュートンの運動方程式を変形しそれぞれ i, j, k について係数比較すると

ẍ = 2nẏ + n2 (R0 + x) −z̈ = −µz + uz
R3µ(R0 + x) + ux ,R3µy + uy ,R3である. 面内運動は c = 2x0 + ẏn0 = 0 のとき周期解となり, この条件を CW 条件という. x = (x, y, ẋ, ẏ, z, ż) とお
くと式 (2) の状態方程式はẋ = Ax + Bu, x(0) = x0
ここでẍ − 2nẏ − 3n2 x = ux ,ÿ + 2nẋ = uy ,]1
3x0 + 2ẏ0 2ẋ0 2 2ẏ0a= () +( ), c = 2x0 + ,nnn2ẋ0ẋ0d = y0 −, sin α = − ,nna
3x0 + 2ẏ0ż0 1cos α = −, b = [z02 + ( )2 ] 2 ,na
nz0ż0cos β = , sin β = −(6)bnb(1)となる. ここで [ux uy uz ]T = u は従衛星に働く推力,R =[(R0 + x)2 + y 2 + z 2 ]1/2 である. この方程式を原点 x =y = z = 0 で線形化すると

z̈ + n2 z = uz

(5)ijÿ = −2nẋ + n2 y −(4)となる. これは周期解となる. 上で求めた解 (3),(4) をパラメータ表現するとFollowery(3)(2)0 0 3n2
A=
 0

 0001000010002n00 −2n 000000000 −n2(7)00
 00 

 10 
, B =  00 

 01 


00 

0 

0 

0 01000
100 0

パラメータ (a, c, d, α; b, β) により表される状態方程式 (7)の解を γ H = (a, c, d, α; b, β) と表す.c = 0 のときこの解

は周期軌道となり γ H = (a, d; b) と表す. フィードバック制御による通常のフォーメーション形成問題とは, 状態方程式 (7) の解を与えられた周期軌道 γfH = (af , df ; bf ) に漸近的に追従させることである. このときの評価関数は,制御に使う推力の絶対積分であり, これは消費燃料に比例する.

3とすればよい. このとき[ẇ1 = S1 w1 ,[0−2ω

S1 =w1 (0) =]ω2]
a0,(10)0

出力レギュレーション理論を用いたフィー となり, この S1 を用いるとレギュレーター方程式 (9) の解はドバックの設計

状態方程式 (7) を無次元化するときのシステムを考える. 軌道面内運動のシステムは
0 00 0 1 0
 0 0  0 0 0 1 ẋ = u
x + 
 1 0  3 0 0 2 0 10 0 −2 0
また, 軌道面外運動のシステムは][ ][00 1uẋ =x+−1 01である. 軌道面内運動のシステムにおいて可制御性行列
Mc は rankMc = 4 であり可制御である. 同様に, 軌道面外運動においても rankMc = 2 となり可制御である.
以下の一般的システムを考える.

ẇ = Sw
とおくと, 上の式の S は S の固有値がすべて ０ 以上のとき反安定である. 出力レギュレーション問題は, ある外生信号に対し, 制御対象のシステムの出力を漸近的に追従する追従制御であり, 評価用出力を 0 に収束させる問題である. この場合, 任意の x(0), w(0) に対して,z を 0 に収束させる. この問題が解けるための必要十分条件は,２ つの行列 P1 , P2 が存在し, 次の等式を満たすことである [4].

AP1 − P1 S + B1 + B2 P2 = 0,
C1 P1 + D1 + D2 P2 = 0[P21 =−3 + 4ω − ω 20(9)(11)]0ω − ω2(12)となる. 次に, 軌道面外運動について考える. 軌道面外運動では z は w2 = (a1 cosωt, −a1 ωsinωt) を目標軌道とする追従問題として考える. ここで a1 は任意の正数である.
このとき
[]
[]
B1 = D2 = 0, C1 = 1 0 0 0 , D1 = − 1 0とすればよい. このとき

S2 =[w2 (0) =[(8)上の u は入力,z は評価用出力,w は外生信号で, 反安定の外部システムによって生成される. 外部システムを
01 ω ,2 01 0
P11 =  0−2ωẇ2 = S2 w2 ,ẋ = Ax + B1 w + B2 u,z = C1 x + D1 w + D2 u]0−ω 210]a10,(13)となり, これを用いるとレギュレーター方程式 (9) の解は[]
1 0P12 =(14)
0 1[P22 =−ω 2 + 1]
0(15)となる. x = [x y ẋ ẏ z ż]T ,u = [ux uy uz ]T とおくと状態方程式 (7) を無次元化したシステムは, 円軌道の半径およびその角速度を用いて
0 0 1 0 0 00 0 0
 0 0 0 1 0 0 
 0 0 0 
 3 0 0 2 0 0  1 0 0 ẋ = 
x +  u, x(0) = x0
 0 0 −2 0 0 0 
 0 1 0 
 0 0 0 0 0 1  0 0 0 

この時 u = Kx + (P1 − KP2 )w により出力レギュレーションが達成される. ここで,K は A + BK が安定となる
0 0 0 0 −1 00 0 1フィードバックゲインである.
軌 道 面 内 運 動 で は (x, y) は そ れ ぞ れ w1= となる. このとき w(0) は w1 (0) と w2 (0), また S は式
(acosωt, −2asinωt) を 目 標 軌 道 と す る 追 従 問 題 に つ (10),(13) によりいて考える. ここで a は任意の正数である. このとき
ω
00a0
[]
[]
2
 0 
 −2ω 01 0 0 01 000 

B1 = D2 = 0, C1 =, D1 = −w(0) = 
S=,0 1 0 00 1
 a1 
 0001 00 −ω 2 00

となり, レギュレーター方程式 (9) の解は式 (11),(14) と(12),(15) によって1
 0 0
P1 = 
 −2ω
 00
−3 + 4ω − ω 2P2 = 0001
ω2000xxf
8

00 
0 

,0 

0 10000100ω − ω2000−ω 2 + 164AT X + XA + Q − XBR−1 B T X = 00-4-605101520253035time

00 0

図 3 x と xf(16)の解を用いて

K = R−1 B T Xとする. また, 出力レギュレーションのフィードバックは

u = Kx + (P2 − KP1 )w

2-2となる. 状態フィードバックゲインはリッカチ方程式(17)となる.

4x and xfx and xf

-5x 1010

通常のフォーメーション形成問題の x 座標の制御軌道と目標軌道 xf の時間変化を描いた結果は図 3 のようになる. 図上の破線が制御軌道, 実線が目標軌道である. 無次元化により, 軌道の周期は 2π となる.
次に第 3 章の出力レギュレーション理論によって得られたフィードバックを用い, 衛星の軌道のシミュレーションを行う. 式 (16) の重み Q,R は Q = I, R = I(I = 単位行列) とおく. このときの制御軌道の初期値と目標軌道の初期値は, 通常のフォーメーション形成問題のシミュレーションと同じ値とし, また 5 周期以内に収束させる.
第 3 章で求めた S,P1 ,P2 の式に ω = 1 を代入すると, 衛星の軌道は図 4 のようになる.

シミュレーション

3D motion
controlled
final
initial-4

x 101z

第 2 章の HCW 方程式を用いた通常のフォーメーシ
ョン形成問題によって得られた衛星の軌道のシミュレーションを行う. 式 (16) の重み Q,R は Q = I, R =I(I = 単位行列) とおく. このときの制御軌道の初期値は,[x y ẋ ẏ z ż]T = [0.0001 0 0 − 0.0002 0.00005 0]とし, 目標軌道を初期軌道の 21 倍の軌道とし,5 周期以内に収束させる. 目標軌道の初期値は,[x y ẋ ẏ z ż]T =[0.00005 0 0 − 0.0001 0.000025 0] とする. 地球の半径を6387Km, 主衛星の地球の中心からの距離を 6787Km とする. 主衛星と従衛星の距離は, 制御軌道, 目標軌道それぞれの初期値に 6787Km を掛けた値となる. したがって, 初期軌道では,x = 0.6787Km,y = 0,z = 0.3362Km となる.
通常のフォーメーション形成問題によって得られた衛星の軌道は図 2 となる. 図の破線は初期軌道で, 実線は制御軌道で, 点線は目標軌道である.

20-1-221-2-40x 10-40x 10-1x-22y

図 4 ω = 1 のときの 3D モーション

ω = 1 の x 座標の制御軌道と目標軌道 w1 の時間変化を描いた結果は図 5 のようになる. この図から x 座標は-5x and w1x 1010xw1864controlled
finalinitial-4x 10x and w13D motion

202-21z-40-60510

1520253035-1time-2
21-2-40x 10-40x 10-1x-22

図 5 x と w1 のグラフy

図 2 3D モーション

周期 2π の軌道を描くことが分かる. また y 座標,z 座標の軌道の周期も 2π となる. よって, 出力レギュレーション理論を使い, 通常のフォーメーション形成問題と同じ軌道を描くことができる.
ω = 2 のときの x 座標の制御軌道と目標軌道 w2 の時間変化を描いた結果は図 6 のようになる.

図 8 は r の値によらず一定の値をとる.
出力レギュレーションのフィードバックである

u = Kx + (P2 − KP1 )w(18)

u = K(x − P1 w) + P2 w

(19)

-5x and w1x 1010を変形して

xw1
86x and w1

420-2-4-605

101520253035time

図 6 x と w1
図 5, 図 6 より, 目標軌道の周期は 2πω となり, 制御軌道はそれに追従するので, 出力レギュレーション理論を用いると, 周期の短い目標相対軌道に乗せることができる.

5とする. この 1 周期の L1 ノルムを求めるシミュレーションでは, 目標軌道と制御軌道は同じ軌道に乗っている. よって, 式 (19) の x − P1 w = 0 となり, P2 w だけが残るので,r の値によらず一定の値をとる.
ω = 1 から ω = 5 に対する L1 ノルムの値は次の表のようになる.

ωL1 ノルム
L1 ノルム (m/s)

10020.000859
0.068230.00239
0.189

40.004910.389

表 1 ω と L1 ノルム

1 周期の燃料消費

次に ω の値を変えた時の 1 周期にかかる L1 ノルムの値を調べる. このとき目標軌道と制御軌道は同じ軌道上を動いている.
このときの制御軌道の初期値は,[x y ẋ ẏ z ż]T =[0.00005 0 0 − 0.0001 ∗ ω 0.000025 0] となる.
ω = 2 のとき, 衛星の軌道は図 7 のようになる.

ω = 1 のとき,L1 ノルムの値が 0 となるのは,P2 の式にω = 1 を代入すると,P2 の値が 0 となるからである. 表 1から ω の値が大きくなるほど, 衛星の 1 周期にかかる燃料の消費が増えることが分かる.

6

おわりに

3D motion

円軌道上の主衛星に対して, 出力レギュレーション理論を用い, 従衛星を周期の短い目標相対軌道に乗せるフィードバックを設計することができた. また通常のフォーメーション形成問題と同じ軌道を描くことができる. また ω の値を大きくしていくと, 衛星の周期は短くなっていく. さらに ω の値を大きくすればするほど, 主衛星の総速度変化の値が大きくなり, より燃料を消費することがわかった.
この研究により, 円軌道上にある衛星や宇宙機などの衝突を回避したり, 宇宙機同士を接近させることができる.

controlled
final
initial-4x 10z

10-1
1-1-40x 10-40x 10x-11y

図 7 ω = 2 の 3D モーション
図 7 から, 目標軌道と制御軌道は同じ軌道上にあることが分かる.
ω = 2 のときの L1 ノルムは, 図 8 のようになる.
-4x 108.59488.59488.59488.5948[1] A. Ichikawa: Recent Developments in Formation
Flying, Lecture Notes ver.1, 2010.
[2] M. Shibata and A. Ichikawa: Orbital Rendezvous
and Flyaround Based on Null controllability withVanishing Energy, Journal of Guidance, Control,and Dynamics, Vol. 30, No. 4, pp. 934-945, 2007.
[3] B. Wie: Space Vehicle Dynamics and Control, AIAA
Education Series, Reston, Virginia, 1998.

8.5948
L1 ノルム参考文献

8.59488.5948[4] Saberi, A.,Stoorvogel, A.A., and Sannuti, P.S.,”Control of Linear Systems with Regulation and Input Constraints,”Springer,2000.

8.59488.59488.59488.5948-2-1.5-1-0.5

00.5

11.5

2r

図 8 ω = 2 のときの L1 ノルム

50.00840
0.667