ロータ軸の回りを回転するロータと、当該ロータ軸に直交するジンバル軸の回りに上記ロータを回転させるジンバルとを有し、人工衛星に搭載された複数台のCMGと、当該CMGの上記ジンバルの回転を制御するジンバル制御部と、人工衛星の姿勢を制御する姿勢制御部と、姿勢制御部の出力である姿勢制御トルクをCMGのジンバル角速度に分解するステアリング制御部と、を備え、当該ステアリング制御部において、上記CMGが特異状態に近い場合には、人工衛星の微小角運動量変化と上記CMGの微小ジンバル角変化の関係から、CMGの微小ジンバル角変化を求め、求めたCMGの微小ジンバル角変化から得られるCMGのジンバル角速度をステアリング制御部の出力とすることを特徴とする人工衛星の姿勢制御装置。ステアリング制御部の出力とする上記CMGのジンバル角速度を求める場合に使用する上記微小ジンバル角変化は、人工衛星の微小角運動量変化をCMGの微小ジンバル角変化の2次項まで展開し、その関係に基づいてCMGの微小ジンバル角変化の近似値として求めたものであることを特徴とする請求項1に記載の人工衛星の姿勢制御装置。求めた微小ジンバル角変化の近似値は、人工衛星の微小角運動量変化をCMGの微小ジンバル角変化の2次項まで展開した関係に基づいて、CMGの微小ジンバル角変化を関係づける連立代数方程式を解析的に解いて得た値としたことを特徴とする請求項2に記載の人工衛星の姿勢制御装置。人工衛星の微小角運動量変化から上記CMGの微小ジンバル角変化を求める場合において、上記CMGの微小ジンバル角変化量の大きさに閾値を設け、当該微小ジンバル角変化量が所定の大きさを超える場合には、その大きさを閾値の値で制限することを特徴とする請求項1に記載の人工衛星の姿勢制御装置。特許文献1に記載のように、人工衛星の姿勢制御を行うにあたっては、CMGが特異状態になって、CMGのジンバル角速度と人工衛星の姿勢制御トルクとを関係づけるヤコビ行列の逆行列を求められないことが問題となる。この問題は、以下のように表すことができる。人工衛星の姿勢制御トルクをτ、CMGのジンバル角をθ、ジンバル角速度をdθdtとするときに、姿勢制御トルクτとジンバル角速度dθdtの間には次の関係が成り立つ。ここでτは3次元のベクトル、dθdtはのCMG数をnとすればn次元のベクトルであり、Aは両者を関係づける3×n次の行列であり、ヤコビ行列と呼ばれる。人工衛星のフィードバック姿勢制御を行うにあたっては、姿勢制御トルクτからCMGのジンバル角速度dθdtを求める必要があり、ヤコビ行列Aのランクが3であれば以下のように求めることができる。ここで右肩添字のTは転置、-1は逆行列をとる操作を示す以下同様。しかしこのヤコビ行列AのランクはCMGのジンバル角によって変化する。とくにCMGが特異と呼ばれる状態ではヤコビ行列Aのランクが2以下になって、上式の逆行列を求めることができない。このことはCMGが特異に近い状態であっても同様であり、具体的にはヤコビ行列Aの条件数が10以上であると、上式の逆行列を求めることはできても、その値が過大になって、ジンバル角速度dθdtの値も過大になり支障を生じる。このような特異に近い状態も含めてCMGの特異状態と呼ぶと、CMGの特異状態はしばしば出現する状態であり、CMGを用いる人工衛星の姿勢制御系の大きな技術的課題であり、その対処法がつねに問題となっている。実際、CMGの特異状態において、ヤコビ行列の逆行列を近似的に求めるために、さまざまな工夫が試みられてきた。その中で代表的な方法のひとつは、上述の特許文献1に記載されている方法である。この方法によれば、姿勢制御トルクτからCMGのジンバル角速度dθdtをつぎのように求める。ここでκはCMGが特異状態になっても逆行列を求められるようにする任意の正定数であり、Iは単位行列、またPはCMGが特異状態を回避できるように微小擾乱を与える行列であり、つぎの形で表される。ここでεii=1、2、3は、時変の微小擾乱を与える項であり、CMGが特異状態にとどまるのを回避する働きをもつ。εiは、例えば次式5で与えられる。この方法では、定数κと行列Pを従来のヤコビ行列に加えることによって、CMGが特異状態になったときに、ジンバル角速度dθdtが発散するのを防ぐとともに、CMGが特異状態にとどまるのを回避している。これらは、CMGの特異状態において一定の効果を発揮するが、計算の中にκやεiのように任意に設定することのできる定数が含まれており、その与え方によって効果も異なったものとなるが、与え方の明確な指針はなく、試行錯誤的に決定される。また、κによってCMGの特異状態におけるジンバル角速度dθdtの発散は防ぐことができるが、実現される姿勢制御トルクτが姿勢制御系から与えられる姿勢制御トルクと大きく異なる可能性がある。さらに、εiによってCMGが特異状態にとどまるのを回避する効果は期待できるが、つねに有効であるという保証はなく、εiを用いても特異状態にとどまる状況もありうる。以上述べたように、このような人工衛星の姿勢制御装置にあっては、CMGが特異状態にあるときに、精度よくかつ確実に特異状態を回避するように、人工衛星の姿勢制御トルクτからCMGのジンバル角速度dθdtを求める手段を提供するものではない。また、任意に設定することのできるパラメータが多く、それらは人工衛星の慣性モーメントやCMGのもつ角運動量、CMGのジンバル角速度の上限値などを考慮して試行錯誤的に決定する必要があり、効果的に適用するためには多くの事前の検証を必要とするなどの問題点もあった。この発明は、上記のような問題点を解決するためになされたものであり、CMGが特異状態にあっても姿勢制御トルクからジンバル角速度を、確実にかつ精度よく求めることのできる人工衛星の姿勢制御装置を得ることを目的としている。本発明に係る人工衛星の姿勢制御装置は、ロータ軸の回りを回転するロータと、当該ロータ軸に直交するジンバル軸の回りに上記ロータを回転させるジンバルとを有し、人工衛星に搭載された複数台のCMGと、当該CMGの上記ジンバルの回転を制御するジンバル制御部と、人工衛星の姿勢を制御する姿勢制御部と、姿勢制御部の出力である姿勢制御トルクをCMGのジンバル角速度に分解するステアリング制御部と、を備え、当該ステアリング制御部において、上記CMGが特異状態に近い場合には、人工衛星の微小角運動量変化と上記CMGの微小ジンバル角変化の関係から、CMGの微小ジンバル角変化を求め、求めたCMGの微小ジンバル角変化から得られるCMGのジンバル角速度をステアリング制御部の出力とするようにしたものである。この発明によれば、複数台のCMGにおいて、CMGが特異状態にあっても、微小角運動量変化と微小ジンバル角変化の間の非線形関係に基づいて微小ジンバル角変化の解を求めるようにしたので、確実に微小ジンバル角変化の解を求めることができる。またCMGが特異状態にあってもジンバル角の変化を最小限に抑えることができ、かつ、CMGが特異状態にとどまることも自然に回避できる、といった従来にない顕著な効果を奏するものである。実施の形態1.図1は本発明において用いる人工衛星の姿勢制御用アクチュエータである単一ジンバルCMGの概念図である。CMGはスピン軸回りに高速で回転するホイールとそのホイールの回転軸を回転させるためのジンバル軸から構成される。ホイールの高速回転によりスピン軸回りに一定量の角運動量を有するが、その角運動量の向きをジンバル軸回りに回転させることにより人工衛星本体に大きな反作用トルクが働き、このトルクを姿勢制御に利用する。図2はこのCMGを人工衛星に複数台配置する場合の一例を示す。CMGは通常、この図のように4台程度が人工衛星に配置され、これらの複数台のCMGを協調して動作させることで人工衛星の姿勢制御を実現する。図2の場合は、各CMGが四角錐の底辺に配置されることからピラミッド配置と呼ばれ、CMGの典型的な配置例の一つである。図3は本発明における人工衛星の姿勢制御装置のブロック図の一例を示すものである。人工衛星の姿勢目標値2が与えられると、姿勢制御部3では、姿勢目標値2と人工衛星の姿勢角との偏差から、人工衛星の姿勢角を姿勢目標値2に近づけるようにフィードバック制御を行う。この姿勢制御部3の出力として人工衛星の姿勢制御トルクが得られるが、この姿勢制御トルクを実現するためにはCMGのジンバル角の動きによって人工衛星本体に所望の反作用トルクが働くようにする必要がある。このためにCMGのジンバル角の角速度、角度の目標値を与えるのがステアリング制御部4であり、本発明はこのステアリング制御部4に係るものである。人工衛星の姿勢制御トルクτ3次のベクトルとCMGのジンバル角速度dθdtn次のベクトルとの間にはつぎの線形関係がある。ただしヤコビ行列はジンバル角θの関数となるので、それを明示するために以下ではAと表す。上式6において、CMGが特異状態にあるときには、ヤコビ行列Aのランクが2以下になるので、姿勢制御トルクτを実現するようなジンバル角速度dθdtを求めることができない。一方、CMG全体のもつ角運動量をhとすれば、hもジンバル角θの関数hとなり、上式を次式で表すこともできる。ここで姿勢制御トルクτが一定と見なされる短い時間を微小時間Δtとする。微小時間Δtの間隔は姿勢制御系の帯域によって異なるが、ほぼ02~03秒程度の値と考えることができる。この微小時間Δtの間のジンバル角変化をΔθとすれば、次式8が得られる。本発明では、CMGが特異状態にあって、姿勢制御トルクτからジンバル角速度dθdtを決定できないときに、上式に基づいてジンバル角度の微小変化分Δθを決定するものである。これを図4のフローチャートに基づいてさらに詳しく説明する。上式はΔθに対する非線形の方程式なので、通常はこれを解析的に解くことはできない。しかし、時間を上式のΔtまで延ばすことによって得られるこの非線形関係を用いると、線形関係では省略されていた2次以上の微小項の影響も含まれることになるので、線形関係では解の得られない場合にも解を得ることが可能になる。この解を繰り返し計算で得ることにすると、まず解Δθの適当な近似値を設定する必要がある。これがS1にあたる。この近似値をΔθiとおくと、次式9によって角運動量誤差Δherrを求めることができる。上式の計算を行うのがS2にあたる。このΔherrの大きさを0と見なすことができれば、このときのΔθiをΔθの解として採用することができる。また、このΔherrの大きさを0と見なすことができない場合には、Δθiの値を更新して、この値をより小さくするようにする必要がある。この判断を行うのがS3である。仮にΔherrが十分に小さくはなく、更新の必要がある場合には、の状態に近づけるようなΔθiを、Δθの更新値Δθi1として用いる。この更新を行うのがS4である。Δherrの大きさが十分に小さくなるまで、この更新を繰り返すことにより、十分な精度を有するジンバル角の変化量Δθが得られる。こうしてΔθが得られれば、ジンバル角速度dθdtはとして求めることができる。このように、ステアリング制御部4において、従来の姿勢制御トルクτとジンバル角速度dθdtの関係ではなく、姿勢制御トルクに微小時間をかけたτΔtと微小ジンバル角変化Δθの間の、Δθの2次以上の項も含む関係を用いるようにしたので、2次以上の項の効果によって、Δθの1次だけを含む線形関係では解の得られない場合にもΔθの解を求めることができ、CMGが特異状態にあっても、CMGのジンバル角が発散することのない、信頼性の高い人工衛星の姿勢制御装置を構成できる。実施の形態2.図5は本発明の実施の形態2に係るジンバル角の微小変化量の近似値を求めるフローチャートである。hをテイラー級数展開することによってτΔtは次式12で近似できる。ここでHはヘッセ行列と呼ばれ、各成分を3次のベクトルとするようなn×n行列となる。上式の右辺第1項は通常の姿勢制御トルクτとジンバル角速度dθdtの線形関係を表しており、右辺第2項が2次の非線形項を表している。この関係が図5のS6にあたる。実施の形態2においては、この関係をもとに、この関係を満たすΔθの値を導く方法を以下に示す。この関係は、n次のベクトルであるΔθの2次項を含む非線形な関係なので、2次ではあるが解を解析的に与えることは容易ではない。そこで、この場合にも解を繰り返し計算で求めることにして、最初にΔθのの近似値を初期値として与える。これが図5のS7にあたる。この場合の近似値はΔθ=0、または、適当な微小量を設定しておけばよい。つぎにこの近似値をΔθiとして、次式によって角運動量誤差Δherrを求める。上式の計算を行うのがS8にあたる。このΔherrの大きさを0と見なすことができれば、このときのΔθiをΔθの解として採用することができる。また、このΔherrの大きさを0と見なすことができない場合には、Δθiの値を更新して、この値をより小さくするようにする必要がある。この判断を行うのがS9である。かりにΔherrが十分に小さくはなく、更新の必要がある場合には、Δθの更新値Δθi1を次式14で求める。ただし記号†は擬似逆行列を表す。この更新を行うのがS10である。Δherrの大きさが十分に小さくなるまで、この更新を繰り返すことにより、Δθの2次の項を含む関係に対して、十分な精度を有するジンバル角の変化量Δθが得られる。この場合、対象とする方程式は2次の非線形項しか含んでいないので、一般の非線形の場合に比べれば解を得ることは容易である。この方法は、Δherr=0を実現するようなΔθを近似的に求めることになるので、図4のフローチャートを実行する過程において用いれば、最終的な、2次以上の非線形を含む方程式の解を効率的に求めることができる。具体的には、図4のフローチャートのS1およびS4においてΔherr=0を実現するΔθの近似値が必要であり、そのステップにおいて図5のフローチャートを実行する。その結果、本来の目的である、τΔt=h-hを満たすΔθを短時間で得ることが可能になる。以上のように、実施の形態2に係る発明では、求める方程式をΔθの2次までの項に限定し、その範囲でΔθの解を求めるようにしたので、この結果を用いれば、本来のΔθの3次以上の項を含む方程式の解も効率的に求めることができ、CMGが特異状態にあってもジンバル角速度の解を短時間で確実に求めることができる。実施の形態3.図6は本発明の実施の形態3に係るジンバル角の微小変化量の近似値を求めるフローチャートである。前述したように、hをテイラー級数展開することによってτΔtは次式15で近似できる。ここで上式はn次元ベクトルΔθに関する3個の連立方程式である。この関係を図6のS12に示すが、これは図5のS6と同じものである。通常、この方程式はΔθの各成分が連立した複雑な2次方程式となるので、容易に解くことはできず、実装性の面からは、図5のように繰り返し計算を用いるのが一般的方法となる。しかし、高次の代数方程式を解くことにすれば、上式を解析的に扱うこともできる。上式は3成分を有するので、上式を満たすのに必要なΔθの成分は3つとなる。これをΔθにおける従属変数ΔθDとする。それに対して、Δθの残りのn-3個の成分は独立に値を設定しうる成分となるので、これを独立変数ΔθIとする。このステップがS13にあたる。この独立変数ΔθIについては任意に値を設定できるが、たとえば0ゼロを設定しておけばよい。この値を設定する部分がS14である。このとき、グレブナー基底例えば、野呂正行、横山和弘著、グレブナー基底の計算基礎篇、東京大学出版会、2003年参照のような解計算を単純化する数学的手段を用いることにより、S12に示す3個の連立方程式の解計算は、以下の3個の代数方程式を解くことに帰着できる。1個はΔθDのある1成分に対する8次の代数方程式である。残りの2個は、ΔθDの残りの2成分に対する1次の代数方程式である。ここで後者の残りの2個の代数方程式の解計算では、前者の8次の代数方程式の解が用いられる。したがって、この3個の代数方程式を解くことにすれば、繰り返し計算を用いないでも解を求めることが可能になる。このステップがS15である。こうして、Δθに関する解は得られるが、ΔθIの値の設定によってΔθDの解も異なるものとなるので、解を評価して、それが満足のいくものであるかどうかを判断する。これを行うステップがS16である。これは、たとえばΔθの解の大きさΔθのノルムで評価すればよい。もしも、解が望ましくない場合には、Δθの解の大きさが小さくなるようにΔθIの値を更新することになるが、これを行うステップがS17になる。ただし、この方法ではS12に示す方程式を満たす解析解が得られるので、ΔθIの値の更新は必ずしも必要ではない。例えば、ΔθIの初期値によって得られるΔθの解の大きさが満足のいくものであれば、ステップS17は省略することができる。以上のように、実施の形態3に係る発明では、求める方程式をΔθの2次までの項に限定し、その範囲でΔθの解析解を求めるようにしたので、繰り返し計算を用いることなく、CMGが特異状態にあってもジンバル角速度の解を短時間で確実に求めることができる。実施の形態4.つぎに実施の形態4に係る姿勢制御装置の動作について説明する。CMGが特異状態にあるときに、人工衛星の姿勢制御トルクを微小時間Δtで積分した微小角運動量変化を実現するような、CMGのジンバル角変化Δθを求めると、CMGが通過困難な特異状態にある場合には、Δθの値が必ずしも微小ではない状況が発生しうる。この場合、Δt間でこのΔθを実現することになるので、ΔθΔtの大きさが過大になり、ジンバル角速度がハードウェアの上限値を超えてしまう可能性がある。このような状況に対処するために、ジンバル角の微小変化Δθの大きさ|Δθ|に閾値δを設けて、|Δθ|>δの場合には、次式16でΔθを修正する。このΔφをΔθの代わりに用いると、結果的に実現される姿勢制御トルクτも修正を受ける。これは特異状態において回避しづらいトルクの向きを避けることになり、結果的に特異状態をスムーズに回避することにつながる。以上のように、実施の形態4に係る発明では、CMGが特異状態にあるときに、人工衛星の微小角運動量変化姿勢制御トルクτを一定と見なせる微小時間Δtにおける、角運動量の変化からCMGのジンバル角変化を求め、その大きさがある閾値を超える場合には、その大きさを閾値で制限するようにしたので、実現される姿勢制御トルクも特異状態を回避するように修正を受け、CMGが通過困難な特異状態にあっても自然に特異状態を回避することが可能になる。なお、ここではCMGのジンバル角変化をベクトルとしてとらえ、その大きさを制限するようにしたが、各CMGのジンバル角変化に対して閾値を設定するようにしても、同様の効果を得ることができる。