JPWO2019044735A1 - 宇宙機制御装置、宇宙機制御方法、およびプログラム - Google Patents

Abstract

角運動量算出部（１２）は、宇宙機の角運動量を算出する。指令値算出部（１３）は、宇宙機の軌道を定める軌道要素の変化量と軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致する第１の制約条件の下で、複数の噴射区間におけるスラスタ（２１，２２，２３，２４）の噴射量の合計が目的関数である非線形計画問題の解を求めることで、スラスタ指令値および角指令値を算出する。スラスタ制御部（１４）は、スラスタ指令値に基づき、スラスタ（２１，２２，２３，２４）の噴射を実行する。ジンバル制御部（１５）は、角指令値に基づき、ジンバル機構（３１，３２，３３，３４）の角度を調節する。

Description

本発明は、天体の回りの軌道を周回する宇宙機の軌道保持制御および角運動量のアンローディングを行う宇宙機制御装置、宇宙機制御方法、およびプログラムに関する。

宇宙機には、宇宙機の姿勢を制御する姿勢制御アクチュエータ、および宇宙機の経度および緯度を所望の範囲に維持するスラスタが搭載される。宇宙機は、例えば、静止衛星、地球観測衛星、宇宙船等である。姿勢制御アクチュエータの一例として、リアクションホイール、ＣＭＧ（Control Moment Gyro）等が宇宙機に搭載される。宇宙機には、太陽輻射トルク、磁気トルク等の種々の外乱トルクが働く。そこで、リアクションホイールを回転させて、外乱トルクを吸収して、宇宙機の姿勢を所望の向きに維持する必要がある。リアクションホイールの回転数には限度があるため、リアクションホイールの回転数の増加を抑え、リアクションホイールの回転数を正常範囲内に維持することが行われる。これを姿勢制御アクチュエータに蓄積された角運動量のアンローディングという。詳細には、スラスタを用いてアンローディングトルクを発生させることで、アンローディングが行われる。

宇宙機の経度および緯度を所望の範囲に維持する軌道保持制御と、角運動量のアンローディングとを行う技術が特許文献１，２に開示されている。特許文献１に開示される衛星に搭載されたスラスタは、軌道の下降ノードの近傍、および上昇ノードの近傍で噴射される。特許文献２に開示される制御システムは、超過運動量をアンロードするように、宇宙船の姿勢およびスラスタへの入力に対する制約を条件とした後退ホライズンにわたるコスト関数の最適化を用いて、宇宙船のスラスタおよび宇宙船の運動量交換デバイスを同時に制御する制御入力コマンドを生成する。

特開平９−２７７９９７号公報 特開２０１６−１２４５３８号公報

上述の軌道保持制御と角運動量のアンローディングとを行うためには、推薬が必要である。推薬が枯渇すると、宇宙機は軌道保持制御および角運動量のアンローディングを行うことができない。そのため、軌道保持制御と角運動量のアンローディングとにおける推薬の使用量の低減が求められている。

一方で、宇宙機の通信能力を増大させるために宇宙機の発生電力の大電力化が求められ、宇宙機の太陽電池パドルが大型化する傾向にある。太陽電池パドルが大型化すると、宇宙機に働く太陽輻射圧の影響が増大する。その結果、太陽輻射圧の影響によって増大する衛星の軌道の離心率を補償するために、衛星の経度を所望の範囲内に維持する東西制御の制御量が増大する。また宇宙機に働く太陽輻射圧トルクも増大するため、宇宙機に搭載された姿勢制御アクチュエータの角運動量が飽和しやすくなる。この場合、例えば、昇交点近傍および降交点近傍の二箇所、または、特許文献１に開示されるように、軌道の下降ノードの近傍および上昇ノードの近傍において、軌道保持制御と角運動量のアンローディングとを行うことは、推薬の使用量の削減の観点から、最適ではないことがある。

特許文献２に開示される制御システムは、後退ホライズンにわたるモデル予測制御（Model Predictive Control：ＭＰＣ）を用いて、宇宙船の動作の最適化を行う。この制御システムは、宇宙船の現在の状態、宇宙船の現在のモデル、および宇宙船への現在の制御入力に基づいて、宇宙船の次の状態を求める。そのため、ＭＰＣの最適化変数の次元が数十から数百と膨大な数となって、処理が複雑化するため、宇宙機に搭載される計算機での処理に適さないという課題がある。

本発明は上述の事情に鑑みてなされたものであり、推薬の使用量が少なく、軌道保持制御と角運動量のアンローディングとを行うための処理が、宇宙機に搭載される計算機で処理が可能となるように簡易化された宇宙機制御装置、宇宙機制御方法、およびプログラムを提供することを目的とする。

上記目的を達成するため、本発明に係る宇宙機制御装置は、角運動量算出部、指令値算出部、スラスタ制御部、およびジンバル制御部を備える。角運動量算出部は、天体の周りの軌道を周回する宇宙機が有する角運動量を算出する。指令値算出部は、それぞれが宇宙機の構体に、ジンバル機構を介して取り付けられた複数のスラスタの内、少なくともいずれかのスラスタの噴射を指示するスラスタ指令値、およびスラスタの噴射時におけるジンバル機構の角度を指示する角指令値を算出する。スラスタ制御部は、スラスタ指令値に基づき、スラスタを制御する。ジンバル制御部は、角指令値に基づき、ジンバル機構を制御する。指令値算出部は、第１の制約条件の下で、複数の噴射区間におけるスラスタの噴射量の合計が目的関数である非線形計画問題の解を求めることで、複数の噴射区間のそれぞれのスラスタ指令値および角指令値を算出する。非線形計画問題の第１の制約条件は、第１の関数および第２の関数に基づく。第１の関数は、軌道における複数の噴射区間でのスラスタ指令値および角指令値を引数として、該スラスタ指令値に基づいてスラスタが制御され、該角指令値に基づいてジンバル機構が制御された場合に、宇宙機が軌道を一周する間の、宇宙機の軌道を表す軌道要素の変化量と角運動量の変化量とを出力する。第２の関数は、軌道要素、軌道要素の時間変化率、および角運動量を引数として、宇宙機が軌道を一周する間の、軌道要素の制御量と角運動量の制御量とを出力する。第１の制約条件は、軌道要素の変化量と軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致することである。

本発明によれば、軌道要素の変化量と軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致する第１の制約条件の下で、複数の噴射区間におけるスラスタの噴射量の合計が目的関数である非線形計画問題の解が求められる。非線形計画問題の解を求めて、スラスタ指令値および角指令値を算出するため、推薬の使用量が少なく、軌道保持制御および角運動量のアンローディングを行うための処理が簡易化される。

本発明の実施の形態１における静止衛星の構成の例を示す図 実施の形態１におけるスラスタの宇宙機への取付例を示す図 実施の形態１におけるスラスタおよびジンバル機構の例を示す図 実施の形態１に係る宇宙機制御装置の構成例を示すブロック図 実施の形態１に係る指令値算出部の構成例を示すブロック図 実施の形態１に係る指令値算出部が行う指令値算出の動作の一例を示すフローチャート 実施の形態１に係る第２関数算出部が行う第２の関数の算出の動作の一例を示すフローチャート 本発明の実施の形態２に係る指令値算出部の構成例を示す図 実施の形態２に係る指令値算出部が行う指令値算出の動作の一例を示すフローチャート 実施の形態に係る宇宙機制御装置のハードウェア構成の一例を示す構成図 実施の形態に係る宇宙機制御装置のハードウェア構成の他の一例を示す構成図

以下、本発明の実施の形態に係る宇宙機制御装置について図面を参照して詳細に説明する。なお図中、同一または同等の部分には同一の符号を付す。

（実施の形態１）
実施の形態１に係る宇宙機制御装置を、宇宙機の一例である地球の周りを周回する静止衛星に搭載される宇宙機制御装置を例に説明する。この宇宙機制御装置は、静止衛星を目標とする軌道上に維持する軌道保持制御と、静止衛星の姿勢を制御する姿勢制御アクチュエータに蓄積された角運動量のアンローディングとを行う。

上述したように実施の形態１に係る宇宙機制御装置の制御対象である静止衛星を図１に示す。図１に示す静止衛星１０は、構体５と、構体５に取り付けられるジンバル機構３１，３２，３３，３４と、ジンバル機構３１，３２，３３，３４のそれぞれに取り付けられるスラスタ２１，２２，２３，２４と、を備える。ジンバル機構３１，３２，３３，３４は、いずれも同じ構造であり、構体５の外面の内、地球に向いている面５０の反対側に位置する面５１に取り付けられる。スラスタ２１，２２，２３，２４は、噴射方向が互いに異なる向きで面５１に取り付けられる。スラスタ２１，２２，２３，２４をそれぞれ、ＮＷスラスタ２１、ＮＥスラスタ２２、ＳＷスラスタ２３、およびＳＥスラスタ２４とも呼ぶこととする。スラスタ２１，２２，２３，２４およびジンバル機構３１，３２，３３，３４は、後述する宇宙機制御装置によって制御される。詳細には、宇宙機制御装置は、図１に示すジンバル機構３１，３２，３３，３４の角度を制御し、ジンバル機構３１，３２，３３，３４を介して構体５に取り付けられるスラスタ２１，２２，２３，２４を噴射させることで、静止衛星１０の軌道保持制御と角運動量のアンローディングとを行う。

図２に示すように、理解を容易にするため、ｘ_Ｈ軸、ｙ_Ｈ軸、ｚ_Ｈ軸を基底とし、静止衛星１０の質量中心Ｃ．Ｍ．（Center of Mass）を原点とする軌道座標系を設定し、適宜参照する。ｘ_Ｈ軸は、軌道の中心に位置する天体である地球４の中心と静止衛星１０の質量中心Ｃ．Ｍ．とを通る。ｚ_Ｈ軸は、静止衛星１０の軌道面の法線方向である。ｙ_Ｈ軸は、ｘ_Ｈ軸およびｚ_Ｈ軸と右手系の座標を構成する。ジンバル機構３１，３２，３３，３４のそれぞれから質量中心Ｃ．Ｍ．までの距離を示す取付距離ｒは互いに同じである。ｘ_Ｈｙ_Ｈ平面に投影された質量中心Ｃ．Ｍ．からスラスタ２１，２２，２３，２４のそれぞれを指すベクトルと、ｘ_Ｈ軸とがなす鋭角である角度ψは互いに同じである。ｘ_Ｈｚ_Ｈ平面に投影された質量中心Ｃ．Ｍ．からスラスタ２１，２２，２３，２４のそれぞれを指すベクトルと、ｚ_Ｈ軸とがなす鋭角である角度γは互いに同じである。

ジンバル機構３１，３２，３３，３４の構造はいずれも同じであるから、ジンバル機構３１について図３を用いて説明する。ジンバル機構３１の回転軸を、互いに直交するｘ_Ｔ軸およびｙ_Ｔ軸とする。ｘＴ軸およびｙＴ軸に直交する軸をｚ_Ｔ軸とする。ｚ_Ｔ軸は、ジンバル機構３１に取り付けられるスラスタ２１の推力軸に平行に位置する。ｘ_Ｔ軸周りのジンバル角をαとし、ｙ_Ｔ軸周りのジンバル角をβとする。ジンバル角α，βが０である場合、ｚ_Ｔ軸は静止衛星１０の質量中心Ｃ．Ｍ．を貫き、スラスタ２１を噴射しても静止衛星１０にトルクが働かない。ジンバル角α，βの少なくとも１つが０でない場合は、スラスタ２１の噴射時に静止衛星１０にトルクが働く。

上述のスラスタ２１，２２，２３，２４およびジンバル機構３１，３２，３３，３４は、図４に示す宇宙機制御装置１によって制御される。宇宙機制御装置１は、ジンバル機構３１，３２，３３，３４の角度を制御し、スラスタ２１，２２，２３，２４の内、少なくとも１つを噴射させることで、静止衛星１０の軌道保持制御と角運動量のアンローディングとを行う。
宇宙機制御装置１は、静止衛星１０の軌道要素を算出する軌道決定部１１と、静止衛星１０の角運動量を算出する角運動量算出部１２と、スラスタ２１，２２，２３，２４に対するスラスタ指令値およびジンバル機構３１，３２，３３，３４に対する角指令値を算出する指令値算出部１３と、スラスタ指令値に基づきスラスタ２１，２２，２３，２４を制御するスラスタ制御部１４と、角指令値に基づきジンバル機構３１，３２，３３，３４を制御するジンバル制御部１５とを備える。

軌道決定部１１は、例えば、静止衛星１０に搭載されたＧＰＳ（Global Positioning System：全地球測位システム）受信機から取得したＧＰＳ情報、地上局から受信したレンジ、レンジレート、静止衛星１０の方位角および仰角の情報等から、静止衛星１０の軌道要素を算出する。なおレンジは、地上局から静止衛星１０までの距離を示す。レンジレートは、地上局から静止衛星１０までの距離の変化の割合を示す。地上局は、地上の光学カメラを用いた観測から静止衛星１０の方位角および仰角を得ることができる。軌道決定部１１は、静止衛星１０の軌道要素の瞬時値から周期的な変動成分を取り除いて平均軌道要素を算出する。軌道要素は、例えば、離心率、傾斜角等を含む。以下の説明において、静止衛星１０の軌道要素として平均軌道要素を用いる。周期的な変動成分とは、地球重力の摂動、例えば太陽、月等の摂動天体の潮汐力等によって発生し、周期が閾値以下の周期的な変動成分である。軌道決定部１１は、算出した平均軌道要素を指令値算出部１３に送る。軌道決定部１１はさらに、平均軌道要素の時間変化率を算出し、平均軌道要素の時間変化率を指令値算出部１３に送る。平均軌道要素の時間変化率は、十分に短い時間における平均軌道要素の変化量を示す。

角運動量算出部１２は、静止衛星１０に搭載された恒星センサ、ジャイロセンサ、磁気センサ、地球センサ、太陽センサ等のセンサから得られる静止衛星１０の姿勢情報を用いて、静止衛星１０の姿勢および姿勢角速度を求め、静止衛星１０の角運動量を算出する。角運動量算出部１２は、算出した角運動量を指令値算出部１３に送る。

指令値算出部１３は、詳細については後述するが、非線形計画問題の解を求めることで、静止衛星１０を目標とする軌道から定められた範囲内に保持するための、スラスタ指令値および角指令値を算出する。なおスラスタ指令値は、スラスタ２１，２２，２３，２４の内、少なくともいずれかのスラスタに対する指令値である。詳細にはスラスタ指令値は、例えば、噴射タイミング、噴射位相、スラスタ２１，２２，２３，２４の噴射量の合計、噴射量の合計に対する噴射量の割合等を含む。角指令値は、スラスタ２１，２２，２３，２４のそれぞれの噴射時におけるジンバル機構３１，３２，３３，３４の少なくともいずれかのジンバル角α，βを示す。

スラスタ制御部１４は、指令値算出部１３が算出したスラスタ指令値に基づき、スラスタ２１，２２，２３，２４を制御する。ジンバル制御部１５は、指令値算出部１３が算出した角指令値に基づきジンバル機構３１，３２，３３，３４を制御する。

指令値算出部１３は、図５に示すように、後述する第１の関数を求める第１関数算出部１３１と、後述する第２の関数を求める第２関数算出部１３２と、第１の関数の出力と第２の関数の出力とが一致する第１の制約条件を決める第１制約条件設定部１３３と、ジンバル機構３１，３２，３３，３４の角度が定められた範囲内にある第２の制約条件を決める第２制約条件設定部１３４と、第１の制約条件と第２の制約条件の下で、複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計が目的関数である非線形計画問題の解を求めることで、スラスタ指令値および角指令値を算出する求解部１３５と、を備える。なお指令値算出部１３は、静止衛星１０の緯度および経度の許容範囲、角運動量の許容範囲、ならびにジンバル角α，βの許容範囲を予め保持している。

第１関数算出部１３１は、複数の噴射区間でのスラスタ指令値および角指令値を引数として、後述する下記（１３）−（１６）式を含む第１の関数を求める。第１の関数は、引数であるスラスタ指令値に基づいてスラスタ２１，２２，２３，２４が制御され、引数である角指令値に基づいてジンバル機構３１，３２，３３，３４が制御された場合に、静止衛星１０が軌道を一周する間の、軌道要素の変化量および角運動量の変化量を出力する。

第２関数算出部１３２は、軌道決定部１１から平均軌道要素および平均軌道要素の時間変化率を取得し、角運動量算出部１２から角運動量を取得する。そして、第２関数算出部１３２は、平均軌道要素、平均軌道要素の時間変化率、および静止衛星１０が有する角運動量を引数として、後述する下記（１７），（２０），（２２），（２４）式を含む第２の関数を求める。第２の関数は、静止衛星１０の平均軌道要素の制御量および静止衛星１０が有する角運動量の制御量を出力する。

第１制約条件設定部１３３は、第１の関数を第１関数算出部１３１から取得し、第２の関数を第２関数算出部１３２から取得する。そして、第１制約条件設定部１３３は、軌道要素の変化量と平均軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致することを示す条件式である第１の制約条件を求める。

第２制約条件設定部１３４は、ジンバル機構３１，３２，３３，３４の角度が定められた範囲内にあることを示す条件式である第２の制約条件を求める。

求解部１３５は、第１の制約条件と第２の制約条件の下で、複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計が目的関数である非線形計画問題について、目的関数を最小化する解を求める。その結果、複数の噴射区間のそれぞれのスラスタ指令値および角指令値が算出される。

指令値算出部１３が予め保持している許容範囲の詳細について説明する。指令値算出部１３は、静止衛星１０の目標位置からのずれの許容範囲を予め保持している。詳細には、静止衛星１０の緯度の目標値φ_ｒｅｆおよび保持精度Δφ_ｍａｘ、ならびに、静止衛星１０の経度の目標値λ_ｒｅｆおよび保持精度Δλ_ｍａｘを予め保持している。緯度の目標値φ_ｒｅｆおよび経度の目標値λ_ｒｅｆは、上述の目標とする軌道を示す。緯度の保持精度Δφ_ｍａｘおよび経度の保持精度Δλ_ｍａｘは、静止衛星１０が保持される、目標とする軌道を含む定められた範囲を示す。
静止衛星１０の緯度φの範囲は下記（１）式で表される。静止衛星１０の経度λの範囲は、下記（２）式で表される。
φ_ｒｅｆ−Δφ_ｍａｘ≦φ≦φ_ｒｅｆ＋Δφ_ｍａｘ （１）
λ_ｒｅｆ−Δλ_ｍａｘ≦λ≦λ_ｒｅｆ＋Δλ_ｍａｘ （２）

静止衛星１０の場合、φ_ｒｅｆは０であり、静止衛星１０の観測対象が日本であれば、λ_ｒｅｆは例えば１４０°である。緯度の保持精度Δφ_ｍａｘおよび経度の保持精度Δλ_ｍａｘを、例えば０．１°とする。この場合、静止衛星１０の緯度の範囲は−０．１°以上、かつ、０．１°以下であり、静止衛星１０の経度の範囲は、１３９．９°以上、かつ、１４０．１°以下である。

さらに、指令値算出部１３は、静止衛星１０の角運動量の目標値からのずれの許容範囲を予め保持している。詳細には、指令値算出部１３は、慣性系における静止衛星１０の角運動量の目標値ｈ_ｘｒｅｆ，ｈ_ｙｒｅｆ，ｈ_ｚｒｅｆ、および角運動量の保持精度Δｈ_ｘｍａｘ，Δｈ_ｙｍａｘ，Δｈ_ｚｍａｘを予め保持している。

静止衛星１０の角運動量の範囲は、下記（３）−（５）式で表される。
ｈ_ｘｒｅｆ−Δｈ_ｘｍａｘ≦ｈ_ｘ≦ｈ_ｘｒｅｆ＋Δｈ_ｘｍａｘ （３）
ｈ_ｙｒｅｆ−Δｈ_ｙｍａｘ≦ｈ_ｙ≦ｈ_ｙｒｅｆ＋Δｈ_ｙｍａｘ （４）
ｈ_ｚｒｅｆ−Δｈ_ｚｍａｘ≦ｈ_ｚ≦ｈ_ｚｒｅｆ＋Δｈ_ｚｍａｘ （５）
例えば、ゼロモーメンタム衛星の場合は、ｈ_ｘｒｅｆ、ｈ_ｙｒｅｆ、および、ｈ_ｚｒｅｆはいずれも０である。バイアスモーメンタム衛星の場合は、ｈ_ｘｒｅｆ、ｈ_ｙｒｅｆ、ｈ_ｚｒｅｆはいずれも０でない。Δｈ_ｘｍａｘ、Δｈ_ｙｍａｘ、および、Δｈ_ｚｍａｘは、リアクションホイール、ＣＭＧ（Control Moment Gyro）等の静止衛星１０に搭載された姿勢制御アクチュエータで補償可能な角運動量の大きさに応じて定められる。Δｈ_ｘｍａｘ、Δｈ_ｙｍａｘ、および、Δｈ_ｚｍａｘを、例えば３０Ｎｍｓとする。

さらに、指令値算出部１３は、ジンバル機構の駆動角度の許容範囲を予め保持している。詳細には、指令値算出部１３は、ジンバル機構３１，３２，３３，３４の駆動角度の上限α_ｍａｘ，β_ｍａｘおよび下限−α_ｍａｘ，−β_ｍａｘを予め保持している。

ジンバル機構３１，３２，３３，３４の駆動角度の上限α_ｍａｘ、β_ｍａｘおよび下限−α_ｍａｘ，−β_ｍａｘは、静止衛星１０の姿勢情報を取得するための上述のセンサ、姿勢制御アクチュエータ、姿勢制御アクチュエータを制御する制御装置等を含む姿勢制御系が補償可能な外乱トルクの大きさに応じて定められる。ジンバル機構３１，３２，３３，３４を回転させて、スラスタ２１，２２，２３，２４を噴射することで、静止衛星１０にトルクが働く。トルクの大きさＴは、下記（６）式で表されるように、ジンバル角の二乗ノルム、スラスタ２１，２２，２３，２４の推力の大きさＦ、および取付距離ｒに比例する。

静止衛星１０の打ち上げ後は、取付距離ｒを変更することはできない。スラスタ２１，２２，２３，２４の推力の大きさＦは、軌道上で一定とみなされる。そのため、トルクの大きさＴの調節は、ジンバル角α，βによって行われる。ジンバル角α，βを大きくすることで、より大きなトルクを静止衛星１０に与えることが可能である。ジンバル角α，βを大きくすることで、静止衛星１０に蓄積された角運動量をアンローディングする能力は向上する。しかしながら、スラスタ２１，２２，２３，２４の噴射によって静止衛星１０に作用するトルクは、静止衛星１０の姿勢制御系にとっては、外乱トルクである。そのため、トルクの大きさＴが大きすぎると、姿勢制御系が不安定になり、静止衛星１０の姿勢を保持できなくなる。そのため、ジンバル角α，βは、姿勢制御系で補償可能な外乱トルクの上限値Ｔ_ｍａｘに応じて定められる。ジンバル角αの最大値α_ｍａｘは、下記（７）式で表される。ジンバル角βの最大値β_ｍａｘは、下記（８）式で表される。下記（７），（８）式におけるκは、１未満の正数である。

ジンバル角α，βの範囲は、下記（９），（１０）式で表される。
−α_ｍａｘ≦α≦α_ｍａｘ （９）
−β_ｍａｘ≦β≦β_ｍａｘ （１０）

指令値算出部１３は、上述の静止衛星１０の緯度および経度の範囲、角運動量の範囲、ならびにジンバル角α，βの範囲を予め保持している。指令値算出部１３は、これらの許容範囲に基づいて、非線形計画問題の解を求めることで、スラスタ指令値および角指令値を算出する。詳細には、指令値算出部１３は、後述する第１の制約条件の下で、目的関数を最小化する解を求める。図６を用いて、指令値算出部１３の全体の流れを説明した後、個々のステップについて順に詳細に説明する。第１関数算出部１３１は、複数の噴射区間でのスラスタ指令値および角指令値を引数として、第１の関数を求める（ステップＳ１１）。第２関数算出部１３２は、軌道決定部１１から平均軌道要素を取得し、角運動量算出部１２から角運動量を取得する（ステップＳ１２）。第２関数算出部１３２は、平均軌道要素、平均軌道要素の時間変化率、および静止衛星１０が有する角運動量を引数として、第２の関数を求める（ステップＳ１３）。第１制約条件設定部１３３は、軌道要素の変化量と平均軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致することを示す条件式である第１の制約条件を求める（ステップＳ１４）。第２制約条件設定部１３４は、ジンバル機構３１，３２，３３，３４の角度が定められた範囲内にあることを示す条件式である第２の制約条件を求める（ステップＳ１５）。求解部１３５は、目的関数を求める（ステップＳ１６）。求解部１３５は、第１の制約条件、第２の制約条件、および複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計である目的関数に基づく非線形計画問題について、目的関数を最小化する解を求めることで、スラスタ指令値および角指令値を算出する（ステップＳ１７）。指令値算出部１３は、スラスタ指令値および角指令値をそれぞれ、スラスタ制御部１４およびジンバル制御部１５に出力する（ステップＳ１８）。

図６の各ステップの詳細について順に説明する。上記ステップＳ１１の処理の詳細について説明する。第１関数算出部１３１は、複数の噴射区間でのスラスタ指令値および角指令値を引数として、第１の関数を求める。軌道の一周あたりＮ個の噴射区間が設けられる場合を例にして説明する。スラスタ指令値が、Ｎ個の噴射区間のそれぞれにおいて、スラスタ２１，２２，２３，２４の内、２本に対する噴射の指示である場合を例にして説明する。２本の組み合わせは、ＮＷスラスタ２１とＮＥスラスタ２２の組、ＳＷスラスタ２３とＳＥスラスタ２４の組、ＳＷスラスタ２３とＮＷスラスタ２１の組、および、ＮＥスラスタ２２とＳＥスラスタ２４の組のいずれかである。１つの噴射区間でスラスタ２１，２２，２３，２４の内、２本のみ噴射することで、推薬の使用量を少なくすることが可能である。

ｊ（１≦ｊ≦Ｎ）番目の噴射区間において噴射するスラスタ２１，２２，２３，２４の内の２本を、ａ_ｊ，ｂ_ｊとする。ａ_ｊとｂ_ｊは、ＮＷスラスタ２１とＮＥスラスタ２２の組、ＳＷスラスタ２３とＳＥスラスタ２４の組、ＳＷスラスタ２３とＮＷスラスタ２１の組、および、ＮＥスラスタ２２とＳＥスラスタ２４の組のいずれかである。ｊ番目の噴射区間において噴射するスラスタａ_ｊ，ｂ_ｊの噴射量の合計をＦ_ｊとすると、ａ_ｊ，ｂ_ｊのそれぞれの噴射量を決定する配分係数ξ_ｊを用いて、スラスタａ_ｊ，ｂ_ｊのそれぞれの噴射量ｆ_ａｊ，ｆ_ｂｊは、下記（１１），（１２）式で表される。配分係数ξ_ｊは、−１≦ξ_ｊ≦１の条件を満たす係数である。ξ_ｊ＝１の場合には、スラスタａ_ｊのみが噴射し、ξ_ｊ＝−１の場合は、スラスタｂ_ｊのみが噴射し、ξ_ｊ＝０の場合は、スラスタａ_ｊ，ｂ_ｊが等しく噴射する。

第１の関数のパラメータは、ｊ番目の噴射区間における噴射量の合計Ｆ_ｊ，ｊ番目の噴射区間における配分係数ξ_ｊ，ｊ番目の噴射区間における平均経度λ_ａｖｇｊ、およびｊ番目の噴射区間における角度α_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊである。角度α_ａｊ，β_ａｊは、スラスタａ_ｊを構体５に取り付けるジンバル機構３１，３２，３３，３４のｊ番目の噴射区間における角度である。角度α_ｂｊ，β_ｂｊは、スラスタｂ_ｊを構体５に取り付けるジンバル機構３１，３２，３３，３４のｊ番目の噴射区間における角度である。

第１の関数の出力である軌道要素の変化量は、線形化されたガウスの惑星方程式に基づいて、噴射量の合計Ｆ_ｊ、配分係数ξ_ｊ、平均経度λ_ａｖｇｊ、および角度α_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊで表される。軌道要素の変化量は、平均離心率ベクトル、平均傾斜角ベクトル、および平均直下点経度のそれぞれの変化量を含む。さらに第１の関数の出力である角運動量の変化量は、線形化されたオイラーの方程式に基づいて、噴射量の合計Ｆ_ｊ、配分係数ξ_ｊ、平均経度λ_ａｖｇｊ、および角度α_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊで表される。

Ｎ個の噴射区間におけるスラスタ２１，２２，２３，２４の噴射によって生じる軌道要素の変化量の合計および角運動量の変化量の合計は、下記（１３）−（１６）式で表される。下記（１３）式において、Ｇ_ｅ（Ｆ_ｊ，ξ_ｊ，λ_ａｖｇｊ，α_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊ）は、離心率ベクトルについてのガウスの惑星方程式である。下記（１４）式において、Ｇ_ｉ（Ｆ_ｊ，ξ_ｊ，λ_ａｖｇｊ，α_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊ）は、傾斜角ベクトルについてのガウスの惑星方程式である。下記（１５）式において、Ｇ_λ（Ｆ_ｊ，ξ_ｊ，λ_ａｖｇｊ，α_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊ）は、平均直下点経度についてのガウスの惑星方程式である。下記（１６）式において、Ｅ_ｈ（Ｆ_ｊ，ξ_ｊ，λ_ａｖｇｊ，α_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊ）は、静止衛星１０が有する角運動量に関するオイラーの方程式である。第１関数算出部１３１は、上記パラメータを確定して、下記（１３）−（１６）式を含む第１の関数を求め、第１制約条件設定部１３３に送る。

図６のステップＳ１１の処理の説明は以上であるから、図６のステップＳ１１と並行して行われるステップＳ１２，Ｓ１３の処理について以下に説明する。ステップＳ１２において、第２関数算出部１３２は、軌道決定部１１から平均軌道要素および平均軌道要素の時間変化率を取得し、角運動量算出部１２から角運動量を取得する。そして、ステップＳ１３において、第２関数算出部１３２は、平均軌道要素、平均軌道要素の時間変化率、および静止衛星１０が有する角運動量を引数として、パラメータを確定し、静止衛星１０の平均軌道要素の制御量および静止衛星１０が有する角運動量の制御量を出力する第２の関数を求める。平均軌道要素の制御量とは、スラスタ２１，２２，２３，２４を噴射することで達成すべき平均軌道要素の制御量である。制御対象となる平均軌道要素は、平均離心率ベクトル、平均傾斜角ベクトル、および平均直下点経度を含む。

第２関数算出部１３２が第２の関数を算出するステップＳ１３の処理の概要について図７を用いて説明した後に、各ステップの詳細について説明する。第２関数算出部１３２は、平均離心率ベクトルの目標値ｅ_{ｔａｒｇｅｔ}を定める軌道座標系の原点を中心とする円の半径ｅ_Ｒを決定する（ステップＳ２１）。第２関数算出部１３２は、軌道座標系の原点を中心とする半径ｅ_Ｒの円と、地球４と太陽とを結ぶ線分の交点に設定される平均離心率ベクトルの目標値ｅ_{ｔａｒｇｅｔ}に平均離心率ベクトルを近づけるための制御量Δｅ_ｃを算出する（ステップＳ２２）。第２関数算出部１３２は、緯度の保持精度Δφ_ｍａｘの大きさに応じて、年平均傾斜角ベクトル、月平均傾斜角ベクトル、および日平均傾斜角ベクトルのいずれかを、制御対象である平均傾斜角ベクトルとして決定する（ステップＳ２３）。第２関数算出部１３２は、ステップＳ２３で決定した平均傾斜角ベクトルに対する制御量Δｉ_ｃを算出する（ステップＳ２４）。第２関数算出部１３２は、平均直下点経度の制御量Δλ_ｃを算出する（ステップＳ２５）。第２関数算出部１３２は、角運動量の制御量Δｈ_ｃを算出する（ステップＳ２６）。図７に示すように、第２関数算出部１３２は、ステップＳ２１，Ｓ２３，Ｓ２５，Ｓ２６の処理を並行に行うことができる。

上記ステップＳ２１，Ｓ２２の処理について説明する。第２関数算出部１３２は、軌道座標系の原点を中心とする半径ｅ_Ｒの円と、地球４と太陽とを結ぶ線分の交点に平均離心率ベクトルの目標値ｅ_{ｔａｒｇｅｔ}を設定し、下記（１７）式で表される平均離心率ベクトルの制御量Δｅ_ｃを算出する。下記（１７）式において、ｅ_ＯＤは、軌道決定部１１が算出した平均離心率ベクトルである。Ｋ_ｅは正の係数である。下記（１７）式の第２項−Ｋ_ｅ（ｅ_ＯＤ−ｅ_{ｔａｒｇｅｔ}）は、フィードバック制御量に相当する。下記（１７）式における第１項−δｅ_ＦＦは、フィードフォワード制御量に相当する。フィードフォワード制御量は、静止衛星１０に働く摂動力、例えば、太陽輻射圧、月の潮汐力、太陽の潮汐力等による、平均離心率ベクトルの変動を補償する。下記（１７）式におえるδｅ_ＦＦは、下記（１８）式で表される。下記（１８）式におけるｅ’_ＯＤは、ｅ_ＯＤの時間変化率であり、Ｔ_ｃは、静止衛星１０が軌道を周回する周期である。
Δｅ_ｃ＝−δｅ_ＦＦ−Ｋ_ｅ（ｅ_ＯＤ−ｅ_{ｔａｒｇｅｔ}） （１７）
δｅ_ＦＦ＝ｅ’_ＯＤＴ_ｃ （１８）

上記（１７）式における平均離心率ベクトルの目標値ｅ_{ｔａｒｇｅｔ}は、軌道座標系の原点を中心とする半径ｅ_Ｒの円と、地球４と太陽とを結ぶ線分の交点に設定される。ｅ_Ｒは、下記（１９）式で表される。Δλ_Ｍは、軌道決定部１１が算出する平均軌道要素の誤差、スラスタ２１，２２，２３，２４の噴射誤差等によって生じる経度の誤差である。Δλ_Ｍは、０＜Δλ_Ｍ＜Δλ_ｍａｘを満たす値であり、例えば、０．０２°である。

図７のステップＳ２１，Ｓ２２の処理については以上であるから、図７においてステップＳ２１，Ｓ２２と並行して行われるステップＳ２３，Ｓ２４の処理について説明する。第２関数算出部１３２は、平均傾斜角ベクトル、例えば、年平均傾斜角ベクトル、月平均傾斜角ベクトル、および日平均傾斜角ベクトル等を制御する。第２関数算出部１３２は、緯度の保持精度Δφ_ｍａｘの大きさに適した平均傾斜角ベクトルを制御する。例えば、０．０４０°×（１＋κ_φ）≦Δφ_ｍａｘの場合、第２関数算出部１３２は、年平均傾斜角ベクトルを制御対象として決定する。０．００４０°×（１＋κ_φ）＜Δφ_ｍａｘ＜０．０４０°×（１＋κ_φ）の場合、第２関数算出部１３２は、月平均傾斜角ベクトルを制御対象として決定する。Δφ_ｍａｘ＜０．００４°×（１＋κ_φ）の場合、第２関数算出部１３２は、日平均傾斜角ベクトルを制御対象として決定する。κ_φは、正の係数であり、例えば、１．１に設定される。

平均傾斜角ベクトルの制御量Δｉ_ｃは、下記（２０）式で表される。下記（２０）式において、ｉ_ＯＤは、軌道決定部１１が算出した平均傾斜角ベクトルである。ｉ_{ｔａｒｇｅｔ}は、定められた平均傾斜角ベクトルの目標値である。Ｋ_ｉは、正の定数である。下記（２０）式における第２項−Ｋ_ｉ（ｉ_ＯＤ−ｉ_{ｔａｒｇｅｔ}）は、フィードバック制御量に相当する。下記（２０）式における第１項−δｉ_ＦＦは、フィードフォワード制御量に相当する。フィードフォワード制御量は、静止衛星１０に働く摂動力による、平均傾斜角ベクトルの変動を補償する。δｉ_ＦＦは、下記（２１）式で表される。下記（２１）式において、ｄｉ_ＦＦ／ｄｔは、平均傾斜角ベクトルの時間変化率である。
Δｉ_ｃ＝−δｉ_ＦＦ−Ｋ_ｉ（ｉ_ＯＤ−ｉ_{ｔａｒｇｅｔ}） （２０）

図７のステップＳ２３，Ｓ２４の説明は以上であるから、図７においてステップＳ２３，Ｓ２４と並行して行われるステップＳ２５の処理について説明する。第２関数算出部１３２は、下記（２２）式で表される平均直下点経度の制御量Δλ_ｃを算出する。下記（２２）式において、λ_ＯＤは、軌道決定部１１で算出された平均直下点経度である。下記（２２）式において、λ_{ｔａｒｇｅｔ}は、平均直下点経度の目標値である。Ｋ_λは正の係数である。下記（２２）式における第２項−Ｋ_λ（λ_ＯＤ−λ_{ｔａｒｇｅｔ}）は、フィードバック制御量に相当する。下記（２２）式における第１項−δλ_ＥＦは、フィードフォワード制御量に相当する。フィードフォワード制御量は、静止衛星１０に働く摂動力による、平均直下点経度の変動を補償する。δλ_ＥＦは、下記（２３）式で表される。下記（２３）式において、λ’_ＯＤは、平均直下点経度の角速度である。下記（２３）式においてλ”_ＯＤは、平均直下点経度の角加速度である。
Δλ_ｃ＝−δλ_ＦＦ−Ｋ_λ（λ_ＯＤ−λ_{ｔａｒｇｅｔ}） （２２）

図７のステップＳ２５の説明は以上であるから、図７においてステップＳ２５と並行して行われるステップＳ２６の処理について説明する。第２関数算出部１３２は、下記（２４）式で表される角運動量の制御量Δｈ_ｃを算出する。下記（２４）式において、ｈ_ＡＤは、角運動量算出部１２で算出された慣性系における静止衛星１０の角運動量である。下記（２４）式の第２項−Ｋ_ｈ（ｈ_ＡＤ−ｈ_ｒｅｆ）は、フィードバック制御量に相当する。Ｋ_ｈは、正の係数である。下記（２４）式の第１項−δｈ_ＦＦは、フィードフォワード制御量に相当する。フィードフォワード制御量は、摂動トルクによって生じる慣性系での角運動量の変動量を補償する。摂動トルクの内、最も支配的なトルクは、太陽輻射圧によるトルクである。静止衛星１０に働く太陽輻射圧トルクをτ_ＳＲＰとし、δｈ_ＦＦを下記（２５）式で近似する。ｈ_ｒｅｆは、慣性系における角運動量の目標値であり、下記（２６）式で表される。角運動量の目標値ｈ_ｒｅｆの各成分は、予め指令値算出部１３が保持している。
Δｈ_ｃ＝−δｈ_ＦＦ−Ｋ_ｈ（ｈ_ＡＤ−ｈ_ｒｅｆ） （２４）
δｈ_ＦＦ＝Ｔ_ｃτ_ＳＲＰ （２５）
ｈ_ｒｅｆ＝［ｈ_ｘｒｅｆ，ｈ_ｙｒｅｆ，ｈ_ｚｒｅｆ］^Ｔ （２６）

第２関数算出部１３２は、上記（１７），（２０），（２２），（２４）式を含む第２の関数を第１制約条件設定部１３３に送る。

図６のステップＳ１２，１３の処理の説明は以上であるから、後続のステップＳ１４の処理について説明する。第１制約条件設定部１３３は、軌道要素の変化量と平均軌道要素の制御量とが一致し、かつ、角運動量の変化量と角運動量の制御量とが一致することを示す条件式である第１の制約条件を求める。第１制約条件設定部１３３は、下記（２７）−（３０）式を含む第１の制約条件を求解部１３５に送る。
Δｅ_ｃ＝δ_ｅ （２７）
Δｉ_ｃ＝δ_ｉ （２８）
Δλ_ｃ＝δ_λ （２９）
Δｈ_ｃ＝δ_ｈ （３０）

図６のステップＳ１４の処理の説明は以上であるから、後続のステップＳ１５の処理について説明する。第２制約条件設定部１３４は、ジンバル機構３１，３２，３３，３４の角度が定められた範囲内にあることを示す条件式である第２の制約条件を求める。第２制約条件設定部１３４は、下記（３１）−（３４）式を含む第２の制約条件を求解部１３５に送る。
−α_ｍａｘ≦α_ａｊ≦α_ｍａｘ （３１）
−β_ｍａｘ≦β_ａｊ≦β_ｍａｘ （３２）
−α_ｍａｘ≦α_ｂｊ≦α_ｍａｘ （３３）
−β_ｍａｘ≦β_ｂｊ≦β_ｍａｘ （３４）

スラスタ２１，２２，２３，２４の噴射量は０より大きいことは明らかであるから、下記（３５）式が成り立つ。また上述のように、噴射区間ｊにおける配分係数については、下記（３６）式が成り立つ。
０≦Ｆ_ｊ （３５）
−１≦ξ_ｊ≦１ （３６）

図６のステップＳ１５の説明は以上であるから、後続のステップＳ１６，Ｓ１７の処理について説明する。求解部１３５は、第１の制約条件、第２の制約条件、および複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計である目的関数に基づく非線形計画問題について、目的関数を最小化する解を求めることで、複数の噴射区間のそれぞれのスラスタ指令値および角指令値を算出する。目的関数は、下記（３７）式で表される。

求解部１３５は、噴射区間のそれぞれにおける、噴射量の合計Ｆ_ｊ、配分係数ξ_ｊ、平均経度λ_ａｖｇｊ、および角度α_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊの解を求める。非線形計画問題を解くにあたっては、逐次二次計画法、内点法等の手法が用いられる。求解部１３５は、非線形計画問題を解くことで得られるλ_ａｖｇｊの解から、スラスタ２１，２２，２３，２４の噴射区間の開始位置および終了位置を含む噴射タイミングを得る。また求解部１３５は、非線形計画問題を解くことで得られるＦ_ｊ，ξ_ｊの解から、噴射区間におけるスラスタ２１，２２，２３，２４の噴射量を得る。さらに求解部１３５は、非線形計画問題を解くことで得られるα_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊの解から、スラスタ２１，２２，２３，２４を噴射する際のジンバル角を得る。

図６のステップＳ１６，Ｓ１７の説明は以上であるから、後続のステップＳ１８の処理について説明する。求解部１３５は、スラスタ２１，２２，２３，２４の噴射タイミングおよび噴射量を示すスラスタ指令値を、スラスタ制御部１４に送り、ジンバル角を示す角指令値を、ジンバル制御部１５に送る。

スラスタ制御部１４は、スラスタ指令値に基づいてスラスタ２１，２２，２３，２４を制御して噴射させる。ジンバル制御部１５は、角指令値に基づいてジンバル機構３１，３２，３３，３４のジンバル角を制御する。上述のように非線形計画問題を解いて、スラスタ指令値および角指令値を算出することで、推薬の使用量を少なくしながら、軌道保持制御および角運動量のアンローディングを行うための処理を簡易化することが可能である。

以上説明したとおり、本実施の形態１に係る宇宙機制御装置１によれば、非線形計画問題について、スラスタの噴射量の合計である目的関数を最小化する解を求め、スラスタ指令値および角指令値を算出することで、推薬の使用量を抑制しながら、軌道保持制御および角運動量のアンローディングを行うための処理を簡易化することが可能である。

実施の形態１に係る宇宙機制御装置１によれば、昇交点および降交点から離隔した位置でスラスタ２１，２２，２３，２４を噴射する必要がある場合、また軌道を一周する間にスラスタ２１，２２，２３，２４を複数回噴射する必要がある場合でも、推薬の使用量を抑制しながら、軌道保持制御および角運動量のアンローディングを行うための処理を行うことができる。

（実施の形態２）
実施の形態１における非線形計画問題の制約条件は、第１の制約条件と第２の制約条件であったが、非線形計画問題は、さらに他の制約条件を含んでもよい。実施の形態２に係る宇宙機制御装置１の構成は、図４に示す実施の形態１に係る宇宙機制御装置１の構成と同様である。図８は、本発明の実施の形態２に係る指令値算出部の構成例を示す図である。実施の形態２に係る指令値算出部１６は、図５に示す実施の形態１に係る指令値算出部１３の構成に加えて、第３制約条件設定部１３６をさらに備える。第３制約条件設定部１３６は、複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、角運動量が定められた範囲内にあることを示す条件式である第３の制約条件を求め、求解部１３５に送る。

指令値算出部１６がスラスタ指令値および角指令値を算出する処理の概略について図９を用いて説明する。ステップＳ１１−Ｓ１６，Ｓ１８の処理は、図６に示す処理と同様である。実施の形態２に係る指令値算出部１６においては、第３制約条件設定部１３６が、複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、角運動量が定められた範囲内にあることを示す条件式である第３の制約条件を求める（ステップＳ１９）。求解部１３５は、目的関数を求める（ステップＳ１６）。求解部１３５は、第１の制約条件、第２の制約条件、第３の制約条件、および複数の噴射区間におけるスラスタ２１，２２，２３，２４の噴射量の合計である目的関数に基づく非線形計画問題について、目的関数を最小化する解を求めることで、スラスタ指令値および角指令値を算出する（ステップＳ２０）。後続の処理は、図６に示す処理と同様である。

上記（３１）−（３４）式における、α_ｍａｘ，β_ｍａｘを小さくすると、スラスタ２１，２２，２３，２４の噴射中に静止衛星１０に作用するトルクは小さくなる。すなわち、スラスタ２１，２２，２３，２４の噴射によって生じる角運動量の変動が抑制される。α_ｍａｘ，β_ｍａｘを小さくすることで、角運動量の変動が抑制されても、噴射区間の噴射直前または噴射直後において、角運動量が、上記（３）−（５）式で規定される範囲内にないことがある。そこで、第３制約条件設定部１３６は、複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、角運動量が定められた範囲内にあることを示す条件式である第３の制約条件を求める。すなわち、第３制約条件設定部１３６は、下記（３８），（３９）式の少なくともいずれかが成立することを示す条件式である第３の制約条件を求める。下記（３８）式におけるｈ_ｊ ^（−）は、ｊ番目の噴射区間の噴射直前における角運動量である。ｊ番目の噴射区間の噴射直前における角運動量とは、ｊ番目の噴射区間におけるスラスタ２１，２２，２３，２４の噴射開始直前の静止衛星１０の角運動量である。下記（３９）式におけるｈ_ｊ ^（＋）は、ｊ番目の噴射区間の噴射直後における角運動量である。ｊ番目の噴射区間の噴射直後における角運動量とは、ｊ番目の噴射区間におけるスラスタ２１，２２，２３，２４の噴射終了直後の静止衛星１０の角運動量である。ｈ_ｊ ^（−）およびｈ_ｊ ^（＋）は、変数Ｆ_ｊ，ξ_ｊ，λ_ａｖｇｊ，α_ａｊ，β_ａｊ，α_ｂｊ，β_ｂｊの関数である。
ｈ_ｒｅｆ−Δｈ_ｍａｘ≦ｈ_ｊ ^（−）≦ｈ_ｒｅｆ＋ｈ_ｍａｘ （３８）
ｈ_ｒｅｆ−Δｈ_ｍａｘ≦ｈ_ｊ ^（＋）≦ｈ_ｒｅｆ＋ｈ_ｍａｘ （３９）

ｈ_ｊ ^（−）の各成分について、上記（３）−（５）式と同様に、下記（４０）−（４２）が成り立つ。またｈ_ｊ ^（＋）の各成分について、上記（３）−（５）式と同様に、下記（４３）−（４５）が成り立つ。
ｈ_ｘｒｅｆ−Δｈ_ｘｍａｘ≦ｈ_ｘｊ ^（−）≦ｈ_ｘｒｅｆ＋Δｈ_ｘｍａｘ （４０）
ｈ_ｙｒｅｆ−Δｈ_ｙｍａｘ≦ｈ_ｙｊ ^（−）≦ｈ_ｙｒｅｆ＋Δｈ_ｙｍａｘ （４１）
ｈ_ｚｒｅｆ−Δｈ_ｚｍａｘ≦ｈ_ｚｊ ^（−）≦ｈ_ｚｒｅｆ＋Δｈ_ｚｍａｘ （４２）
ｈ_ｘｒｅｆ−Δｈ_ｘｍａｘ≦ｈ_ｘｊ ^（＋）≦ｈ_ｘｒｅｆ＋Δｈ_ｘｍａｘ （４３）
ｈ_ｙｒｅｆ−Δｈ_ｙｍａｘ≦ｈ_ｙｊ ^（＋）≦ｈ_ｙｒｅｆ＋Δｈ_ｙｍａｘ （４４）
ｈ_ｚｒｅｆ−Δｈ_ｚｍａｘ≦ｈ_ｚｊ ^（＋）≦ｈ_ｚｒｅｆ＋Δｈ_ｚｍａｘ （４５）

以上説明したとおり、本実施の形態２に係る宇宙機制御装置１によれば、複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、角運動量が定められた範囲内にあることを、非線形計画問題の第３の制約条件とすることで、静止衛星１０の角運動量を保持範囲内に維持することができる。そのため、姿勢制御アクチュエータの角運動量を飽和させることなく、軌道保持制御および角運動量のアンローディングを行うことが可能である。

次に、本実施の形態のハードウェア構成について説明する。図１０は、ハードウェア構成の例を示す図である。図１０は、宇宙機制御装置１が処理回路１０１によって実現される例を示している。図１０のように、静止衛星１０は、衛星の位置、速度、姿勢角、角速度を観測するセンサ１１０、宇宙機制御装置１、スラスタ２１，２２，２３，２４、およびジンバル機構３１，３２，３３，３４を有する。

処理回路１０１が、専用のハードウェアである場合、処理回路１０１は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Application Specific Integrated Circuit）、ＦＰＧＡ（Field Programmable Gate Array）、またはこれらを組み合わせたものが該当する。軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の各部の機能それぞれを処理回路１０１で実現してもよいし、各部の機能をまとめて処理回路１０１で実現してもよい。

図１１は、本実施の形態のハードウェア構成の別の例を示す図である。図１１は、宇宙機制御装置１がプロセッサ（ＣＰＵ：Central Processing Unit）１０２とメモリ１０３によって実現される例を示している。図１１のプロセッサ１０２とメモリ１０３は、図１０の処理回路１０１に相当する。

宇宙機制御装置１がプロセッサ１０２とメモリ１０３によって実現される場合、軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の機能は、ソフトウェア、ファームウェア、またはソフトウェアとファームウェアとの組み合わせにより実現される。ソフトウェアやファームウェアはプログラムとして記述され、メモリ１０３に格納される。プロセッサ１０２は、メモリ１０３に記憶されたプログラムを読み出して実行することにより、各部の機能を実現する。すなわち、メモリ１０３には、軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の処理を実行するためのプログラムが格納される。メモリ１０３は、例えば、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read-Only Memory）、フラッシュメモリー、ＥＰＲＯＭ（Erasable Programmable Read Only Memory）、ＥＥＰＲＯＭ（Electrically Erasable and Programmable Read-Only Memory）等の、不揮発性または揮発性の半導体メモリ、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤ（Digital Versatile Disc）等を含む。

メモリ１０３は、上述のプログラムの他、軌道決定部１１で決定される軌道決定値、角運動量算出部１２で設定される衛星角運動量値、指令値算出部１３で算出されるスラスタ指令値およびジンバル角指令値を記憶する。上述のメモリ１０３に記憶された情報は適宜読みだされて、軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の処理に用いられ、修正される。

なお、軌道決定部１１、角運動量算出部１２、指令値算出部１３、スラスタ制御部１４、およびジンバル制御部１５の各機能について、一部を専用のハードウェアで実現し、一部をソフトウェアまたはファームウェアで実現するようにしてもよい。例えば、軌道決定部１１については専用のハードウェアとしての処理回路１０１でその機能を実現し、指令値算出部１３については、プロセッサ１０２がメモリ１０３に格納されたプログラムを読み出して実行することによってその機能を実現することが可能である。

本発明の実施の形態は上述の実施の形態に限られない。宇宙機は、静止衛星１０に限られず、天体の周りの軌道を周回する任意の飛行体である。飛行体の一例は、地球観測衛星、宇宙船等である。スラスタ２１，２２，２３，２４の本数は任意であり、それぞれの噴射区間において噴射するスラスタ２１，２２，２３，２４の本数および組み合わせも任意である。またジンバル機構３１，３２，３３，３４の自由度は１でもよい。ジンバル機構３１，３２，３３，３４の自由度が１である場合、上記数式のβ＝０となるため、第１の関数は、変数Ｆ_ｊ，ξ_ｊ，λ_ａｖｇｊ，α_ａｊ，α_ｂｊの関数である。β＝０の場合、第１の関数と同様に、ｈ_ｊ ^（−）およびｈ_ｊ ^（＋）は、変数Ｆ_ｊ，ξ_ｊ，λ_ａｖｇｊ，α_ａｊ，α_ｂｊの関数である。図１の例では、宇宙機の構体５に、スラスタ２１，２２，２３，２４が取り付けられているが、宇宙機の構体５に取り付けられるスラスタ２１，２２，２３，２４の数は任意である。スラスタ２１，２２，２３，２４のそれぞれの角度ψおよび角度γは互いに同じでもよいし、異なってもよい。

緯度および経度の範囲、ならびに、緯度および経度の保持精度は、上述の例に限られず、任意に定めることができる。上述のように、軌道決定部１１が平均軌道要素の時間変化率を算出する代わりに、第２関数算出部１３２が平均軌道要素の時間変化率を算出してもよい。軌道決定部１１は、軌道要素を算出する代わりに、外部装置から起動要素を取得してもよい。軌道決定部１１は、軌道要素の瞬時値を指令値算出部１３に送ってもよい。第２関数算出部１３２は、図７に示す、平均離心率ベクトルの制御量Δｅ_ｃの算出、平均傾斜角ベクトルΔｉ_ｃの算出、平均直下点経度の制御量Δλ_ｃの算出、および、角運動量の制御量Δｈ_ｃの算出を順次行ってもよい。

求解部１３５は、第１の制約条件の下で、非線形計画問題の解を求めることで、スラスタ指令値およびジンバル指令値を算出してもよい。あるいは、求解部１３５は、第１の制約条件と第３の制約条件の下で、非線形計画問題の解を求めてもよい。また、求解部１３５は、第１の制約条件の下で、目的関数を小さくする準最適解を求めることで、スラスト指令値およびジンバル指令値を算出してもよい。

本発明は、本発明の広義の精神と範囲を逸脱することなく、様々な実施の形態及び変形が可能とされるものである。また、上述した実施の形態は、この発明を説明するためのものであり、本発明の範囲を限定するものではない。すなわち、本発明の範囲は、実施の形態ではなく、特許請求の範囲によって示される。そして、特許請求の範囲内及びそれと同等の発明の意義の範囲内で施される様々な変形が、この発明の範囲内とみなされる。

本出願は、２０１７年９月１日に出願された、日本国特許出願特願２０１７−１６８６４９号に基づく。本明細書中に日本国特許出願特願２０１７−１６８６４９号の明細書、特許請求の範囲、図面全体を参照として取り込むものとする。

１ 宇宙機制御装置、４ 地球、５ 構体、１０ 静止衛星、１１ 軌道決定部、１２ 角運動量算出部、１３，１６ 指令値算出部、１４ スラスタ制御部、１５ ジンバル制御部、２１，２２，２３，２４ スラスタ、３１，３２，３３，３４ ジンバル機構、５０，５１ 面、１０１ 処理回路、１０２ プロセッサ、１０３ メモリ、１１０ センサ、１３１ 第１関数算出部、１３２ 第２関数算出部、１３３ 第１制約条件設定部、１３４ 第２制約条件設定部、１３５ 求解部、１３６ 第３制約条件設定部。

Claims (14)

1. 天体の周りの軌道を周回する宇宙機が有する角運動量を算出する角運動量算出部と、
   それぞれが前記宇宙機の構体に、ジンバル機構を介して取り付けられた複数のスラスタの内、少なくともいずれかのスラスタの噴射を指示するスラスタ指令値、および前記スラスタの噴射時における前記ジンバル機構の角度を指示する角指令値を算出する指令値算出部と、
   前記スラスタ指令値に基づき、前記スラスタを制御するスラスタ制御部と、
   前記角指令値に基づき、前記ジンバル機構を制御するジンバル制御部と、
   を備え、
   前記指令値算出部は、前記軌道における複数の噴射区間での前記スラスタ指令値および前記角指令値を引数として、該スラスタ指令値に基づいて前記スラスタが制御され、該角指令値に基づいて前記ジンバル機構が制御された場合に、前記宇宙機が前記軌道を一周する間の、前記宇宙機の軌道を表す軌道要素の変化量と前記角運動量の変化量とを出力する第１の関数、ならびに、前記軌道要素、前記軌道要素の時間変化率、および前記角運動量を引数として、前記宇宙機が前記軌道を一周する間の、前記軌道要素の制御量および前記角運動量の制御量を出力する第２の関数に基づき、前記軌道要素の変化量と前記軌道要素の制御量とが一致し、かつ、前記角運動量の変化量と前記角運動量の制御量とが一致する第１の制約条件の下で、前記複数の噴射区間における前記スラスタの噴射量の合計が目的関数である非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
   宇宙機制御装置。
2. 前記指令値算出部は、前記目的関数を最小化する前記非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
   請求項１に記載の宇宙機制御装置。
3. 前記指令値算出部は、前記ジンバル機構の角度が定められた範囲内にある第２の制約条件の下で、前記非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
   請求項１または２に記載の宇宙機制御装置。
4. 前記指令値算出部は、前記ジンバル機構の角度が定められた範囲内にある第２の制約条件と、前記複数の噴射区間のそれぞれの噴射直前および噴射直後の少なくともいずれかにおいて、前記角運動量が定められた範囲内にある第３の制約条件との下で、前記非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
   請求項１または２に記載の宇宙機制御装置。
5. 前記第２の関数が算出する前記軌道要素の制御量は、前記宇宙機が前記軌道を一周する間に、前記宇宙機に作用する摂動力によって生じる前記軌道要素の変動を補正するフィードフォワード制御量、および、前記軌道要素と前記軌道要素の目標値とに基づくフィードバック制御量を含む、
   請求項１から４のいずれか１項に記載の宇宙機制御装置。
6. 前記第２の関数が算出する前記角運動量の制御量は、前記宇宙機が前記軌道を一周する間に、前記宇宙機に作用する摂動トルクによって生じる前記角運動量の変動を補正するフィードフォワード制御量、および、前記角運動量と前記角運動量の目標値とに基づくフィードバック制御量を含む、
   請求項１から５のいずれか１項に記載の宇宙機制御装置。
7. 前記軌道要素は、離心率ベクトル、傾斜角ベクトル、および直下点経度を含む、
   請求項１から６のいずれか１項に記載の宇宙機制御装置。
8. 前記第２の関数は、前記宇宙機の軌道保持制御において求められる緯度の保持精度に応じた、年平均傾斜角ベクトル、月平均傾斜角ベクトル、および日平均傾斜角ベクトルのいずれかに対する制御量を含む前記軌道要素の制御量を出力する、
   請求項７に記載の宇宙機制御装置。
9. 前記構体の外面の内、いずれかの面は、常に前記天体に向いており、
   前記複数のスラスタは、噴射方向が互いに異なる向きで、前記構体の外面の内、常に前記天体に向いている面と反対側に位置する面に取り付けられる、
   請求項１から８のいずれか１項に記載の宇宙機制御装置。
10. 前記ジンバル機構の自由度は２であり、
    前記複数のスラスタの個数は４つであり、
    前記スラスタは、前記宇宙機の質量中心から離れる方向に伸び、
    前記噴射区間において４つの前記スラスタの内、２つの前記スラスタが噴射される、
    請求項１から９のいずれか１項に記載の宇宙機制御装置。
11. 前記ジンバル機構の角度は、前記スラスタの噴射によって前記宇宙機に作用するトルクを、前記宇宙機が有する姿勢制御系で補償可能な外乱トルクの上限値以下にする範囲内の値である、
    請求項１から１０のいずれか１項に記載の宇宙機制御装置。
12. 前記スラスタ指令値は、前記スラスタの噴射位相、前記複数のスラスタの噴射量の合計、前記合計において前記スラスタ指令値によって噴射が指示される前記スラスタの噴射量の割合である、
    請求項１から１１のいずれか１項に記載の宇宙機制御装置。
13. 軌道における複数の噴射区間でのスラスタ指令値および角指令値を引数として、該スラスタ指令値に基づいてスラスタが制御され、該角指令値に基づいてジンバル機構が制御された場合に、軌道要素の変化量および角運動量の変化量を出力する第１の関数、ならびに、軌道要素、前記軌道要素の時間変化率、および前記角運動量を引数として、前記軌道要素の制御量および前記角運動量の制御量を出力する第２の関数に基づき、予め定めた制約条件を満たしつつ、前記複数の噴射区間における前記スラスタの噴射量の合計が目的関数である非線形計画問題について、前記目的関数の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
    宇宙機制御方法。
14. コンピュータを、
    天体の周りの軌道を周回する宇宙機が有する角運動量を算出する角運動量算出部、
    それぞれが前記宇宙機の構体に、ジンバル機構を介して取り付けられた複数のスラスタの内、少なくともいずれかのスラスタの噴射を指示するスラスタ指令値、および前記スラスタの噴射時における前記ジンバル機構の角度を指示する角指令値を算出する指令値算出部、
    前記スラスタ指令値に基づき、前記スラスタを制御するスラスタ制御部、および、
    前記角指令値に基づき、前記ジンバル機構を制御するジンバル制御部、
    として機能させるためのプログラムであって、
    前記指令値算出部は、前記軌道における複数の噴射区間での前記スラスタ指令値および前記角指令値を引数として、該スラスタ指令値に基づいて前記スラスタが制御され、該角指令値に基づいて前記ジンバル機構が制御された場合に、前記宇宙機が前記軌道を一周する間の、前記宇宙機の軌道を表す軌道要素の変化量と前記角運動量の変化量とを出力する第１の関数、ならびに、前記軌道要素、前記軌道要素の時間変化率、および前記角運動量を引数として、前記宇宙機が前記軌道を一周する間の、前記軌道要素の制御量および前記角運動量の制御量を出力する第２の関数に基づき、前記軌道要素の変化量と前記軌道要素の制御量とが一致し、かつ、前記角運動量の変化量と前記角運動量の制御量とが一致する第１の制約条件の下で、前記複数の噴射区間における前記スラスタの噴射量の合計が目的関数である非線形計画問題の解を求めることで、前記複数の噴射区間のそれぞれの前記スラスタ指令値および前記角指令値を算出する、
    プログラム。

Images