2 つの代数リッカチ方程式とフォーメーションフライト
M2012MM051 吉田賢司
指導教員：市川1朗が得られる．ここで u = [ux uy uz ]T は従衛星に働く制御加速度，R = [(R0 + x)2 + y 2 + z 2 ]1/2 である [2, 3]．
この方程式 (1) を原点 x = y = z = 0 で線形化するとはじめに円軌道上のフォーメーションフライトは多くの研究者によって考察されている．フィードバック制御を設計するとき，通常は線形 2 次レギュレータ理論のリッカチ方ẍ − 2nẏ − 3n2 x = ux
程式をもとに設計する．一方で，指数関数の重みを課す，
ÿ + 2nẋ = uy(2)異なる種類のレギュレータ問題が知られており，最適解z̈ + n2 z = uz
は特異リッカチ方程式によって与えられる [1]．しかし，2つのリッカチ方程式によるフィードバック制御の比較はが 得 ら れ る ．こ の 方 程 式 は Hill-Clohessy-Wiltshire
検討されていない．そこで，本研究では円軌道上のフォー(HCW) 方程式とよばれる．推力を u = 0，初期値をメーション再構成問題について，2 つの代数リッカチ方程[x0 y0 ẋ0 ẏ0 z0 ż0 ]T として解き (2) 式をパラメータ表現
式によるフィードバック制御の設計を比較する．フィーするとドバック制御の燃費は入力の絶対積分 (L1 ノルム) で評価する．線形 2 次レギュレータ理論のリッカチ方程式では
x(t) = 2c + a cos(nt + α)
制御入力の重みをパラメータとして変化させ，入力に指y(t) = d − 3nct − 2a sin(nt + α)(3)数関数の重みを課すレギュレータ理論の特異リッカチ方z(t) = b cos(nt + β)
程式では指数関数のパラメータを変化させる．この 2 つの設計法について，整定時間 (目標軌道到達時間) と L1となる. ここでノルムのグラフによる比較を行う．このグラフをもとに[]1/2
フィードバックを設計して制御入力の大きさやのフィーa = (3x0 + 2ẏ0 /n)2 + (ẋ0 /n)2, c = 2x0 + ẏ0 /n
ドバックゲインの比較も行う．
d = y0 − 2ẋ0 /n, sin α = −ẋ0 /na

2

フォーメーション問題の方程式

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

LeaderR0

Center of Earth

Circular Orbitcos α = − (3x0 + 2ẏ0 /n) /a, b = [z02 + (ż0 /n)2 ]1/2cos β = z0 /b, sin β = −ż0 /nbである. 面外運動 (z, ż) は周期解となっている．面内運動(x, y) は (3) 式において c = 0 のとき周期解となる (CW条件)．
x = [x y ẋ ẏ z ż]T とおくと (2) 式の状態方程式はẋ = Ax + Bu, x(0) = x0と表せる．ここで
001
 000
 3n2 00
A= 00 −2n 00000000102n000000 −n2(4)00
 00  10 
, B = 
 00 01 00 
0 0 
0 

01000100 0

図 1 円軌道上の主衛星である．
このとき相対位置ベクトルを r = xi + yj + zk として， パラメータ (面内運動; 面外運動) = (a, c, d, α; b, β) に運動方程式を変形しそれぞれi,j,kについて係数を比較す より表される (4) 式の解を γ H = (a, c, d, α; b, β) と表す.
c = 0 のときこの解は周期軌道となり γ H = (a, d; b) と表るとす．フィードバック制御によるフォーメーション形成問題µẍ = 2nẏ + n2 (R0 + x) − 3 (R0 + x) + uxとは，(4) 式の解を与えられた周期軌道 γfH = (af , df ; bf )Rµに漸近的に追従させることである．このときの評価関数(1)ÿ = −2nẋ + n2 y − 3 y + uy
Rは，制御に使う推力の絶対積分 (L1 ノルム) であり，これµz̈ = − 3 z + uz
は消費燃料に比例する．特に従衛星の初期軌道が周期軌R道であるときは，フォーメーション再構成問題となる [2]．

と変形する．X が逆行列をもつとすると Y = X −1 はフィードバックの設計

3

HCW システム (4) の周期解を目標軌道 xf とし，初期軌道から，この軌道にフィードバック制御で移行させることを考える．目標軌道の方程式を
ẋf = Axf , xf (0) = xf 0とおくと，この軌道との誤差 e = x − xf はė = Ae + Bu, e(0) = e0

u = −Ke(6)で与えられる．ここで K は，A − BK が安定となる任意のフィードバックゲインである．K の設計には，2 次形式評価関数を最小にする次の 2 つの方法を用いる．
線形 2 次レギュレータ理論のリッカチ方程式

評価関数

∫ ∞
J(u; x0 ) =(x′ (t)Qx(t) + u′ (t)Ru(t))dt(7)0を最小化するフィードバックゲイン K は

K = R−1 B ′ X(8)A′ X + XA + Q − XBR−1 B ′ X = 0(9)を満たす解である．ここで Q は半正定，R は正定行列で√ある．( Q, A) が可検出であるとき，リッカチ方程式は安定化解 (A − BR−1 B ′ X が安定となる解）をもち，可観測であるときその解は正定となる．HCW システムはエネルギー零収束原点可制御 (NCVE)[4] であるので，Q を0 に近づけるとき（または R を大きくしていくとき）Xは 0 に収束し，フィードバック制御の 2 乗積分 (L2 ノルム) は 0 に収束する．
入力に指数関数の重みを課した特異リッカチ方程式

安定化の条件の下で，評価関数∫ ∞Jγ (u; x0 ) =exp(γt)u′ (t)Ru(t)dt(A − BR−1 B ′ X)′ X + X(A − BR−1 B ′ X)
= −γX − XBR−1 B ′ X < 0であるから A−BR−1 B ′ X は安定となる．このときフィードバック Kγ = BR−1 B ′ X は減衰率 γ/2 をもつことを示す．リッカチ方程式

A′ X + XA + Q − XBR−1 B ′ X + γX = 0
は，A + γ/2I − BR−1 B ′ X が安定となる解をもつ．このとき A − BR−1 B ′ X は減衰率 γ/2 となる．Q を 0 に近づけると X はリッカチ方程式 (12) の解に収束する．従って
Kγ = BR−1 B ′ X は減衰率 γ/2 をもつ [1]．HCW システ
ムは NCVE であるので γ を 0 に近づけるとき，リッカチ方程式 (12) の解は 0 に収束し，フィードバック制御の 2乗積分 (L2 ノルム) も 0 に収束する．

シミュレーション

4

ここでは円軌道の半径および角速度を用いて状態方程式 (4) を無次元化したシステム ẋ = Ax + Bu,で与えられる．ただし X はリッカチ方程式

3.2を満たす．A の固有値は虚軸上にあるので (−A−γ/2 I) は安定となる．従って，正定解 Y が存在し，(−A−γ/2 I, B)の可制御グラミアンという．以上よりリッカチ方程式 (12)は正定解をもつ．さらに(5)となり，フィードバック制御は3.1−Y (A + γ/2 I)′ − (A + γ/2 I)Y + BR−1 B ′ = 00 0 3A=
 0

 00

A X + XA − XBR′B X + γX = 0

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


00 

0 

0 

0 1(13)

xf = [0.005 0 0 − 0.01 0 0]T(10)(11)で与えられる．ただし Xγ はリッカチ方程式−1012000x0 = [0.01 0 0 − 0.02 0.01 0]T ,を最小化するフィードバックゲイン Kγ は′100−200

により，2 つの設計法について L1 ノルム，L2 ノルムと整定時間の違いを，リッカチ方程式の重み R，γ を変えながらシミュレーションで考察する．初期値と目標値はそれぞれ0Kγ = R−1 B ′ Xγ000000とする．整定条件である目標軌道までの軌道誤差 e は1.00 × 10−5 以内とする．このときのパラメータは γ0H =(a0 , d0 ; b0 ) = (0.01, 0, 0.01)，γfH = (af , df ; bf ) =(0.005, 0, 0) である．

4.1(12)を満たす解である．この方程式を解くには (12) 式を(A + γ/2 I)′ X + X(A + γ/2 I) − XBR−1 B ′ X = 0(14)線形 2 次レギュレータ理論のリッカチ方程式まず，線形 2 次レギュレータ理論のリッカチ方程式 (9) のフィードバックについて考える．フィードバックは (6),(8)で与えられる．このときの評価関数は (7) の x を e に置き換えたものである．このシステムは NCVE であるので入力 u を小さくして L1 ノルムを小さくするために，状態の0.25
in-planeout-of-plane

in-planeout-of-plane

0.08

0.2
0.07
0.06

0.15
0.05

L2-norm

L1-norm

重み Q は小さな値に固定する．リッカチ方程式 (9) のパラメータをそれぞれ Q = diag(qi ), qi = 1.00 × 10−3 (i =1, 2, 5), qi = 0.00(i = 3, 4, 6)，R = 10r I3×3 として r の値を大きくして制御入力の重みを変化させる．

0.04

0.1

0.03
0.02

0.05
0.0350.035
in-planeout-of-plane

0.03

0.01

in-planeout-of-plane

0.03

0-1

0.02
L2-norm

0.02
L1-norm

010.025

0.025

0.015

0-1201ρρ(a) ρ − L1 ノルム(b) ρ − L2 ノルム

20.015

0.01

0.01

0.005

0.005

図 4 ρ − L1 , L2 ノルム

00-4-20-42-202rr

700(a) r − L1 ノルム(b) r − L2 ノルム

SARE controller
600500settling time

図 2 r − L1 , L2 ノルム
図 2(b) から，r の値が大きくなると L2 ノルムの極限が 0 に向かうことが確認できる (NCVE である)．また図2(a) から，r の値が大きくなると入力 u が小さく抑えられるため L1 ノルムが小さくなることがわかる．L1 ノルムは一定値 (L1 = 0.00308, L1 z = 0.0127) まで小さくなる．

400
300
200
1000-1012ρ

図 5 ρ − 整定時間4500
LQR controller
4000

3000

settling time

フィードバック設計法の比較

4.3

3500

2 つのリッカチ方程式による設計法の違いを比較したものが図 6 である．

2500
2000
1500-3x 101000

LQR controller
SARE controller

14500

LQR controller
SARE controller

0.028
0.026

120-4-2020.024r

図 3 r − 整定時間

L1z-norm

L1-norm

1080.022
0.020.018

60.0160.014
40.012

図 3 から，r の値が大きくなると入力 u が小さく抑えられるため整定時間が急激に大きくなることがわかる．

4.20100200settling time300(a) 面内運動050settling time100(b) 面外運動

入力に指数関数の重みを課した特異リッカチ方程式

図 6 整定時間 − L1 ノルム
次に，入力に指数関数の重みを課したときのレギュレータ理論の特異リッカチ方程式 (12) のフィードバックについ図 6 から，リッカチ方程式 (LQR) の重み Q の速度成分て考える．フィードバックは (6)，(11) となる．このときの を 0 としたときは，整定時間が小さいとき，特異リッカチ評価関数は (10) である．ここでもシステムは NCVE であ 方程式 (SARE) を用いたほうが短時間・低燃費なフィーるので入力 u を小さくして L1 ノルムを小さくする．特異 ドバックを設計できることがわかる．
リッカチ方程式 (12) のパラメータをそれぞれ R = I3×3 ， ここで，L1 ノルムに差があるときについて，整定時間
γ = 10−ρ として ρ の値を大きくして γ を変化させる．
を 4π(2 周期) と指定したときのフィードバックを設計して 2 つの代数リッカチ方程式による設計の違いを調べる．
リッカチ方程式 (9) のときと同様に，図 4(b) から，ρ の値が大きくなると L2 ノルムの極限が 0 に向かうことが確認できる (NCVE である)．また図 4(a) から，ρ の値が大きくなると入力 u が急激に抑えられるため L1 ノルムが小さくなることがわかる．L1 ノルムは最適レギュレータと同じ値まで小さくなる．
図 5 から，ρ の値が大きくなると入力 u が小さく抑えられるため整定時間が大きくなることがわかる．

4.3.1

整定時間を指定する場合

整定時間が 4π であるときのパラメータを図 3 と図 5 から求めてフィードバックを設計する．そのときの L1 ノル
ムは図 2 と図 4 や図 6 から求める．
表 1 から，整定時間 (ST) が同じとき L1 ノルムは SAREのほうが小さくなる．SARE のほうが L1 ノルムが小さい理由は，図 7 から入力 u が小さく抑えられているからである．そのため入力の絶対積分である L1 ノルムも小さくなるのである．

異リッカチ方程式 (SARE) のグラフの間に表れている．

-3x 10

表 1 シミュレーション結果rρ
L1
L1 z-3.24

0.00976
0.003270.273-30.0137
0.0129LQR controller(q=10)-212LQR controller(q=1010LQR controller(q=10LQR controller(q=1))-14π4π

L1-norm

LQR
SARE

LQR controller(q=0)14ST)
1

LQR controller(q=10SARE controller)

86-3-3-3x 10x 10x 10844LQR controller
SARE controller

2LQR controller
SARE controller

3LQR controller
SARE controller

64042y-4-61002input u

input ux

input uz-2

0100
200
settling time300

-2-8-1-4-10-2-6

図 9 整定時間 − L1 ノルム (面内運動)-12-8-3-14
0510015510015510timetimetime(a) 入力 ux(b) 入力 uy(c) 入力 uz

15

図 7 時間 − 入力u
このとき設計したフィードバックゲインは次のようになった．

3.92 −1.27 2.14 0.691 0.00 0.00

K =  3.15 0.371 0.691
1.74 0.00 0.00  ，
0.00
0.00 0.00 0.00 0.659 1.15

リッカチ方程式 (LQR) の速度成分の重みを少しずつ大きくすると特異リッカチ方程式 (SARE) のグラフに傾きが近づいていくことがわかる．しかし，q = 1 と q = 101のグラフはほとんど重なっている．そのため，速度成分に重みを課しただけでは特異リッカチ方程式 (SARE) の燃費性能までは近づけられないことがわかった．

5

考察

特異リッカチ方程式 (SARE) のほうが燃費性能が良い理
由の 1 つに，重みの課し方の違いがある．特異リッカチ方
1.19 −0.0732 0.550 0.525 0.00 0.00
 程式 (SARE) はリッカチ方程式 (LQR) の重みを Q = γX
Kγ =  2.98 −0.0814 0.525 1.58 0.00 0.00  ．としたものである．物理的に意味のある対角成分以外に0.000.00 0.00 0.00 0.284 1.07
も重みを課すことによって燃費性能を向上させることがゲインの各要素の絶対値は Kγ のほうが小さくなるため， できると考えられる．
SARE のほうが入力 u が小さいことがわかる．また，こ6 おわりにのときの制御軌道は図 8 のようになった．
整定時間―L1 ノルムのグラフ (図 6) から，2 つの代数リッカチ方程式によるフィードバックの設計法を比較した．このグラフから，特異リッカチ方程式を用いたほうが整定時間と L1 ノルムが小さくなり，短時間・低燃費なフィードバックを設計できることがわかった．この理由は，整定時間を指定したときの入力のグラフ (図 7) を比較することで，入力が小さく抑えられていることによって確認できた．

controlled
initialfinal

controlled
initialfinal0.015

0.015

z0.010.005

z0.010.005

00-0.005-0.005-0.01-0.01-0.015-0.015

0.01

0.01

-0.02

-0.02

0000-0.01-0.01

xx0.020.02

y(a) LQR controllery(b) SARE controller

図 8 制御軌道
図 7 からリッカチ方程式 (LQR) のほうがはじめに大きな入力を与えるため，LQR(図 8(a)) のほうが早く目標軌道に近づくことがわかる．

4.4

速度に重みを付加するとき

次に，重み Q の速度成分に重みを課して整定時間が小さいときのリッカチ方程式 (LQR) の L1 ノルムの改善を図る．
図 6(b) から，面外運動の L1 ノルムの差はなかったため，リッカチ方程式 (LQR) の重み Q の面内運動の速度成分 (3,3)，(4,4) に重みを課して L1 ノルムの改善を図った．その結果が図 9 である．重みを 0 としたグラフと特参考文献[1] B. Zhou, G. Duan and Z. Lin:A Parametric Lyapunov Equation Approach to the Design of Low Gain
Feedback, IEEE Transaction on Automatic Control,Vol. 53, No. 6, pp. 1548-1554, 2008.
[2] A. Ichikawa:Recent Developments in Formation Flying, Lecture Notes ver.1, 2010.
[3] 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.
[4] A. Ichikawa:Null Controllability with Vanishing Energy for Discrete-Time Systems, Systems & ControlLetters, Vol. 57, No. 1, pp. 34-38, 2008.