(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2023-11-22
(45)【発行日】2023-12-01
(54)【発明の名称】誘導プログラム、誘導方法、及び誘導装置
(51)【国際特許分類】
G01C 21/24 20060101AFI20231124BHJP
B64G 1/24 20060101ALI20231124BHJP
G05B 11/36 20060101ALI20231124BHJP
G05B 13/02 20060101ALI20231124BHJP
G05D 1/10 20060101ALI20231124BHJP
【ＦＩ】
G01C21/24
B64G1/24 Z
G05B11/36 G
G05B13/02 A
G05D1/10
【請求項の数】
10
(21)【出願番号】P 2019128447
(22)【出願日】2019-07-10
(65)【公開番号】P2021015351
(43)【公開日】2021-02-12
【審査請求日】2022-06-16
(73)【特許権者】
【識別番号】503361400
【氏名又は名称】国立研究開発法人宇宙航空研究開発機構
(74)【代理人】
【識別番号】100106909
【弁理士】
【氏名又は名称】棚井 澄雄
(74)【代理人】
【識別番号】100161702
【弁理士】
【氏名又は名称】橋本 宏之
(74)【代理人】
【識別番号】100188592
【弁理士】
【氏名又は名称】山口 洋
(74)【代理人】
【識別番号】100181124
【弁理士】
【氏名又は名称】沖田 壮男
(74)【代理人】
【識別番号】100163496
【弁理士】
【氏名又は名称】荒 則彦
(72)【発明者】
【氏名】伊藤 琢博
(72)【発明者】
【氏名】坂井 真一郎
【審査官】松本 泰典
(56)【参考文献】
【文献】中国特許出願公開第１０６３７９５５５（ＣＮ，Ａ）
【文献】特開２０１４－１９１７３６（ＪＰ，Ａ）
【文献】特開平１１－３０１５９８（ＪＰ，Ａ）
【文献】実開平３－７０６０９（ＪＰ，Ｕ）
【文献】特開２０２０－４１８５８（ＪＰ，Ａ）
【文献】特開昭６１－７７９０５（ＪＰ，Ａ）
【文献】Ping Lu，Propellant-Optimal Powered Descent Guidance，Journal of Guidance, Control, and Dynamics，米国，American Institute of Aeronautics and Astronautics，2018年，Vol.14 No.4，pp.813-826
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０５Ｄ １／１０
Ｂ６４Ｇ １／２４
Ｇ０５Ｂ １１／３６
Ｇ０５Ｂ １３／０２
Ｇ０１Ｃ ２１／２４
(57)【特許請求の範囲】
【請求項1】
移動体に搭載されるコンピュータに、スイッチング関数理論に基づいて、前記移動体が推力を出力する方向の時間履歴を示す推力方向時間履歴と、前記移動体の推力の大きさを切り替えるタイミングを示すスイッチング時刻と、前記移動体を航行させる軌道の終端における時刻を示す終端時刻とを含む制御量を説明するパラメータの暫定的な初期値、または前記終端時刻における前記移動体の状態を示す終端状態を説明するパラメータの暫定的な初期値から、前記スイッチング時刻を導出させ、
前記導出させたスイッチング時刻と前記初期値とに基づいて、前記推力方向時間履歴、前記終端時刻、および前記終端状態のうち一部または全部を予測させ、
前記パラメータを説明変数としたスイッチング関数がゼロとなる方程式を解かせることで、前記スイッチング時刻を目的変数として導出させる、
誘導プログラム。
【請求項2】
前記パラメータには、前記移動体の運動方程式の随伴変数に関する定数と、前記終端時刻とのそれぞれを未知数としたときの暫定的な初期値が含まれる、請求項１に記載の誘導プログラム。
【請求項3】
前記コンピュータに、前記推力の加速度の上限および下限を一定とした条件下において、
２次方程式の前記スイッチング関数がゼロとなる方程式を解かせることで、前記スイッチング時刻を導出させる、
請求項１に記載の誘導プログラム。
【請求項4】
前記コンピュータに、前記推力の上限および下限を一定とした条件下において、
２次方程式に近似した前記スイッチング関数がゼロとなる方程式を解かせることで、前記スイッチング時刻を導出させる、
請求項１または３に記載の誘導プログラム。
【請求項5】
前記コンピュータに、前記スイッチング関数がゼロとなる方程式の解である前記スイッチング時刻が第１時刻の前の時刻である場合、前記スイッチング時刻を、前記第１時刻と同じ時刻に変更させ、
前記スイッチング関数がゼロとなる方程式の解である前記スイッチング時刻が第２時刻の後の時刻である場合、前記スイッチング時刻を、前記第２時刻と同じ時刻に変更させ、
前記スイッチング関数がゼロとなる方程式の解である前記スイッチング時刻が前記第１時刻以降であり、且つ前記第２時刻以前である場合、前記スイッチング時刻を変更させない、
請求項１から４のうちいずれか一項に記載の誘導プログラム。
【請求項6】
前記コンピュータに、前記予測させた終端状態と目標とする終端状態との誤差を導出させ、
前記導出させた誤差に基づいて、前記パラメータを更新させる、
請求項１から５のうちいずれか一項に記載の誘導プログラム。
【請求項7】
前記コンピュータに、前記誤差が閾値以上である場合、前記誤差に基づく感度行列を用いて補正量を導出させ、
前記導出させた補正量と前記パラメータとの和に前記パラメータを更新させる、
請求項６に記載の誘導プログラム。
【請求項8】
前記コンピュータに、前記誤差を前記パラメータで偏微分することで、前記感度行列を導出させ、
前記導出させた感度行列を用いて前記補正量を導出させる、
請求項７に記載の誘導プログラム。
【請求項9】
移動体に搭載されるコンピュータが、スイッチング関数理論に基づいて、前記移動体が推力を出力する方向の時間履歴を示す推力方向時間履歴と、前記移動体の推力の大きさを切り替えるタイミングを示すスイッチング時刻と、前記移動体を航行させる軌道の終端における時刻を示す終端時刻とを含む制御量を説明するパラメータの暫定的な初期値、または前記終端時刻における前記移動体の状態を示す終端状態を説明するパラメータの暫定的な初期値から、前記スイッチング時刻を導出し、
前記導出したスイッチング時刻と前記初期値とに基づいて、前記推力方向時間履歴、前記終端時刻、および前記終端状態のうち一部または全部を予測し、
前記パラメータを説明変数としたスイッチング関数がゼロとなる方程式を解かせることで、前記スイッチング時刻を目的変数として導出する、
誘導方法。
【請求項10】
移動体に搭載される誘導装置であって、スイッチング関数理論に基づいて、前記移動体が推力を出力する方向の時間履歴を示す推力方向時間履歴と、前記移動体の推力の大きさを切り替えるタイミングを示すスイッチング時刻と、前記移動体を航行させる軌道の終端における時刻を示す終端時刻とを含む制御量を説明するパラメータの暫定的な初期値、または前記終端時刻における前記移動体の状態を示す終端状態を説明するパラメータの暫定的な初期値から、前記スイッチング時刻を導出する導出部と、
前記導出部によって導出された前記スイッチング時刻と前記初期値とに基づいて、前記推力方向時間履歴、前記終端時刻、および前記終端状態のうち一部または全部を予測する予測部と、を備え、
前記導出部は、前記パラメータを説明変数としたスイッチング関数がゼロとなる方程式を解かせることで、前記スイッチング時刻を目的変数として導出する、
誘導装置。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、誘導プログラム、誘導方法、及び誘導装置に関する。
【背景技術】
【0002】
スペースシャトル等に起源をもつ「エクスプリシット誘導」と呼ばれる飛翔体の誘導方式は、目標へ精密に誘導するための軌道や制御コマンド、つまり最適解を実時間内でオンライン出力するアルゴリズムである。エクスプリシット誘導は、既にフライト実績も多数ある方式である。エクスプリシット誘導は、（１）解析的または数値的に終端状態量を予測する予測プロセスと、（２）予測した終端状態量が所望の値になるよう変分原理を用いて想定解を補正する補正プロセスという２つの動作を繰り返し、想定解を最適解に漸近させていくことで実現される。具体的には、例えば、制御量や終端状態を説明する未知数について、未知数と同数の連立非線形方程式を数値的に解くことで求められる。
【0003】
初期の位置および速度と、終端の位置および速度とを指定したうえで、重力一定、推力または推力加速度の上下限制約一定、大気の影響が無視できる場合において、消費燃料の最小化を目的関数とした時の最適解が、最適制御理論分野では理論的に知られている。特殊なケースを除いて、推力の大きさは、制約の上下限範囲内で、最大値－最小値－最大値のように最大２回変更する。推力の大きさの変更タイミングは、スイッチング関数と呼ばれる関数の値に基づいて決定される。推力方向は双線形タンジェント則と呼ばれる方法に基づいて決定される。
【0004】
実用的には、消費燃料最小となる最適な推力方向および終端時刻のみを求めるエクスプリシット誘導が利用されてきた。例えばスペースシャトルにおいては、推力の大きさは既知の時間関数として扱われ、エクスプリシット誘導は最適な推力方向と終端時刻のみを求める方法として利用されている。これに推力の大きさの変更タイミングを加えたエクスプリシット誘導方式は実用化されるには至っていない。
【0005】
従来は、前述のスイッチング関数の直接的な計算により誘導則を構築することは、数値安定性などの観点から困難と考えられてきた。これに関連し、スイッチング関数を直接計算するのではなく、推力の大きさの変更タイミングを最適化パラメータに加えて最適化問題をオンラインで解くことで、最適なスイッチング時刻のほか、推力方向、終端時刻を導出する誘導則が知られている（非特許文献１参照）。
【先行技術文献】
【非特許文献】
【0006】
【文献】P. Lu, "Propellant-Optimal Powered Descent Guidance," Journal of Guidance, Control, and Dynamics, vol. 41, no. 4, pp. 813--826, 2018.
【発明の概要】
【発明が解決しようとする課題】
【0007】
例えば、月や惑星などの重力天体への高精度な着陸のように、短時間内に大きな位置修正が求められる誘導問題においては、推力スイッチングを取り入れることが、消費燃料最小性及び大移動性能（高精度誘導性能）に効果を発揮する。従って、最適な推力方向、推力スイッチング、終端時刻を同時に実時間で計算するエクスプリシット誘導アルゴリズムが求められている。
【0008】
しかしながら、従来の技術では、一部の推力スイッチング時刻をオフラインで求め、そのスイッチング時刻を固定パラメータとし、残りのスイッチング時刻を他制御変数とは独立な最適化パラメータとしてオンライン最適化を実行するため、適用範囲が限定的となる仮定をする必要があった。この結果、スイッチング時刻に関する適切な固定パラメータが事前に見つけられない場合、燃料をより多く消費することになり、移動性能が低下する場合があった。また、このような課題は、スペースシャトルのような飛翔体に限られず、誘導則を適用する移動体全般に共通するところである。
【0009】
本発明は、このような事情を考慮してなされたものであり、移動体のエネルギーの消費効率を向上させつつ、移動性能を向上させることができる誘導プログラム、誘導方法、及び誘導装置を提供することを目的の一つとする。
【課題を解決するための手段】
【0010】
本発明の一態様は、移動体に搭載されるコンピュータに、スイッチング関数理論に基づいて、前記移動体が推力を出力する方向の時間履歴を示す推力方向時間履歴と、前記移動体の推力の大きさを切り替えるタイミングを示すスイッチング時刻と、前記移動体を航行させる軌道の終端における時刻を示す終端時刻とを含む制御量を説明するパラメータの暫定的な初期値、または前記終端時刻における前記移動体の状態を示す終端状態を説明するパラメータの暫定的な初期値から、前記スイッチング時刻を導出させ、前記導出させたスイッチング時刻と前記初期値とに基づいて、前記推力方向時間履歴、前記終端時刻、および前記終端状態のうち一部または全部を予測させるための誘導プログラムである。
【発明の効果】
【0011】
本発明の一態様によれば、移動体のエネルギーの消費効率を向上させつつ、移動性能を向上させることができる。
【図面の簡単な説明】
【0012】
【図1】実施形態の移動体Ｍの一例を示す図である。
【図2】実施形態の航法装置１００の構成の一例を示す図である。
【図3】天体ＰＬの目標地点に移動体Ｍを着陸させる際の降下シーケンスを模式的に示す図である。
【図4】実施形態に係る誘導演算部１２０の構成の一例を示す図である。
【図5】実施形態に係る誘導演算部１２０によるＴＥＧの一連の処理の流れの一例を示すフローチャートである。
【図6】実施形態に係る誘導演算部１２０による終端状態量の計算処理の一連の処理の流れを示すフローチャートである。
【図7】パターンＡのスイッチング関数κ（ｔ）の一例を示す図である。
【図8】パターンＡのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。
【図9】パターンＢのスイッチング関数κ（ｔ）の一例を示す図である。
【図10】パターンＢのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。
【図11】パターンＣのスイッチング関数κ（ｔ）の一例を示す図である。
【図12】パターンＣのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。
【図13】パターンＤのスイッチング関数κ（ｔ）の一例を示す図である。
【図14】パターンＤのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。
【図15】パターンＥのスイッチング関数κ（ｔ）の一例を示す図である。
【図16】パターンＥのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。
【発明を実施するための形態】
【0013】
以下、図面を参照し、本発明の誘導プログラム、誘導方法、及び誘導装置の実施形態について説明する。本実施形態における誘導装置は、移動体に搭載されたコンピュータが誘導プログラムを実行することで各種処理を行う装置である。
【0014】
図１は、実施形態の移動体Ｍの一例を示す図である。移動体Ｍは、例えば、月や惑星などの重力天体（以下、単に天体ＰＬと称する）に着陸し、その天体ＰＬの探査を行うような宇宙機（宇宙探査機）であってよい。宇宙機である移動体Ｍには、誘導装置を利用した航法装置１００が搭載される。また、宇宙機である移動体Ｍには、航法装置１００に加えて、カメラ１０と、慣性計測装置（Inertial Measurement Unit；ＩＭＵ）２０と、レーダ３０と、推進力出力装置４０とが搭載される。
【0015】
カメラ１０は、例えば、移動体Ｍの下側に設けられ、天体ＰＬに移動体Ｍが着陸する際に、天体ＰＬの地表を撮像する。「下側」とは、例えば、天体ＰＬの地表と接する脚が設けられた筐体側である。言い換えれば、カメラ１０は、天体ＰＬの重力が移動体Ｍに作用する方向（鉛直方向下向き）側に設けられる。カメラ１０は、撮像した画像を航法装置１００に出力する。
【0016】
慣性計測装置２０は、例えば、三軸式加速度センサと、三軸式ジャイロセンサとを含む。慣性計測装置２０は、これらのセンサによって検出された検出値を航法装置１００に出力する。慣性計測装置２０による検出値には、例えば、水平方向、垂直方向、奥行き方向の各加速度及び／又は角速度や、ピッチ、ロール、ヨーの各軸の速度（レート）などが含まれる。
【0017】
レーダ３０は、例えば、移動体Ｍの下側に設けられ、天体ＰＬに移動体Ｍが着陸する際に、天体ＰＬの地表に電波を照射し、照射した電波の反射波を受信する。レーダ３０による検出値には、例えば、地表に対する相対的な移動体Ｍの高度や速度などが含まれる。レーダ３０は、例えば、ドップラーレーダである。
【0018】
推進力出力装置４０は、例えば、リアクションホイールやコントロールモーメントジャイロ、軌道制御用および姿勢制御用のスラスタ（例えば化学燃料を使用したスラスタ）、推力偏向装置（Thrust Vector Controller；ＴＶＣ）などを含む。例えば、推進力出力装置４０は、リアクションホイールやコントロールモーメントジャイロなどのアクチュエータを駆動させることで地表に対する機体の向きを変化させることによって、あるいはＴＶＣによってスラスタ噴射方向を機体に対して変化させることによってスラスタが発生させた推力（推進力）の出力方向を、任意の方向に変更する。
【0019】
航法装置１００は、誘導プログラムを実行することで、例えば、宇宙機である移動体Ｍを、カメラ１０によって撮像された画像と、慣性計測装置２０の検出結果と、レーダ３０の検出結果とを利用して航行させる。以下、カメラ１０によって撮像された画像を利用した航法を「画像航法」と称し、慣性計測装置２０の検出結果を利用した航法を「慣性航法」と称し、レーダ３０の検出結果を利用した航法を「レーダ航法」と称して説明する。なお、移動体Ｍは、宇宙機に限られず、地球上において飛行するＵＡＶ（Unmanned Aerial Vehicle）といった他の移動体であってもよい。以下、一例として、移動体Ｍが宇宙機であるものとして説明する。
【0020】
図２は、実施形態の航法装置１００の構成の一例を示す図である。航法装置１００は、例えば、通信部１０２と、制御部１１０と、記憶部１３０とを備える。
【0021】
通信部１０２は、例えば、テレメータ回線などで利用されるような周波数帯の電波を用いて地球上の監視装置と無線通信する。地球上の監視装置は、例えば、通信部１０２を介して航法装置１００を遠隔制御（遠隔誘導）することで、天体ＰＬの目標位置に移動体Ｍの位置を制御することもできる。
【0022】
制御部１１０は、例えば、取得部１１２と、画像航法演算部１１４と、慣性航法演算部１１６と、レーダ航法演算部１１８と、誘導演算部１２０と、制御演算部１２２とを備える。
【0023】
制御部１１０の構成要素は、例えば、ＣＰＵ（Central Processing Unit）やＧＰＵ（Graphics Processing Unit）などのプロセッサが記憶部１３０に格納された誘導プログラムを実行することにより実現される。また、制御部１１０の構成要素のうち一部または全部は、ＬＳＩ（Large Scale Integration）、ＡＳＩＣ（Application Specific Integrated Circuit）、またはＦＰＧＡ（Field-Programmable Gate Array）などのハードウェアにより実現されてもよいし、ソフトウェアとハードウェアの協働によって実現されてもよい。また、プロセッサにより参照される誘導プログラムは、予め記憶部１３０に格納されていてもよいし、ＤＶＤやＣＤ－ＲＯＭなどの着脱可能な記憶媒体に格納されており、その記憶媒体が航法装置１００のドライブ装置に装着されることで記憶媒体から記憶部１３０にインストールされてもよい。
【0024】
記憶部１３０は、例えば、ＨＤＤ（Hard Disc Drive）、フラッシュメモリ、ＥＥＰＲＯＭ（Electrically Erasable Programmable Read Only Memory）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）などの記憶装置により実現される。記憶部１３０には、ファームウェアやアプリケーションプログラム（誘導プログラムを含む）などの各種プログラムのほかに、地形データ１３２や参照軌道情報１３４、スイッチングパターン情報１３６などが格納される。
【0025】
地形データ１３２は、例えば、移動体Ｍに着陸させる天体ＰＬの地表面の三次元形状がモデル化されたデータである。天体ＰＬの地表面の三次元形状を示すモデルには、天体ＰＬの地表面に凹凸を形成するクレータや岩などの視覚的に目立つ特徴（特徴点）が含まれる。
【0026】
参照軌道情報１３４は、参照軌道（ノミナル軌道）を定義した情報である。参照軌道とは、移動体Ｍを天体ＰＬに誘導する際に、移動体Ｍがとるべき位置、速度、姿勢などの状態量や、推力を出力する方向（推力方向）、推力を出力するタイミングを示す時刻（推力スイッチング時刻）などの制御量を、例えば、ある時間間隔ごとまたは距離間隔ごとに指定した情報、あるいは解析的に状態量や制御量を表現する数式及び参照軌道を表現するパラメータ値のセットである。
【0027】
スイッチングパターン情報１３６は、後述する複数のスイッチング関数のパターンを定義した情報である。
【0028】
取得部１１２は、カメラ１０から画像を取得したり、慣性計測装置２０から加速度や速度といった検出値を取得したり、レーダ３０から天体ＰＬの地表までの距離といった検出値を取得したりする。また、取得部１１２は、通信部１０２を介して、地上の監視装置から参照軌道情報１３４を取得し、これを記憶部１３０に記憶させることで、記憶部１３０に格納された参照軌道情報１３４を更新することができる。
【0029】
画像航法演算部１１４は、取得部１１２がカメラ１０から取得した画像と、記憶部１３０に格納された地形データ１３２とを照合（比較）して、移動体Ｍの状態量を導出する。例えば、画像航法演算部１１４は、カメラ１０の画像に対して画像処理を行って、岩やクレータなどの特徴点を抽出し、抽出した特徴点と、地形データ１３２が示す三次元モデルに含まれる特徴点とのパターンマッチングによって、天体ＰＬに対する移動体Ｍの３次元的な位置を、状態量として導出する。
【0030】
慣性航法演算部１１６は、例えば、地上の監視装置からの開始信号を受け取った後、あるいは記憶部１３０にプリセットされた開始時刻に基づいて、慣性計測装置２０の検出結果を利用した慣性航法を開始し、移動体Ｍの状態を定量的に示す状態量を導出する。「状態量」には、例えば、位置、速度、機体座標系から慣性座標系へのクォータニオンなどが含まれる。機体座標系は、例えば、移動体Ｍの重心を座標の原点とし、移動体Ｍに固定された座標系である。慣性座標系は、例えば、宇宙空間や天体に固定された座標系である。また、「状態量」には、慣性計測装置２０によって検出された加速度や角速度に対するセンサバイアス値が含まれていてもよい。
【0031】
レーダ航法演算部１１８は、例えば、天体ＰＬの地表面を基準とした移動体Ｍの高度が所定高度以下となった場合、レーダ３０の検出結果を利用したレーダ航法によって、移動体Ｍの状態を定量的に示す状態量を導出する。この状態量には、例えば、天体ＰＬの地表面からの距離、すなわち、天体ＰＬに対する移動体Ｍの相対的な位置や、相対的な速度などが含まれる。
【0032】
一般的に、地球上の監視装置との無線通信によって遠隔制御する場合、地球と天体ＰＬとの距離のオーダーが億キロメートル以上である場合があり、電波の往復伝搬時間に換算して数十分以上掛かる場合がある。そのため、天体ＰＬの地表面に対する移動体Ｍの高度を低下させ、移動体Ｍを天体ＰＬに着陸させていく際、地上からの指令に基づく遠隔制御では対応しきれないことがある。従って、制御部１１０は、自律航法を行うために、カメラ１０の画像を利用した画像航法や、慣性計測装置２０の検出値を利用した慣性航法、レーダ３０の検出値を利用したレーダ航法等を複合的に利用することで、移動体Ｍの状態量を導出する。
【0033】
誘導演算部１２０は、エクスプリシット誘導方式（エクスプリシット誘導則ともいう）を用いて、天体ＰＬの地表面の目標地点（例えば、科学的関心の強い地点や、居住に優れた地点）に到達するための移動体Ｍの軌道を予測する。具体的には、誘導演算部１２０は、現在の移動体Ｍの状態量（以下、初期状態量と称する）から、軌道終端における移動体Ｍの状態量（以下、終端状態量と称する）を、消費燃料最小で飛行する条件を満たす解を逐次的に（オンラインで）探索することで導出する。軌道終端とは、例えば、天体ＰＬの目標地点から数キロメートル上空の位置であってもよいし、天体ＰＬの目標地点そのものであってもよく、任意に与えることができる。以下、一例として、軌道終端が、天体ＰＬの目標地点の上空であるものとして説明する。
【0034】
制御演算部１２２は、画像航法演算部１１４によって導出された状態量と誘導演算部１２０によって予測された軌道（状態量）との差分量等に基づいて推進力出力装置４０を制御することで、移動体Ｍの位置、速度、加速度、姿勢などの状態を制御する画像航法を行う。また、制御演算部１２２は、慣性航法演算部１１６によって導出された状態量と誘導演算部１２０によって予測された軌道（状態量）との差分量等に基づいて推進力出力装置４０を制御することで、移動体Ｍの位置、速度、姿勢などの状態を制御する慣性航法を行う。また、制御演算部１２２は、レーダ航法演算部１１８によって導出された状態量と誘導演算部１２０によって予測された軌道（状態量）との差分量等に基づいて推進力出力装置４０を制御することで、移動体Ｍの位置、速度、姿勢などの状態を制御するレーダ航法を行う。これによって、制御演算部１２２は、移動体Ｍを、天体ＰＬの目標地点にピンポイントで着陸させる。
【0035】
図３は、天体ＰＬの目標地点に移動体Ｍを着陸させる際の降下シーケンスを模式的に示す図である。降下シーケンスには、例えば、ＰＤ（Powered Descent）フェーズとＶＤ（Vertical Descent）フェーズとが含まれる。ＰＤフェーズでは、カメラ１０を使用した画像マッチングによる画像航法と慣性航法とを複合的に組み合わせて移動体Ｍを降下させる。ＰＤフェーズ下において、制御演算部１２２は、画像マッチング時にスロットルを一時的に遮断し、推力の出力を停止してもよい。ＶＤフェーズでは、カメラ１０を使用した画像マッチングによる画像航法と、慣性計測装置２０の検出値を利用した慣性航法と、レーダ３０の検出値を利用したレーダ航法とを複合的に組み合わせて行い、ＰＤフェーズにおいて生じた鉛直方向の位置誤差を小さくしつつ、移動体Ｍを天体ＰＬの目標地点に垂直降下させる。
【0036】
例えば、地点Ｐ１から地点Ｐ２の第１区間と、地点Ｐ２から地点Ｐ３の第２区間と、地点Ｐ３から地点Ｐ４の第３区間では、ＰＤフェーズに従って画像航法が行われる。地点Ｐ４において移動体Ｍの高度が所定高度以下となった場合、ＶＤフェーズに従って画像航法、慣性航法、およびレーダ航法が行われる。
【0037】
ＰＤフェーズの第１区間から第３区間のうち、特に第３区間では、ＶＤフェーズにスムーズに移行するため、短時間（例えば百秒程度）で、数キロメートルの位置誤差を修正する必要がある。そのため、本実施形態において誘導演算部１２０は、エクスプリシット誘導方式の一種であるＴＥＧ（Throttled Explicit Guidance）を用いて、ＶＤフェーズ移行前の移動体Ｍの状態量、すなわち終端状態量を導出する。
【0038】
ＴＥＧとは、推力を調整するスロットルの制御タイミングの最適解を陽（明示的）に求めることができ、その最適解に基づいて推力の出力タイミングを切り替えることで燃料の最適化を行うことが可能なエクスプリシット誘導方式である。ＴＥＧには、後述する７つの未知数を初期化する初期化プロセスと、終端状態誤差が小さくなるように（例えばゼロに近づくように）７つの未知数を収束させる予測・補正プロセスとが含まれる。
【0039】
図４は、実施形態に係る誘導演算部１２０の構成の一例を示す図である。誘導演算部１２０は、例えば、初期化部１２０Ａと、導出部１２０Ｂと、補正部１２０Ｃと、予測部１２０Ｄとを備える。以下、これらの構成要素のＴＥＧの処理についてフローチャートを用いて詳細に説明する。
【0040】
図５は、実施形態に係る誘導演算部１２０によるＴＥＧの一連の処理の流れの一例を示すフローチャートである。本フローチャートの処理は、所定の周期で繰り返し行われてよいし、任意のタイミングで実行処理されてもよい。
【0041】
［初期化プロセス］
まず、初期化部１２０Ａは、制御量または終端状態を説明する未知のパラメータ（以下、未知数と称する）の暫定的な初期値を設定するとともに、そのパラメータから求めることが可能な終端状態量と目標終端状態量との誤差（以下、終端状態誤差と称する）を計算する（ステップＳ１００）。数式（１）は未知数の初期値（以下、初期未知数と称する）を表し、数式（２）は、終端状態誤差を表し、数式（３）は、初期化した終端状態誤差を表している。後述するように、終端状態誤差は、多次元ベクトル間の距離、すなわちユークリッドノルムによって表される。以下、アルファベットに（→）を付けたものは、ベクトルを表すものとする。また、文字の右上に記載される“Ｔ”の記号は、ベクトルまたは行列の転置を意味する。
まず、初期化部１２０Ａは、制御量または終端状態を説明する未知のパラメータ（以下、未知数と称する）の暫定的な初期値を設定するとともに、そのパラメータから求めることが可能な終端状態量と目標終端状態量との誤差（以下、終端状態誤差と称する）を計算する（ステップＳ１００）。数式（１）は未知数の初期値（以下、初期未知数と称する）を表し、数式（２）は、終端状態誤差を表し、数式（３）は、初期化した終端状態誤差を表している。後述するように、終端状態誤差は、多次元ベクトル間の距離、すなわちユークリッドノルムによって表される。以下、アルファベットに（→）を付けたものは、ベクトルを表すものとする。また、文字の右上に記載される“Ｔ”の記号は、ベクトルまたは行列の転置を意味する。
【0042】
【数1】
【0043】
【数2】
【0044】
【数3】
【0045】
初期未知数ｕ０（→）は、移動体Ｍの位置の随伴変数に関する定数νＲ（→）と、その移動体Ｍの速度の随伴変数に関する定数νＶ（→）と、軌道終端における時刻（以下、終端時刻と称する）ｔｆとを含む多次元ベクトルによって表される。定数νＲ（→）およびνＶ（→）は、それぞれ３次元ベクトルによって表される。すなわち、初期未知数ｕ０（→）は、７つの未知数を要素として含む７次元のベクトルである。
【0046】
例えば、初期化部１２０Ａは、初期未知数ｕ０（→）を予め決められた数値（例えば参照軌道を表現する値）に置き換えることで初期化し、終端状態誤差ｈ（ｕ０（→））（→）を計算する。すなわち、初期化部１２０Ａは、７つの未知数（３次元の定数νＲ（→）と、３次元の定数νＶ（→）と、１次元の終端時刻ｔｆ）の値を暫定的な値（いわゆるプリセット値）として、終端状態誤差ｈ（ｕ０（→））（→）を算出する。あるいは、初期化部１２０Ａは、当ＴＥＧの実行が２回目以降である場合、初期化として、初期状態量ｕ０（→）を前回のＴＥＧ実行時の最終的な解ｕｋ（→）とし、終端状態誤差ｈ（ｕ０（→））（→）を計算してもよい。このように初期未知数ｕ０（→）および終端状態誤差ｈ（ｕ０（→））（→）は、必ずしも最適値である保証はなく、最適値かどうか確かでない任意の値となる。初期化した７つの未知数、すなわち初期未知数ｕ０（→）は、「パラメータの暫定的な初期値」の一例である。
【0047】
数式（２）におけるｗＲ（→）、ｗＶ（→）、およびｗＨは、それぞれ位置、速度、およびハミルトニアンの終端状態誤差の重み係数を表している。終端状態誤差ｈ（ｕ（→））（→）は、終端時刻ｔｆにおける目標位置ベクトルＲｆ（→）と、未知数ｕ（→）およびスイッチング時刻（ｔ１およびｔ２）を説明変数とした終端時刻ｔｆにおける位置ベクトルＲ（ｔｆ）（→）との差に、重み係数ベクトルｗＲ（→）を掛けたものと、終端時刻ｔｆにおける目標速度ベクトルＶｆ（→）と、未知数ｕ（→）およびスイッチング時刻（ｔ１およびｔ２）を説明変数とした終端時刻ｔｆにおける速度ベクトルＶ（ｔｆ）（→）との差に、重み係数ベクトルｗＶ（→）を掛けたものと、終端時刻のハミルトニアンに重み係数ｗＨをかけたもの、の合計７次元ベクトルとして定義される。
【0048】
終端状態誤差ｈ（ｕ（→））（→）は、上述した定数νＲ
Ｔ（→）及びνＶ
Ｔ（→）と、終端時刻ｔｆと、後述するスイッチング時刻ｔ１およびｔ２の関数である。しかしながら、後述するステップＳ２００の処理を行うことで、スイッチング時刻ｔ１およびｔ２は未知数ｕ（→）の関数と見なすことができる。したがって、Ｓ２００を用いて終端状態誤差ｈ（ｕ（→））（→）を計算するとき、終端状態誤差ｈ（ｕ（→））（→）は定数νＲ
Ｔ（→）とνＶ
Ｔ（→）と終端時刻ｔｆの関数、つまりｕ（→）の関数と見なすこともできる。
【0049】
また、未知数ｕ（→）によって説明される制御量が消費燃料最小解であるための必要条件として、ハミルトニアンＨ（ｔ）が時刻ｔに沿ってゼロでなければいけないことが知られている。従って、終端時刻ｔｆにおけるハミルトニアンの目標値もゼロであり、数式（２）のハミルトニアンに関する部分は、目標値がゼロであることを考慮した表現となっている。
【0050】
次に、初期化部１２０Ａは、終端状態誤差ｈ（ｕ０（→））（→）が終了条件を満たすか否かを判定する（ステップＳ１０２）。終了条件とは、終端状態誤差ｈ（ｕ０（→））（→）が十分に小さいことであり、例えば、ある閾値εｈ（例えばゼロより大きくかつゼロに十分近い値）未満であることである。
【0051】
終端状態誤差ｈ（ｕ０（→））（→）が終了条件を満たす場合、本フローチャートの処理が終了する。終端状態誤差ｈ（ｕ０（→））（→）が終了条件を満たさない場合、後述の予測・補正プロセスの反復に移行する。
【0052】
［予測・補正プロセス］
導出部１２０Ｂは、終端状態誤差ｈ（ｕｋ（→））（→）が終了条件を満たさない場合、当反復処理ｋにおける、ｕｋ（→）に対する終端状態誤差ｈ（ｕｋ（→））（→）の感度行列Ｓｋを導出し、終端状態誤差ｈ（ｕｋ（→））（→）を減らすために、修正すべきｕk（→）の最良の方向を探索する（ステップＳ１０４）。数式（４）は、感度行列Ｓｋを表している。終端状態誤差ｈ（ｕｋ（→））（→）の予測（導出）方法については、別途サブフローチャートを用いて説明する。
導出部１２０Ｂは、終端状態誤差ｈ（ｕｋ（→））（→）が終了条件を満たさない場合、当反復処理ｋにおける、ｕｋ（→）に対する終端状態誤差ｈ（ｕｋ（→））（→）の感度行列Ｓｋを導出し、終端状態誤差ｈ（ｕｋ（→））（→）を減らすために、修正すべきｕk（→）の最良の方向を探索する（ステップＳ１０４）。数式（４）は、感度行列Ｓｋを表している。終端状態誤差ｈ（ｕｋ（→））（→）の予測（導出）方法については、別途サブフローチャートを用いて説明する。
【0053】
【数4】
【0054】
感度行列Ｓｋ（→）は、ヤコビアン行列によって表される。ｕｋ（→）は、ｋ回目の反復処理における７つの未知数を表している。例えば、導出部１２０Ｂは、数式（５）に示す線形方程式を解くことで、終端状態誤差ｈ（ｕｋ（→））（→）を減らすために、修正すべきｕｋ（→）の最良の方向を探索する。
【0055】
【数5】
【0056】
例えば、導出部１２０Ｂは、数式（５）に示す線形方程式を解く際の数値的な安定性を確保するために、７つの未知数と探索方向を正規化する。数式（６）の左辺は、正規化された７つの未知数を表し、数式（７）の左辺は、正規化された探索方向を表している。
【0057】
【数6】
【0058】
【数7】
【0059】
数式（６）におけるｕｋ，ｉおよび数式（７）におけるｄｉは、それぞれｕｋ（→）とｄｋ（→）のｉ番目の要素を表している。同様に、ｕｉ
＊とｄｉ
＊は正規化されたベクトルのｉ番目の要素を表している。また、ｕ＾ｉは、正規化に用いる定数ベクトルｕ＾（→）のｉ番目の要素を表している。ここで引数ｉは、１，２，…，７である。正規化された未知数を使用することにより、正規化された感度行列Ｓｋ
＊は、数式（８）および（９）のように表すことができる。
【0060】
【数8】
【0061】
【数9】
【0062】
まず、導出部１２０Ｂは、数式（８）に従って、終端状態誤差ｈ（ｕｋ（→））（→）を、正規化された終端状態量ｕ＊で偏微分することで、正規化された感度行列Ｓｋ
＊を計算する。次に、導出部１２０Ｂは、一次方程式である数式（９）を解くことで、正規化した探索方向ｄｋ，ｉ
＊を算出する。そして、導出部１２０Ｂは、上述した数式（７）に従って、正規化していない探索方向ｄｋ，ｉを算出する。
【0063】
具体的には、導出部１２０Ｂは、正規化された感度行列Ｓｋ，ｉ
＊の各ｉ番目の列を、数式（１０）に基づいて計算する。
【0064】
【数10】
【0065】
数式（１０）におけるεｕは、十分に小さいスカラ値であり、δｕｉ（→）は、十分に小さい偏差ベクトルである。偏差ベクトルδｕｉ（→）に含まれるそれぞれの要素δｕｉ，ｊは、数式（１１）によって定義される。引数ｊは、１，２，…，７である。
【0066】
【数11】
【0067】
導出部１２０Ｂによって数式（１０）が計算されるまえに、予測部１２０Ｄは、後述するスイッチング時刻ｔ１およびｔ２を、（ｕｋ（→）＋δｕｉ（→））と、ｕｋ（→）のそれぞれについて算出しておき、その後で、ｈ（ｕｋ（→）＋δｕｉ（→））（→）と、ｈ（ｕｋ（→））（→）とを算出する。スイッチング時刻ｔ１およびｔ２の算出方法については後述する。
【0068】
次に、補正部１２０Ｃは、（ｋ＋１）回目の処理のために、ｋ回目の反復処理における終端状態量ｕｋ（→）に含まれる７つの未知数を、導出部１２０Ｂによって算出された探索方向ｄｋ（→）を基に補正し、補正した終端状態量ｕｋ（→）を、（ｋ＋１）回目の反復処理における終端状態量ｕｋ＋１（→）に更新する（ステップＳ１０６）。数式（１２）は、状態量の更新式を表している。
【0069】
【数12】
【0070】
数式（１２）におけるαｋは、ｋ回目の反復処理におけるスケーリング係数を表している。スケーリング係数αｋは、０＜αｋ≦１の値をとる。終端状態量ｕｋが最適解に十分に近いとき、αｋ＝１であることが知られている。補正部１２０Ｃは、例えば、数式（１３）に示すように、減衰ニュートン法の要領で、終端状態誤差のユークリッドノルムが減少するようにスケーリング係数αｋを決定する。減衰ニュートン法は、連立非線形方程式を解くための一般的な方法であり、大域的に収束した解を見つけることができる。スケーリング係数αｋと、探索方向ｄｋとの積は、「補正量」の一例である。
【0071】
【数13】
【0072】
例えば、補正部１２０Ｃは、減衰ニュートン法において、αｋ，０＝１から計算を開始して、数式（１４）が満たされるまで、αｋ，０の値を小さくする。数式（１４）におけるρは、０＜ρ＜１の値をとるスカラ値である。これによって、終端状態誤差のユークリッドノルムが減少するようなスケーリング係数αｋをより高速に求めることができる。
【0073】
【数14】
【0074】
補正部１２０Ｃは、数式（１５）に示す関係を満たす場合、スケーリング係数αｋの更新を終了する。数式（１５）の終了条件は、算出された探索方向ｄｋが終端状態誤差を減少させるための良い方向ではない場合に特に有効である。
【0075】
【数15】
【0076】
補正部１２０Ｃは、スケーリング係数αｋを使用して、状態量の更新式である数式（１２）の終端状態量ｕｋ＋１（→）を更新し、状態量の更新処理を終了する。
【0077】
次に、補正部１２０Ｃは、終端状態誤差ｈ（ｕｋ＋１（→））（→）が終了条件を満たすか否か、或いは、反復回数ｋ＋１が終了条件を満たすか否かを判定する（ステップＳ１０８）。
【0078】
例えば、補正部１２０Ｃは、終端状態誤差ｈ（ｕｋ＋１（→））（→）のユークリッドノルムが閾値εｈ未満である場合、制御演算部１２２に対して、更新した未知数ｕｋ＋１（→）を出力し、本フローチャートの処理を終了する。また、補正部１２０Ｃは、反復回数ｋ＋１が所定回数ｋｍａｘを超えた場合、制御演算部１２２に対して、更新した未知数ｕｋ＋１（→）を出力し、本フローチャートの処理を終了してもよい。
【0079】
一方、補正部１２０Ｃは、終端状態誤差ｈ（ｕｋ＋１（→））（→）のユークリッドノルムが閾値εｈ以上である場合、或いは、反復回数ｋ＋１が所定回数ｋｍａｘ以下である場合、Ｓ１０４に処理を戻し、探索方向の算出処理と、終端状態量の更新処理とを繰り返す。
【0080】
［終端状態量の計算プロセス］
以下、誘導演算部１２０による終端状態量の計算処理の一連の処理の流れについてフローチャートを用いて説明する。図６は、実施形態に係る誘導演算部１２０による終端状態量の計算処理の一連の処理の流れを示すフローチャートである。
以下、誘導演算部１２０による終端状態量の計算処理の一連の処理の流れについてフローチャートを用いて説明する。図６は、実施形態に係る誘導演算部１２０による終端状態量の計算処理の一連の処理の流れを示すフローチャートである。
【0081】
まず、導出部１２０Ｂは、初期状態量ｕ０（→）が初期化されると、スイッチング関数理論に基づいて、初期時刻０から終端時刻ｔｆまでの間に、移動体Ｍが推力を出力すべきタイミングを表すスイッチング時刻を導出する（ステップＳ２００）。第２スイッチング時刻ｔ２は、第１スイッチング時刻ｔ１よりも将来の時刻を表している。
【0082】
例えば、導出部１２０Ｂは、数式（１６）および（１７）と、後述するスイッチング関数κ（ｔ）とに基づいて、一時的な第１スイッチング時刻ｔ１
＊と、一時的な第２スイッチング時刻ｔ２
＊とを導出する。
【0083】
【数16】
【0084】
【数17】
【0085】
例えば、導出部１２０Ｂは、時刻ｔ１
＊が０未満である場合、第１スイッチング時刻ｔ１を０とし、時刻ｔ１
＊が０以上、且つ最終時刻ｔｆ以下である場合、第１スイッチング時刻ｔ１をｔ１
＊とし、時刻ｔ１
＊が最終時刻ｔｆを超える場合、第１スイッチング時刻ｔ１をｔｆとする。同様に、導出部１２０Ｂは、時刻ｔ２
＊が０未満である場合、第２スイッチング時刻ｔ２を０とし、時刻ｔ２
＊が０以上、且つ最終時刻ｔｆ以下である場合、第２スイッチング時刻ｔ２をｔ２
＊とし、時刻ｔ２
＊が最終時刻ｔｆを超える場合、第２スイッチング時刻ｔ２をｔｆとする。
【0086】
［推力加速度の上限および下限を一定とした条件］
スイッチング関数κ（ｔ）は、推力加速度の上限および下限を一定とした場合、数式（１８）によって表される。
スイッチング関数κ（ｔ）は、推力加速度の上限および下限を一定とした場合、数式（１８）によって表される。
【0087】
【数18】
【0088】
数式（１８）におけるνＲ（→）は、上述したように移動体Ｍの位置の随伴変数に関する定数を表し、νＶ（→）は、移動体Ｍの速度の随伴変数に関する定数を表している。また、ｔｇｏは、終端時刻ｔｆから現時刻ｔを減算した時間（以下、残フライト時間と称する）を表している。すなわち、残フライト時間ｔｇｏは、以下の数式（１９）に示す通りである。
【0089】
【数19】
【0090】
数式（１８）に示すスイッチング関数κ（ｔ）の値をゼロとしたときの方程式は、時刻ｔに関する２次方程式となっている。スイッチング関数κ（ｔ）の形状によって、Ａ、Ｂ、Ｃ、Ｄ、Ｅの５つのパターンが存在する。これらの各パターンは、スイッチングパターン情報１３６として記憶部１３０に格納される。
【0091】
例えば、導出部１２０Ｂは、κ（ｔ）＝０という２次方程式を解き、その方程式の解を時刻ｔ１
＊およびｔ２
＊として導出する。２次方程式の解は重解であってもよい。例えば、導出部１２０Ｂは、２次方程式の解がない場合、ｔ１
＊＝ｔ２
＊を満たす値を返せばよく、例えば、数式（２０）のようにκ（ｔ）の軸を返すと好適である。スイッチング関数κ（ｔ）が０となる２次方程式を解くことには、例えば、上記のように、解がある場合にスイッチング時刻を目的変数として導出することに加えて、更に、解がない場合であってもスイッチング時刻を目的変数として導出することが含まれる。
【0092】
【数20】
【0093】
図７は、パターンＡのスイッチング関数κ（ｔ）の一例を示す図であり、図８は、パターンＡのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。図９は、パターンＢのスイッチング関数κ（ｔ）の一例を示す図であり、図１０は、パターンＢのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。図１１は、パターンＣのスイッチング関数κ（ｔ）の一例を示す図であり、図１２は、パターンＣのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。図１３は、パターンＤのスイッチング関数κ（ｔ）の一例を示す図であり、図１４は、パターンＤのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。図１５は、パターンＥのスイッチング関数κ（ｔ）の一例を示す図であり、図１６は、パターンＥのスイッチング関数κ（ｔ）に基づく推力の最適なスイッチングの一例を示す図である。いずれのパターンのスイッチング関数κ（ｔ）は、下に凸の関数である。例えば、導出部１２０Ｂは、上記のスイッチング関数κ（ｔ）から導出される推力の出力タイミングを示すスイッチング時刻を、目的変数として導出する。
【0094】
例えば、パターンＥのスイッチング関数κ（ｔ）である場合、ｔ１
＊≠ｔ２
＊となり、第１スイッチング時刻ｔ１および第２スイッチング時刻ｔ２が互いに異なる時刻となる。この場合、推力は、時刻０から最大推力または最大推力加速度で飛行し、時刻ｔ１において最小推力または最小推力加速度に切り替わり、その後時刻ｔ２のタイミングで最大推力または最大推力加速度となる。
【0095】
［推力の上限および下限を一定とした条件］
また、スイッチング関数κ（ｔ）は、推力加速度の上限および下限を一定とした条件以外でも定義できる。例えば、スイッチング関数κ（ｔ）は、推力の上限および下限を一定とした関数によって表されてもよい。数式（２１）は、推力加速度の上限および下限を一定とした場合のスイッチング関数κ（ｔ）を表している。
また、スイッチング関数κ（ｔ）は、推力加速度の上限および下限を一定とした条件以外でも定義できる。例えば、スイッチング関数κ（ｔ）は、推力の上限および下限を一定とした関数によって表されてもよい。数式（２１）は、推力加速度の上限および下限を一定とした場合のスイッチング関数κ（ｔ）を表している。
【0096】
【数21】
【0097】
式中のｃは、スラスタ（エンジン）の有効排気速度を表し、Ｍ（ｔ）は、時刻ｔにおける移動体Ｍの質量を表し、λＭ（ｔ）は、質量の随伴変数を表している。数式（２１）に示すスイッチング関数κ（ｔ）は、時間微分することにより、その変曲点の数は最大１つであることがわかる。そのため、ｔに関する方程式κ（ｔ）＝０の解は、最大でも２つである。しかしながら、数式（２１）に示すスイッチング関数κ（ｔ）は、数式（１８）に示すスイッチング関数κ（ｔ）と異なり、方程式κ（ｔ）＝０が時刻ｔに関する２次方程式となっていない。
【0098】
そこで、導出部１２０Ｂは、Ｍ（ｔ）λＭ（ｔ）＝Ｋ（定数）と近似することで、κ（ｔ）＝０を２次方程式の根の求解問題として扱う。導出部１２０Ｂは、ｔについて方程式κ（ｔ）＝０を解き、その方程式の解を時刻ｔ１
＊およびｔ２
＊として導出する。２次方程式の解は重解であってもよい。例えば、導出部１２０Ｂは、２次方程式の解がない場合、ｔ１
＊＝ｔ２
＊を満たす値を返せばよく、例えば、上述した数式（２０）のようにκ（ｔ）の軸を返すと好適である。スイッチング関数κ（ｔ）が０となる２次方程式を解くことには、例えば、上記のように、解がある場合にスイッチング時刻を目的変数として導出することに加えて、更に、解がない場合であってもスイッチング時刻を目的変数として導出することが含まれる。
【0099】
定数Ｋの値は、λＭ（ｔｆ）＝１であることが知られており、例えば、これを利用して、Ｋ＝Ｍ（ｔｆ）λＭ（ｔｆ）＝Ｍ（ｔｆ）のように、ノミナル軌道の終端質量を目安に決定してよいし、ノミナル終端質量以外の値に調整することも可能である。
【0100】
予測部１２０Ｄは、導出部１２０Ｂによって第１スイッチング時刻ｔ１および第２スイッチング時刻ｔ２が導出されると、それら時刻ｔ１、ｔ２と、７つの未知数を含む未知数ｕｋ（→）とに基づいて、終端状態誤差ｈ（ｕｋ（→））（→）を予測する（ステップＳ２０２）。
【0101】
［推力加速度の上限および下限を一定とした条件］
数式（２２）および（２３）は、推力加速度の上限および下限を一定とした条件下において、位置速度に関する終端状態を予測するための予測式を表している。
数式（２２）および（２３）は、推力加速度の上限および下限を一定とした条件下において、位置速度に関する終端状態を予測するための予測式を表している。
【0102】
【数22】
【0103】
【数23】
【0104】
式中のＶ（ｔｆ）（→）は、終端時刻ｔｆにおける移動体Ｍの速度の状態量を表し、Ｒ（ｔｆ）は、終端時刻ｔｆにおける移動体Ｍの位置の状態量を表している。Ｖ（ｔｆ）（→）およびＲ（ｔｆ）（→）は、それぞれ３次元ベクトルによって表される。そのため、終端状態量は、３次元の速度の状態量と、３次元の位置の状態量と、１次元のハミルトニアンを含む７次元のベクトルによって表される。また、ｇｍ（→）は重力ベクトルを表し、ｌ（ｔ）（→）は推力方向単位ベクトルを表す。また、ａｍａｘは最大加速度を表す。推力方向単位ベクトルｌ（ｔ）（→）を時系列に集めたものは、「推力方向時間履歴」の一例である。
【0105】
上記のＶ（ｔｆ）（→）およびＲ（ｔｆ）（→）は、数式（２４）から（２７）に従うものとする。
【0106】
【数24】
【0107】
【数25】
【0108】
【数26】
【0109】
【数27】
【0110】
なお、Ｒ０（→）は初期の位置ベクトルを表し、Ｖ０（→）は初期の速度ベクトルを表す。またａｍｉｎは、最小加速度を表す。
【0111】
推力方向単位ベクトルｌ（ｔ）（→）は、最適制御理論から導かれる最適解である、双線形タンジェント則に従う。すなわち、推力方向単位ベクトルｌ（ｔ）（→）は、以下の数式（２８）に示す通りである。
【0112】
【数28】
【0113】
数式（２８）から分かるように、未知数νＲ（→）、νＶ（→）、ｔｆを与えると、任意の時刻ｔに対して推力方向単位ベクトルｌ（ｔ）（→）が定まる。
【0114】
数式（２２）から数式（２７）にあらわれる、推力方向単位ベクトルｌ（ｔ）（→）に関する積分項は、解析的に計算してもよいし、数値的に計算してもよい。
【0115】
数式（２９）は、推力加速度の上限および下限を一定とした条件下において、ハミルトニアンに関する終端状態を予測するための予測式を表している。
【0116】
【数29】
【0117】
ただし、ａ（ｔｆ）は、終端時刻ｔｆにおける推力加速度であり、終端時刻ｔｆにおけるスイッチング関数の符号によりａｍｉｎもしくはａｍａｘのどちらかに決定される。
【0118】
［推力の上限および下限を一定とした条件］
数式（３０）および（３１）は、推力加速度の上限および下限を一定とした条件下において、終端状態量を予測するための予測式を表している。
数式（３０）および（３１）は、推力加速度の上限および下限を一定とした条件下において、終端状態量を予測するための予測式を表している。
【0119】
【数30】
【0120】
【数31】
【0121】
ただしＴｍａｘは、最大推力である。また上記のＶ（ｔｆ）（→）およびＲ（ｔｆ）（→）は、数式（３２）から（３５）に従うものとする。
【0122】
【数32】
【0123】
【数33】
【0124】
【数34】
【0125】
【数35】
【0126】
ただしＴｍｉｎは、最小推力である。推力方向単位ベクトルｌ（ｔ）（→）および移動体Ｍの質量Ｍ（ｔ）に関する積分項は、解析的に計算してもよいし、数値的に計算してもよい。
【0127】
数式（３６）は、推力の上限および下限を一定とした条件下において、ハミルトニアンに関する終端状態を予測するための予測式を表している。
【0128】
【数36】
【0129】
ただし、ｍ（ｔｆ）（・）は、終端時刻ｔｆにおける、単位時間あたりの燃料消費量である。なお、時刻ｔにおける単位時間あたりの燃料消費量ｍ（ｔ）（・）と推力Ｔ（ｔ）の間には、有効排気速度ｃを用いて数式（３７）に示す関係がある。
【0130】
【数37】
【0131】
終端時刻ｔｆにおける燃料消費量ｍ（ｔｆ）（・）は、終端時刻ｔｆにおけるスイッチング関数の符号により、推力下限もしくは上限に相当する単位時間あたりの燃料消費量ｍｍｉｎ（・）またはｍｍａｘ（・）のどちらかに決定される。
【0132】
以上説明した実施形態によれば、移動体Ｍに搭載される航法装置１００が、スイッチング関数理論に基づいて、初期化した初期未知数ｕ０（→）から、第１スイッチング時刻ｔ１と第２スイッチング時刻ｔ２を求めたうえで、更に、その第１スイッチング時刻ｔ１および第２スイッチング時刻ｔ２と初期未知数ｕ０（→）の関数と見做して終端状態量誤差ｈ（ｕ０（→））（→）を予測する。終端状態誤差ｈ（ｕ０（→））（→）が終了条件を満たす場合、誘導処理が終了する。満たさない場合、予測・補正プロセスの反復に移行する。
【0133】
航法装置１００は、ｋ回目の反復処理における終端状態誤差ｈ（ｕｋ（→））（→）を予測し、その終端状態誤差ｈ（ｕｋ（→））（→）の感度行列Ｓｋを導出し、終端状態誤差ｈ（ｕｋ（→））（→）を減らすために修正すべきｕｋ（→）の最良の方向ｄｋ，ｉを探索する。
【0134】
航法装置１００は、ｋ回目の反復処理における未知数ｕｋ（→）を探索方向ｄｋ，ｉを基に補正し、補正された未知数を、（ｋ＋１）回目の反復処理における終端状態量ｕｋ＋１（→）として更新する。
【0135】
航法装置１００は、数式（１５）に示す関係が満たされるまでスケーリング係数αｋの更新を繰り返し、数式（１５）に示す関係が満たされるとスケーリング係数αｋの更新を終了する。航法装置１００は、更新したスケーリング係数αｋを使用して、終端状態量ｕｋ＋１（→）を更新する。
【0136】
そして、航法装置１００は、更新した未知数ｕｋ＋１（→）が終了条件を満たす場合、その未知数ｕｋ＋１（→）から導出される、スイッチング時刻ｔ１、ｔ２、推力方向ｌ（ｔ）（→）、終端時刻ｔｆなどの制御量と、位置、速度、姿勢などの状態量とに基づいて、推進力出力装置４０を制御することで、目標とする軌道終端に移動体Ｍを精度よく移動させることができる。この結果、移動体Ｍのエネルギーの消費効率を向上させつつ、移動性能を向上させることができる。
【0137】
例えば、着陸重量２００ｋｇ程度の着陸機を仮定した場合、上述したＶＤフェーズのように、天体ＰＬへの降下終盤フェーズでは、およそ２０ｋｇ程度の燃料を消費する場合がある。ＴＥＧを利用した場合、２０ｋｇの燃料のうち、約３．６％程度の燃料の消費効率が向上する。すなわち、０．７２ｋｇ程度の燃料を節約することができる。着陸重量２００ｋｇの着陸機に搭載可能なペイロードの重量は数ｋｇ程度以下であることが典型的であり、０．７２ｋｇは、１、２割程度であることから、積載重量を増やすという点においても非常に効果的である。
【0138】
［変形例］
以下、上述した実施形態の変形例について説明する。上述した実施形態では、移動体Ｍが宇宙機であり、その移動体Ｍに搭載されるセンサが、カメラ１０や慣性計測装置２０、レーダ３０であるものとして説明したがこれに限られない。例えば、移動体Ｍが上述したＵＡＶである場合、その移動体Ｍに搭載されるセンサには、カメラ１０やレーダ３０に代えて、あるいは加えて、ＬＩＤＡＲ（Laser(Light) Imaging Detection and Ranging）が含まれてよい。
以下、上述した実施形態の変形例について説明する。上述した実施形態では、移動体Ｍが宇宙機であり、その移動体Ｍに搭載されるセンサが、カメラ１０や慣性計測装置２０、レーダ３０であるものとして説明したがこれに限られない。例えば、移動体Ｍが上述したＵＡＶである場合、その移動体Ｍに搭載されるセンサには、カメラ１０やレーダ３０に代えて、あるいは加えて、ＬＩＤＡＲ（Laser(Light) Imaging Detection and Ranging）が含まれてよい。
【0139】
また、位置および速度に関する終端状態量の拘束の一部を外してもよい。例えば、２点間の軌道が物理的に困難なために条件緩和した軌道を生成したい場合がある。この場合、一部の終端状態量の拘束が必ずしも必要でない場合（例えば、円軌道投入はダウンレンジの拘束が必要なくなる）、終端状態量の拘束の一部を外してよい。また、終端状態量の拘束を外す場合、除外する終端状態量に付随する随伴変数に関する定数ν（未知数）をゼロとしたうえで、該当する未知数と条件式を非線形連立方程式から外してよい。
【0140】
また、上述した実施形態では、推力加速度または推力の上下限を一定としたがこれに限られない。例えば、推力加速度または推力の上下限値を、「制約」ではなく「目標」としてＴＥＧを行ってもよい。これにより、例えば、制約として扱うと物理的に軌道生成が不可能な場合でも、推力（加速度）の上下限制約を緩和することで物理的に成立する軌道を生成する、といったより実用的な誘導ができるようになる。
【0141】
終端状態誤差ｈ（ｕ（→））（→）の連立非線形方程式を解くのを「制約」として扱う場合、ｈ（ｕ（→））（→）＝０の制約条件のもと、ｆ＝ρｍｉｎ（Ｔｍｉｎ－Ｔｍｉｎ，ｔ）２＋ρｍａｘ（Ｔｍａｘ－Ｔｍａｘ，ｔ）２を最小化（推力の上下限制約の場合）、または、ｆ＝ρｍｉｎ（aｍｉｎ－aｍｉｎ，ｔ）２＋ρｍａｘ（aｍａｘ－aｍａｘ，ｔ）２を最小化（推力加速度の上下限制約の場合）するのが「目標」として扱ってよく、これを非線形計画問題として求解する。ρは重みパラメータ、Ｔｔは推力の、aｔは推力加速度の目標値である。
【0142】
この際、「制約」から「目標」に扱いを変えるときの変更箇所は、推力（加速度）の上下限値に関するパラメータ（２つ）を未知数に追加したうえで、ＴＥＧの処理フローにおいて、例えば「非線形連立方程式の解法（例えば減衰ニュートン法）」を「等式制約つき最適化問題の解法（例えばラグランジュの未定乗数法）」に置き換える点である。
【0143】
以上、本発明を実施するための形態について実施形態を用いて説明したが、本発明はこうした実施形態に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々の変形及び置換を加えることができる。
【符号の説明】
【0144】
１０…カメラ、２０…慣性計測装置、３０…レーダ、４０…推進力出力装置、１００…航法装置、１０２…通信部、１１０…制御部、１１２…取得部、１１４…画像航法演算部、１１６…慣性航法演算部、１１８…レーダ航法演算部、１２０…誘導演算部、１２０Ａ
…初期化部、１２０Ｂ…導出部、１２０Ｃ…補正部、１２０Ｄ…予測部、１２２…制御演算部、１３０…記憶部、Ｍ…移動体
…初期化部、１２０Ｂ…導出部、１２０Ｃ…補正部、１２０Ｄ…予測部、１２２…制御演算部、１３０…記憶部、Ｍ…移動体
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】