パルス制御による円軌道上のフォーメーション
2011SE054 堀尾 一孝

2011SE133 小玉 純乃

指導教員：市川朗
る. このとき, 主衛星に対する従衛星の相対位置ベクトル

はじめに1を,r = xi + yj + zk とする.

フォーメーションフライト (編隊飛行) により, 単一の人

o：主衛星の質量中心i：半径方向j：飛行方向k：軌道面の外向き
nk：角速度工衛星だけでは, 困難なミッションや, 高精度の観測が可能となる. その例として, レーダーによる大気観測や高精度重力測定などがある. その基礎研究として, 地球周回軌道上の衛星とその近くを飛行する従衛星の相対運動が注目されている. このとき, 従衛星は適切な相対位置関係を維持することが望ましいと考えられる. 本研究では, フォー メーショ運動方程式より,ンフライトで一般的なパルス制御を用いる. パルス制御とは, 離散時間型の一定時間入力を与える制御方法である. まẍ = 2nẏ + n2 (R0 + x) −た, インパルス制御とは瞬間的に入力を与える離散時間型

ÿ = −2nẋ + n2 y −の制御方法である. 研究課題として, 燃費消費量の削減のたµz̈ = − 3 z + uz
Rめに, 総速度変化を抑える制御を行うことである. ここで,インパルス制御とパルス制御での燃費消費量を比較検討すµy + uyR3(1)Tが得られる. このとき,u = [ux uy uz ]は従衛星に働くẍ−2nẏ − 3n2 x = ux
ÿ+2nẋ = uy

(2)

[]1
加速度であり,R = (R0 + x)2 + y 2 + z 2 2 である. ここで,(1) を原点 x = y = z = 0 で線形化すると,る.

2µ(R0 + x) + ux
R3

円軌道上の相対運動方程式
ニュートンの運動方程式を衛星の質量で割ると,R̈ = −µR+u
R3

R：地球から衛星への位置ベクトルµ：GMe は地球の重力パラメーター
R：∥R∥u：衛星に働く推力

半径 R0 の円軌道上にある主衛星とその近傍の従衛星との相対運動を考える. まず, 図 1 のように主衛星の重心z̈+n2 z = uzとなり, これを,Hill-Clohessy-Wiltshire(HCW) 方程式という. HCW 方程式の状態方程式は,ẋ = Ax + Bu, x(0) = x0(3)

となり,

0
 0
 3n2A= 0
 0001000010002n00 −2n 000000000 −n2 000  0 10 ,B =  00 1 0000001000 0 0 
0 01と な る. 軌 道 面 内 と 軌 道 面 外 の 解 は 独 立 し て い る . 加
速 度 を u = [0 0 0]

Tと し, 面 内 の 初 期 値 を, x0 in =[x0 y0 ẋ0 ẏ0 ] とすると面内運動の解は,Tx(t) = 4x0 + 2y˙0 /n − (3x0 + 2y˙0 /n)cosnt + (x˙0 /n)sinnt
y(t) = y0 − 2x˙0 /n + (2x˙0 /n)cosnt +( 6x0 + 4y˙0 /n)sinnt − (6nx0 + 3y˙0 )t(4)ẋ(t) = x˙0 cosnt + (3nx0 + 2y˙0 )sinntẏ(t) = (6nx0 + 4y˙0 )cosnt − 2x˙0 sinnt − (6nx0 + 3y˙0 )Tとなる. また面外の初期値を,x0out = [z0 z˙0 ] とする面外
図 1 円軌道の図を原点とする右手系の回転座標系 o − {i, j, k} を考え運動の解は,

z(t) = z0 cosnt + (z˙0 /n)sinntż(t) = −nz0 sinnt + z˙0 cosnt(5)となる.(4)(5) の式をパラメータ表現すると,x(t) = 2c + acos(nt + α)y(t) = d − 3nct − 2asin(nt + α)z(t) = bcos(nt + β)となり,t = τ で飛びが生ずる. 次に一定時間毎にインパルス入力を入れる. 時間 kτ にインパルス uk−1 を加えると(6)し,xk = x(kτ −) とおくと,x1 = x(τ +)= eAt x0 + Bu0
xk+1 = lim x(t)となる, ここで[]1/2
a = (3x0 + 2y˙0 /n)2 + (x˙0 /n)2[]1/2
b = z0 2 + (z˙0 /n)2(13)t↘(k+1τ )=[eA(t−kτ ) x(kτ +) + eA(t−(k+1)τ ) Buk ]limt↘(k+1)τ

= eAτ xk + Buk

c = 2x0 + y˙0 /nd = y0 − 2x˙0 /nsinα = −x˙0 /nαsinβ = −z˙0 /nb
cosα = −(3x0 + 2y˙0 /n)/acosβ = z0 /b(7)(14)となる. 本研究では,1 周期ごとに λ = 4 回のインパルス入2πA λとおくと,力を考える. このとき τ = 2πλ である.Ad = exk+1 = Ad xk + Bukとなり, 入力をした直後を表す漸化式が得られる. ここで,である. 面内運動は c = 2x0 + y˙0 /n = 0 を満たす時周期解

目標軌跡を自由解,

xf = Axf , xf (0) = xf0となり, この条件を CW 条件と呼ぶ, フォーメーションフ
ライトでは従衛星をこの周期解にのせる.

fとし, インパルス入力によりこの目標を追従させる. xk =xf (kτ ) とおく

3xfk+1 = Ad xfk

インパルス制御
状態方程式,fであるから, 誤差 ek = xk − xk はẋ = Ax + Bu, x(0) = x0

ek+1 = Ad ek + Bukに時間 s をおいて, インパルス入力 us を入れた時の応答を求める. まず区間 [τ, τ + h] において一定入力 uhτ を与えた時の応答を考える. t < τ のとき,x(t) = eAt x0(8)t > τ + h のとき,4パルス制御える. 区間 [kτ − h2 , kτ + h2 ] において一定入力 uk−1 を加eA(t−r) Bu(r)drえるとし,xk = x(kτ + h2 ) とおく. このとき,0∫ τ +h= eAt x0 +eA(t−r) Bτ= eAt x0 + eAtを満たす. 従って安定化フィードバック uk = −Kd ek により誤差 ek は,0 に収束し, 目標の軌跡に追従する.

パルス制御では, インパルス入力をパルス入力に置き換∫ tx(t) = eAt x0 +(15)1h∫ τ +huτ
drhe−Ar Buτ drx0 = xh(9)τとなる. 区間 [τ, τ + h] での入力の積分は一定値 uτ となる. ここで, h を 0 に近付けたときの極限を τ におけるインパルス uτ の応答といい, これを求める.

x(t) = eAt x0 + eAt e−Aτ Buτ
At= e x0 + eA(t−τ )Buτh2= eA 2 x(0)(16)
hxk+1 = x((k + 1)τ + )2∫ (k+1)τ + h2hhAτ
= e x(kτ + ) +eA((k+1)τ + 2 −r) Buk dr
2(k+1)τ − h2
= Ad xk + Bd uk(10)これは, u(t) = uτ δ(t − τ ) を代入し形式的に t = τ の時の値を取り出した時の値に等しい. また,τ における左右の極限値は,となる. ここで,Ad = eAτ∫ (k+1)τ + h2
heA((k+1)τ + 2 −r) Bdr
Bd =(k+1)τ − h2

∫ hx(τ −) ≡ lim x(t)t↗τ= eAτ x0x(τ +) ≡ lim x(t)eA(h−r) Bdr=(11)0である. 従って, パルス入力直後の状態の漸化式は,t↘τ= eAτ x0 + Buτ(17)(12)

xk+1 = Ad xk + Bd uk(18)

fとなる. 目標軌道に関しても,xk = xf (kτ + h2 ) とおくとpulse in-plane-motion
0.01

xfk+1 = Ad xfk

0.008
0.006

hxf0 = eA 2 xf0(19)0.004

f0.002xであるから, 誤差 ek = xk − xk は,

ek+1 = Ad ek + Buk(20)

0-0.002-0.004を満たす. これより, 安定化フィードバック uk = −Kd ek-0.006

により誤差 ek は 0 に収束し, 目標の軌跡に追従する.

-0.008-0.01
0.020.015

0.010.005

最適レギュレータ

50-0.005-0.01-0.015-0.02

y図 2 パルス制御での面内運動の軌道図

本研究では, 軌道半径 R0 及び, 円軌道の角速度 n を用いtに置き換え,(x̄, ȳ, z̄) =て無次元化する. 時間 t を τ = (1/n)
yxz
1( R0 , R0 , R0 ),(ūx , ūy , ūz ) = ( R0 n2 )(ux , uy , uz ) とし,′ は τ0.01

による微分を表すと, 状態方程式は0.008

pulse out-of-plane-motion′0.006x̄ = Āx̄ + B ū, x̄(0) = x̄00.004

(21)
dz/dt

0.002

となり, ここで 0 0 3Ā =  0
000000000 0 0 1 0 0 2 0 0 
0 0 0 0 0 10 −1 0[]= diag Āin, Āout

0-0.002100−200-0.004-0.006-0.008-0.01-0.01-0.00500.0050.01z

図 3 パルス制御での面外運動の軌道図

このシステムを安定化するために最適レギュレータによる離散リッカチ方程式,

X = Q + Ad T XAd − Ad T XB(R + B T XB)−1 B T XAd
であり, 目標軌道の自由運動, 目標軌道の初期値を xf0 とすると

xfk+1 = Ad xfk , xf0 = xf0(22)となり, 漸化式 xk+1 = Ad xk + Bd uk を用いて (14) への移項を考える. 誤差 ek = xk − xfk とおいて,で, 燃料消費と比例関係にある ∆V と収束するまでの時間,整定時間 ST を用いる. 本研究では, 整定時間は, 制御軌道と目標軌道の誤差が 10−4 以下になる時間と設定した. 総速度変化 ∆V を面内運動 ∆Vin , 面外運動 ∆Vout に分けて表すと,∆Vin =(23)ūin =となり,uk = −Kek を用いて誤差 ek は 0 に収束し目標の∆Vout =軌跡に追従する.

√ū2x + ū2yST∑∥ūzk ∥k=1となる.∆V の変化を見るため, 今回最適レギュレータでシミュレーション
本 研 究 で は,

ūink

k=1

ek+1 = Ad ek − Buk , x( 0) = x0

6ST
∑従 衛 星 の 初 期 値 をT用いた重み行列 Q を Q = 10I6 で固定し, 重み行列 R を

xs0=[0.01 0 0 − 0.02 0.010]
目 標 軌 道 の 初 期 値 をT
xf 0 = [0.005 0 0 − 0.01 0.005 0] と 設 定 し た. 本 節
では,∆V の変化を見るため, パルス入力回数を 1 周期に 4回とする. 最適レギュレーターによって設計された制御系の面内運動と, 面外運動の軌道を図 2, 図 3 に表す. ここ

R = 10r I3 の r を変化させる. このときの ∆V のグラフと整定時間のグラフを図 4, 図 5 に表す. 図 4 より,r の増加に伴って操作量が制限され,∆Vin を抑えて, 燃料の消費を抑えていることが分かる. 一方, ∆Vout の変化は小さい. これは, 面外運動が r に依存せず, 早く収束してしまうからだと考えられる.

pulse r-dV-3x 10pulse ST-dV-3x 105.54.4dVin
dVout
4.2544.5dV

dV

3.8

3.643.4
3.53.23-1

3-0.500.5

11.522.530r200

400

600

800

1000

1200

1400

ST

図 4 パルス制御による r の変化に伴う ∆V の推移

図 6 パルス制御による整定時間の変化に伴う ∆V の推移

pulse r-ST

次に,r の変化に伴う整定時間の推移を比較する.

1400

ここで, 整定時間の推移をみるとパルス制御のほうが約 101200倍程度長くかかることが分かった. これは, インパルス制御
1000

では, 入力時間を 0, パルス制御では, 入力時間を与えているので, このような差ができたと考えた.

ST

800

インパルス制御とパルス制御を比較したところ, 同じ重み600だと, インパルス制御の方が良い結果が得られ,Q を大きくすると, 同等の性能フィードバックが得られることが分

400かった.

2000-1-0.5

00.5

11.5

22.5

3r

図 5 パルス制御による r の変化に伴う整定時間終りに8

本研究では, パルス制御によるフォーメーションについて研究を行った. フォーメーションフライトにより, 高精度な観測や数々のミッションを行うことが可能になる. 人工

図 5 より,r の値が大きくなると入力 u が小さく抑えら

衛星には, 燃料補給ができないため, 低燃費の運用が必要とれるため, 整定時間が大きくなり, 収束時間が長くなることなる. その問題に対して, 最適レギュレーターの離散型リッがわかる.

カチ方程式により, フィードバック制御を設計し, 最適な推力, 整定時間を制御対象に得られた. そして, これらをパル

7

性能評価比較
本章では, 前章のインパルス制御とパルス制御の総速度変化の図を比較する. まず,r の変化に伴う ∆V の推移を比較する.
面外運動の推移を見てみると, インパルス制御とパルス制
御ともに変化があまり見られないことから,r に依存していス制御とインパルス制御の 2 通りの制御方法で制御し, それらの, 性能評価を行い比較検討し, どちらの制御方法の方が燃費を向上させ, 制御できるかの比較を行った. 結果として, 同じ重みだと, インパルス制御の方が良い結果が得られ,Q を大きくすると, 同等の性能フィードバックが得られた.

ないことが見て取れる. これは, 早く収束してしまうためこのような推移になると考えた.
次に, 面外運動の推移を見てみると,r が −1 のとき ∆V はインパルス制御では,9.708 × 10−3 パルス制御では,5.04 ×10−3 とインパルス制御の方が約 2 倍大きくなっていることがわかる. これは, パルス制御の方が燃費がよいことが分かった. しかし,r が 3 の時の ∆V はインパルス制御では,3.07 × 10−3 パルス制御では,3.175 × 10−3 とほぼ近い値になっていることから, 操作量を制限することにより rが 3 以降はほぼ燃費消費量が変わらないと考えた.

参考文献[1] A. Ichikawa:Recent Developments in Formation Flying, Lecture Notes ver.1, 2010.
[2] T.Wakazono,Y.Washino:インパルス制御による円軌
道上のフォーメーション,2014(南山大学 2013 年卒業
研究)