惑星間航行・地球周回編隊飛行のための最適軌道制御理論の構築

京都大学生存圏研究所 ミッション専攻研究員
坂東 麻衣

Mar. 9. 2010

１．スペースデブリ・地球接近小惑星探査に関する多数回ランデヴーのための軌
道計画
２．燃料消費の少ないフォーメーションフライトや軌道間遷移に関する基礎理論

Mar. 9. 2010

はじめに
•

地球接近小惑星やスペースデブリを観測するための軌道設計を行うためには，従来の軌道理論では燃料消費の意味で最適な軌道を設計することが難しい．

新しい軌道理論

Mar. 9. 2010

はじめに
•

新しい軌道理論– 非線形性を直接考慮し，最適軌道を求める理論•

複数の対象へのフライバイ軌道最適化– 地球接近小惑星(Near Earth Asteroids)や地球に接近して衝突する可能性のある小惑星(Potentially Hazardous Asteroids, PHA)などを複数の対象を観測するための軌道計画
直接，運動方程式から最適な軌道を決定

軌道を生成する母関数を決定

非線形ダイナミクスを直接扱い，３体問題，摂動ポテンシャルなど様々な問題
へもポテンシャル関数の形を変えることでそのまま適用できる．
膨大な情報を効率よく扱うための軌道理論
Mar. 9. 2010

Orbital Control Theory- Classical Approach•

Lambert’s Problem- Transfer of a spacecraft from one point to another

Impulse ΔVcImpulse ΔV

決められた時間でPからQへ移動するために必要な速度を求めよ．

Mar. 9. 2010

Orbital Control Theory- Classical ApproachTwo-point boundary value problem for Keplerian motion•

Classical approach
Lambert’s problem

Lambert’s equation

Impulse ΔV

cImpulse ΔV

Lambert方程式を解くことにより，必要な速度を求めることができるが，終端位置，初期位置が変わるたびに解き直す必要がある．
Mar. 9. 2010

ダイナミクスが変わると解くことができない．

Lambert’s Problem- Direct ApproachTwo-point boundary value problem for Keplerian motion•

Direct approach (Bando and Yamakawa, 2010)
Lambert’sproblem
Hamilton-Jacobiequation

Impulse ΔV

cImpulse ΔV

Initial and final velocity is determined as
an optimal control.

作用積分を直接最小化することにより，初速度，終端速度を得る.

Mai Bando and Hiroshi Yamakawa, “A New Lambert Algorithm Using the Hamilton-Jacobi-Bellman Equations,” AIAA
Journal of Guidance, Control, and Dynamics, (to appear)

Mar. 9. 2010

Lambert’s Problem- Direct Approach
PからQに飛行時間 tf - t0 で移動するための速度を求めよ•

Hamiltonの原理：運動は，作用積分の停留値として決定される．

I[ q] = ∫ tt0f L (q, qÝ,t )dt•

最適制御問題：評価関数を最小とする制御入力を決定する．

J [ u ] = ∫ tt0f Q(x) + R (u )dt
xÝ= f (x) + g (x)utf - t0 の間にx からx へ移す評価関数を最小にする入力を求めよ0Mar. 9. 2010fLambert’s Problem- Direct ApproachTwo-point boundary value problem for Keplerian motion•Direct approach (Bando and Yamakawa, 2010)
Lambert’s
problemHamilton-Jacobi
equation

Impulse ΔV

cImpulse ΔV

Initial and final velocity is determined as
an optimal control.

作用積分を直接最小化することにより，初速度，終端速度を得る.
非線形性を直接扱い，３体問題，摂動ポテンシャルなど様々な問題へもポテンシャル関数の形を変えることでそのまま適用できる．
Mai Bando and Hiroshi Yamakawa, “A New Lambert Algorithm Using the Hamilton-Jacobi-Bellman Equations,” AIAA
Journal of Guidance, Control, and Dynamics, (to appear)

Mar. 9. 2010

New Orbital Control Theory- Generating Functions ApproachTwo-point boundary value problem for Keplerian motion•Generating Functions approach (Guibout, Scheeres, 2004)
Lambert’sproblem
Hamilton-Jacobiequation

Impulse ΔV

cImpulse ΔV

F1 (q f ,q0 ,t f ) Generating
function
F1母関数は，初期位置，終端位置の関数であるため，関数の評価を行うことで必要な速度が求まる．
Mar. 9. 2010

New Orbital Control Theory- Generating Functions Approach
Generating Functions Approach（2004, Guibout, Scheeres ）••Canonical transformation
Local solution about nominal trajectory
Local Solution about nominal trajectory

Global Solution
非線形ダイナミクスに対して直接，母関数を求める手法を提案．
Bando, M. and Yamakawa, H.: A New Optimal Orbit Control for Two-Point Boundary-Value Problem Using
Generating Functions, 19th AAS/AIAA Space Flight Mechanics Meeting, 09-117, Savannah, Georgia, 2009.

Mar. 9. 2010

運動方程式••

Newton’s equation

Hamilton’s equation: canonical variables•• Hamiltonの原理：運動は，作
Euler-Lagrange’s equation
用積分の停留値として決定さ: Hamiltonsystemれる．
: generalized
coordinate t

I[ q] = ∫ t0f L (q, qÝ,t )dt: Hamiltonian: Lagrangian••

すべての運動方程式は等価.
以下ではHamiltonの方程式を用いる．

Mar. 9. 2010

: generalized
momentumpq

正準変換
•

２つのハミルトン系の間の座標変換

New variables (Q, P)

Old variables (q, p)

Transformation equationsp?

Pq

Hamiltonian : H
Mar. 9. 2010

QHamiltonian : K

Lambert問題の解法
V. M. Guibout,D. J. Scheeres, Journal of Guidance, Control, and Dynamics 2004, vol.27 no.4 (693-704)•

固定した終端状態を変数とするハミルトン系への正準変換(Q,P) = (q f , p f )K≡0Transformation equations（Final velocity）

固定した終端状態（velocity at t）
時刻tでの状態（HJ equation）(q, p) = (q(t ), p(t ))H = T +U

•F1母関数により初期速度，終端速度が初期値，終端位置の関数として得られる．
Mar. 9. 2010

Proposed Method
Lambert’sLambert’s Problem
ProblemCanonical Transformation
Hamilton-Jacobi
Hamilton-Jacobi eq.
eq.
(Nonlinear(Nonlinear PDE)
PDE)
Step 1

Nonlinear
Nonlinear Optimal
Optimal Control
Control
Step 2The key observation we made is that F3
generating funtion solves the associated nonlinear
optimal control problem.

Successive ApproximationLinear
Linear PDE
PDE
Step 3

LinearLinear ODEODEMar. 9. 2010

Spectral method(Chebyshev Galerkin approximation)

4 types of Generating Functions
HJ equation

Boundary condition• F1, F4母関数は終端時刻(t=tf)で特異となるが，F2, F3母関数は恒等変換となり，境界条件をもつ.
Mar. 9. 2010

Hamilton-Jacobi equation forF3-generating function•

F3母関数のHamilton-Jacobi方程式
母関数のみたす偏微分方程式

Hamilton-Jacobi方程式

評価関数を最小化する問題

最適制御問題システム

運動エネルギー+ポテンシャルエネルギー
評価関数

Lagrangian• F3, F4母関数は作用積分を評価関数を最小化するpを求める最適化問題とみなすことができる．
Mar. 9. 2010

Successive Approximation
Saridis, G. N. and Lee, C.-S. G.,, IEEE Transactions on Systems, Man and Cybernetics 1979, vol.9 no.3 (152-159)•

逐次近似の評価関数

GHJB:

Linear PDE

HJB:

Nonlinear PDE•

ハミルトニアンが力学的エネルギーで表される場合，逐次近似により最適解への収束性が保証される．

Mar. 9. 2010

Spectral Method•The spectral method is based on the idea to assume the solution can beapproximated by a sum of N+1 basis functions.

Taylor series solution••

Spectral method

In the limit N →∞, the approximation must converge exponentially fast tothe exact solution.
By using spectral method we can obtain the solution in the arbitary region.
And approximation error is expect to be uniform in the region.

Mar. 9. 2010

Chebyshev Polynomials•

n-th order Chebyshev polynomial•

The Chebyshev polynomials form a complete orthogonal system on the interval[ -1, 1] with respect to the the weighting function．

Mar. 9. 2010

Proposed Method
Lambert’sLambert’s Problem
Problem
Step 3

終端速度 p fと初期位置 qの関数
0

LinearLinear ODEODEGenerating
Generating Function
Function FF

F1 (q f ,q,T )F2 (q f , p,T )F3 ( p f ,q,T )F4 ( p f , p,T )Mar. 9. 2010

F3 ( p f ,q,T )

∂F3
p0 =
( p f , q 0 ,T )∂qqf =∂F3( p f ,q 0 ,T )∂p f

p 0 = p 0 ( p f , q 0 ,T )q f = q f ( p f ,q 0 ,T )
初速度，終端位置を関数として表すことができる．

Simulation: 1-dimensional problem•

Consider 1-dimensional motion described by
Hamiltonianqf•

The equations of motion•

We solve two-point boundary-value
problem by finding F3 generating
function.

（Final position）

q0（HJ eq.）

（initial momentum）（GHJB eq.）F3(q,pf) for 1-dimensional motion)• F3 ( p f ,q,Tから求めた軌道q0 ∈ [−1, 1]vf=0•

境界条件 (終端速度vf=0)を満たすすべての軌道が，ひとつの母関数から求めることができた．

Mar. 9. 2010

F3 ( p f ,q,T )Ex) F3(q,pf) for 2-dimensional motion•

終端状態での速度誤差

Taylor series solution about circular nominaltrajectory.

Proposed method with Chebyshev
polynomials with N=2.

• 提案手法ではよりより広い領域で有効な解が得られた．
Mar. 9. 2010

Multiple Flyby Design•

従来までの問題設定との相違点?
従来の軌道設計で対象とする問題ではあらかじめ対象がきまっていた．
地球接近小惑星のように膨大な数の対象を観測するミッションを考えた場合，一度のフライトでできるだけ多くの観測を行うことが望ましい．
Mar. 9. 2010

Multiple Flyby Design•

従来までの問題設定との相違点?
フライバイごとに異なる境界条件であるため，単一のフライバイ問題を複数回解く必要がある．
膨大な数のターゲットから最適な組み合わせと順番を選ぶのは多くの計算量を必要とす
る.
Mar. 9. 2010

Assumptions•
•
•
•••
•

各時刻でのターゲットの位置は正確にわかるものとする.
ターゲットは質点としてあつかい，衝突等は考慮しない．
各軌道移行は一度のインパルスで行うものとする．
インパルスにより瞬時に速度ベクトルが変化するものとする．
最終目的地へ行くための速度変化は考慮しない．
特定のターゲットを優先して選ぶ必要はないものとする．
関心のある領域内に，十分多くのターゲットが存在するものとする．
t = t3t = t2Asteroid n3

ΔVn3 (t 3 )Asteroid n2

ΔVn2 (t 2 )Mar. 9. 2010

t = t1
ΔVn1 (t1 )Asteroid n1

N 0 : すべてのターゲットの数
N : フライバイするターゲットの数

Problem Statement•= r0
時刻 t 0に r(t 0 ) を出発した宇宙機が個の位置が既知のターゲットのうN0
ち 個のターゲットに固定した時間間隔Tでフライバイを行うとする．
N
総速度変化量：ΔV = ∑ | ΔVni (ti ) |=∑ | vni (ti+1 ) − vni (ti ) |•を最小とするターゲットの組み合わせと各フライバイ時に必要な速度変化量を求めよ．

t = t3
ΔVn3 (t 3 )t = t2T (固定した時間)

Asteroid n3

ΔVn2 (t 2 )Mar. 9. 2010

t = t1
T (固定した時間)

Asteroid n2

ΔVn1 (t1 )Asteroid n1

N 0 : すべてのターゲットの数
N : フライバイするターゲットの数

Multiple Flyby Design•

Analytical expression for total velocity change

F1 generating function:
Initial and final velocity:

Total velocity:

ΔV

F1 (q f ,q,T )∂F1v0 =
(q f ,q0 ,T )∂q
∂F1vf =(q f ,q0 ,T )∂q f

ΔV = ∑ | ΔVni (ti ) |=∑ | vni (ti+1 ) − vni (ti ) |母関数を用いることにより，必要な速度変化量は，フライバイの時刻での位置の関数として表される.
一度，母関数 F1 (q f ,q,Tをみつけておけば，すべての組み合わせで必要な速度量が直)ちに計算できる．
Mar. 9. 2010

Simulation Results•

地球近傍の軌道をもつ仮想的な１００個の小惑星に対してアルゴリズムをテスト

An exmaple of multiple flyby
ひとつの母関数により，複数フライバイを行うことに成功．
Mar. 9. 2010

Simulation Results•

Optimal sequence

100 asteroids

9,034,502,400 possibilities

Selection of 5 asteroids
９０億通りの組み合わせの中から最適な複数フライバイ軌道を決定．

Mar. 9. 2010

Data of Asteroids•

The MPC Orbit (MPCORB) Database(June 8, 2009)

離心率

0.9 AU < a < 1.1
AU

Mar. 9. 2010

長半径0< i < 5.0 [deg]軌道傾斜角Data of Asteroids•

The MPC Orbit (MPCORB) Database(June 8, 2009)t=0
Sunz˜y˜x˜

113 asteroids

16,843,743,840 possibilities

５つの小惑星を観測する場合
Mar. 9. 2010

地球円軌道に固定した回転座標系 ( x˜ , y˜, z˜ )Optimal Multiple Flyby•

１年で５回フライバイ(5T=2year)•

２年で５回フライバイ(5T=2year)• 実際には，問題に応じて時間間隔Tを適切に設定する必要がある．

Mar. 9. 2010

Conclusions•

運動の母関数を用いた新たな軌道理論の解法を提案した．

•

膨大な数の対象から，少ない燃料消費で複数フライバイを達成するためにどのターゲットにいつフライバイするのがよいのか決定するための設計方法を提案した．

•応用例– Low-Thrust– Formation flight– Space situational awareness

Mai Bando and Hiroshi Yamakawa, "A New Optimal Orbit Control for Two-Point Boundary-Value Problem Using
Generating Functions”, Advances in the Astronautical Sciences, Vol. 134, pp. 245-260, (2009).
Mai Bando and Hiroshi Yamakawa, “Orbital Design for Multiple Flyby Mission,” 27th International Symposium onSpace Technology and Science, Tsukuba, 7・2009.
Mai Bando and Hiroshi Yamakawa, “A New Lambert Algorithm Using the Hamilton-Jacobi-Bellman Equations,” AIAA
Journal of Guidance, Control, and Dynamics, (accepted).
Mar. 9. 2010

１．スペースデブリ・地球接近小惑星探査に関する多数回ランデヴーのための軌道計画
２．燃料消費の少ないフォーメーションフライトや軌道間遷移に関する基礎理論

Mar. 9. 2010

はじめに•

フォーメーションフライト– 近年，軌道上サービスや高精度の科学観測を目的として，複数の宇宙機が相互の位置などの決められた条件を保ったまま周回軌道を航行するフォーメーションフライトの研究が注目を集めている．

円軌道上

楕円軌道上
LeaderLeaderFollower

Follower
EarthEarthLeader-follower formation flightMar. 9. 2010

はじめに
•

円軌道上のフォーメーションフライト

1:2の楕円– Hill-Clohessy-Wiltshire方程式（HCW方程式）– 線形時不変システム– 周期解の条件•

楕円軌道上のフォーメーションフライト– Tschauner-Hempel方程式（TH方程式）– 線形時変システム– 周期解をもつ．

周期軌道を利用することでフライアラウンドを行うことができる
非線形方程式の周期軌道となる条件を導出し，フォーメーションフライトへの応用を考える．
Mar. 9. 2010

Hill座標系•

地球周回軌道を回っている宇宙機(Leader)からのfollowerの相対運動について考える．

リーダーに固定された座標系

x:楕円の動径方向
y:進行方向Mar. 9. 2010

Followerの運動方程式•

相対運動の方程式（非線形）Mar. 9. 2010

非線形周期軌道のための条件•

周期解となる条件
慣性系の楕円軌道

回転座標系の非線形周期軌道

軌道の周期は長半径によって決まる
回転座標系での初期値

リーダー衛星の長半径と同じ長半径をもつ楕円Mar. 9. 2010

非線形周期軌道を決めるパラメータ•周期解となる条件
慣性系の楕円軌道
離心率Mar. 9. 2010

回転座標系の非線形周期軌道

相対軌道の大きさ(形状)を決める．

慣性系

回転座標系

非線形周期軌道を決めるパラメータ•

周期解となる条件
慣性系の楕円軌道
離心ベクトルの向き

Mar. 9. 2010

慣性系

回転座標系の非線形周期軌道

相対軌道上での位置を決める．

回転座標系

非線形周期軌道を決めるパラメータ•

周期解となる条件
慣性系の楕円軌道

回転座標系の非線形周期軌道

Leader衛星がPerigee方向を通過した時刻で位置相対軌道の位置を決める．

Perigee方向
LeaderMar. 9. 2010

慣性系

Leader

回転座標系

非線形周期軌道を決めるパラメータ•周期解となる条件
慣性系の楕円軌道

回転座標系の非線形周期軌道

Target衛星がPerigee方向を通過した時刻で位置Mar. 9. 2010

慣性系

相対軌道の位置を決める．

回転座標系

非線形周期軌道を決めるパラメータ•

非線形周期軌道を決めるパラメータ

離心率

離心ベクトルの向き

Target衛星がPerigee方向を通過した時刻で位置

Mar. 9. 2010

面外運動を含めた非線形周期軌道•

３次元の場合には，軌道面での３つのパラメータに加え，相対的な軌道傾斜角がパラメータとなる．

相対的な軌道傾斜角
Mar. 9. 2010

楕円軌道上の非線形周期軌道•

楕円軌道上のはリーダー衛星に対する周期軌道

Mar. 9. 2010

回転座標系でみた非線形周期軌道の例

楕円軌道上の非線形周期軌道•

面外運動はFollower衛星の慣性座標系での軌道の２つの軸まわりの回転で表される．

回転座標系でみた３次元非線形周期軌道の例

Mar. 9. 2010

非線形周期軌道間の移行問題•

周期軌道を決めるパラメータ
離心率
離心ベクトルの向き

Target衛星がPerigee方向を通過した時刻で位置

相対的な軌道傾斜角•

周期軌道間の軌道移行問題

re1 , e1 , ϑ1 , i1e 2 , e 2 , ϑ 2 , i2

目標軌道上のVirtual Satelliteへのランデブー問題として扱う．
Mar. 9. 2010

Example 1•

リーダー衛星近傍への移行問題

３次元周期軌道を利用した移行

線形の周期軌道を利用した移行

慣性系
非線形の周期軌道を利用した移行Mar. 9. 2010

回転座標系

Example 2
初期時刻の異なる周期軌道を利用した移行•

リーダー衛星への同期問題

慣性系
Follower

Leader

Virtual Satellite

慣性系
Mar. 9. 2010

回転座標系

Conclusions•

Hill座標系での非線形方程式のすべての周期解を４つのパラメータにより生成できることを示した．

•

非線形周期軌道の初期値を近似を行うことなく完全に特徴づける方法を示した．

•応用例– 非線形周期軌道間の移行問題(reconfiguration問題)– Initialization– Satellite Constellation

Mai Bando and Akira Ichikawa, “Periodic Orbits of Nonlinear Relative Dynamics and Satellite Formation,” AIAA
Journal of Guidance, Control, and Dynamics, Vol. 32, No. 4, pp. 1200-1208, 2009.
Mai Bando and Akira Ichikawa, “Periodic Orbits of Nonlinear Relative Dynamics along an Eccentric Orbit,” AIAA
Journal of Guidance, Control, and Dynamics , Vol. 33 No. 2, pp. 385--395, 2010.
Mar. 9. 2010