出力レギュレーションによる複数の人工衛星のフォーメーション形成
M2014SC005 池山拓弥
指導教員：大石泰章

1

はじめにクリッドノルムを R0 ，R とすると Newton の運動方程式
より，µẍ = 2nẏ + n2 (R0 + x) − 3 (R0 + x) + ux ,Rµÿ = −2nẋ + n2 y − 3 y + uy ,(1)Rµz̈ = − 3 z + uz
R

現在，数多くの宇宙機が打ち上げられそれらによる様々なミッションが実行されている．その中でも，複数の宇宙機の編隊飛行によりミッションを達成する方法をフォーメーションフライトという．フォーメーションフライトでは主衛星に対して従衛星が適切な相対位置を維持することが求められる．また，衛星運用を考えるとき地球周回T上の主衛星とその付近を飛行する従衛星の相対運動が重 が得られる．ここで，u= [ux uy uz ] は宇宙機の制御加要である．フォーメーションフライトに関する研究の例 速度である．式 (1) を原点 x = y = z = 0 で線形化すとして文献 [2] では，地球の周りの円軌道上の主衛星近傍 ると,ẍ − 2nẏ − 3n2 x= ux ,にある 1 機の従衛星に対して出力レギュレーション理論を適用しフォーメーション形成を行っている．また文献
ÿ + 2nẋ= uy ,(2)2[3] では，出力レギュレーション理論を適用し楕円軌道上z̈ + n z= uzの主衛星近傍にある複数の従衛星のフォーメーション形成を行っている．本研究では，出力レギュレーション理 が得られる．式 (2) は HCW 方程式と呼ばれる [1][2]．ま論を用い，円軌道上の主衛星近傍にある複数の従衛星の た，式 (2) より z の運動は x，y と独立しているのがわ
フォーメーション形成をインパルス制御で行う．そして， かる．
本研究では 式 (2) を軌道半径 R0 および円軌道の角速制御器の性能評価を行うことを目的とする．なお，制御器の性能評価には総速度変化 (∆V ) とフォーメーション 度 n を用いて R0 が 1 になり，√ 1/n が 1 になるように精度を用いる．また，前述の文献に対し，本研究では特にフォーメーション精度に着目して制御器の設計 (追加入力の導入) を行う．

2

円軌道上の相対運動無次元化する．ここで n =µ/R0 3 であり，µ は地球の重力定数である. また時間 t を τ = t/(1/n) により置
き換え，(x, y, z) = (x/R0 , y/R0 , z/R0 ), (ux , uy , uz ) =(1/R0 n2 )(ux , uy , uz ) とすると, 式 (2) を無次元化した状態方程式は地球を中心とした半径 R0 ，軌道周期 T = 2π/n (n はx′ = Ax + Bu,x(0) = x0(3)円軌道の角速度) の円軌道上の主衛星とその近傍を運動する従衛星を考える，主衛星に対する従衛星の相対運動を となる．ここで x = [x y ẋ ẏ z ż]T ,u = [u u u ]T であx y z
考えるため，主衛星の重心を原点とする図 1 の回転座標 り，′ は τ による微分を表している．なお，
系 o - {i,j,k} を考える．ここで，o は主衛星の質量中心，
0 0 1 0 0 00 0 0i は動径方向の単位ベクトル，j は主衛星の飛行方向の単 0 0 0 1 0 0 
 0 0 0 位ベクトルであり，k はそれらに直交する単位ベクトル
 3 0 0 2 0 0  1 0 0 である．
A=,B = 
 0 0 −2 0 0 0 
 0 1 0 
 0 0 0 0 0 1  0 0 0 0 000−1 00 01である．また，等間隔の離散的な時刻 τ1 , τ2 · · · にインパルス入力を行ったときの (5) の離散時間の動特性は

xj+1 = Ad xj + Buj(4)となる．ここで，x(τj ) = xj であり u(τj ) = uj である．
また，Ad は expA(τj+1 − τj ) であり，τj+1 −τj が j によらず一定なので Ad も定数行列になる．
図 1 円軌道上の主衛星と従衛星
このとき主衛星から見た従衛星の位置ベクトルを r =xi + yj + zk とおく. 従衛星の地球の質量中心からの位置ベクトルは R = R0 + r であり，ベクトル R0 ，R のユー

3

出力レギュレーション理論を用いたフィードバック設計

従衛星のフォーメーション形成を行うために，出力レギュレーション理論を用いる [2]. 出力レギュレーション

とは制御対象のシステムの出力を外生信号に追従させる制御系の設計手法である. ここで外部システムを′

w = Sw, w(0) = w0(5)とする．w は外生信号であり，w= [w1 w2 w3 w4 ]T とす
る．式 (5) の離散化は

wj+1 = Sd wj , w(0) = w0(6)となる．ここで，Sd は expS(τj+1 − τj ) である．また, 評
価用出力をz j = Cxj + Dwj

1 0
D = − 0 10 0001

0
0 0とする．評価用出力 z j = [z1j z2j z3j ]T とすると式 (7) は

z1j = xj − w1 ,z2j = y j − w2 ,(8)z3j = z j − w3と表せ，z j が零に近づくとき，xj と wj の差が小さくなる．これは制御対象が外生信号に追従していることになる [2]．式 (4)，式 (7) に対する出力レギュレーション問題の可解条件は以下のレギュレーション方程式

Ad P1 − P1 Sd + BP2 = 0,
CP1 + D = 0(9)が解 P1 ,P2 を持つことである. この解を用いたフィードバック制御器

uj = F xj + (P2 − F P1 ) wj

(10)

により，出力レギュレーションが達成される. すなわち,評価用出力 z j が 0 に収束し，式 (8) より制御対象の x 座標，y 座標，z 座標と外生信号の第 1 成分，第 2 成分，第3 成分の差が 0 になることを表す．また，フィードバック
ゲイン F は離散時間リッカチ代数方程式 (DARE)X = Q + Ad T XAd − Ad T XB(R + B T XB)−1 B T XAdの解 X によりF = (R + B T XB)−1 B T XAd
で求める [2]．ここで Q，R は重み行列である．

シミュレーション

4.1

フォーメーション形成

本研究では，フォーメーションの例として従衛星 3 機が正三角形のフォーメーションを形成する場合について考える．すなわち，各従衛星が外部システムにより生成された周期軌道上を運動し，時刻 2π/3 間隔ごとに各従衛星が周期軌道上の 3 個の参照点を相異なる位相で通り，一定の距離を保つことでフォーメーションを形成することを考える．なお，インパルス入力も時刻 2π/3 間隔ごとに行うこととする．まず，従衛星が運動する周期軌道を以下のように設定する：(7)と定義し，これを零に収束させることを考える．ここで本研究では C ，D を

1 0 0 0 0 0
C =  0 1 0 0 0 0 ,0 0 0 0 1 0



4xr (τ ) = 0.0005 cos ωτ,yr (τ ) = −0.001 sin ωτ,√zr (τ ) = 0.0005 3 cos ωτ.

(11)

式 (11) の周期軌道は式 (5) の外部システムにより生成される. このとき S は


0
ω/200
 −2ω000 
S= 0001 00−ω 2 0とする．本研究では ω = 1 の場合についてシミュレーションを行った. これより S は
0 1/2 0 0
 −200 0 S= 000 1 00−1 0(1)(2)(3)となる．参照点 xref ，xref ，xref を以下に示し，各時刻に各衛星がそれぞれどの参照点を通るかを表 1 に示す：

0.0005 cos 2π/3
(1)xref =  −0.001 sin 2π/3  ,√0.0005 3 cos 2π/30.0005 cos 4π/3(2)
xref =  −0.001 sin 4π/3  ,√0.0005 3 cos 4π/3


0.0005 cos 2π(3)xref =  −0.001 sin 2π  ．
√
0.0005 3 cos 2π

表 1 各時刻における衛星の位置
従衛星 2π/3 4π/3
2π8π/3 ・・・衛星 1
衛星 2
衛星 3(1)xref
(2)
xref
(3)
xref

(2)

xref
(3)
xref
(1)
xref

(3)

xref
(1)
xref
(2)
xref

(1)

xref
(2)
xref
(3)
xref

・・・・・・・・・

外部システムにより生成された周期軌道を図 2 に示す．

ときの衛星 1 の ∆V と誤差のグラフを図 3 に，r = −2 における衛星 3 機の軌道の x 座標と外生信号 w の第 1 成分
w1 の時間変化のグラフを図 4 に示す．また，r = −2 における衛星 3 機の 3 次元軌道を図 5 に示す．

×10-3
1z0.5
0-3-4×10

2×10

4-0.5-1
150×10-3

y0-1-5x×10-41.52

図 2 外部システムの周期軌道と参照点
また衛星 1，衛星 2，衛星 3 に対する外部システムの初期値を w01 ,w02 ,w03 とし
0.00050√ 
w01 = , 0.0005 3 00.0005 cos 2π/3
 −0.001 sin 2π/3 √
w02 = ,
 0.0005 3 cos 2π/3 √−0.0005 3 sin 2π/3
0.0005 cos 4π/3
 −0.001 sin 4π/3 √
w03 =  0.0005 3 cos 4π/3 √−0.0005 3 sin 4π/3とする．以上のように設定することで 3 機の従衛星のフォーメーション形成を行う [3]．

4.2

シミュレーション条件

本研究で行うシミュレーションの条件を以下に示す．

1-2• r をパラメータとして変化させ，5 周期分のシミュレーションを行う．
• 制御器の性能評価には総速度変化 (∆V ) とフォーメー
ション精度 (以下では誤差と表記する) を用いる．
• 誤差は，5 周期後の各衛星とそれぞれに対応する参照点との距離とする．
また，∆V の計算式を以下に示す．
ni∑∆V =(|uxk | + |uyk | + |uzk |)k=1

ここで，ni はインパルスの入力回数である．なお，式 (3)の線形システムで制御器を設計する．その制御器を無次元化した非線形システム (1) に適用する．

4.3

シミュレーション結果 (追加入力なし)

r をパラメータとして変化させながらシミュレーション行い，制御器の性能について検証を行う．r を変化させた-1-0.500.5

11.520r

図 3 r を変化させたときの衛星 1 の ∆V と誤差
1×10-30-1
10×10-351015202530350-3×105101520253035051015202530350-110-1

図 4 r = −2 における衛星 3 機の x と w1 (x：青線 w1 ：赤線)
×10-31• 各従衛星の初期値は原点 (主衛星近傍) とする．

0.5z

• 重み行列 R は 10r I3 とし， 重み行列 Q は I6 と固定する．ただし，In は n 次元単位行列である．

-1.50-0.5-121×10-30y0-2-1x×10-3

図 5 r = −2 における衛星 3 機の 3 次元軌道
図 3 より r の値が小さくなることで入力にかかる重みが小さくなり，その結果 ∆V が増加し，誤差は減少していることがわかる．衛星 2，衛星 3 についてもグラフの概形はほぼ同様であった．図 4 より各従衛星の x 座標が外生信号 w1 に追従していることがわかる．y 座標，z 座標についてもグラフの概形はほぼ同様である．なお，入力の重みを大きくして r = 2 とした場合は r = −2 の場合と比べて外生信号に追従するまでに時間がかかっていた．なお，図 4 の ∗ はそれぞれの時刻における参照点を表している．また，図 5 より各従衛星が相異なる位相で参照点を通り，三角形のフォーメーションを形成しているのがわかる．次に r = −2 の場合の各従衛星の ∆V と誤差を表 2 に示し，r = 2 の場合の各従衛星の ∆V と誤差を表 3 に示す．
表 2 r = −2 の場合の各従衛星の ∆V と誤差
従衛星
∆V
誤差
衛星 1 1.84×10−3 1.30×10−5
衛星 2 2.43×10−3 1.43×10−5
衛星 3 2.62×10−3 5.49×10−6

表 4 r = −2 の場合の各従衛星の ∆V と誤差 (追加入力あり)
従衛星
∆V
誤差
衛星 1 2.54×10−3 1.17×10−5
衛星 2 4.03×10−3 1.46×10−5
衛星 3 4.25×10−3 3.85×10−6
表 5 r = 2 の場合の各従衛星の ∆V と誤差 (追加入力あり)
従衛星
∆V
誤差−3
衛星 1 1.15×10
2.38×10−4−3
衛星 2 1.14×10
1.21×10−4
衛星 3 1.15×10−3 1.19×10−4

表 3 r = 2 の場合の各従衛星の ∆V と誤差
従衛星
∆V
誤差−3
衛星 1 1.08×10
3.00×10−4−3
衛星 2 1.06×10
1.50×10−4
衛星 3 1.07×10−3 1.51×10−4
表 2，表 3 より r = −2 の場合では r = 2 より ∆V は大きいが誤差は小さいことがわかる．また，r を r = −2 より小さくしても ∆V ，誤差ともに変化はほとんど見られなかった。
4.4 シミュレーション結果 (追加入力の導入)

表 2，表 3 と表 4，表 5 を比較すると追加入力を行った
場合はほとんどの場合で誤差が減少していることが確認できる．また，追加入力を行っているので ∆V は増加傾向にある．次に，各従衛星の誤差についての増減率を表6 に示す．

フォーメーション形成において追加入力を導入したシミュレーションを行う．先ほどの結果より，一定の誤差が残ってしまうことがわかった．本研究では，誤差を抑える方法として，追加入力を行う手法を提案する．式 (10)に以下で定義する制御入力 (追加入力) を加える．

表 6 各従衛星の誤差についての増減率[a1 ux a2 uy a3 uz ]T(12)

従衛星
衛星 1
衛星 2
衛星 3

増減率 (r = −2 の場合)
約 10% 減約 2% 増約 30% 減増減率 (r = 2 の場合)
約 21% 減約 19% 減約 21% 減ここで，ux ,uy ,uz は式 (10) の uj の第 1 成分，第 2 成分， 表 6 より追加入力を行うことでほとんどの場合で誤差
第 3 成分であり a1 ，a2 ，a3 は設計パラメータである．本 が 10%〜30% 程度減少していることがわかる．しかし，研究では a1 ，a2 ，a3 を以下のように定める．
衛星 2 の r = −2 の場合では誤差が 2% 程度増加してい• 手順 1：誤差の基準値を定める．(本研究では1.0 × 10−6 とする．)

ることがわかる．これは，追加入力が余分な入力であったと考えられる．以上の結果より追加入力を行うことで• 手順 2：インパルス入力時刻において x 軸，y 軸， 誤差を減少させることができたが，追加入力を行わないz 軸ごとに誤差 (各軸とそれに対応する外生 方が良い場合もあることがわかった．
信号とのずれ) を計算する．

• 手順 3：x 軸の誤差が 1.0 × 10−6 以下の場合は a1 = 0とする．(y 軸の場合は a2 = 0 とし，z 軸の場合は a3 = 0 とする．)• 手順 4：x 軸の誤差が 1.0×10−6 より大きく 1.0×10−4より小さい場合は a1 = 0.1 とする．(y 軸の場合 a2 = 0.1，z 軸の場合 a3 = 0.1 とする．)• 手順 5：x 軸の誤差が 1.0 × 10−4 以上で 1.0 × 10−3より小さい場合は a1 = 0.2 とする．(y 軸の場合 a2 = 0.2，z 軸の場合 a3 = 0.2 とする．)• 手順 6： x 軸の誤差が 1.0 × 10−3 以上場合は a1 = 0.3とする．(y 軸の場合 a2 = 0.3，z 軸の場合a3 = 0.3 とする．)
このように各軸の誤差に応じて設計パラメータを定めることでフォーメーション精度の向上を図る．以上の手順を各従衛星について行い追加入力を決め，シミュレーションを行った．次に r = −2 の場合の各従衛星の ∆V と誤差を表 4 に示し，r = 2 の場合の各従衛星の ∆V と誤差を表 5 に示す．

5おわりに

本研究では出力レギュレーションを用いて円軌道上の主衛星近傍にある 3 機の従衛星のフォーメーション形成を行った．その結果，フォーメーションを形成する際に誤差が生じてしまうことが分かった．その対策として本研究では追加入力を行う手法を提案した．その結果，誤差が微増してしまう場合もあったが，ほとんどの場合で誤差を減少させていることが確認できた．

参考文献[1] R. Jifuku, A. Ichikawa, M. Bando: Optimal pulse
strategies for relative orbit transfer along a circularorbit, Journal of Guidance, Control, and Dynamics,
Vol. 34, No. 5, pp. 1329–1341, 2011.
[2] M. Bando, A. Ichikawa: Active formation along acircular orbit by pulse control, Proceeding of theSICE Annual Conference 2013, pp. 2197–2203.
[3] M. Bando, A. Ichikawa: Precise formation at discrete points along elliptic orbit, preprint.