楕円制限三体問題におけるハロー軌道の維持制御
2015SC013

坂井祐介

指導教員：大石泰章

1

通重心 Ob に原点をとり，地球から月に向かう方向に x 軸はじめにを，地球-月の軌道面の中で x 軸に垂直な方向に y 軸を，軌

ハロー軌道とは，制限三体問題の力学的平衡点の周りを道面に対して垂直な方向に z 軸をとる．楕円制限三体問題周回する周期軌道であり，この軌道上に宇宙機を周回させでは，二天体は共通重心 Ob の周りに楕円運動を行うと仮ることで，2 天体の同時通信や観測が行えることから注目

定し，宇宙機は二天体から影響を受けて運動する．独立変を集めている．しかし，不安定周期軌道であることから， 数としては，時間の代わりに真近点離角 ν すなわち，共通宇宙機を軌道上に維持するには制御が必要である．本研究 重心を焦点とし，楕円軌道上で運動する月が，近地点からでは月，地球と，宇宙機による三体問題を扱う．月や地球どれだけ離れているかを表す角度を用いる (図１を参照)．の質量に対して宇宙機の質量は無視できるほど微小であ

また月が近地点にいる場合の地球と月の間の距離 a(km)り，宇宙機は 2 天体により生成される重力場内を運動してを基準とした無次元化を行う．また二天体の総質量に対いると考えることができる．しかしこの 2 天体の運動は楕し，月の質量比を ρ と表す．さらに原点を共通重心から力円運動であり，離心率が 0.0542 と比較的大きいため，こ学的平衡点である L2 点 (l２ , 0, 0) に移す．そのときの宇宙れを考慮した楕円制限三体問題を考える必要がある．この機の座標を (x, y, z) とすると，運動方程式は以下のように

楕円制限三体問題における重力場の力学的平衡点はラグラ表される：

ンジュ点と呼ばれ，L1 から L5 まで存在し，それぞれにハ

x′′ − 2y ′ − ( 1+e1cos ν )x1
 y ′′ + 2x′ − (=1+e cos ν )ye cos ν′′z − ( 1+e cos ν )zl2 − 1−ρ(x + l2 + ρ)r13ρ− r3 (x + l2 − 1 + ρ) 1
2
.
1−ρρ
−y−
y1 + e cos ν 
r13r23
ρ
− 1−ρ
z−
z33rr

ロー軌道をもつが，楕円運動により 2 天体の距離間が周期的に変動するため，ラグランジュ点の位置も同様の周期で変動する．よってそのラグランジュ点の変動周期にあったハロー軌道しか存在せず，軌道生成自体が難しい．
本研究では, 楕円制限三体問題における L2 点周りのハロー軌道に宇宙機を維持することを考える．制御器設計
は, 最適レギュレータ理論に基づいて行う．制御性能の評

1価には, 燃料消費に比例する加速度変化の総和を用い, 最2

ただし,()′ は ν に関する微分をあらわし，適レギュレータの入力の重みを変化させることで, 燃料消√
r1 = (x + l2 + ρ)2 + y 2 + z 2 ,√r2 = (x + l2 − 1 + ρ)2 + y 2 + z 2

費の最も少ないフィードバックゲインを求める．連続時間制御の場合とパルス制御の場合の二つを扱うが，スペースの都合で本稿ではパルス制御の場合のみを述べる．先行研究である文献 [2] のパルス制御の手法と比較するとともである．e は離心率である．原点を L2 点に移動させてに，入力の重み行列だけでなく，パルス回数，入力幅も変いるため，この回転座標系においてはラグランジュ点は化させることで，最小の燃料消費を示すシステムの設計を不動であり，平衡点として定義できる．また微分方程式行った．

中に cos ν といった項が存在することから，システムは2周期 2π の周期性を持つ．以降での扱いのため状態量を運動方程式

q = [x y z x′ y ′ z ′ ]T として以下のように表す：
q ′ = f (q).

3(1)

ハロー軌道
円制限三体問題のハロー軌道は Ricardson らによって近似的な周期解を求めることができた [3]．しかし，独立変数が方程式中に存在する楕円制限三体問題では周期解が 2πの整数倍でなければならないという制約があることから直接生成することが難しい [4]．よって，周期 π の円制限三
図 1 楕円制限三体問題の位置関係体問題のハロー軌道 2 周分に離心率を付加する形で楕円制限三体問題のハロー軌道を生成した [5]．このハロー軌道

図 1 に地球，月，宇宙機からなる系を示す．ただし，M1は，2 天体が近地点に存在する時を初期位置としたハローは地球，M2 は月，m は宇宙機を表している．地球と月の共軌道 (left halo orbit) と 2 天体が遠地点にいる時を初期位1

置としたハロー軌道 (right halo orbit) の 2 種類のみが存である．

在する．実際に運動方程式中に cos ν といった項が含まれ4.2ることから，この項の変動周期に合わせた軌道，つまり周離散化表現

前章により得られた宇宙機の軌道方程式を N 点に離散期が 2π で，かつ xz 平面に対して対称になるハロー軌道化する．離散化の幅は ν に関して等分であり，したがっしか存在しない [6]．
本研究では，近地点を初期位置としたハロー軌道を用いて w = 2π/N である．初期時刻 ν0 を起点として離散時刻て宇宙機の軌道維持制御システムの設計を目指す．その軌

νk = ν0 + kw を定める．各離散時刻 νk から入力幅 l の間は入力を入れ，その後次の離散時刻 νk+1 までの自由運動
幅 s = w − l の間は入力を入れない (図 3 参照)

道を図 2 に示す．

0.1
L2 point
ν=0left halo orbitz方向0図 3 離散化の定義-0.1-0.2-0.2

N 点に離散化したうちの k 番目の時刻 νk から，w 進
0-0.1
0.2-0.05

んだ時のシステムの変化を考える．まず，次の式を満たす

x方向

6 × 6 の行列値関数 Φk (τ ) を考える：

0y方向

Φ′ (τ ) = A(νk + τ )Φ(τ )図 2 近地点を初期時刻としたハロー軌道 (ν0 = 0)(5)(0 ≤ τ ≤ w)
k = 1, 2, 3...N − 1, N

4

制御器設計

4.1

時変システム行列

τ = w におけるこの関数の値を離散化後のシステム行列とみなす：

動特性を考慮したシステム行列を生成するために，前章で得られた周期軌道を qf (ν) とし，そこからの微小な誤差

Ad (k) = Φk (w)qe (ν) の振る舞いを考える．qf (ν) も qf (ν) + qe (ν) もともに式 (1) を満たすので

次に，離散化後の入力行列を求める．次の式を満たす(6)

6 × 3 の行列値関数 Ψk (τ ) を考える：

d(qe ) = f (qf + qe ) − f (qf )dνが成り立つ．この式をハロー軌道上の点 qf (ν) の周りで線

Ψ′k (τ ) = A(νk + τ )Ψk (τ ) + B
Ψ′k (τ ) = A(νk + τ )Ψk (τ )形近似することで，以下のような時変システムを得ること(0 ≤ τ ≤ l)(7)(l ≤ τ ≤ l + s)(8)

ただし，ができる:Φk (0) = O6×3 ．
d(qe ) = A(ν)qe ．
dν(2)よって τ = w におけるこの関数の値を離散化後の入力行列とみなす．

ただし,Bd (k) = Ψk (w)．

∂f= A(ν)．
∂q q=qf (ν)

以上により，離散化した宇宙機の軌道方程式

qf の周期性より，A(ν + 2π) = A(ν) であることから, このシステム (2) は周期性を持つ．この時変システム行列
A(ν) を用いて宇宙機の軌道方程式を表すと；
qe′ = A(ν)qe + N (qe , ν) + Bu

qe (k + 1) = Ad (k)qe (k) + Bd (k)u(k)4.3

力 u は 3 方向すべてに入力を入れることができると仮定しB=O3×3I3×3

差分方程式の安定化

4.2 節で得た宇宙機の方程式ていない非線形項であり，Bu は入力項である．ここで入[(10)を得る．ただし，qe (νk ) を qe (k) と書いている．

(3)と表すことができる．ただし，N (qe , ν) は A(ν)qe で表せている．したがって入力行列 B は(9)

qe (k + 1) = Ad (k)qe (k) + Bd (k)u(k)の制御を考える．ここで入力 u(k) を](4)

u(k) = −Kd (k)qe (k)2(11)

とすることで，差分方程式はqe (k + 1) = (Ad (k) − Bd (k)Kd (k))qe (k)4.5

本研究の提案手法によるシミュレーション結果(12)

パルス制御では, 加速度の総和である ∆V で性能を評価となる．ここで注意しなくてはならないのは，連続時間シする．∆V は，消費燃料に比例関係であることが分かっており，u = [ux uy uz ]T とすると，∆V は以下のように表ステム同様，システム行列の動特性が周期的に変化するたされる．

め，フィードバックゲインも時間とともに変化することである．よって適切な時変のフィードバックゲイン Kd (k)J∆V =を考えることが自然である．

N−1∑(uxj 2 + uyj 2 + uzj 2 )1/2 l.

j=0

離散時間行列リッカチ方程式

宇宙機の初期状態は，ハロー軌道上 (qe (0) = 0) にある評価関数を，次のように定める:ものとし，軌道を 50 周させたうちの１周期あたりの平均4.4

JN (u) =N−1∑

[wqeT (k)Qqe (k) + luT (k)Ru(k)]の ∆V で性能を評価する．
本研究では，離散時間の枠組みで制御機設計を行ってい(13)るため，より高い性能が期待できる．また，パルス回数とk=1

ここで Q は状態にかかる重み，R は入力にかかる重みで入力の重み行列を同時に変化させ，最小の ∆V を示す制御ある．各離散幅 w ごとのシステム行列 Ad (k) と入力行列器の生成を行うこと，およびパルス回数を固定し，パルスBd (k) を用いて, 軌道上での動特性を考慮した設計法を考える．ハロー軌道の離散回数が N なので Ad (N ) = Ad (0)幅と入力の重み行列を同時に変化し，最小の ∆V を示す制御器の生成を行うことを考える．
評価関数 (13) 式の入力の重み行列をである．よって，評価関数を最小化する時変フィードバッ
クゲイン Kd (k) も離散回数 N の周期関数となる．Kd (k)R = 10r × I3

は，離散時間行列リッカチ方程式とおき，r を ‐ 3 から 1 まで変化させ，1 周期あたりのパPd (0) = O6×6 ，ルス回数を 3, 4, 5, 6, 8, 10, 20, 30, 40, 50, 60 回にし，Pd (k) = Q + ATd (k)Pd (k + 1)Ad (k) − ATd (k)Pd (k + 1)
Bd (k)(R + Bd (k)T Pd (k + 1)Bd (k))−1
BdT Pd (k + 1)Ad (k)1
入力幅はそれぞれの離散幅の 10の長さとしてシミュレーションを行った．図 5 と表 2 にそのシミュレーション結果を示す．

(14)の周期解 Pd (k) を用いることで Kd (k) = (R + BdT Pd (k +−1-4x 10

BdT (k)Pd (k + 1)Ad として得られる [1]．ただ1)Bd (k))
し, 離散時間行列リッカチ方程式は，解析的に周期解を求めることが難しい [1]．そこで，終端条件 Pd (0) = O ∈ R6×6から出発して，逆時間方向に数値解 P̂ (k) を求め，周期関数に収束したとき，これを Pd (k) とする．計算の模式図を図 4 に示す．

∆ V151050604010 20
frequency N-2values of-40

the exponent rof pulse inputs

図 5 入力の重み行列の指数とパルス回数を同時に変化させたときの ∆V の推移
図 5 と表 2 からわかるように，最小の ∆V を示したのはパルス回数が 3 回のときであった．これは入力を入れる回数が少ないからだと考えられる．しかし，パルス回数が
図 4 離散時間行列リッカチ方程式

3 回の場合では，安定化できる入力の重み行列の範囲が狭く，r = ‐ 0.5 以上になると，リッカチ差分方程式の解がある正整数 n に対して，P̂ (−nN ) − P̂ (−(n − 1)N ) の全収束しなかった．これは，離散幅が長くなると，もともと成分の絶対値が 10−7 以下となった場合に，0 ≦ k ≦ N − 1の連続時間システムの離散化誤差がおおきくなることが原に対しては，Pd (k) = Pˆd (−nN + k) となるように，それ以外の k に対しては，Pd (k) が周期 N の周期関数になるように，関数 Pd (k) を定める．

因であると考えられる．パルス回数を増やしていくと ∆Vの値は大きくなる．60 回以上にすると ∆V は単調減少になり，連続時間制御と同様の傾向を示す．

3

システムから得られたフィードバックゲインをパルス制御

表 1 各パルス回数における最小の ∆V
パルス回数

入力重み

入力幅

∆V

3回
4回
5回
6回
8回
10 回
20 回
30 回
40 回
50 回
60 回

r = −0.875
r = −0.625
r = −0.5
r = −1r = −0.875
r = −0.875
r = −1.625
r = −1.625
r = −1.5r ≥ −0.875r ≥ −0.875

0.2[rad]0.15[rad]0.12[rad]0.1[rad]0.08[rad]0.06[rad]0.03[rad]0.02[rad]0.015[rad]0.013[rad]0.01[rad]3.032 × 10−4
3.499 × 10−4
4.147 × 10−4
4.921 × 10−4
4.989 × 10−4
4.972 × 10−4
7.443 × 10−48.636 × 10−48.700 × 10−4
8.185 × 10−4
6.780 × 10−4に応用し，連続時間制御と同等の効果を出すため，入力を開始する時刻のフィードバックゲインに一つの離散化幅に対して入力の割合を決めるデューティー比の逆数をかけることでパルス制御を行っていた．しかしこの手法では，パルス回数は 1 周期（約 27 日）に 60 回のパルス入力を必要とし，一定以上の入力回数と一定以上の入力幅がなければ安定化が困難であった．文献 [2] の手法におけるシミュ
レーション結果を以下の表 3 に示す．
表 3 文献 [2] の手法によるシミュレーション結果パルス回数

入力幅

入力の重み

∆V60 回

0.05radr = 0.2

1.42 × 10−3

次にパルス回数を 1 周期あたり 30 回（1 日に 1 回の入力）に固定し，各離散幅における入力幅を 0.025[rad] から0.3[rad] まで変化させ，入力の重み行列の指数 r を −4 から 1 まで変化させた時の ∆V の変化を図 6 に示す．

5

おわりに
本研究では, 楕円制限三体問題のハロー軌道を生成し, 近地点を初期時刻としたハロー軌道への宇宙機の維持制御を行った．宇宙機の運動方程式を離散化することで，文献 [2]-3x 10より性能のいいパルス制御を行うことができた．チューニ
ングパラメータが入力の重み，パルス回数，入力幅と 3 つあるため，消費燃料を小さくすることが困難である．しか∆ V2し，傾向として，少ないパルス回数でかつ，短い入力幅でも入力の重みを適切に調節すれば良好な制御が可能である

1ことが確認できた．
0参考文献

212Width of 60-2values of0[1] Akira Ichikawa，Hitoshi Katayama: Linear timevarying systems and sampled-data systems: Lecturenotes in Control and Infomation Science 265-4pulse inputs
exponent r

図 6 入力の重み行列の指数と入力幅を同時に変化させた
ときの ∆V の推移[2] 坂井 祐介， 宇佐美元啓，大石 泰章: 発散度合いを考慮したハロー軌道の維持制御:第 3 回マルチシンポジウム
予稿集，（2016）
最小の ∆V を示したのは，パルス幅が 0.025[rad] で，入力の重み行列の指数 r = −1.5 の時であり，∆V =少なくなるが，入力の最大値は大きくなるため，入力の重[3] David L.Richardson：Analytic Construction of periodic orbit about the collinear points：Celestial Mechanics 22 (1980) 241-253

み行列をある程度小さくしなければ安定化することが難し[4] 幅口雄太, 長田圭介, 坂本登： 楕円制限三体問題におけいといえる．逆に入力幅を長くすると，入力の最大値は小る準ハロー軌道計算に対する中心多様体法の適用：第さくて済むが，長い時間入力を入れなければならないこと

58 回自動制御連合講演会予稿集, (2014)−48.583 × 10であった．入力幅が短いほど，必要な ∆V はから，∆V が増加する．

[5] 幅口雄太, 坂井祐介, 石田翔也, 坂本登： 楕円制限三体問題における準ハロー軌道計算：第 60 回自動制御連合
講演会予稿集, (2016)

表 2 パルス回数 30 回における最小の ∆V
パルス回数

入力幅 [時間]

入力重み∆V30 回

2.8 時間r = −1.58.583 × 10−44.6[6] 幅口雄太：制限三体問題における不変多様体解析と最適軌道計画：名古屋大学大学院 2015 年度修士論文

連続時間制御器の近似によるシミュレーション結果

制御対象を近地点を初期位置としたハロー軌道 (lefthalo orbit) とし，文献 [2] で提案した手法では，連続時間4