出力フィードバックによるフォーメーション
2009SE075 井上陽介

2009SE116 桔梗祐太
指導教員：市川1

はじめに2009SE256 下出拓弥朗

m を消去し，地球の重力定数 µ = GM を代入すると，µR̈ = − 3 R(2)Rとなる．x，y ，z 方向の単位ベクトルをそれぞれ i，j，k，主衛星から従衛星への位置ベクトルを r ，地球の中心から主衛星への位置ベクトルを R0 = R0 i，とすると，図 1 より

本研究では出力フィードバックを用いて円軌道上の主衛星とその近傍の従衛星のフォーメーションについて考える．円軌道上の主衛星とその近傍の従衛星の相対運動の方程式は，時不変非線形微分方程式で与えられる．原点で線形化した方程式は Hill-Clohessy-Wiltshire 方程式とr = xi + yj + zkよばれ，軌道面内の運動は楕円で表される周期解をもつ．
軌道面外の運動は単振動の方程式により表される．状態
R = R0 + r = (R0 + x)i + yj + zk
フィードバック制御を用いた人工衛星の制御は，全てのが導出できる．角速度を n とすると，単位ベクトル i，j，状態が観測可能であるという条件があるが，状態が一部kは不可能な場合もある．この場合，出力フィードバックをi̇ = nj
用いて制御を行う．
j̇ = −ni
人工衛星の操作量と観測量から状態変数を推定するオブザーバーを使った出力フィードバック制御を用いて，フォーメーション形成・再構成問題を定式化する．ここでの評価関数は燃料消費を表す総速度変化とする．10 周期程度で収束するフィードバックゲインとオブザーバゲインを求めて，状態フィードバックに近い制御を行うことを目標とする．なお，人工衛星が目標軌道に収束する条件をストッピングルールとする．

2

相対運動の方程式

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

FollowerryxijoRLeaderk̇ = 0と表される．R を時間 t で 2 回微分した式はR̈ = [ẍ − 2nẏ − n2 (R0 + x)]i + (ÿ + 2nẋ − n2 y)j + z̈k(3)となる．
ここで，(2) 式に従衛星に働く推力 u = [ux , uy , uz ]T を含めた次式µR̈ = − 3 R + u(4)Rを考える．推力のベクトル u = ux i + uy j + uz k として，式 (3) と式 (4) を単位ベクトル i，j，k で係数比較すると，µẍ = 2nẏ + n2 (R0 + x) − 3 (R0 + x) + ux
Rµÿ = −2nẋ + n2 y − 3 y + uy
(5)Rµz̈ = − 3 z + uz
Rが得られる．(5) 式を原点 x = y = z = 0 で線形化すると

R0ẍ − 2nẏ − 3n2 x = ux

Center of Earthÿ + 2nẋ = uy

(6)2z̈ + n z = uz

Circular Orbit

図 1 円軌道上の主衛星

地球と人工衛星の間に働く万有引力はmR̈ = −GM mR
R3(1)と表される．ここで，m は人工衛星の質量，M は地球の質量，G は万有引力定数，R は地球の中心から人工衛星(図 1 での従衛星) までの位置ベクトルを示す．(1) 式からとなる．(6) 式は Hill-Clohessy-Wiltshire(HCW) 方程式とよばれる．推力を u = 0, 初期値を [x0 y0 ẋ0 ẏ0 ]T とする面内運動の解は()
2ẏ02ẏ0ẋ0x(t) = 4x0 +− 3x0 +cos nt +sin ntnnn()
2ẋ02ẋ04ẏ0y(t) = y0 −+cos nt + 6x0 +sin ntnnn− (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 )(7)で与えられる [3][5]．初期値を [z0 ż0 ]T とする面外運動の解はż0sin ntnż(t) = −nz0 sin nt + ż0 cos ntz(t) = z0 cos nt +(8)であり，周期解となる．(7) 式，(8) 式の解をパラメータ表現するととする．このときの推定誤差を e = x − x̂ とすると

ė = (A − HC)e(12)となる．A − HC が安定であるとき，e(t) は漸近的に 0に収束する．このとき推定システムをオブザーバー (状態観測器) といい，H をオブザーバーゲインという．
ここで，u = −K x̂ として，式 (10) に代入すると

ẋ = Ax − BK x̂.

x(t) = 2c + a cos(nt + α),y(t) = d − 3nct − 2a sin(nt + α),(9)(13)x̂ = x − e より，式 (13) に代入して変形すると

z(t) = b cos(nt + β).

ẋ = (A − BK)x + BKe.

ここで，[]12ẏ0 2 ẋ0 2 2ż0 2 1a = (3x0 +) +, b = [z02 +]2nnnẏ02ẋ0c = 2x0 + , d = y0 −nn2ẏ03xẋ00+ n
sin α = − , cos α = −naaż0z0
sin β = − , cos β =nbb(14)式 (12) と式 (14) を拡大行列系で表すと))() ((xA − BKBKẋ=ėO
A − HCe

この系の固有値は σ(A − BK) ∪ σ(A − HC) であるから安定となる．出力フィードバックによって安定化を行
うには状態フィードバック K とオブザーバーゲイン H を独立に設計すればよい．これを安定化補償器 (制御器) 設
ẏ0より，面内運動は c = 2x0 + n = 0 のとき周期解となり， 計の分離原理という．オブザーバーと元のシステムの次元が同じであるとき同一次元オブザーバーという．
この条件を CW 条件という．次に状態フィードバックの安定化µ = 398600[N]，Re = 6.39 × 103 [km]3R0 = Re + 400[km]10 周期前後で収束するフィードバックゲイン K を求める．(9) 式にとおくと，n は√n=µ= 1.1 × 10−3 [km/s2 ]R03

a = 5, b = 1, c = 0, d = 0と表される [1][2]．
一般に状態空間表現はẋ = Ax + Bu, x(0) = x0(10)y = Cxで表される．(6) 式を，(10) 式で表すと行列 A,B,C は001000 000100  3n2 002n00 A= 00 −2n 000  000001 



00
 0
 1

B=
 0

 000000100

00−n20

0 

1 0

0 

,C =  0 10 
0 00 を代入して，システムの初期値 x0 を

5 cos α
−10 sin α
−3−5.50 × 10 sin α
x0 = 

−3−2
 −5.50 × 10 − 1.10 × 10 cos α 
−1.10 × 10−2 cos β−3−1.10 × 10 sin βとする．ここで，α は初期値を定める角度であり，このα を π2 ごとに変えていく．β は面外運動の座標を決める任意の角度である．目標軌道 xf の初期値は

a = 0.5, b = 00000000001

00 01となる．次に出力フィードバック制御について考える. 観
測値y に基づく状態の推定を考える．推定システムを元のシステムである (10) 式を模倣した動的システム

˙ = Ax̂ + Bu + H(y − C x̂)
x̂

(11)

とおき

12 cos α− sin α −5.50 × 10−4 sin α xf = 
 −1.10 × 10−3 cos α 00に設定する．ẋf の方程式は

ẋf = Axf , xf (0) = xf 0(15)となる．目標軌道とシステムの誤差を e1 = x − xf とすると ė1 = ẋ − ẋf となり，式 (10) の ẋ と式 (15) からė1 = Ae1 + Buとなる．制御入力uをu = −K(x − xf ) = −Ke1となった．α の角度がちょうど π だけずれており，反対側から目標軌道までの楕円運動を同じ距離で行うので，r, st, L1 ノルムが同じ値になった．
また，α = π2 及び α = 3π2 のときの L1 ノルムのシミュ
レーション結果を図 6,7 に，整定時間 st のシミュレーショ
ン結果を図 8,9 に表示する．
-3-3x 10

AT X + XA + Q − XBR−1 B T X = 0
K=R−111101099887(16)TB X

L1 ノルム

1211L1 ノルムとし，フォーメーションを行う．
A − BK を安定化させるため，最適レギュレータのリッカチ方程式を解く．

x 10127665544332244.555.5644.55r(16) 式はフィードバックゲイン K を与えるリッカチ方程式である．(A,B) が可制御のとき，安定となるフィードバックゲイン K を求めることができる．この時のフィードバックゲイン K の設計を Q = 10−7 I に固定し，操作量の重み R = 10r I の r を 4 から 8 まで 0.125 刻みで動かして設計する．ストッピングルールを「目標軌道と人工衛星との誤差が時刻刻みで 3 回連続 10−2 以下に収まったら収束した」と決める．α = 0 及び α = π のときの L1
ノルムのシミュレーション結果を図 2,3 に，整定時間 stのシミュレーション結果を図 4,5 に表示する．
-3121010L1 ノルム

L1 ノルム12856442244.55

5.5

66.5

7

7.5

8

1.5

1.5
st

2st

2

1

1

0.5

0.5

0

4.5

5

5.5

66.5

77.5

図 3 L1 ノルム (α = π)

55x 102

1.5

1.5
st

2st

2.5

1

1

0.5

0.5

00
4.5

5

5.5

6

6.5

7

7.5

84

4.5

r

5

5.5

66.5

77.58r

図 4 st(α = 0)

4.5

5

5.5

66.5

7

7.5

844.555.5r66.577.58r

図 9 st(α = 3π2 )8rx 1004st = 5.73 × 104 [s]42.5x 102.5

α = π2 と α = 3π2 のときのグラフもほぼ一致している．
10 周期前後で収束する st の値は図 8,9 からr

図 2 L1 ノルム (α = 0)5x 10

図 8 st(α = π2 )866

図 7 L1 ノルム (α = 3π2 )2.5x 10144

図 6 L1 ノルム (α = π2 )-3x 10145.5r

図 5 st(α = π)となり，r = 6.75 のときに得られた．図 6,7 から L1 ノルムは
L1 = 3.70[m/s]となった．楕円運動で π2 ずれると目標軌道までの距離が変わり，α = 0, π のときと比べて L1 ノルムの値も小さいので，α = π2 , 3π2 のときの方が距離が小さいと分かる．
よって，初期軌道から目標軌道までの面内運動の移動距
離は α = 0, π のときが最も長く，α = π2 , 3π2 のときが最も短いと言える．角度 α の各値から r を決定したことにより一意に Q と R が得られたので，(16) 式からフィードバックゲインは[L1 ノルムは燃料消費を表す評価関数であり，単調に減少していく [2][5]．r を大きくすることで，操作量が過大でないことを表す．4 つの図を見ると，α = 0 と α = πのときのグラフはほぼ一致していることがわかる．図 4,5から，10 周期近くで収束する st の値
st = 5.73 × 104 [s]となり r = 6.75 のときに得られた．図 2,3 から L1 ノル
ムは
L1 = 3.36[m/s]K =9.29e − 7
2.50e − 6
5.31e − 21−8.72e − 8−1.09e − 7−2.80e − 22

3.03e − 4
3.62e − 4−6.85e − 19

3.62e − 4
1.10
2.47e − 18

6.10e − 22
3.28e − 21
6.96e − 9−6.85e − 19
2.47e − 18
1.18e − 4と定まった．フィードバックゲインの各成分は非常に小さく 0 に近いものとなり，入力も非常に小さい値になっていることが分かる．角度 α が α = 0, π2 , π, 3π2 以外で値をとっても上記のフィードバックゲインをとる．

4

オブザーバーゲインの決定

10 周期近くで収束するオブザーバーゲイン H を決める．目標値は前章の値を使い，オブザーバーの初期値は]

5x 1022

1.5

1.5
st2.5

110.5

0.500234567892345r1−1.21 × 10−3 sin βe2 = x̂ − xfとすると，-3ė2 = Ae2 + Bu + H(y − C x̂)

-35

4.5

4.5

L1 ノルム

4L1 ノルムu = −K(x̂ − xf ) = −Ke2

3.5

3.5

33

2.5

2.5
2

345

67892

345r1として，フォーメーションを行う．A − HC を安定化させるためには，AT − C T H T を安定化させる最適レギュ
レータのリッカチ方程式を解けばよい．

図 14 L1 ノルム (α = π2 )CY = 0(17) 式は Y を解とする H T のリッカチ方程式で，安定となる H T から H を得られる．オブザーバーゲイン H の設計を Q1 = 10−7 I に固定し，操作量の重み R1 = 10r1 I のr1 を 2 から 9 まで動かして設計する．ストッピングルー
ルを「人工衛星と目標軌道との差が 10−2 以下に各時刻ごとに 3 回連続で入れば収束した」とみなす．L1 ノルムが4.00[m/s] 以下でかつ状態フィードバックの性能に近づけられたら良いオブザーバーと評価する．初期値の速度成分が 1 割増の場合において，α = 0, π のときの L1 ノルムのシミュレーション結果を図 10,11 に，整定時間 st のシミュレーション結果を図 12,13 に表示する．
-3-3x 101098877L1 ノルム

L1 ノルムx 109665544332223456789r1

図 10 L1 ノルム (α = 0)2349521.51.556789r1

図 11 L1 ノルム (α = π)st(17)CY108x 102st

H = R1−17

図 15 L1 ノルム (α = 3π2 )5T6r1x 10

AY + Y A + Q1 − Y C R19x 1054となる．出力フィードバック制御の入力を−18

図 13 st(α = π)x 10T7

α = 0 と α = π のグラフはほぼ一致する．図 12, 図 13より，r1 が 2.5 ≦ r1 ≦ 7.25 の範囲なら状態フィードバック制御を用いた場合の st と L1 ノルムに近い値をとる．図10,11 より，L1 ノルムは α = 0, π で誤差が+10% のとき，図 12, 図 13 より,L1 ノルムが低い値の r1 を選べばよい．
同様に α = π2 , 3π2 のときの L1 ノルムのシミュレーショ
ン結果を図 14,15 に，整定時間 st のシミュレーション結果を図 16,17 に表示する．

とおく．出力フィードバックと目標軌道の誤差をT

6r1

図 12 st(α = 0)

5 cos α−10 sin α −4.95 × 10−3 sin α x̂−10% =  −9.90 × 10−3 cos α cos β−3−9.90 × 10 sin β

5x 102.5st

元のシステムの速度成分を± 10% 変化させたものとして，それぞれ
5 cos α
−10 sin α −6.05 × 10−3 sin α x̂+10% = 
 −1.21 × 10−2 cos α cos β110.50.523456789234r1

図 16 st(α = π2 )56789r1

図 17 st(α = 3π2 )α = π2 や α = 3π2 のときのグラフもほぼ同じ形をしている．図 14,15 より，α = π2 , 3π2 で誤差が-10% のとき，
L1 ノルムの最小値は最小値 : L1 = 3.66[m/s]となる．r1 = 6.75 のとき L1 = 4.00[m/s] より大きくなる．上記のことと図 16,17 より 2.5 ≦ r1 < 6.75 の範囲の値を選べばよい.
次に，初期値の速度成分を 1 割減にした場合において，
α = 0 及び α = π のときの L1 ノルムのシミュレーション結果を図 18,19 に表示し，整定時間 st のシミュレーション結果を図 20,21 に表示する．

-3x 1076.5665.55.5
L1 ノルムL1 ノルムα = π2 と α = 3π2 のグラフも形はほぼ同じになった．図22,23 より，α1 = π2 , 3π2 で誤差が+10% のとき，L1 ノル
ムの最小値, 最大値は-3x 1076.554.554.5443.53.53

最小値 : L1 = 3.45[m/s]32.534567892345r16789r1

図 18 L1 ノルム (α = 0)5となり，L1 ノルムの誤差の範囲内なので図 24,25 より,，
図 19 L1 ノルム (α = π) 2.5 ≦ r1 ≦ 7.75 の範囲の値を選べばよい．角度 α による L1 ノルムの値から，目標までの楕円軌道の運動をする人工衛星の移動距離は，α = 0, π のときが最短になり，α = π2 , 3π2 のときが最長になることが分かった．
のどの点からも収束する
H の設計をした
α = 0, π2 , π, 3π2いので r1 = 2.5 のときにどの点でも L1 = 4.07[m/s] を下回り，状態フィードバック制御を用いたときの st に近い値をとる r1 の値は
r1 = 2.55x 10x 102.5

2.522st1.5st

1.5110.50.50023456789r123456789r1

図 20 st(α = 0)となる．r1 = 2.5 のときの L1 ノルムと整定時間 st は以下の表のようになる．

図 21 st(α = π)α = 0 のときと α = π のときのグラフはほとんど一致していることが見られた．図 18,19 より，α = 0, π で誤差が-10% のとき，L1 ノルムが低い値の r1 を選べばよい．図 20,21 より，r1 が 2.12 ≦ r1 ≦ 7.37 の範囲なら状態フィードバック制御を用いた場合の st と L1 ノルムに近い値をとる．
同様に α = π2 , 3π2 の場合における L1 ノルムのシミュ
レーション結果を図 22,23 に，整定時間 st のシミュレー
ション結果を図 24,25 に表示する．
-3-3x 10443.8

3.8

3.6

3.6

3.4

3.4
L1 ノルムL1 ノルムx 10

3.2

3.2

33

2.8

2.8

2.6

2.6

2.4

2.4234567892345r16789r1

図 22 L1 ノルム (α = π2 )図 23 L1 ノルム (α = 3π2 )55x 10x 1022

1.8

1.8

1.6

1.6

1.4

1.4stst

最大値 : L1 = 3.79[m/s]

2.52

1.2

1.2

11

0.8

0.8

0.6

0.6

0.4

0.4
2

345678r1

図 24 st(α = π2 )92345678r1

図 25 st(α = 3π2 )9

表 1 ±10% の比較

α0, π0, ππ 3π2, 2
π 3π2, 2

誤差-10%+10%-10%+10%

r1
2.5
2.5
2.5
2.5

L1
3.99[m/s]
3.22[m/s]
3.69[m/s]
3.71[m/s]

st5.74 × 104 [s]5.42 × 104 [s]5.74 × 104 [s]5.46 × 104 [s]以上の表より α = 0, π2 , π, 3π
2 のどの点からも収束する
r1 は r1 = 2.5 となる. そのときのオブザーバーゲイン Hは
0.64e − 1−1.24e − 4 −9.28e − 19
 −1.24e − 4
5.80e − 3
1.17e − 20  2.03e − 5
6.09e − 6
2.55e − 22 H= −7.60e − 6
1.67e − 5
1.72e − 21 
 −9.28e − 19 1.73e − 20
5.80e − 3 −5.52e − 22−9.80e − 22

1.66e − 5となった．オブザーバーゲインの各値も非常に小さい値が多く，出力フィードバックの操作量が少ないことが分かった．ただし，L1 ノルムは 1 割増の方が小さくなったので操作量自体は 1 割増の方が小さくなる．
目標までの楕円軌道運動による移動距離は，L1 ノルムπから α = 0, π のときが最長で，α = 3π2 , 2 のときが最短3πになる．よって，誤差を含めた α = 0, 2 , π, π2 以外のすべての角度 α でも上記のオブザーバーゲインを得られる．
オブザーバーゲインを求める際に初期速度成分を 1 割減しても 1 割増しても L1 ノルムは 4.00[m/s] 以内に収まり，かつ収束時間も 10 周期前後になったため，どの角度 α でも有効なオブザーバーを用いることができると分かった．

出力フィードバックによる軌道制御

こちらも同様に状態フィードバック制御とほぼ同じ軌
道を構成することができた．各軌道を見比べると，角度 α
前章で得られたフィードバックゲイン K とオブザーバーで初期位置は異なるものの，目標軌道までの人工衛星のゲイン H を用いてシミュレーションし軌道制御を行う．
動きも x, y 軸の平面において角度 α だけ異なるものだとπ1 割増にした場合の角度 α = 0, 3π2 , π, 2 のときの制御軌分かる．また，初期値の速度成分を誤差の範囲が ±10%
道をそれぞれ図 26,27.28,26 に表示する．
以内なら初期値を変えても r と r1 の値は変わらず, 収束時間が変化するだけで軌道そのものに変化はないと見られた．
誤差がシステムの初期値の速度成分成分に ±10% なら初期値がどこからでも制御できる K と H が決まり, 状態フィードバックと同じようにオブザーバーを用いても制御できることが分かった.

53D motion

3D motionsystemfeedback
final
observer10.5z0.5z

systemfeedback
final
observer100-0.5-0.5

-1
4-1
62

4-10-1502-10-2

00-5-4-2

0xx-4

5-6

10yy

6

終わりに誤差が ±10% 以内ならどこからでも L1 ノルムが L1 =4.000[m/s] 以下で収束するオブザーバーゲインを求め, 状態フィードバック制御とオブザーバーを用いた出力フィーバック制御の比較を行った. システムの初期値 x0 の α を
α = 0, π2 , π, 3π2 の 4 点を与え, オブザーバーを用いたシス
テムの初期値を速度成分に誤差を ±10% に設定し,L1 ノルムが L1 = 4.00[m/s] 以下で収束するオブザーバゲイン
H を求めた. 初期値や誤差に関わらず収束させるフィー
ドバックゲイン K の r の値は r = 6.75 になった. 初期値や誤差に関わらず収束させるオブザーバーゲイン H の r1の値は r1 = 2.5 になった.
制御軌道図より, 出力フィードバック制御を用いても状
図 29 制御軌道 (α = 3π
図 28 制御軌道 (α = π)2 ) 態フィードバック制御と同じような軌道を描ているので4 点の初期値の全てから同じ r1 を用いて制御でき,α の値
図 26,27,28,29 から，状態フィードバック制御と出力フに関わらず誤差が速度成分の ±10% ならオブザーバーをィードバック制御の軌道がほぼ同じであることが分かる．
用いた出力フィードバックでも状態フィードバックと同L1 ノルムや整定時間の条件を満たしたオブザーバーで軌じように制御できる.
道修正も状態フィードバックに近似できた．
次に，速度成分を 1 割減にしたものをそれぞれ図 参考文献30,31,32,33 に表示する．
[1] M. Bando and A. Ichikawa:Formation Flying Along
制御軌道 (α = π2 )図 27

図 26 制御軌道 (α = 0)

3D motion

3D motionsystemfeedback
final
observersystemfeedback
final
observer1

10.5zz0.5

00-0.5-0.5-1
10-65-1-5-4-20000521015-104-5xxyy

3D motion

3D motionsystemfeedback
final
observer

1

10.5z

0.5z

systemfeedback
final
observer0[2] A. Ichikawa:Null Controllability with Vanishing Energy for Discrete-Time Systems, Systems & Control
Letters, Vol. 57, No. 1, pp. 34-38, 2008.

0-0.5-0.5-16-1424-10-1002-5-2000-4-2xx-45-610yy

図 31 制御軌道 (α = π2 )図 30 制御軌道 (α = 0)

3D motion

3D motionsystemfeedback
final
observersystemfeedback
final
observer1

1[3] A. Ichikawa:Recent Developments in Formation Flying, Lecture Notes ver.1, 2010.
[4] R. Jifuku, A. Ichikawa and M. Bando:Satellite Formation by Pulse Control Along a Circular Orbit,
Journal of Guidance, Control, and Dynamics, Vol.
34, No. 5, pp. 1329-1341, 2011.

0.5z0.5za Circular Orbit with Control Constraints, Proceedings of AAS/AIAA Astrodynamics Specialist Conference, Alaska, AAS 11-491, 2011.

00-0.5-0.5-1
1065-1-54-502000-55-210-105x-4xyy

図 32 制御軌道 (α = π)

図 33 制御軌道 (α = 3π2 )[5] M. Shibata and A. Ichikawa:Orbital Rendezvous andFlyaround Based on Null controllability with Vanishing Energy, Journal of Guidance, Control, andDynamics, Vol. 30, No. 4, pp. 934-945, 2007.