JP4156544B2 - 出力分配算出装置 - Google Patents

Description

この発明は、例えば、人工衛星の姿勢制御装置や、４輪独立駆動車両の姿勢制御装置等において、システム全体で所望の出力を実現するための、各アクチュエータ出力を決定する場合に用いられる出力分配算出装置に関するものである。

機械システム等において、与えられた力、あるいはトルク等をある個数のアクチュエータで実現する際に、与えられた力、あるいはトルクの自由度に対して、アクチュエータの自由度（個数）が多ければ、その力、あるいはトルクを実現するためのアクチュエータへの力、あるいはトルクの分配は一般には一意ではない。ここでは、与えられた力、あるいはトルク等をＭ次元ベクトルｆで表し、決定すべき個々のアクチュエータの力あるいはトルクなどを並べてＮ次元ベクトルｘで表したときに、ｆとｘの関係が以下のように線形であるようなものを考える。

ｆ＝Ａｘ （１０１）
ここで、Ａはｘからｆを求めるＭ×Ｎの係数行列である。上述のようにＮ＞Ｍの場合を考えると、与えられたｆに対して式（１０１）を満たすｘは一意ではなく、ある程度設定の自由度がある。そこで、この冗長自由度を利用して与えられた目的において最適となるようにｘを決定することが考えられる。ここでは、個々のアクチュエータ出力を極力抑えることを目的とし、以下のような目的関数を設定する。

ｍｉｎ ｍａｘ_ｉ｜ｘ_ｉ｜ （１０２）
ここで、式（１０２）は式（１０１）を満たすｘのうちで、各要素ｘ_ｉの絶対値の最大値が最小となるようなｘを求めることを意味する。

上述のように定式化される問題は、現実問題として様々な分野で現れる。例えば、アクチュエータとして４個以上のリアクションホイールを持つ人工衛星の姿勢制御装置において、人工衛星に所望の姿勢変更を与えるための姿勢駆動トルクコマンドｆから、それを実現するためのリアクションホイールの駆動トルクｘを決定する場合に、ホイールの数が姿勢駆動トルクｆの次元である３よりも多いために、ｆからｘの決定は一意ではない。

そのため、ホイールトルクの飽和を抑えつつ、より効率的に姿勢を駆動するために、ホイールトルク最大値を最小にするようにホイールトルクｘを決定することが考えられる。例えば、特許文献１に記載の人工衛星の姿勢制御装置においては、ホイールトルクの飽和を抑えつつ衛星の姿勢角加速度を最大にするために、与えられた姿勢角加速度目標値からホイールトルクを決定する際に、冗長自由度を表す自由パラメータを１軸あるいは２軸まわりの姿勢変更時の近似最適値に比例するように設定し、ホイールトルクの分配を求めていた。

特開２００３−２５２２９９号公報

しかしながら、従来の出力分配方法では、計算量を少なくするためにいくつかの近似が用いられるため、任意の姿勢駆動トルクに対して、式（１０１）、（１０２）で表される問題の最適解としてのホイールトルクを与えるとは限らず、ホイールトルクの飽和を抑え、姿勢駆動の効率を向上するという意味で、確実に最適な分配を与えるわけではない。また、近似的に導出されたアルゴリズムであることから、実際に人工衛星に搭載して実施する上で、信頼性と性能に不安があった。

この発明は上記のような課題を解決するためになされたもので、各アクチュエータの出力を極力抑えつつ、システム全体として所望の出力を得ることのできる出力分配算出装置を得ることを目的とする。

この発明に係る出力分配算出装置は、複数のアクチュエータの出力であるベクトルｘと、制御対象の出力であるベクトルｆが、ｆ＝Ａｘによって関係付けられ、かつベクトルｘの次元Ｎがベクトルｆの次元Ｍよりも大きいモデルに対し、ベクトルｆが与えられたときにそれを実現するためのベクトルｘを求める場合、行列Ａの第ｉ列のＭ次元ベクトルをａ_ｉ（ｉ＝１，・・・Ｎ）として、全てのベクトルａ_ｉ（ｉ＝１，・・・Ｎ）との内積の絶対値の総和がある正の定数ＣとなるようなＭ次元ベクトルのうちで、ベクトルｆとの内積の値が最大となるようなベクトルｙ^＊を求めるｙ^＊演算手段と、ベクトルｙ^＊とベクトルｆとの内積の値を定数Ｃで割った値ｘ_ｍを求めるｘ_ｍ演算手段と、ベクトルａ_ｉとベクトルｙ^＊との内積が０でないｉに対して、ｘ_ｉの値を、ｘ_ｍにａ_ｉとｙ^＊との内積の符号を付した値とする第１のｘ_ｉ演算手段と、ベクトルａ_ｉとベクトルｙ^＊との内積が０となるｉに対して、ｘ_ｉの値を、その絶対値がｘ_ｍを超えないように求める第２のｘ_ｉ演算手段とを備えたものである。

この発明の出力分配算出装置は、行列Ａの第ｉ列のＭ次元ベクトルをａ_ｉとして、全てのベクトルａ_ｉとの内積の絶対値の総和がある正の定数ＣとなるようなＭ次元ベクトルのうちで、ベクトルｆとの内積の値が最大となるようなベクトルｙ^＊を求め、ベクトルｙ^＊とベクトルｆとの内積の値を定数Ｃで割った値ｘ_ｍを求め、ベクトルａ_ｉとベクトルｙ^＊との内積が０でないｉに対しては、ｘ_ｉの値を、ｘ_ｍにａ_ｉとｙ^＊との内積の符号を付した値とし、また、ベクトルａ_ｉとベクトルｙ^＊との内積が０となるｉに対しては、ｘ_ｉの値を、その絶対値がｘ_ｍを超えないように求めたので、与えられたベクトルｆを実現するベクトルｘのうちで、その各要素絶対値の最大値が最小となるような分配を常に求めることができ、アクチュエータの出力を極力抑えつつ全体として所望の出力が得られるような分配を行うことができる。

実施の形態１．
図１は、この発明の実施の形態１による出力分配算出装置のブロック構成図である。
図の装置は、ｙ^＊演算手段１、ｘ_ｍ演算手段２、第１のｘ_ｉ演算手段３、第２のｘ_ｉ演算手段４からなる。ｙ^＊演算手段１は、行列Ａの第ｉ列のＭ次元ベクトルをａ_ｉ（ｉ＝１，・・・Ｎ）として、全てのベクトルａ_ｉ（ｉ＝１，・・・Ｎ）との内積の絶対値の総和がある正の定数ＣとなるようなＭ次元ベクトルのうちで、ベクトルｆとの内積の値が最大となるようなベクトルｙ^＊を求める機能部である。

ｘ_ｍ演算手段２は、 ベクトルｙ^＊とベクトルｆとの内積の値を前記定数Ｃで割った値ｘ_ｍを求める機能部である。第１のｘ_ｉ演算手段３は、ベクトルａ_ｉとベクトルｙ^＊との内積が０でないｉに対して、ｘ_ｉの値を、ｘ_ｍにａ_ｉとｙ^＊との内積の符号を付した値とする機能部である。第２のｘ_ｉ演算手段４は、ベクトルａ_ｉとベクトルｙ^＊との内積が０となるｉに対して、ｘ_ｉの値を、その絶対値がｘ_ｍを超えないように求める機能部である。

本発明の出力分配算出装置は、コンピュータを用いて実現され、上記ｙ^＊演算手段１〜第２のｘ_ｉ演算手段４は、それぞれの機能に対応したソフトウェアと、これらのソフトウェアを実行するためのＣＰＵやメモリ等のハードウェアから構成されている。

次に、本実施の形態の出力分配算出装置を用いた出力分配方法を説明する。
本発明では、以下のようにモデル化された力あるいはトルク分配問題を考える。
ｆ＝Ａｘ ・・・（１）
ｍｉｎ_ｘ ｍａｘ_ｉ｜ｘ_ｉ｜ ・・・（２）
ここで、ｆはＭ次元ベクトル、ｘはＮ次元ベクトル、ＡはＭ×Ｎ行列で、Ｎ＞Ｍとする。ｘの各要素ｘ_ｉ（ｉ＝１，・・・，Ｎ）は各アクチュエータの駆動力、あるいは駆動トルク等の出力を表し、ｆの各要素ｆ_ｉ（ｉ＝１，・・・，Ｍ）はこれらアクチュエータの駆動力あるいはトルクの合成により実現される力あるいはトルク等のシステム全体としての出力を表す。

式（２）は、式（１）を満たすｘのうちで、各要素の絶対値｜ｘ_ｉ｜の最大値を最小にするようなものを求めることを意味する。
全体として実現すべき力あるいはトルクｆが与えられたときに、問題（１）、（２）の最適解となるようなアクチュエータの出力ｘを求めることができれば、全てのアクチュエータの出力を極力抑えることができ、ｆを実現する上で最も効率的と言える。

問題（１）、（２）は一般にミニマックス問題と呼ばれる最適化問題であり、変数の次元Ｎがそれほど大きくなければ既存の解法を用いて市販のコンピュータ等で解くことができる。しかし、一般に、そのような解法は収束計算等が必要になるため、制御装置等に組み込んで制御ループの中で用いる際には、制御の時間刻みに対して十分に短い時間で確実に解を求めることが必要であり、計算量の観点で現実的ではない。

本発明においては、問題（１）、（２）を直接解くことなく、より簡単な計算により最適な分配ｘを求める。以下、本実施の形態の出力分配算出装置の動作としての分配方法をフローチャートに従って説明する。

図２は、この発明の実施の形態１による出力分配算出装置における出力分配方法を示すフローチャートである。
先ず、与えられたｆに対し、ｙ^＊演算手段１は、次のような問題の解を求める（ステップＳＴ１）。
ｍａｘ_ｙ ｆ・ｙ ・・・（３）
制約条件：Σ_ｉ＝_１ ^Ｎ｜ａ_ｉ・ｙ｜＝Ｃ ・・・（４）
ここで、ｆ・ｙは、ベクトルｆとｙとの内積を表す。ａ_ｉ（ｉ＝１，・・・，Ｎ）は、行列Ａの第ｉ列であるＭ次元ベクトル、Σ_ｉ＝_１ ^Ｎはｉ＝１からｉ＝Ｎまでの和を表す。Ｃは適当な正の定数、例えばＣ＝１とする。問題（３）、（４）は、式（４）の制約条件を満たすｙのうちで、ｆとの内積ｆ・ｙの値が最大となるものを求めることを意味する。以下では、問題（３）、（４）をもとのミニマックス問題（１）、（２）に対して、双対問題と呼ぶ。双対問題の最適解を求めるための具体的な手法については、後に説明する。

双対問題の最適解ｙ＝ｙ^＊が得られると、次に、ｘ_ｍ演算手段２は、ｘ_ｍを以下のように与える（ステップＳＴ２）。
ｘ_ｍ＝ｆ・ｙ^＊／Ｃ ・・・（５）
ｘ_ｍは、求められるアクチュエータへの力、あるいはトルクの分配において、力、あるいはトルクの絶対値の最大値となる。式（５）において、内積ｆ・ｙ^＊の値は双対問題（３）、（４）の目的関数値であるため、双対問題を解く過程で既に計算されている。式（５）の計算は実際にはそれをＣで割る計算だけである。Ｃ＝１としておけば、その計算も不要であり、ステップＳＴ２の計算はステップＳＴ１の計算と同時に行うことができる。

次に、第１のｘ_ｉ演算手段３は、ｙ^＊およびｘ_ｍを用いて、ベクトルａ_ｉとｙ^＊の内積ａ_ｉ・ｙ^＊が０とならないｉに対するｘ_ｉの値を、ｘ_ｍにａ_ｉ・ｙ^＊の符号を付した値として与え、それ以外のｉに対するｘ_ｉの値は、一旦ｘ_ｉ＝０としておく。即ち、ｘ_ｉを以下のように与える（ステップＳＴ３）。

ｘ_ｉ＝ｓｉｇｎ（ａ_ｉ・ｙ^＊）ｘ_ｍ （ｉ＝１，・・・，Ｎ） ・・・（６）
ここで、ｓｉｇｎ（ｔ）はｔ＞０のとき＋１、ｔ＜０のとき−１、ｔ＝０のとき０となる符号関数である。ａ_ｉ・ｙ^＊＝０となるｉに対してｘ_ｉ＝０とするのは、単に記述の簡単化のためであり、実際はここで値を与える必要はない。０とする手順を省いた場合には、後の式（８）の右辺第２項において、Ａとｘとから対応する列、要素をそれぞれ除いておけばよい。

次に、第２のｘ_ｉ演算手段４によって、残りのｘ_ｉ、即ち、ａ_ｉ・ｙ^＊＝０となるｘ_ｉの値を求める。ここで、ｙ^＊との内積が０となるａ_ｉの本数をＮ^（１）とおき、これらを列とするＭ×Ｎ^（１）行列をＡ^（１）、Ａ^（１）の第ｉ列をＭ次元ベクトルａ_ｉ ^（１）（ｉ＝１，・・・Ｎ^（１））、対応するｘ_ｉをまとめてＮ^（１）次元ベクトルｘ^（１）とする。ｘ^（１）の値は、次式を満たし、かつ、各要素ｘ_ｉ ^（１）が｜ｘ_ｉ ^（１）｜≦ｘ_ｍを満たすように与えられる（ステップＳＴ４）。
Ａ^（１）ｘ^（１）＝ｆ^（１） ・・・（７）
ここで、ｆ^（１）は次式で与えられる。
ｆ^（１）＝ｆ−Ａｘ ・・・（８）

但し、式（８）の右辺第２項のｘにおいて、ｘ^（１）に対応する要素は０となっている。上述のように、式（６）でこれらの要素を０にする手順を省いた場合には、式（８）の右辺第２項は、ＡとｘとからＡ^（１）、ｘ^（１）に対応する列、要素を除いて計算する。先に求められていたｘ_ｉと、ここで求められたｘ_ｉ ^（１）をあわせて、ｘの全ての値が決まる。

次に、実施の形態１に係る、双対問題（３）、（４）の具体的な解法について図３を用いて説明する。
図３は、実施の形態１に係る図２のフローチャートのステップＳＴ１のより詳細な動作を説明するものである。

双対問題（３）、（４）の最適解として、行列Ａを構成するベクトルａ_ｉ（ｉ＝１，・・・，Ｎ）のうちの少なくともＭ−１本との内積が０となるようなものが必ず存在する。そこで、最適解ｙ^＊の候補としては、これらを考えれば十分である。これは以下の理由からである。

今、相対問題（３）、（４）の最適解ｙ＝ｙ^＊が求められたとする。ｙ^＊はＭ次元ベクトルであるので、ａ_ｉのうちのＭ本以上に直交することはできない（それらａ_ｉのうちの任意の１本がその他のａ_ｉの線形結合で表される場合を除く）。そこで、一般にｙ^＊は高々Ｍ−１本以下のａ_ｉに直交する。残りのＮ−Ｍ＋１本以上の（ｙ^＊が直交しない）ａ_ｉに対しては、そのｉに対するｘ_ｉの値は、上記の式（６）のように一意的に与えられる。一方、ｘ_ｉは式（１）を満たす必要がある。即ち、ＡのＮ本の列ベクトルａ_ｉにＮ個のパラメータｘ_ｉをかけた形としてベクトルｆを表す必要がある。

ここで、例えばｙ^＊がＭ−２本のａ_ｉに直交する場合、Ｎ−Ｍ＋２個のｘ_ｉは式（６）で表されるように一つのパラメータｘ_ｍで決められている。よって、パラメータの数は、式（６）で値の決まらないＭ−２個のｘ_ｉとｘ_ｍの、合計Ｍ−１個しかなく、一般にはＭ次元ベクトルｆを表すことはできない。よって、ｙ^＊はＭ−１本以上のａ_ｉに直交する必要がある。以上により、相対問題の最適解ｙ^＊は一般にＭ−１本のａ_ｉに直交する。

上記で、例えばＭ−２個のｘ_ｉとｘ_ｍの合計Ｍ−１個のパラメータでＭ次元ベクトルｆを表すことができるような場合も考えられる。しかしそのような場合にも、それらｘ_ｍで値を与えているＮ−Ｍ＋２個のｘ_ｉのうちの一つをフリーパラメータとしてパラメータをＭ個にしても、同じ結果になることは容易にわかるので、そのような特別の場合も含めてｙ^＊としては、一般にＭ−１本のａ_ｉに直交するものを考えればよい。

以上のような理由から、ステップＳＴ１０では、準備段階として、Ｍ次元ベクトルａ_ｉ（ｉ＝１，・・・，Ｎ）の内の任意のＭ−１本に直交し、全てのａ_ｉとの内積の絶対値の総和がある正の定数Ｃとなるようなベクトルｙの集合Ｙを求める。

集合Ｙに含まれるベクトルｙとしては、方向が正反対で大きさが同じものをそれぞれ２本ずつ考えることができる。例えば、あるベクトルｙがＹに含まれていれば、符号を反転した−ｙもＹに含まれる。そこで、そのうち１本ずつを省けば、Ｙに含まれるベクトルの本数は、ホイールといったアクチュエータの個数Ｎに対して、_ＮＣ_Ｍ−１となる。

ここで、_ＮＣ_Ｍ−１は異なるＮ個の要素からＭ−１個を選ぶ組み合わせの数であり、Ｍが大きくなければそれほど大きな値ではない。例えば、Ｍ＝３の場合には_ＮＣ_Ｍ−１＝Ｎ（Ｎ−１）／２である。これらＭ次元ベクトルの集合Ｙは、行列Ａが定数として与えられる場合には、予め計算して、メモリ等の記憶部に記憶しておくことができる。そして、次のステップ以降が、ｆが与えられるごとに計算される。

ｆが与えられると、記憶されている集合Ｙを呼び出し（ステップＳＴ１１）、ｆとの内積の絶対値が最大となるようなベクトルｙ^＊を集合Ｙの中から選び、更にｙ^＊とｆとの内積が負であれば、ｙ^＊の符号を反転しておく（ステップＳＴ１２）。これを式で表すと以下のようになる。
ｙ^＊＝ａｒｇｍａｘ_ｙ∈Ｙ｜ｆ・ｙ｜ ・・・（９）
ｙ^＊＝ｓｉｇｎ（ｆ・ｙ^＊）ｙ^＊ ・・・（１０）
ここで、ａｒｇｍａｘ_ｙ∈Ｙ｜ｆ・ｙ｜は集合Ｙに含まれるｙのうち、｜ｆ・ｙ｜が最大となるものを意味する。

また、ステップＳＴ１０で集合Ｙを求める際、方向が正反対で大きさが同じ２本のｙのうち１本を省くことでＹに含まれるベクトルの本数を半分にしたが、これら２本を別のものとしてＹに含めておけば、式（９）、（１０）の計算は、式（９）で絶対値を外した以下の式にまとめることができる。
ｙ^＊＝ａｒｇｍａｘ_ｙ∈Ｙｆ・ｙ ・・・（９’）
これにより、Ｙを記憶するためのメモリ量は２倍に増え、内積ｆ・ｙの計算回数も２倍になるが、プログラムの記述を簡単にすることはできる。

また、図３のフローチャートでは、ｙ^＊の候補Ｙを予め計算しておき、記憶部に蓄えておいたが、これをｆが与えられる毎に計算することもできる。この場合計算量は増加するが、その計算量は予め見積もることができるため、計算量が許容される程度である限り、悪影響を及ぼすことはない。

次に、実施の形態１に係る、図２のフローチャートにおけるステップＳＴ４の具体的計算方法について、図面を用いて詳細に説明する。
図４は、実施の形態１に係る図２のフローチャートのステップＳＴ４のより詳細な動作を説明するものである。

ステップＳＴ４では、式（７）を満たし、かつ全ての要素が｜ｘ_ｉ ^（１）｜≦ｘ_ｍを満たすようなＮ^（１）次元ベクトルｘ^（１）を求めればよい。ここで、Ｎ^（１）はベクトルｙ^＊に直交するベクトルａ_ｉの本数であった。Ｎ^（１）≦Ｍ−１であれば、式（７）を満たすｘ^（１）は一意に求められるが、Ｎ^（１）≧Ｍであれば、ｘ^（１）は一意には求められない。そこで、ステップＳＴ４０として、ベクトルｙ^＊に直交するベクトルａ_ｉの本数の判定部分を設ける。Ｎ^（１）≦Ｍ−１である場合に、式（７）を満たすｘ^（１）の求め方としては様々な手法が考えられ、どのような手法を用いてもよいが、例えば以下のように求める（ステップＳＴ４１）。

ｘ^（１）＝（（Ａ^（１））^ＴＡ^（１））^−１（Ａ^（１））^Ｔｆ^（１）
・・・（１１）
ここで、上付きのＴは行列の転置、−１は逆行列を表す。
Ｎ^（１）≧Ｍであった場合、上述のように式（７）を満たすｘ^（１）は一意ではないため、ｍａｘ_ｉ｜ｘ_ｉ ^（１）｜≦ｘ_ｍを満たす範囲内で、何らかの基準に従ってその値を決定する必要がある。アクチュエータ出力を極力抑えることを考えると、ｍａｘ_ｉ｜ｘ_ｉ ^（１）｜が最小になるように求めることが最適である。そこで、上記同様、次のような問題の最適解となるベクトルｙ＝ｙ^（１）＊を求める。

ｍａｘ_ｙ ｆ^（１）・ｙ ・・・（１２）
制約条件：Σ_ｉ＝_１ ^Ｎ（１）｜ａ_ｉ ^（１）・ｙ｜＝Ｃ^（１） ・・・（１３）
ｙ^＊・ｙ＝０ ・・・（１４）
ここで、Ｃ^（１）は適当な正の定数、例えばＣ^（１）＝１とする。

問題（１２）、（１３）、（１４）は、形式上、問題（３）、（４）をｙ^＊に直交する部分空間に制限した形になっている。そのため、問題（３）、（４）と同様にして解くことができ、求められたｙ＝ｙ^（１）＊を用いて、ｘ^（１）の値が決定される。その際、ａ_ｉ ^（１）・ｙ^（１）＊＝０となるａ_ｉ ^（１）がＭ−１本以上存在する場合、対応するｘ_ｉ ^（１）の値はやはり一意には求められず、更にｙ^＊とｙ^（１）＊の二つに直交する部分空間での上記と同様の問題を考える必要がある。
そこで、次のような問題を考える（これをk-部分問題と呼ぶことにする）。

k-部分問題：
ｆ^（ｋ）＝Ａ^（ｋ）ｘ^（ｋ） ・・・（１５）
ｍｉｎ ｍａｘ_ｉ｜ｘ_ｉ ^（ｋ）｜ ・・・（１６）
最初のミニマックス問題（１）、（２）は形式的に0-部分問題と呼ぶことにする。即ち、ｆ^（０）＝ｆ、Ａ^（０）＝Ａ、ｘ^（０）＝ｘとする。また、1-部分問題が上記で考えたｘ^（１）を求める問題である。以下、k-部分問題の解ｘ^（ｋ）を求める手順を、図面を用いて説明する。

図５は、k-部分問題の解ｘ^（ｋ）を求める手順を示すフローチャートである。
k-部分問題（１５）、（１６）に対して、次のようにk-双対問題を考える。
k-双対問題：
ｍａｘ_ｙ ｆ^（ｋ）・ｙ ・・・（１７）
制約条件：Σ_ｉ＝_１ ^Ｎ（ｋ）｜ａ_ｉ ^（ｋ）・ｙ｜＝Ｃ^（ｋ） ・・・（１８）
ｙ^（ｊ）＊・ｙ＝０ （ｊ＝０、・・・、ｋ−１） ・・・（１９）

ここで、Ｃ^（ｋ）は適当な正の定数、例えばＣ^（ｋ）＝１とする。また、ａ_ｉ ^（ｋ）（ｉ＝１，・・・，Ｎ^（ｋ））は行列Ａ^（ｋ）の第ｉ列であり、ｙ^（ｊ）＊（ｊ＝０、・・・、ｋ−１）はj-双対問題の最適解である。最初の双対問題（３）、（４）は形式的に0-双対問題と呼ぶものとする。即ち、Ｎ^（０）＝Ｎ、Ｃ^（０）＝Ｃ、ｙ^（０）＊＝ｙ^＊とする。0-双対問題（３）、（４）と同様、k-双対問題の最適解として、ａ_ｉ ^（ｋ）（ｉ＝１，・・・，Ｎ^（ｋ））のうちのＭ−ｋ−１本との内積が０となる（かつ、制約条件（１９）よりｙ^（ｊ）＊（ｊ＝０、・・・、ｋ−１）との内積も０となる）ようなものが必ず存在する。これら最適解の候補からｙ^（ｋ）＊を求める手続きは、上記0-双対問題（３）、（４）の場合と同様である。

k-双対問題の最適解ｙ^（ｋ）＊が求められる（ステップＳＴ１０１）と、ｘ_ｉ ^（ｋ）の絶対値の最大値ｘ_ｍ ^（ｋ）は、式（５）と同様に以下の式により求められる（ステップＳＴ１０２）。
ｘ_ｍ ^（ｋ）＝ｆ^（ｋ）・ｙ^（ｋ）＊／Ｃ^（ｋ）・・・（２０）
また、内積ａ_ｉ ^（ｋ）・ｙ^（ｋ）＊が０でないｉに対しては、ｘ_ｉ ^（ｋ）の値は、式（６）と同様に以下の式により求められる（ステップＳＴ１０３）。
ｘ_ｉ ^（ｋ）＝ｓｉｇｎ（ａ_ｉ ^（ｋ）・ｙ^（ｋ）＊）ｘ_ｍ ^（ｋ）
（ｉ＝１，・・・，Ｎ^（ｋ）） ・・・（２１）

また、ａ_ｉ ^（ｋ）・ｙ^（ｋ）＊＝０となるｉに対するｘ_ｉ ^（ｋ）の値は、ａ_ｉ ^（ｋ）・ｙ^（ｋ）＊＝０となるベクトルａ_ｉ ^（ｋ）の本数をＮ^{（ｋ＋１）}として、Ｎ^{（ｋ＋１）}≦Ｍ−ｋ−１であれば、式（１１）と同様に、以下の式により一意に求められる（ステップＳＴ１０５）。
ｘ^{（ｋ＋１）}＝（（Ａ^{（ｋ＋１）}）^ＴＡ^{（ｋ＋１）}）^−１
（Ａ^{（ｋ＋１）}）^Ｔｆ^{（ｋ＋１）} ・・・（２２）
ここで、ａ_ｉ ^（ｋ）・ｙ^（ｋ）＊＝０となるａ_ｉ ^（ｋ）を列とするＭ×Ｎ^{（ｋ＋１）}行列をＡ^{（ｋ＋１）}、対応するｘ_ｉ ^（ｋ）をまとめてＮ^{（ｋ＋１）}次元ベクトルｘ^{（ｋ＋１）}とする。また、ｆ^{（ｋ＋１）}は式（８）と同様、式（２１）で求められたｘ^（ｋ）を用いて以下のように定義される。
ｆ^{（ｋ＋１）}＝ｆ^（ｋ）−Ａ^（ｋ）ｘ^（ｋ） ・・・（２３）

Ｎ^{（ｋ＋１）}≧Ｍ−ｋであれば、ｘ^{（ｋ＋１）}の値は一意ではなくなるため、更にk+1-部分問題を考え、その値を決定する（ステップＳＴ１０６）。最後に、上記ステップＳＴ１０５あるいはステップＳＴ１０６により求められたｘ^{（ｋ＋１）}の値を、ｘ^（ｋ）の対応する位置に代入し、全てのｘ^（ｋ）の値を得る（ステップＳＴ１０７）。

再び図４に戻って説明する。ステップＳＴ４０においてＮ^（１）≧Ｍであった場合には、1-部分問題の手続きにより、ｘ^（１）を求める（ステップＳＴ４２）。その際、上述のように、必要であれば、更に2-部分問題、3-部分問題と、k-部分問題をｋに関して必要な回数だけ逐次計算する（最大ｋ＝Ｍ−２まで）。それにより最終的にｘ^（１）が求められ、結局全てのｘの値が求められることになる。

尚、与えられたモデルにおいて、行列Ａの構成によってはこのようなk-部分問題を逐次解くことも必要ない場合がある。即ち、ａ_ｉのうちの任意のＭ本が１次独立となっていれば、最初の双対問題の最適解ｙ^＊に対して内積が０となるようなａ_ｉの本数は常にＭ−１本以下である。よってそのような場合にはk-部分問題（ｋ≧１）の計算は不要である。また、Ｍ≦２の場合にも当然、k-部分問題（ｋ≧１）の計算は不要である。

このように計算を行うことで、与えられたｆに対して常にｘの各要素絶対値の最大値が最小となることが保証されるようなｘを簡単な演算により求められる。それにより、各アクチュエータの出力を極力抑えつつ、所望の力、あるいはトルクを出力できるため、非常に効率的である。また、本装置によると、ｆの値が連続的に変化したときに、求められるｘの値も連続的に変化することが保証される。更に、これらの計算において、特にk-双対問題の最適解を有限個の候補の中から選択するだけで求められるため、収束計算が全く必要なくなり、計算上の信頼性が高い。

このような手法が必要となる例として、姿勢制御用のアクチュエータとしてホイールを備えた人工衛星の姿勢制御装置に適用した場合について、図面を用いて説明する。
図６は、ホイールによる人工衛星の姿勢制御装置の概要を示すブロック図である。
図６において、姿勢制御トルク生成部１１は、人工衛星が所望の姿勢駆動を行うための姿勢制御トルクコマンドを生成する機能部である。ホイールトルク分配生成部１２は、姿勢制御トルク生成部１１の出力である姿勢制御トルクコマンドからホイール１３を駆動するためのホイールトルクコマンドを生成する機能部である。ホイール１３は人工衛星本体１４の姿勢制御を行うためのアクチュエータである。以下、ホイールトルク分配生成部１２の働きについて説明する。

３軸の姿勢制御トルクを並べた３次元ベクトルをｆとし、ホイール３のＮ個のホイールの各駆動トルクを並べたＮ次元ベクトルをｘとすると、ｆとｘの関係は以下のように表される。

ｆ＝Ａｘ ・・・（２４）
ここで、Ａはホイール１３の幾何学的配置から一意に定まる３×Ｎ行列である。即ち、Ａの第ｉ列が、ｉ番目のホイールの回転軸方向ベクトルａ_ｉ（ｉ＝１，・・・，Ｎ）であるような行列である。式（２４）は式（１）において、Ｍ＝３の場合に相当する。ホイールの個数Ｎが４以上であれば、式（２４）を満たすホイール駆動トルクｘは一般には一意ではなくなる。通常、ホイールトルク分配生成部１２において、姿勢制御トルクｆからそれを実現するホイールトルクｘを求めるために、次式を用いることが多い。

ｘ＝Ｂｆ ・・・（２５）
ここで、ＢはＡＢ＝Ｉ_３を満たすＮ×３の定数行列であり、典型的にはＡの擬似逆行列、即ち、Ｂ＝Ａ^Ｔ（ＡＡ^Ｔ）^−１が用いられる。また、Ｉ_３は３×３の単位行列である。式（２５）によるトルク分配は計算が簡単なため、衛星に搭載のコンピュータで実現する上での問題はないが、ホイールのトルクを有効に利用するという点では効率がよいとはいえない。そこで、ホイールトルク分配生成部１２に、本発明に係る出力分配算出装置を適用する。その具体的な構成および動作については、上述の通りである。この場合、Ｍ＝３であるから、最大でも1-部分問題までを考えればよく、また任意の３本以上のホイール軸が同一平面内に含まれないようなホイール配置であれば、1-部分問題の計算は不要となる。また、双対問題（および1-双対問題）の最適解の候補Ｙは、ホイールの幾何学的配置から決まるベクトルの集合となる。

このように、ホイールによる人工衛星の姿勢制御装置に本発明に係る出力分配算出装置を適用することで、与えられた３軸姿勢制御トルクｆに対して、簡単な演算によりホイールトルク最大値を最小にするようなトルク分配ｘが求められ、効率的な姿勢駆動が実現できる。本手法においては収束計算が不要なため、衛星搭載のコンピュータにおいて軌道上で用いる上で、十分に信頼性が高い。

また、他の適用例として、本発明に係る出力分配算出装置を、４輪独立駆動車両の姿勢制御装置へ応用する場合について説明する。この種の問題の例が、坂井、佐渡、堀、“４輪独立駆動電気自動車における動的な制駆動力配分法”、電気学会論文誌Ｄ、ｖｏｌ．１２０−Ｄ、Ｎｏ．６、ｐｐ．７６１−７６８（２０００年）において述べられている。ここでは、駆動力と制動力をまとめて制駆動力と呼んでいる。

上記論文において述べられているように、電気自動車等においては分散配置されたモータにより、各駆動輪に独立にトルクを分配して、動的に車体姿勢の制御を行うことが可能である。

図７は、このような制御を行う際の車両制御系のブロック図である（上記論文７６２ページの図２に相当している）。
図中の、車両運動制御部１５において、車両の姿勢を制御するための総駆動力（制動力）およびヨーモーメントのコマンドを生成し、制駆動力分配生成部１６において、その総駆動力およびヨーモーメントを実現するための各車輪の駆動力（制動力）を求める。制駆動力制御装置１７は、制駆動力分配生成部１６の出力に従って各車輪に制駆動力を与え、車両本体１８の運動が実現される。以下、制駆動力分配生成部１６の働きについてより詳細に説明する。

各車輪の制駆動力をｘ_ｉ（ｉ＝１，・・・，４）とし、実現すべき車両全体の総駆動力をＦ_ｄ、ヨーモーメントをＭ_ｚとすると、これらは以下の式を満たす必要がある。
ｘ_１＋ｘ_２＋ｘ_３＋ｘ_４＝Ｆ_ｄ ・・・（２６）
−ｄ_ｆｘ_１＋ｄ_ｆｘ_２−ｄ_ｒｘ_３＋ｄ_ｒｘ_４＝Ｍ_ｚ・・・（２７）
ここで、ｄ_ｆ、ｄ_ｒはそれぞれ前輪、後輪における車輪間隔を２で割った値である。尚、車輪の番号は左前輪、右前輪、左後輪、右後輪の順としている。式（２６）、（２７）は上記論文中でも用いられている簡単化された車両のモデルである。車輪の数４に対して、制御したい量はＦ_ｄ、Ｍ_ｚの二つであるから、式（２６）、（２７）は冗長である。明らかに、ｆ＝［Ｆ_ｄ、Ｍ_ｚ］^Ｔ、ｘ＝［ｘ_１、ｘ_２、ｘ_３、ｘ_４］^Ｔとおけば、式（２６）、（２７）は式（１）でＭ＝２、Ｎ＝４の場合に相当する。

車輪の駆動力が大きくなると車輪の空転が起こりやすくなることから、各車輪の制駆動力の大きさ｜ｘ_ｉ｜は極力小さい値であることが好ましい。そこで、｜ｘ_ｉ｜の最大値を最小にするように各車輪の制駆動力の分配を求めることを考えると、やはりこの場合もミニマックス問題（１）、（２）の形にモデル化でき、本発明に係る出力分配算出装置を適用可能である。この場合、Ｍ＝２であるから、k-部分問題（ｋ≧１）は考える必要はなく、また、双対問題の最適解の候補は４本の定数ベクトルであるから、駆動力分配に要する計算量は非常に少ない。

このように、本発明に係る出力分配算出装置を４輪独立駆動車両の姿勢制御装置に用いることで、非常に簡単な演算により最適な制駆動力分配を求めることができるため、このような動的な車両姿勢の制御を実装する上での信頼性を高めつつ、性能を向上することができる。

尚、上記の例はＭ≦３の場合であるが、これに限定されるものではない。例えば、多関節のロボットアームの手先を物体表面に押し付けるような制御として、物体表面に加える力とトルクを全て所望の値にしたい場合が考えられる。このような場合、Ｍは、３次元空間における力とトルクの自由度を表すので、Ｍ＝３＋３＝６となる。Ｎはアームの各関節に取り付けられたモータの個数で、関節が多数あってＮ＞６である場合であれば本発明が同様に適用可能である。このような場合でも、特定の関節（モータ）に負荷が偏ることなく制御できるという上記の例と同様の効果を得ることができる。

以上のように、実施の形態１によれば、複数のアクチュエータの出力で制御対象の出力が実現され、前記複数のアクチュエータの個々の出力をＮ次元ベクトルｘで表し、制御対象の出力をＭ次元ベクトルｆで表した場合、ベクトルｘとベクトルｆが線形の関係式ｆ＝Ａｘによって関係付けられ、かつベクトルｘの次元Ｎがベクトルｆの次元Ｍよりも大きいモデルに対し、ベクトルｆが与えられたときにそれを実現するためのベクトルｘを求める出力分配算出装置において、行列Ａの第ｉ列のＭ次元ベクトルをａ_ｉ（ｉ＝１，・・・Ｎ）として、全てのベクトルａ_ｉ（ｉ＝１，・・・Ｎ）との内積の絶対値の総和がある正の定数ＣとなるようなＭ次元ベクトルのうちで、ベクトルｆとの内積の値が最大となるようなベクトルｙ^＊を求めるｙ^＊演算手段と、ベクトルｙ^＊とベクトルｆとの内積の値を定数Ｃで割った値ｘ_ｍを求めるｘ_ｍ演算手段と、ベクトルａ_ｉとベクトルｙ^＊との内積が０でないｉに対して、ｘ_ｉの値を、ｘ_ｍにａ_ｉとｙ^＊との内積の符号を付した値とする第１のｘ_ｉ演算手段と、ベクトルａ_ｉとベクトルｙ^＊との内積が０となるｉに対して、ｘ_ｉの値を、その絶対値がｘ_ｍを超えないように求める第２のｘ_ｉ演算手段とを備えたので、与えられた制御対象の出力ｆを実現するアクチュエータの出力ｘのうちで、その各要素絶対値の最大値が最小となるような分配を常に求めることができ、各アクチュエータの出力を極力抑えつつ、システム全体として所望の出力を得ることができる。

また、実施の形態１によれば、ｙ^＊演算手段を、ベクトルｙ^＊を求める場合、定数ベクトルの集合Ｙに含まれるベクトルが、Ｎ本のａ_ｉのうちのＭ−１本との内積が０であり、かつ、全てのａ_ｉ（ｉ＝１，・・・Ｎ）との内積の絶対値の総和がある正の定数Ｃとなるようなベクトルであることを制約条件として、ある本数の前記定数ベクトルの集合Ｙの中からベクトルｆとの内積の値が最大となるようなものを選択するよう構成したので、次のような効果がある。即ち、このように、ベクトルｙ^＊を求める際に、有限の本数の定数ベクトルの集合の中から選択することで、ベクトルｙ^＊を求めるための計算が非常に簡単になり、計算量を大幅に減少させることができる。また、収束計算が不要なため、計算の信頼性を向上させることができる。

また、実施の形態１によれば、第２のｘ_ｉ演算手段を、ａ_ｉとｙ^＊との内積が０となるｉに対してｘ_ｉの値を求める際、ｙ^＊との内積が０となるａ_ｉの数がＭ以上であった場合、上述したように定義されるk-部分問題の手続きを、ｋ＝１から順に、最大ｋ＝Ｍ−２まで逐次的に実行し、それらｘ_ｉの値を求めるよう構成したので、次のような効果がある。即ち、ａ_ｉとｙ^＊との内積が０となるｉに対するｘ_ｉの値を、このように求めることで、簡単な計算によりそれらの値がｘ_ｍを超えないように求めることができ、かつ、それらに対応したアクチュエータ出力を極力小さく抑えることができるため、効率向上を図ることができる。

実施の形態２．
図８は、この発明の実施の形態２による出力分配算出装置のブロック構成図である。
図の装置は、ｙ^＊演算手段１、ｘ_ｍ演算手段２、第１のｘ_ｉ演算手段３、第２のｘ_ｉ演算手段４、置き換え手段５からなる。ここで、ｙ^＊演算手段１〜第２のｘ_ｉ演算手段４は、図１に示した実施の形態１の構成と同様であるため、対応する部分に同一符号を付してその説明を省略する。

置き換え手段５は、ｘ_ｍ演算手段２で求められた値に対して、ベクトルｆをｆ（ｘ_ｌｉｍ／ｘ_ｍ）と置き換え、かつ、ｘ_ｍの値をｘ_ｌｉｍに置き換える機能部である。また、置き換え手段５においても、その機能に対応したソフトウェアとこれを実行するためのハードウェアから構成されている。

図９は、実施の形態２の出力分配算出装置の動作を説明するフローチャートである。
図９に示す動作は、図１のフローチャートに示す動作と一部が異なるのみであり、同一の動作をする部分については、同一ステップ符号を付してある。

本実施の形態２において、実施の形態１と異なる点は、ステップＳＴ２とステップＳＴ３の間に、新たにステップＳＴ２１を挿入した点である。ステップＳＴ２１では、置き換え手段５により、ステップＳＴ２で求められたｘ_ｍに対し、ｆをｆ＝（ｘ_ｌｉｍ／ｘ_ｍ）ｆと置き換え、更にｘ_ｍをｘ_ｍ＝ｘ_ｌｉｍと置き換える。ここで、ｘ_ｌｉｍは予め定められているアクチュエータ出力の許容値、あるいはそれより少し小さい値とする。その後の計算は実施の形態１と同様である。

このように、アクチュエータ出力の許容値に応じてｆおよびｘの値を与えることで、アクチュエータ出力許容値の範囲内で、システム全体としての出力が最も大きくなり、非常に効率的である。

例えば、実施の形態１で示した人工衛星の姿勢制御装置に対して、本実施の形態に係る出力分配算出装置を適用することが考えられる。その場合の姿勢制御装置全体のブロック図は図６と同様である。図６において、ホイールトルク分配生成部１２に本実施の形態の出力分配算出装置による計算手順を適用する。この場合、ｘ_ｌｉｍはホイールの最大出力可能トルクである。また、姿勢制御トルク生成部１１から出力される３軸姿勢制御トルクｆは、ホイールトルク分配生成部１２においてホイールの最大出力可能トルクに応じてその大きさを変更される。そして、出力されるホイールトルクコマンドは、ホイールの最大出力可能トルクの範囲内で姿勢駆動トルクの大きさが最大となるようなものとなる。

このように、本実施の形態に係るトルク分配方法を人工衛星の姿勢制御装置に適用することで、高速な姿勢駆動や大角度の姿勢変更が必要となる場合に、ホイールの能力の範囲内で全体として最大限の姿勢制御トルクを出力でき、人工衛星の機動性を大幅に向上できる。

以上のように、実施の形態２によれば、実施の形態１の構成に加えて、ベクトルｆをｆ（ｘ_ｌｉｍ／ｘ_ｍ）と置き換え、かつ、ｘ_ｍの値をｘ_ｌｉｍに置き換える置き換え手段を設けたので、アクチュエータの出力が限られているような場合でも、その範囲内でシステム全体としての最も大きな出力を算出することができる。

この発明の実施の形態１による出力分配算出装置を示す構成図である。 この発明の実施の形態１による出力分配算出装置における出力分配方法を示すフローチャートである。 実施の形態１に係る図２のフローチャートのステップＳＴ１のより詳細な動作を説明するフローチャートである。 実施の形態１に係る図２のフローチャートのステップＳＴ４のより詳細な動作を説明するフローチャートである。 k-部分問題の解ｘ^（ｋ）を求める手順を示すフローチャートである。 実施の形態１を適用する人工衛星の姿勢制御装置の概要を示すブロック図である。 実施の形態１を適用する車両制御系のブロック図である。 この発明の実施の形態２による出力分配算出装置のブロック構成図である。 実施の形態２の出力分配算出装置の動作を説明するフローチャートである

符号の説明

１ ｙ^＊演算手段、２ ｘ_ｍ演算手段、３ 第１のｘ_ｉ演算手段、４ 第２のｘ_ｉ演算手段、５ 置き換え手段。

Claims (4)

1. 複数のアクチュエータの出力で制御対象の出力が実現され、前記複数のアクチュエータの個々の出力をＮ次元ベクトルｘで表し、前記制御対象の出力をＭ次元ベクトルｆで表した場合、前記ベクトルｘと前記ベクトルｆが線形の関係式ｆ＝Ａｘによって関係付けられ、かつベクトルｘの次元Ｎがベクトルｆの次元Ｍよりも大きいモデルに対し、ベクトルｆが与えられたときにそれを実現するためのベクトルｘを求める出力分配算出装置において、
   行列Aの第ｉ列のＭ次元ベクトルをａ_ｉ（ｉ＝１，・・・Ｎ）として、全てのベクトルａ_ｉ（ｉ＝１，・・・Ｎ）との内積の絶対値の総和がある正の定数ＣとなるようなＭ次元ベクトルのうちで、ベクトルｆとの内積の値が最大となるようなベクトルｙ^＊を求めるｙ^＊演算手段と、
   ベクトルｙ^＊とベクトルｆとの内積の値を前記定数Cで割った値ｘ_ｍを求めるｘ_ｍ演算手段と、
   ベクトルａ_ｉとベクトルｙ^＊との内積が０でないｉに対して、ｘ_ｉの値を、ｘ_ｍにａ_ｉとｙ^＊との内積の符号を付した値とする第１のｘ_ｉ演算手段と、
   ベクトルａ_ｉとベクトルｙ^＊との内積が０となるｉに対して、ｘ_ｉの値を、その絶対値がｘ_ｍを超えないように求める第２のｘ_ｉ演算手段とを備えた出力分配算出装置。
2. 複数のアクチュエータの出力で制御対象出力が実現され、前記複数のアクチュエータの個々の出力をＮ次元ベクトルｘで表し、前記制御対象の出力をＭ次元ベクトルｆで表した場合、前記ベクトルｘと前記ベクトルｆが線形の関係式ｆ＝Axによって関係付けられ、かつベクトルｘの次元Ｎがベクトルｆの次元Ｍよりも大きいモデルに対し、ベクトルｆが与えられたときにそれを実現するためのベクトルｘを求める出力分配算出装置において、
   行列Ａの第ｉ列のＭ次元ベクトルをａ_ｉ（ｉ＝１，・・・Ｎ）として、全てのベクトルａ_ｉ（ｉ＝１，・・・Ｎ）との内積の絶対値の総和がある正の定数ＣとなるようなＭ次元ベクトルのうちで、ベクトルｆとの内積の値が最大となるようなベクトルｙ^＊を求めるｙ^＊演算手段と、
   ベクトルｙ^＊とベクトルｆとの内積の値を前記定数Cで割った値ｘ_ｍを求めるｘ_ｍ演算手段と、
   ベクトルｆをｆ（ｘ_ｌｉｍ／ｘ_ｍ）と置き換え、かつ、ｘ_ｍの値をｘ_ｌｉｍに置き換える置き換え手段と、
   ベクトルａ_ｉとベクトルｙ^＊との内積が０でないｉに対して、ｘ_ｉの値を、ｘ_ｍにａ_ｉとｙ^＊との内積の符号を付した値とする第１のｘ_ｉ演算手段と、
   ベクトルａ_ｉとベクトルｙ^＊との内積が０となるｉに対して、ｘ_ｉの値を、その絶対値がｘ_ｍを超えないように求める第２のｘ_ｉ演算手段とを備えた出力分配算出装置。
3. ｙ^＊演算手段は、
   ベクトルｙ^＊を求める場合、定数ベクトルの集合Ｙに含まれるベクトルが、Ｎ本のａ_ｉのうちの少なくともＭ−１本との内積が０であり、かつ、全てのａ_ｉ（ｉ＝１，・・・Ｎ）との内積の絶対値の総和がある正の定数Ｃとなるようなベクトルであることを制約条件として、ある本数の前記定数ベクトルの集合Ｙの中からベクトルｆとの内積の値が最大となるようなものを選択するよう構成されていることを特徴とする請求項１または請求項２記載の出力分配算出装置。
4. 第２のｘ_ｉ演算手段は、
   ａ_ｉとｙ^＊との内積が０となるｉに対してｘ_ｉの値を求める際、ｙ^＊との内積が０となるａ_ｉの数がＭ以上であった場合、次のように定義されるk-部分問題の手続きを、ｋ＝１から順に、最大ｋ＝Ｍ−２まで逐次的に実行し、それらｘ_ｉの値を求めるよう構成されたことを特徴とする請求項１または請求項２記載の出力分配算出装置。
   k-部分問題：
   （１）全てのａ_ｉ ^（ｋ）（ｉ＝１，・・・Ｎ^（ｋ））との内積の絶対値の総和がある正の定数Ｃ^（ｋ）となり、かつ全てのｙ^（ｊ）＊（ｊ＝０，・・・ｋ−１）との内積が０となるようなベクトルのうちで、ｆ^（ｋ）との内積の値が最大となるようなＭ次元ベクトルｙ^（ｋ）＊を求め、
   （２）ｙ^（ｋ）＊とｆ^（ｋ）との内積の値をC^（ｋ）で割った値ｘ_ｍ ^（ｋ）を求め、
   （３）ａ_ｉ ^（ｋ）とｙ^（ｋ）＊との内積が０でないｉに対して、ｘ_ｉ ^（ｋ）の値を、ｘ_ｍ ^（ｋ）にａ_ｉ ^（ｋ）とｙ^（ｋ）＊との内積の符号を付した値とし、
   （４）ａ_ｉ ^（ｋ）とｙ^（ｋ）＊との内積が０となるａ_ｉ ^（ｋ）がＭ−ｋ−１本以下であった場合には、ｘ_ｉ ^{（ｋ＋１）}の値を、一意的な関係式Ａ^{（ｋ＋１）}ｘ^{（ｋ＋１）}＝ｆ^{（ｋ＋１）}から求め、
   （５）ａ_ｉ ^（ｋ）とｙ^（ｋ）＊との内積が０となるａ_ｉ ^（ｋ）がＭ−ｋ本以上であった場合には、更にｋ＋1-部分問題の手続きによりｘ^{（ｋ＋１）}を求め、
   （６）ｘ^{（ｋ＋１）}をｘ^（ｋ）の対応する位置に代入する。
   ここで、形式的にＮ^（０）＝Ｎ、ａ_ｉ ^（０）＝ａ_ｉ（ｉ＝１，・・・Ｎ）、Ａ^（０）＝Ａ、ｆ^（０）＝ｆ、ｙ^（０）＊＝ｙ^＊、ｘ^（０）＝ｘと表したときに、ｋ≧１に対して、
   全てのｙ^（ｊ）＊（ｊ＝０，・・・ｋ−１）との内積が０となるａ_ｉ ^{（ｋ−１）}の数をＮ^（ｋ）、およびそれらを列とするＭ×Ｎ^（ｋ）行列をＡ^（ｋ）、およびＡ^（ｋ）の第ｉ列をａ_ｉ ^（ｋ）（ｉ＝１，・・・Ｎ^（ｋ））、および対応するｘ_ｉ ^{（ｋ−１）}をまとめたＮ^（ｋ）次元ベクトルをｘ^（ｋ）、およびそのｉ番目の要素をｘ_ｉ ^（ｋ）（ｉ＝１，・・・Ｎ^（ｋ））とする。また、k-部分問題において現れるｆ^{（ｋ＋１）}は、その時点で値の与えられていないｘ_ｉ ^（ｋ）を一旦０とおいたときに、ｆ^{（ｋ＋１）}＝ｆ^（ｋ）−Ａ^（ｋ）ｘ^（ｋ）と与える。

Images