最大推力を考慮した円軌道上のフォーメーション
2009SE023 藤原基志

2009SE260 杉山慎太郎

指導教員：市川 朗

1

はじめに運動方程式より,複数の人工衛星が編隊飛行する技術をフォーメーション
フライトという. この技術は, 高精度かつ立体的な観測を行うことができるが, 総燃料の増加や軌道制御の複雑化といった問題がある. また, それぞれの人工衛星には出力できる推力の限界があるため, 最大推力を考慮する必要がある. そのような問題の解決のために, 本研究では燃費や整定時間, 最大推力に条件を与え, その条件を満たすようなシステムの設計を行う. 制御対象を収束させる目標の軌道を目標軌道と呼ぶ. ここでは, 周期解を持ち制御対象に入力を与えずとも軌道上を周り続けることが可能な軌道を目標軌道として Hill-Clohessy-Wiltshire(HCW) 方程式から導出する. この軌道上であれば, 宇宙の観測や通信など人工衛星本来の目的を, 消費燃料を気にすることなく遂行できる.HCW 方程式は主衛星とその近傍の従衛星を原点で線形化した方程式である.
更に目標軌道に制御対象を収束させるために, 状態フィードバックを設計しなければならない. その状態フィードバッ
クの設計方法として, 入力が過大とならないようにシステムを安定化させる最適レギュレータのリッカチ方程式を用いる.
以上の手法を用いて, 制御対象を条件を満たしつつ初期軌道から目標軌道へ収束させるシミュレーションを行う.

2相対運動の方程式

主衛星の軌道を半径 R0 の円軌道とし, その近傍の従衛星の相対運動を考えるために, 主衛星の重心を原点とする.
図 1 の回転座標系 o − {i, j, k} を考える.

ẍ = 2nẏ + n2 (R0 + x) −ÿ = −2nẋ + n2 y −z̈ = −µz + uz
R3µ(R0 + x) + ux
R3µy + uy
R3(1)が得られる [1, 2]. ここで [ux uy uz ]T = u は単位質量あたりの従衛星に働く推力,R = [(R0 + x)2 + y 2 + z 2 ]1/2 である. この方程式を原点 x = y = z = 0 で線形化すると

ẍ − 2nẏ − 3n2 x = uxÿ + 2nẋ = uy(2)z̈ + n2 z = uzが得られる. この方程式が HCW 方程式とよばれる. 推力を u = 0, 初期値を [x0 y0 ẋ0 ẏ0 ]T とする面内運動の解は

x(t) = 4x0 + 2ẏ0 /n − (3x0 + 2ẏ0 /n) cos nt + (ẋ0 /n) sin nty(t) = y0 − 2ẋ0 /n + (2ẋ0 /n) cos nt + (6x0 + 4ẏ0 /n) sin nt− (6nx0 + 3ẏ0 )tẋ(t) = ẋ0 cos nt + (3nx0 + 2ẏ0 ) sin ntẏ(t) = (6nx0 + 4ẏ0 ) cos nt − 2ẋ0 sin nt − (6nx0 + 3ẏ0 )(3)で与えられる. 初期値を [z0 ż0 ]T とする面外運動の解はz(t) = z0 cos nt + (ż0 /n) sin nt(4)ż(t) = −nz0 sin nt + ż0 cos nt(5)Followerryであり, 周期解となる. 上の解をパラメータ表現すると

xijoRx(t) = 2c + a cos(nt + α)
Leadery(t) = d − 3nct − 2a sin(nt + α)R0(6)z(t) = b cos(nt + β)Center of Earth

ここで[]1/2
a = (3x0 + 2ẏ0 /n)2 + (ẋ0 /n)2
Circular Orbitc = 2x0 + ẏ0 /n
d = y0 − 2ẋ0 /n, sin α = −ẋ0 /na

図 1 円軌道上の主衛星
このとき n は角速度,µ =GMとし, 相対位置ベクトルをr=xi+yj+zkとすると,

cos α = − (3x0 + 2ẏ0 /n) /a,(7)b = [z02 + (ż0 /n)2 ]1/2

cos β = z0 /b, sin β = −ż0 /nb
面内運動は c = 2x0 + ẏ0 /n = 0 のとき周期解となり, この条件を CW 条件という.
状態方程式はẋ = Ax + Bu, x(0) = x0(8)ここで0 0 3n2A= 0 00010001002n0 −2n 0000000

00000−n2

00
 00  10 ,B =  00  01 00であり,
0001000st = 5.4375 ∗ 104 [sec] = 15.104[h]0 
L1 = 46.636[m/sec]
0 
0 となり,10 周期の整定時間である 55653[sec] に近い値とな0 ることから, 約 10 周期で収束していることがわかる.
1 また r の増加に伴い,u max,u max,u max は単調減少
xyzとする. このシステムを安定化させる最適レギュレータの している.
リッカチ方程式は,
0.1

AT X + XA + Q − XBR−1 B T X = 0L1(9)0.09

によってフィードバック

u = −R−1 B T Xx

L1ノルム

0.08(10)0.07

0.06を求められる.

0.05

0.04

3

状態フィードバックの評価

0.03
5

5.5

6

6.5

−1

7

7.5

8

8.5

9r

初期軌道を [60.0 0 0 -1.3547*1030.0 0], 目標軌道−2を [6.0 0 0 − 1.3547 ∗ 10 1.0 0] と設定する.
従 衛 星 の 単 位 質 量 あ た り の 推 力 [ux uy uz ]T =u[km/sec2 ] の 絶 対 値 の 最 大 値 を 最 大 推 力
u max[km/sec2 ] とし, 推力の上限を ulim [km/sec2 ], 面内運動の燃費 L1 − ノルム を L1[km/sec] と表現し, その上限を L1lim [km/sec] として, ulim = 2.0 ∗ 10−5 ,L1lim = 5.0 ∗ 10−2 とする. 従衛星が目標軌道に収束するまでの時間を整定時間とし,st と表現する. u max と L1に上限を与え, 目標軌道に収束するまでの st は 15 周期である 8.3480 ∗ 104 [sec] = 23.189[h] 以内を目標とする.
ここで最適レギュレータのリッカチ方程式 (9) により, 設計したフィードバック (10) の重み行列を,R=10r ∗I,Q=10−7 ∗ I と設計する. r の値を変化させながら, それに伴う st と L1 と u max を計算し, 重み行列 R のパラメータを変化させ,st と L1 と u max を指定した値以内にする.

図 3 r の変化に伴う,L1 の推移
図 3 より,r が 6.1250 以上であれば,L1lim を下回る. rの増加に伴い,L1 が単調減少していることがわかる. また、同様に L2 − ノルム [m/sec] も単調減少し 0 に収束することが知られている [3, 4]. r = 6.150,q = −7.0 ときのu max はそれぞれ以下のようになる.

ux = 2.1234 ∗ 10−2 [m/sec2 ]uy = 1.0410 ∗ 10−2 [m/sec2 ]uz = 6.7012 ∗ 10−3 [m/sec2 ]であり,st = 5.1419 ∗ 104 [sec] = 14.283[h]L1 = 49.352[m/sec]-4

history of max ubsolute inputs
x 10となる.

1.2
mabs ux
mabs uymabs uzmax u

15x 1012st
st=83480
100.680.4st

max ubsolute input u

0.8

0.2

6

4
05

5.5

6

6.5

7

7.5

8

8.5

9

r
20

図 2 r の変化に伴う,u max の推移
図 2 より,r = 6.250 以上であれば,ulim を超えない. r =6.250,q = −7.0 ときの u max はそれぞれ以下のようになる.

ux = 1.6889 ∗ 10−2−3uy = 9.6336 ∗ 102[m/sec ][m/sec2 ]uz = 5.9420 ∗ 10−3 [m/sec2 ]

55.56

6.5

7

7.588.59r

図 4 r の変化に伴う,st の推移
図 4 より,r が 6.750 以下であれば, 人工衛星は 15 周期である 8.3480 ∗ 104 [sec] 以内に収束する. また,r の増加に伴い st は単調増加している.

r = 6.750,q = −7.0 のときの u max はそれぞれ
ux = 6.0175 ∗ 10−3 [m/sec2 ]uy = 7.2191 ∗ 10−3 [m/sec2 ]uz = 3.5584 ∗ 10−3 [m/sec2 ]であり,

st = 7.6454 ∗ 104 [sec] = 21.237[h]L1 = 40.288[m/sec]となる.
r の値を 5.0 から 0.1250 刻みで 9.0 まで増加させると, 増加に伴い L1 は小さくなり,st は大きくなった. 図 2, 図 3,
図 4 から r の値が,6.250 ≤ r ≤ 6.750 の範囲内であれば,u max,L1,st を全て与えた条件下で実現できる.
3D motion
controlled
initial
final

6040z

200-20-40-60

50-50
00-5050xy

図 5 制御軌道 (r=6.750)

3D motion
controlled
initial
final

60

衛星が目標軌道に収束できた事も, 図 6 から見てとれる.
このようにして, 図 2, 図 3, 図 4 を見比べ, r の範囲を導出することができ, 定めた条件下でフィードバック制御を設計することが可能となる.
また, 本節では, 初期軌道から目標軌道に直接向かう制御について議論してきた. しかし, 直接的に目標軌道に向かうのではなく, 更に初期軌道と目標軌道の間の軌道を中継し段階的に目標軌道を目指した場合, 最大推力が抑えられるのではないかと考えたため, 次節にてそのシミュレーションを行う.

4

中継軌道の設定

今節では, 初期軌道と目標軌道の間に中継軌道を設ける.
初期軌道 [60.0 0 0 -0.13547 30.0 0], 中継軌道 [30.0 0 0-0.067736 1.0 0], 目標軌道 [6.0 0 0 -0.013547 1.0 0] と設定した.
前節で st や L1,u max をある程度, 与えられた範囲内で設計できた. しかし,st や L1,u max に与える条件によっては, 設計できない場合がある. 例えば, 前節の図 6 の制御設計にて st をそのままにして,u max を小さくする設計はできない. よって, 本節では初期軌道から目標軌道に人工衛星を収束させるまでの間に中継軌道を設定し, それを経由させることで推力を低くする制御を行う. 制御設計
は r = 6.25,q = −7.0 とし, 前節の図 6 の制御設計と比較して評価する.
フィードバック制御の入力は u = −Ke となる. 制御軌
道と目標軌道の誤差である e が大きいほど, 推力 u が大きくなることがわかる. 即ち誤差 e を小さくすれば, 重み行列を変化させなくても推力 u が小さくできる.
上記の方法を利用して, 推力を下げる制御を行う.
初期軌道から制御を開始して, 中継軌道に収束するまでのシミュレーション結果を図 7 に表す. 制御軌道と中継軌
道の誤差が 1000m を下回ったとき制御対象が収束したとみなす.

40z200
3D motion-20
controlled
initialfinal-40-606050
40-50020-50z0500xy-20-40-60

図 6 制御軌道 (r=6.250)

50-5000

図 5 は r=6.750 のとき, 図 6 は r=6.250 のときのシミュレーションである. 図 6 より,r=6.250 のとき約 10 周期で目標軌道に収束する. また,r = 6.250 と r = 6.750 の整定時間を比較すると r = 6.250 の時の方が r = 6.750 の時より, およそ 5 周期分短く収束していることがわかる.
L1 と st,u max に制限を与えた上で制御対象をフィードバック制御で目標軌道に乗せる事ができた. 以上で, 収束までの整定時間は 15 周期以内となり,ulim , L1lim を, それぞれ下回るようなフィードバック制御を設計できた.従-5050xy

図 7 制御軌道 (r=6.250,q=-7.0)このときの u max はそれぞれ,

ux = 9.3828 ∗ 10−3 [m/sec2 ]uy = 5.3326 ∗ 10−3 [m/sec2 ]uz = 5.9473 ∗ 10−3 [m/sec2 ]であり,st = 15202[sec] = 4.2227[h]L1 = 25.903[m/sec]となり, 上記の設計であれば ulim を下回る.
中継軌道に収束したとみなす制御軌道は,制御軌道 [-3.2828 59.521 0.033502 0.0074109 -0.90602
0.0078705] となり, この制御軌道が目標軌道に収束するまでを制御する. そのシミュレーション結果を図 8 に示す.
目標軌道に収束したとみなす条件は, 制御軌道と目標軌道の距離の誤差が 10m を下回った時とする.

比較して, 本節で行った段階を踏んで目標軌道に到達させる制御は u max を小さくでき,st はほとんど変わらなかった.L1 は図 6 の制御設計と比較して大きくなった.
初期軌道から中継軌道を経て目標軌道に到達したときに図 6 の制御設計と st が変わらないことや,L1 が大きくなることの理由は, 制御対象が複数の軌道を経由しているからである. よって L1 は大きくなるが, 中継軌道への収束条件を緩く設定しているため,st は短くなる. 例えば,
図 7 の制御軌道が中継軌道に収束したとみなす条件を,誤差が 10m 以下になったときとする. この場合の st は
st = 48613[sec],L1 = 25.903[m/sec] となり, 収束条件を厳しくすると整定時間は長くなる. 中継軌道は収束条件を厳密に設ける必要がないので条件を緩くできる.
以上の方法より, 段階的に目標軌道へ収束すれば,u maxを抑えつつ, 直接目標軌道に収束させた場合とほぼ整定時間を同じにできる設計が可能となった.st を長くすることなく, 最大推力を抑えたい場合にはこの設計を用いるとよい.

5

おわりに

3D motion
controlled
initialfinal
20z100-10-2050-200020x-50y

図 8 制御軌道 (r=6.250,q=-7.0)
この時の u max は,

ux = 5.6698 ∗ 10−3 [m/sec2 ]uy = 1.0263 ∗ 10−3 [m/sec2 ]uz = 1.6319 ∗ 10−3 [m/sec2 ]フォーメーションフライトは, 複数の人工衛星をフライトさせるので燃料消費量が増大する. そのためにも人工衛星の推力や収束するまでの時間を設計することで消費燃料を軽減させなければならない.
その問題に対して, 最適レギュレータのリッカチ方程式
よりフィードバック制御を設計し, そこで重み行列を変化させることによって, 設計者が望む推力や整定時間を制御対象に与えることができた. それも今回は r と st のグラフ,r と u のグラフを用いたので, 与えた条件を守るような制御を図からわかりやすく明示した.
また, 複数の軌道を経由させることによって, 推力を抑える設計も可能とし, 制約の多い制御設計の幅を広げることができた.
以上の方法をもってすれば, 制御対象の推力, ノルム,整定時間をある程度は設計者の意に沿った設計が可能となる.

参考文献[1] A. Ichikawa: Recent Developments in Formation
Flying, Lecture Notes ver.1, 2010.

であり,

st = 38984[sec] = 10.828[h]L1 = 26.920[m/sec]となり, 上記の設計であれば,ulim を下回る.
a から a1 へ収束する制御,a1 から af に収束する制御,いずれの最大推力も ulim を下回った. そして L1 と st の総和は,

st = 54186[sec] = 15.051[h]L1 = 52.823[m/sec]となった.
前節で r = 6.250,q = −7.0 で制御設計を行った場合と[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] Y. Ichimura and A. Ichikawa:Optimal Impulsive Relative Orbit Transfer Along a Circular Orbit, Journalof Guidance, Control, and Dynamics, Vol. 31, No. 4,
pp. 1014-1027, 2008.
[4] A. Ichikawa:Null Controllability with Vanishing Energy for Discrete-Time Systems, Systems & Control
Letters, Vol. 57, No. 1, pp. 34-38, 2008.