JP2004280792A - リアルタイムモデル予測制御における二次のプログラミングに対する加速されたアクティブセット探索のシステムおよび方法 - Google Patents

Abstract

【課題】 リアルタイムに動的システムを制御するためのアルゴリズムおよび装置を提供すること。
【解決手段】 アクティブセットアルゴリズムは、問題解決への急速な進歩を達成するために、拘束している制約の組に対する「ホットスタート」を、有効な線形代数を伴って、最適に活用する。該線形代数は、要求された因数分解が実行されるとかつ退歩が現れると、最も必要のない予備処理ステップを介さずに退歩制約を処理するように設計される。これらの新規なアプローチは、共に組み合わされることによって、リアルタイムでの制御問題の問題解決を可能にする。
【選択図】 図１

Description

この発明は、米国政府契約N00421-01-2-0131の下での作業の実行において理解される。

本発明は、オンボード制御に対するオンボード最適化技術に関し、かつ、特に、動的システムのモデル予測制御(Model Predictive Control)に対して有用である。

この発明は、オンボード制御に対する既存の最適化技術に対して、重大な増強を記述する。非線形動的システムのオンボードモデル予測制御は、非線形力学を線形化することと、二次のプログラム（Quadratic Program）が所望のプロフィールの出力の近くにとどまることを引き起こすこととをしばしば含む。オンボード制御は、リアルタイムでエラーに強い方法で信頼性のある問題解決（solution）を実行することを要求するので、二次のプログラミングアルゴリズムは、線形代数を最適に使用することを必要とする。

モデル予測制御は、「プラント」のモデルまたは動的システムに基づいて動的プロセスの最適なオペレーティングパラメータを決定する手続きを参照する。このプラントモデルは、ガスタービンエンジンからヘリコプター内の音響学までにわたる全ての物理学ベースのモデルであり得る。システムの物理的な制約を尊重しながら、そのようなプラントを最適に操作すること、則ち、操作中にある一定の目的を満たしまたは超えることは、エンジニアにとって興味深い。この目的のために、プラントの操作中に制約された最適化問題を解決すること、および、システムが時と共に発展しまたは将来の要求の予測が変化すると最適化プログラムのパラメータを更新すること、および、該問題を再解決することは、一般的なことである。

この手続きにおける重大な困難性は、複雑な最適化問題に対する理にかなった問題解決をリアルタイムで取得することができる必要があることである。このことは、この発明によって言及される論点である。

従来は、リアルタイム制御は、化学工学用途において、特に、リアルタイムプロセス制御において試みられていた。それはまた、航空宇宙用途におけるリアルタイム軌道計画において適用されてきた。しかしながら、そこで生じる最適化問題は、数分、しばしば数時間で解決される必要があった。ここで関心がある問題によって許される「リアルタイム」スケールは、ミリ秒のオーダーに基づく。ここで記述される新奇なアルゴリズムは、問題解決の忠実度を犠牲にすることなくこの時間要求に取り組むことができ、故に、既存の方法に対する重大な増強である。

同時継続の出願である米国出願第10/308,285号（２００２年１２月２日出願、題名“Real-Time Quadratic Programming for Control of Dynamical Systems”）は、共通に譲渡され、かつ、発明者のうちの１人はまた、本願の発明者である。本発明が主に取り組む問題は、たった１つのタイムステップにわたりかつ計算に対して非常に小さな予算を伴う「動的反転（Dynamic Inversion ）」制御問題であり、故に、アクティブセットホットスタート（active set hot start）を実行する際に、異なる戦略が使用された。更に、異なる手続きが、アクティブセットに対する探索の間に制約を取り除くために使用された。その結果、異なる線形代数が要求された。全く異なるが、この発明において記述される方法はまた（それは、動的反転に対するより良い方法ではないかも知れないが）動的反転に対する二次のプログラミングアルゴリズムを解決するための代わりの方法であり得る。

本発明は、リアルタイムに動的システムを制御するためのアルゴリズムおよび装置を提供する。そのような動的システムは、ガスタービンエンジンと航空機飛行制御システムとを含むが、これらに限定されない。ここで記述される最適のアクティブセットまたは拘束する制約の組（set of binding constraints）に対する探索のメカニズムは、如何なる凸状の（convex）最適化問題にも適用されることができる。これは、凸状の二次のプログラミングを含む。

本発明は「アクティブセットホットスタート」、則ち、最適に拘束するために推測される制約の組を提供しかつ活用する。これは、既存の方法とは異なる。該既存の方法は、問題解決に対する開始推測（starting guess）を使用し、かつ、アクティブセットを使用しない。

本発明は、一致する（consistent）アクティブセット方法を提供する。拘束していない制約は、乗数の符号に関係なく、アクティブセットから落とされる（dropped from the active set）。これは、目的機能を単調に改善する反復を結果として生じる。

本発明は、また、一致しない（inconsistent）アクティブセット方法を提供する。アルゴリズムの反復は、アクティブセットと一致しない。これは、初期ホットスタート推測における大抵の制約が最適のアクティブセットにおいて終了し故にたとえ不一致でも落とされないことを予期されるからである。反復は目的を単調に改善しないかも知れないが、このことは、開始推測が中程度に良くても、より速い収束とより良いパフォーマンスとをリアルタイムに結果として生じる。

本発明は、現在の因数分解におけるステップを使用して先行する予備処理無しに反復におけるアクティブ制約の組における線形依存を検出するためのメカニズムを提供する。本発明は、また、初期実行可能点の情報無しに実行可能かつ最適な問題解決を同時に達成するためにＭ≧０かつＫ＞０の値に対して二次のプログラムの「ビッグＫ（big-K ）」公式化を提供する。上記技法の何れかまたは全ては、最適化の希薄な（sparse）バージョンに適用されることができる。

本発明の他の利点は、それが、添付図面に関連して考慮されるとき、以下の詳細な記述に対する参照によってより良く理解されると共に、直ちに認識されるであろう。

図面は、本発明の二次のプログラミング方法を使用する１つの型の制御システムを図解する。

〔一般的な問題〕
図１は、モデル予測制御を使用する制御システム１０の一般的なモデルであり、かつ、本発明から利益を受けるタイプの一般的なモデルである。制御システム１０は、所望の軌道発生器１２を具備する。軌道発生器１２は、システム１０の所望のプロフィールの出力を生成する。線形化モジュール１４は、所望の軌道に関する線形化されたモデルを、所望の軌道発生器１２から得る。二次のプログラミング公式化モジュール１６は、如何なる制約をも尊重しながら所望の軌道を最良に達成するための制御プロフィールを決定するために、二次のプログラムを公式化する。二次のプログラミング解決器１８は、最適な制御のプロフィールを生成するために、公式化モジュール１６によって確立された最適化問題を解決する。二次のプログラミング解決器１８は、この発明の焦点である。最適な制御のプロフィールは、作動システム２０に送られる。作動システム２０は、システム１０のプラント２２上で動作する。センサシステム２４は、プラント２２から所望の軌道発生器２２へ、フィードバックを提供する。

ＭＰＣに対する最適化問題の以下の公式は、完全性のためだけに、ここに含まれる。制御変数ｕと状態変数ｘと応答（または出力）ｙとを伴う非線形動的システムを考慮すると、これらは以下の通りに関係付けられる。

一様な時間間隔を伴う上記の離散時間バージョンは、以下の通りに記述される。

非線形関数φおよびｈは、一般に、ベースポイントに関して線形化される。ベースポイントは、安定した状態ポイント、則ち、ｘが消滅するポイントである。そのような安定した状態ベースポイントｘ_ｓ ，ｕ_ｓ 、則ち、φ（ｘ_ｓ ，ｕ_ｓ ）＝０であるポイントが与えられると、離散時間システムは、以下の通りに線形化されることができる。

ここで、φ_ｘ はｘに関するφのヤコビアンを示し、他も同様である。

制御エンジニアは、一般的に、上記システムを以下の通りに表現する。

ここでは、以下の式が成り立つ。

線形化された離散時間システムを定義する上記パラメータの時間依存は、暗黙のうちに（ｘ_ｓ ，ｕ_ｓ ）内に隠される。それは、線形化が実行される点であり、それは、各時間点に対して、新たに選択される。「モデルを変更するこの準線形パラメータ（ｑ−ＬＰＶ）において、線形化のこの点は、いくつかの安定した状態点の凸状の組合せであり得る」ということに注意されたい。ｑ−ＬＰＶモデルは、真の非線形モデルの代わりにはならないが、フィードバック制御枠組において要求される正確さのレベルに対して、しばしば十分である。もし真の非線形モデルが利用可能になるならば、ここに記述される技法は、真の線形化から取得される二次のプログラムに対してまさに適用可能である。

〔目的〕 システムの上記記述が与えられると、目的は、特定の時間窓にわたってシステムを最適に制御することである。それを行う際における成功の程度は「システムの変数が、ある一定の参照軌道に対して、どれだけ近接して追跡するか」ということによって測定される。もし

がｙ_ｔ ，ｘ_ｔ ，ｕ_ｔ に対する参照軌道を示すならば、所望の目的は、以下の通りに表現されることができる。

重み付け行列のダイアゴナル（diagonals ）は、種々の目的の相対的な重みを示す。我々は、通常、たった１つまたは２つの主たる目的を有する。一方、残りは、二次的な目的である。二次的な目的に関する重みは小さな数に設定されるが、それらは、常に非ゼロであり、かつ、問題のヘッシアン（Hessian ）が非常に不十分に条件付けられるわけでないほど十分に大きい。我々は、一般的に、倍精度計算のために、該条件数を１０^７ より下に保持するように試みる。「制約が無いときでさえも、最適制御に対して唯一の問題解決が存在する」ということを保証することもまた必要である。もしそうでないならば、システムが安定した状態であるが、未だ最適制御が方々に跳び回り、エンジニアが「チャタリング」と呼ぶものを引き起こす場合を想像されたい。問題公式化を計算し直すための理由は、凸状の二次のプログラムを解決することを正当化するためであり、かつ、この特許の下での公式化をカバーするためではない。

〔制約〕 拘束された制約は、一般的に、ｘ，ｙ，ｕに課せられる。これは、その変化の比率に課せられるのと同様である。加えて、これらの変数の組合せを含む他の線形不等制約（linear inequality constraint）が存在する。故に、該不等制約は、以下の通りに表現されることができる。
Ａ_０Ｕ≦ｂ_０
ここで、Ｕは、以下に示される通りであり、最適化変数のベクトルを示す。

「ｘ_１ は最適化変数ではない。なぜならば、初期状態および先行する時間からの制御変数が、第１の時間点での状態ｘ_１ を決定するからである」ということに注意されたい。

〔二次のプログラム〕 我々は、上記のことを、厳密な凸状の二次のプログラム（ＱＰ）として、以下の式において示すことができる。

同等および不等の制約を伴う上記ＱＰは、大きくかつ散在（aparse）する。リアルタイム制御に対するアクティブセットアルゴリズムを伴う経験は、より少ない変数を伴う緻密な（dense）公式化を考慮することの利点を示す。以下では、我々は、元のＱＰに対して上記緻密な公式化を取得する方法の概略を述べる。

Ｚを行列とする。その列は、Ｅ_０ のヌル（null）空間に対して、基底を形成する。Ｅ_ｕＵ_ｆ＝ｒ_０を満足する何れかの特定のＵ_ｆ が与えられると、同等制約に対する問題解決のファミリは、以下の通りに表現されることができる。
Ｕ＝Ｕ_ｆ＋Ｚ_ｐ （４）
ここで、

であり、かつ、ｍはＥ_０のヌル空間の次元（dimension）である。

元の問題において（４）を置換することは、ｐにおいて、以下の縮小された緻密なＱＰを生ずる。

次のセクションは、（３）における同等制約を定義する状態式（１）および（２）を活用することによって、かつ、実際の因数分解を実行することなく、行列Ｚを計算するための方法を記述する。

〔ヌル空間に対する基底を計算すること〕 関心のあるヌル空間は、Ｅ_０Ｕ＝０を満足する全てのＵを含む。

変数ＵをＵ_ＢとＵ_ＮＢとに分割することを考える。Ｕ_Ｂに対応するＥ_０の列に部分行列Ｅ_Ｂ形成させ、Ｕ_ＮＢに対応するそれらに部分行列Ｅ_ＮＢ形成させる。我々は「Ｅ_Ｂが正方でありかつ正則であるような方法で分割が行われる」ということを仮定する。そして、ヌル空間における全てのＵに対して、以下の式が成り立つ。

故に、Ｅ_０のヌル空間における全てのベクトルＵは、以下の通りに記述されることができる。

「行列

が、関心のあるヌル空間に対する基底を形成する列を含み、かつ、ベクトルＵ_ＮＢがｐの役割を仮定する」ということを論ずることは困難でない。例えばＱＲ因数分解を使用して計算されたヌル空間と対照的に、その一部は、希薄なアイデンティティブロックである。一層重要なこととして、もし次に示すように賢明に分割する我々の変数を我々が選択するならば、Ｚの残りを計算するために行列Ｅ_Ｂを反転することは決して必要ではない。

時間と共に積み重ねられる制御変数のベクトルとしてＵ_ＮＢを選択する。則ち、以下に示される通りである。

そして、Ｕ_Ｂは、時間と共に積み重ねられる出力および状態変数の対応するベクトルになる。

を調べることによって、「

のｉ番目の列がＵ_Ｂの値によって与えられ、Ｕ_ＮＢ＝ｅ_ｉであるとき、ベクトルはそのｉ番目の位置に１を伴い、他の場合はゼロである」ということが明白である。これは、状態の値を計算することを減少させ、かつ、全ての設定に対してＵ_ＮＢのｅ_ｉを出力する。これは、１つの特定の時間点における制御変数を１に設定することによって、かつ、残りを０に設定することによって、状態式（１）および（２）を通じてシミュレーションを実行することに等しい。故に、実際の因数分解は必要ない。

これは、制御理論において標準基底（standard basis）として知られていることを計算することに本質的に等しい。しかしながら、我々は、それを、分離数学的観点から取得した。

〔縮小された緻密な問題〕 先のサブセクションにおけるＺとｐ＝Ｕ_ＮＢ＝ｕとの選択に基づいて、（５）におけるＱＰは、以下の式を仮定する。

ここで、Ｈ＝Ｚ^ＴＨ_０Ｚ、ｃ＝Ｚ^Ｔ（ｃ_０＋Ｈ_０Ｕ_ｆ）、Ａ＝Ａ_０Ｚ、ｂ＝ｂ_０−Ａ_０Ｕ_ｆである。

本論文の残りにおいて我々がアルゴリズムを開発したものは、不等制約を伴うこの縮小された緻密かつ厳密な凸状のＱＰである。「同様に、ＱＰアルゴリズムにおける反復と基礎的な数学とが、希薄な問題に等しく適用可能である」ということに注目されたい。

〔最適なアクティブセットを探索するためのアルゴリズム〕 アクティブセットアルゴリズムは、拘束している制約の組を、最適に、繰り返し探索する。この組は、通常、アクティブセットとして参照される。モデル予測制御問題に対しては、現在の時間ステップにおけるＭＰＣ問題に対する問題解決が、次の時間ステップにおける最適なアクティブセットに対する推測を提供する。この個々の推測は、ＱＰを収束に向かって解決するために要求される反復の数を縮小するための大きな利点に対して使用され、リアルタイム制御において事実上欠くことのできない特徴である。アクティブセットアルゴリズムの広範なステップは以下に概略を述べられ、かつ、関連する線形代数は以下のセクションにおいて詳述される。

・アクティブセットに対する推測を伴って開始する。我々は、アクティブセットにおける制約のインデックスの組をＷによって示し、また、推測されたアクティブセットに対応する制約行列Ａにおける行を、Ｅによって特徴付ける。「実行可能な点ｕ_ｆが知られている、則ち、Ａｕ_ｆ≦ｂである」ということを仮定する。

・反復ｋにおいて、アクティブ制約に対する推測Ｅが与えられると、同等制約されたＱＰ（Equality-Constrained QP：ＥＱＰ）を解く。

ここで、ｒは、アクティブ制約に対する右手側のサブベクトルを示す。最適なｕ^＊とラグランジュ乗数λ^＊とは、以下の式によって与えられる。

比率テスト（Ratio Test）は、ｕ^ｋ−１＋αｓが実行可能であるような最大のα∈〔０，１〕を決定する。ここで、ｓ＝ｕ^＊−ｕ^ｋ−１ であり、かつ、ｕ^ｋ−１は先行する反復を示す。ｕ^０＝ｕ_ｆであることに注目されたい。言い換えると、Ａｕ^ｋ−１＋αＡｓ≦ｂであるように、〔０，１〕における最も大きな値であるために、αがピックアップされる。このテストは、アクティブセット内にある制約インデックスｉに対して実行される必要がない。なぜならば、それらは、Ａｓ_ｉ＝０（または、反復１においてＡｓ_ｉ≦０）を有しており、かつ、先行する反復ｕ^ｋ−１が実行可能であるからである（脚注：Ａｓ_ｉ＝０を有するアクティブセット内にない退歩した制約（degenerate constarint）を有することは可能であり、かつ、そのような制約に対する比率テストはまたスキップされる）。故に、αは、以下の式によって与えられる。

・次の反復においてアクティブセットに対する推測を更新する。

「どのαが比率テストにおいてその最小値を達成するのか」ということに対応する第１のインデックスｉ_Ｂ をＷに加える。アクティブに対する推測が次の反復において指し示す。対応する制約は「最も厳しい制約」としてゆるく参照される。

もしＥＱＰ（７）内のラグランジュ乗数λ^＊がネガティブである制約が存在するならば、最もネガティブな乗数を伴う制約をアクティブセットから落とす。これは、所謂、最も切り立った端規則（steepest edge rule）である。ブランドの規則（Bland's Rule）は、退歩（degeneracy）がアクティブセットにおいて検出されるときはいつでも展開され、かつ、これらの例では、ネガティブ乗数を伴う第１の制約が落とされる。退歩の処理および検出に関する詳細は、この明細書において、後に現れる。

〔第１の反復に対する特別な場合〕 初期アクティブセットにおける制約を、もし制約がｕ_ｆにおいて拘束しているならば、一致するように定義する。もし第１の反復においてα＜１であるならば、一致しない全ての制約が、それらの乗数の符号に関係なく、落とされる。その理由は、もしいくつかの要素に対する厳密な不等を伴ってＥｕ_ｆ≦ｒであるならば、更新された反復ｕ^１＝ｕ_ｆ＋α（ｕ^＊−ｕ_ｆ）は、次の反復に対するアクティブセットにおける全ての制約に対してＥｕ^１＝ｒを満足せず、矛盾に至るからである。これは、一致するアクティブセット方法である。例えば、もしｕ_ｆが実行可能な組の厳密な内部（strict interior）にあるならば、これは、全ての制約を落としかつ空のアクティブセットを伴って再スタートすることに至り、故に「ホットスタート」のわずかな活用を許可することに至る。しかしながら、もしｕ_ｆが先行する時間ステップにおいてＭＰＣ問題を解くことによって取得された最適な制御の組であるならば、それは、通常、アクティブセットに対する初期の推測と一致する。

〔不一致なアクティブセット方法〕 いくつかの場合において、既知の実行可能な点は、アクティブセットに対する開始推測と主として一致しない。このことは、多数の不一致の制約が第１反復の終わりにアクティブセットから落とされることに至る。もし初期ホットスタート推測が良いものであるならば、これらの制約のバックアップを選び取る際に多数の反復が浪費される。故に、他のアプローチは、アクティブセットにおける不一致の制約を保有し、かつ、たとえ反復およびアクティブセットが一致しなくとも、最適条件に対する探索を続けることである。もし如何なる反復においてもα＝１であるならば、この不一致は消失する。その上、不一致の程度は、全ての反復において縮小する。この不一致方法は、通常、一致方法よりも急速な最適条件への収束を可能にする。

反復を更新する。
ｕ^ｋ＝ｕ^ｋ−１＋αｓ

α＝１でありかつλ^＊≧０であるとき、則ち、アクティブセットが変化しないとき、全体的な最適条件への収束が達成される。

もし収束が達成される前にリアルタイム間隔内の許可された時間が使い果たされるならば、最後の更新された反復が問題解決として戻される。

〔反復を単調に改善すること〕 「一致するアクティブセットアルゴリズムにおける反復が（非ゼロαに対する厳密な改善を伴って）二次の目的を単調に改善する」ということは、凸状（convexity）を使用することによって証明される。

〔線形代数の詳細〕
このセクションでは、我々は、退歩の場合に我々を保護するために必要な、かつ、できるだけ少ない計算を伴ってＥＱＰへの正確な問題解決を取得するために必要な線形代数計算を記述する。

我々の問題においてＨは厳密にポジティブに限定されるので、ｕは（８）において除去されることができ、かつ、問題解決ＥＱＰ（７）は、以下の通りの閉じられた形式で記述されることができる。

我々は、開始時にＨ（Ｈ＝ＬＬ^Ｔ、Ｌはより低い三角形（lower triangular）である）のコレスキー因数分解（Cholesky factorization）を実行でき、かつ、それを全ての反復において再使用できるので、これは、全ての反復においてより大きなシステム（８）を解くよりも手間がかからない。その上、「二次の目的への制約されない問題解決ｕ_ｕｎｃ＝−Ｈ^−１ｃがまた、開始時にコレスキー因数分解を使用して計算されることができ、かつ、（９）および（１０）の計算における全ての反復において再使用されることができる」ということに注目されたい。

〔退歩〕
物理的な等価物が他の制約によって既に示されている冗長な制約を導入することによって、退歩は引き起こされる。それは、他の制約の線形的な組合せとして交互に表現される。しばしば、そのような従属の制約は、同時にアクティブになることがあり、かつ、作業中の組において終わることがある。このことは、アクティブ制約の行列Ｅにおいて行ランクの損失（loss of row rank）を結果として生じ、（９）毎にλ^＊を計算する際における係数行列ＥＨ^−１Ｅ^Ｔに対して相互のゼロまたは非ゼロ条件数に至る。

標準の慣例／論文における示唆されたアプローチは「冗長な制約が、問題公式化において現れない」ということを確実にすることである。これは可能でないかも知れない。なぜならば、我々は「『作動中に』制約が退歩するかも知れない」という可能性を認めなければならないからである。「制約に違反しないことに関するエラー強さを保証するために、多数のセンサ読み出しに基づく殆ど冗長な制約がしばしば課せられる」ということが与えられるならば、これは更にそうである。

大抵の商業的な符号において採用されている第２の保護手段は、予備処理の保護手段であり、または、制約行列の先行する検査によって冗長な制約を除去することである。これは、ガウス除去（Gaussian Elimination）に類似したステップを含み、かつ、退歩が論点でない最も有力な大多数の場合において浪費される特別な計算である。

これに対して、我々は、アクティブセットにおける退歩を、それが現れると、ＥＱＰをとにかく解くために必要とするもの以外の如何なる特別な線形代数をも行うことなく、検出しかつ処理することを目的とする。発明Accelerated Quadratic Programming for On-board Real-Time Controlでは、これは、（９）を決定する際にＥＨ^−１Ｅ^ＴのＱＲ因数分解を実行することによって、かつ、Ｒにおける近似ゼロ軸（near-zero pivot）に対応するλ^＊の要素をゼロに設定することによって達成される。これは「冗長な制約が最適化における役割を果たさない」という事実と一致する。

ここで、我々は「アクティブセットがより伝統的に変化する」という事実を利用する。則ち、我々は、１つの拘束している制約をピックアップし、かつ、ネガティブ乗数を伴う全ての制約を落とすこととは対照的にたった１つの制約のみを落とす。以下に示される計算は「我々が、因数分解の最中に退歩を検出しながら、関連する因数分解は勿論のこと、積ＥＨ^−１Ｅ^Ｔをも有効に更新する」ということを可能にする。

Ｌが与えられると、Ｈのコレスキー係数は、ＥＨ^−１Ｅ^Ｔ＝Ｍ^ＴＭである。ここで、Ｍ＝Ｌ^−１Ｅ^Ｔである。

Ｍ^ＴＭのＱＲ因数分解を実行する代わりに、我々は、Ｍ自身のＱＲ因数分解を実行する。その利点は「これは、Ｍ^ＴＭに関するＱＲの他の選択枝と同様に、Ｍ^ＴＭの条件数（condition number）の平方根のみである条件数を有するＭに関する因数分解を実行しながら、退歩によるＥにおける如何なるランク損失をも示す」ということである。
ＭのＱＲ分解に対して必要とされる一連のハウスホルダー変換（Householder transformation）Ｑ_ｉを計算することを考える。
Ｑ^{（ｍ−１）}...Ｑ^（２）Ｑ^（１）Ｍ＝Ｒ

ｎａをアクティブ制約の数、則ち、ＥまたはＭにおける行の数であるとする。Ｍの列ｊに関して作動するｊ番目のハウスホルダー変換は、列ｊ、行ｊ＋１〜ｎａにおける全ての要素をゼロにする。故に、上部の三角行列（upper triangular matrix）Ｒが形成される。ｖを列ベクトルとする。該列ベクトルにおいて、第１のｊ−１要素は０であり、残りにおいて、Ｑ^{（ｊ−１）}...Ｑ^（１） Ｍ（中間の縮小された行列）の列ｊの要素に釣り合う。そして、このハウスホルダー変換は、以下の通りに表現されることができる。

もしＥの列ｊが退歩制約（degenerate constraint）を示し、かつ、Ｅの列１〜ｊ−１の線形組合せとして形成されることができるならば、Ｍの列ｊについても、同じことが正しい。そして、Ｍに関する上記ハウスホルダー変換Ｑ^{（ｊ−１）}〜Ｑ^（１）は、中間の縮小された行列の列ｊにおける要素ｊ〜ｎａを既にゼロにした。そして、そのサブ列（sub-column）に基づくハウスホルダー変換に対するベクトルｖは全てゼロである。それによって、線形的な従属制約の存在を示す。我々がそのような退歩を検出するときはいつも、我々は、Ｍの対応する列をスキップし、かつ、対応する

を０に設定し、制約を退歩として有効に落とす。そして、我々は、非ゼロｖを提供する次の列を発見し、かつ、対応するハウスホルダー変換を計算する。結果としてのハウスホルダー変換の積は直交であるので、結果としてのＲは、以下の式を満足する。

ここで、

は、退歩行が除外されたアクティブ制約の直交行列Ｅである。故に、退歩の検出と因数分解とは、同時に達成される。

次の式が成り立つとする。

ここで、λ^＊は、前方解決（forward solve）および後方解決（back solve）によって、則ち、システム上の２つのＯ（ｎ^２）処理によって計算されることができる。

ｕ^＊を以下の通りに計算できる。

〔行列積を更新すること〕
スクラッチから計算されるよりむしろ、先行する反復からの行列が与えられると、上記計算において出現する行列ＭおよびＧ＝Ｌ^−ＴＭは更新されることができる。このことは、手間のかかる行列積処理を２つのＯ（ｎ^２ ）処理に縮小する。

Ｅに追加された特別な行ａを有する行列

を考察する。一般性（generality）の損失無しに、これは以下の通りに表現されることができる。

そして、以下の通りである。

Ｍの情報が与えられると、計算されることが要求される全てはＬ^−１ａ^Ｔである。これは、Ｏ（ｎ^２）のみを費やさせるシステム

上における前方解決である。同様に、

であり、かつ、類似のＯ（ｎ^２）処理を伴って更新されることができる。

１つのアクティブセット反復から次までにおいて、行がＥから落とされるか（更新は、単に、対応する列をＭおよびＧから落とすことである）、または、行が追加されるか、または、両方であり、これら全ての場合において、行列積がどのようにして更新されることができるかは明確である。その上、我々の方法では、先行する反復から、全てのベクトル

を記憶する。そのため、もしかつてアクティブであった制約が後の反復でのアクティブセットにおいて再出現するならば、対応する

は、再計算される必要がない。このことは、Ａのサイズの２倍に等しい特別な記憶装置を必要とする。それは、もしＭＰＣ問題の我々の縮小された緻密な公式化と利用可能なメモリが与えられるならば、極微な関心事である。

〔因数分解を更新すること〕
全ての反復においてＭのＱＲ因数分解を更新することも可能である。それは、本質的には、残余のＯ（ｎ^３）処理のみをより手間のかからない計算と置換する。

Ｍを上部の三角形式に入れるハウスホルダ変換のシーケンスを考える。
Ｑ^ＴＭ＝Ｒ

始めに、たった１つの行ａがＥ（および、故にＭ）に追加され、何も落とされない簡単な場合を考える。

列よりも多くの行を有しかつ上部の三角形であるＲに

を付加する。その結果、以下に類似した構成を生じる。

ここで、＊はＲの要素を示し、かつ、ｘはベクトル

の要素を示す。最後の列における１つのハウスホルダ変換、所謂

は、以下に示す上部の三角形式に、行列を置くことができる。

これは、Ｏ（ｎ）のみのコストを要し、最後の列における関連するサブベクトルの単純な標準計算を含む。これは、Ｒ行列を更新し、かつ、更新されたＱは

として示される。

行がＥから落とされるとき、更新はわずかにより複雑である。同様に、中間の列を落とすことによって得られるＭを考慮し、かつ、「Ｑ^Ｔが、Ｍを三角形式に置くハウスホルダ変換のシーケンスを示す」ということを仮定する。そして、以下の式が成り立つ。

ここで、

は、中間の列を落とすことによって第１のより低い三角帯（first lower triangular band）に落ちたＲの対角線上の以前の要素を示す。もしｊ番目の列が落とされたならば、より低い三角帯をゼロにするために、一連の所与の回転（Given Rotations）が

に適用される。このことは、多くとも、２つのハウスホルダ変換のコストを含み、これは、落とされた列の位置に依存する。そのような一連の所与の回転が

則ち以下の通りであるとする。

故に、更新されたＱが

として表現される一方、上記の上位三角化された行列（upper-triangularized matrix）は、更新されたＲを提供する。

「ハウスホルダまたは所与の回転行列の何れも記憶することは、決して必要ない」ということに注意されたい。ハウスホルダ変換に関連するベクトルｖと、計算における各々の所与の変換に関連するパラメータｃｏｓθとを記憶することは容易である。これは、列との乗算を実行するからである。変換の特性を使用して直交変換のシーケンスによって乗算することは、直交変換の積を形成しかつ結果としての緻密な行列によって乗算することより労力がかからない。

〔実行可能性を復元すること〕
我々が解決している問題は非線形力学を有するので、「新たな線形化によって生じる差異によって、先行する時間期間からの最適な制御ｕ_ｆが現在の時間期間において実行可能でない」という場合がしばしばある。そのような実行不可能性（infeasibility）はまた、非線形ＭＰＣに対するBettsのＳＱＰアプローチにおいて生じるＱＰｓを解決するにつれて、Betts によってこうむられ、かつ、彼は、実行不可能性を最小化するために、「フェイズＩ」問題を解くことを示唆している。これは、線形制約を伴う問題における一般的なアプローチである。しかしながら、このアプローチでは、最適化への進歩を成さない一方で、許容された計算時間を消耗させかつ実行可能な問題解決を未だ有しない。これらの理由のために、我々は、所謂「ビッグＭ（big-M）」問題を解決することを選択する。これは、元の目的の重み付けられた合計と実行不可能性の測定とを最小化する。元のＱＰ（６）を解く代わりに、伝統的なビッグＭ方法は、特別な変数τを導入し、かつ、以下の式を解くことを提案する。

ここで、もしａ_ｉｕ_ｆ＞ｂならば（ａ_ｉはＡのｉ番目の行を示す）、ｄは、位置Ｉに１を伴うベクトルであり、他の場合に０を伴う。「ｕ＝ｕ_ｆであり、τ＝ｍａｘ_ｉ｛ａ_ｉｕ_ｆ−ｂ_ｉ｝は、ビッグＭのＱＰに対する実行可能な点であり、それによって、我々は実行可能な点（feasible point）の情報を有する」ということが注目される。

正確なペナルティ機能（exact penalty functions）の理論から、「全てのＭ≧Ｍ＿に対して、ビッグＭのＱＰに対する問題解決が、元のＱＰ（６）に対して実行可能でありかつ最適であるような、即ち、τ＝０を有するような有限のＭ＿が存在する」ということが推定される。

この方法を我々の問題に適応させるために、我々は、ｕおよびτに関するヘッシアン量（Hessian quantity）が逆数であることを必要とする。これを達成するために、項（１／２）×Ｋτ^２を目標に追加する。数学上の実験は、更に、我々にＭ＝０を設定させる。故に、τ≧０制約の必要を軽減する。但し、非ゼロであるが正であるＭと制約τ≧０とを残すことが可能である。我々が「ビッグＭ」問題と呼ぶ我々の問題は、以下の通りである。

「ビッグＭ問題に対する正確なペナルティ機能の理論によって意味される結果が未だ持続する」ということを理解することは容易である。

〔初期アクティブセットを推測すること〕
初期アクティブセットに対する推測は、モデル予測制御問題公式化に密接に関する。Ｔ＝｛ｔ_０，ｔ_０＋Δｔ，ｔ_０＋２Δｔ，……ｔ_０＋（Ｎ−１）Δｔ｝が時間ｔ０における個々のＭＰＣ問題の時間窓内の考察の下における時間点の組を示すとする。動的システムに関する変数はｘ（Ｔ），ｕ（Ｔ），ｙ（Ｔ）として表現される。ここで、例えば、ｕ（Ｔ）は、ｕ（Ｔ）＝｛ｕ（ｔ_０），ｕ（ｔ_０＋Δｔ），ｕ（ｔ_０＋２Δｔ），……ｕ（ｔ_０＋（Ｎ−１）Δｔ）｝という考察の下での時間点における制御変数の組を示す。

更に、わずかに異なる表記法において、（ｕは最適化変数をなお示す）、我々は、ＭＰＣ問題の如何なる個々の時間依存の制約ｇをも以下の通りに記述できる。

「ｔ∈Ｔなので、各々のそのような制約は、最適化問題におけるＮ個の代数的な制約として、実際に現れる」ということに注目されたい。我々の二次のプログラム（quadratic program）に対して、ｇは実際に線形であり、かつ、ｇ_ｔは単にｕ_ｔ−１，ｕ_ｔ，ｔ_ｔ＋１の関数である。しかしながら、アクティブセットに対する初期推測に到達する方法は、非常に多くの一般的な設定において受容可能である。

「

が、時間ｔ_０におけるＭＰＣ問題に対応する最適化問題に関する反復的なアルゴリズムを実行することによって到達される最適な制御の値を示す」ということを仮定する。この問題に対するアクティブセットアルゴリズムの最終的な反復は「どの制約がアクティブセット内にあるのか」ということを我々に教える。「ｇ_ｔが

に対してアクティブである（Ｔ＿は、組Ｔにおけるいくつかのまたは全ての点を示す）」ということを仮定する。これは、則ち、以下の通りである。

そして、次の時間点ｔ_０＋ΔｔにおけるＭＰＣ問題は、シフトされた時間窓Ｔ’＝｛ｔ_０＋Δｔ，ｔ_０＋２Δｔ，……ｔ_０＋ＮΔｔ｝にわたって定義される。

同様に、制約ｇ_ｔ は、この以下のＭＰＣ問題において、全てのｔ∈Ｔ’に対して０≦ｇ_ｔ（ｕ）≦β_ｔとして現れる。これは、我々がアクティブセット推測することを望む最適化問題である。

そして、時間ｔ_０＋ΔｔにおけるＭＰＣ問題に対する初期推測として、制約ｇ_ｔが、全ての

に対してアクティブである（故に、初期アクティブセットＷ_０を入れる）と推測される。更に、もしｇ_ｔ が時間ｔ_０ におけるＭＰＣ問題の最後の時間点ｔ_０＋（Ｎ−１）Δｔにおいてもアクティブであったならば、ｇ_ｔ は、時間ｔ_０＋ΔｔにおけるＭＰＣ問題の最後の時間点ｔ_０＋ＮΔｔにおいてアクティブであると推測される。

この手続きは、アクティブセットに対する初期推測に到達するために、ＭＰＣ問題における全ての制約に対して繰り返される。

特許の法律および法律学の規定によって、上述された例示的な構成は、本発明の好ましい実施形態を示すために考察される。しかしながら、「本発明は、本発明の趣旨または範囲からはずれることなく、具体的に図解されかつ記述される以上の別の方法で実行されることができる」ということが注目されるべきである。本方法の請求項におけるステップに対する英数字の識別子は、従属請求項による参照の容易さのためであり、かつ、もし別の方法で示されないならば、要求されるシーケンスを示さない。

モデル予測制御を使用する制御システム１０の一般的なモデルを示す図面である。

符号の説明

１０ 制御システム
１２ 所望の軌道発生器
１４ 線形化モジュール
１６ 二次のプログラミング公式化モジュール
１８ 二次のプログラミング解決器
２０ 作動システム
２２ プラント
２４ センサシステム

Claims (37)

1. 二次のプログラミング問題を公式化しかつ最適化するための方法であって、
   本方法は、
   ａ）複数の時間ステップの各々において、多数の時間ステップにわたる窓に対して所望の動的な応答を達成するという問題を、二次のプログラミング問題に対する問題解決として公式化する段階と、
   ｂ）最適なアクティブセットを探索する反復的なアルゴリズムを使用して各時間ステップにおける二次のプログラミング問題を解く段階と、
   ｃ）複数の時間ステップのその後の各々の時間ステップにおいて、複数の時間ステップの先行する時間ステップの最終のアクティブセットに基づいて、最良のアクティブセットに対する探索を初期化する段階と
   を具備する
   ことを特徴とする方法。
2. アクティブセットは、最適化された問題解決において拘束している１組の制約を具備する
   ことを特徴とする請求項１記載の方法。
3. 前記段階ｅ）は、負のラグランジュ乗数を伴うたった１つの制約を落とすという手続きを更に具備する
   ことを特徴とする請求項２記載の方法。
4. 前記段階ｅ）は、制約のラグランジュ乗数の符号に関係なく拘束していない全ての制約を落とすという手続きを更に具備する
   ことを特徴とする請求項２記載の方法。
5. 先行する時間ステップの最終のアクティブセットが現在の時間ステップにおいて実行可能でない、ということを決定する段階と、
   元のＱＰ問題の重み付けられた合計と実行不可能性の測定値とを最小化する段階と
   を更に具備する
   ことを特徴とする請求項２記載の方法。
6. 実行可能なポイントがアクティブセットと大きく相反する、ということを決定する段階と、
   最適化を継続しながら、アクティブセットにおける不一致の制約を保持する段階と
   を更に具備する
   ことを特徴とする請求項２記載の方法。
7. 前記段階ｅ）の反復的なアルゴリズムの各反復においてＱＲ因数分解を実行する段階
   を更に具備する
   ことを特徴とする請求項１記載の方法。
8. 反復的なアルゴリズムの第１の反復は、完全なＱＲ因数分解を実行する段階を具備し、かつ、続く反復は、ＱＲ因数分解更新を実行する段階を具備する
   ことを特徴とする請求項７記載の方法。
9. 各ＱＲ因数分解の間に退歩を検出する段階
   を更に具備する
   ことを特徴とする請求項７記載の方法。
10. 退歩制約を落とす段階
    を更に具備する
    ことを特徴とする請求項９記載の方法。
11. 予備処理無しに各々の反復においてアクティブ制約における線形依存を検出する段階
    を更に具備する
    ことを特徴とする請求項１０記載の方法。
12. 請求項１の方法を使用して多変数のシステムを制御する方法であって、
    本方法は、
    ｄ）システムの現在の状態を示す複数のセンサ信号を受信する段階と、
    ｅ）複数のコマンドを受信する段階と、
    ｆ）コマンドおよびセンサ信号に基づいてシステムの所望の動的な応答を決定する段階と、
    ｇ）前記段階ａ）において、所望の動的な応答を達成するという問題を公式化する段階と
    を更に具備する
    ことを特徴とする方法。
13. 二次のプログラミング問題を反復的に解くための方法であって、
    本方法は、
    第１の反復においてＨのコレスキー因数分解を実行する段階と、
    第１の反復に続く反復においてコレスキー因数分解を再使用する段階と、
    反復のうちの１つにおいて制約されない問題解決を計算する段階と、
    反復のうちの前記１つに続く反復において前記制約されない問題解決を再使用する段階と
    を具備する
    ことを特徴とする方法。
14. 反復的なアルゴリズムの各反復においてＱＲ因数分解を実行する段階
    を更に具備する
    ことを特徴とする請求項１３記載の方法。
15. 反復的なアルゴリズムの初期反復は、完全なＱＲ因数分解を実行する段階を具備し、かつ、続く反復は、ＱＲ因数分解更新を実行する段階を具備する
    ことを特徴とする請求項１４記載の方法。
16. 各ＱＲ因数分解の間に退歩を検出する段階
    を更に具備する
    ことを特徴とする請求項１４記載の方法。
17. 退歩制約を落とす段階
    を更に具備する
    ことを特徴とする請求項１６記載の方法。
18. 予備処理無しに各々の反復においてアクティブ制約における線形退歩を検出する段階
    を更に具備する
    ことを特徴とする請求項１７記載の方法。
19. 二次のプログラミング問題を反復的に解くための方法であって、アクティブセットにおける制約のインデックスの組がＷによって示され、アクティブセットに対応する制約行列Ａにおける行がＥによって特徴付けられ、実行可能な点が既知であり、Ａｕ_ｆ≦ｂであり、
    本方法は、
    ａ）初期推測アクティブセットを選択する段階と、
    ｂ）反復ｋにおいて、アクティブ制約に対して推測Ｅが与えられると、以下の同等制約されたＱＰ（ＥＱＰ）を解く段階と
    

    

    を具備し、
    ここで、ｒはアクティブ制約に対する右手側のサブベクトルを示し、
    本方法は、
    ｃ）最適なｕ^＊とラグランジュ乗数λ^＊とを以下の通りに決定する段階
    

    

    を具備する
    ことを特徴とする方法。
20. ｕ^ｋ−１＋αｓが実行可能であるような最大のα∈〔０，１〕を決定する段階
    を更に具備し、
    ｓ＝ｕ^＊−ｕ^ｋ−１であり、ｕ^ｋ−１は先行する反復を示し、ｕ^０＝ｕ_ｆである
    ことを特徴とする請求項１９記載の方法。
21. αは以下の式によって与えられる
    ことを特徴とする請求項２０記載の方法。
    

    
22. 次の反復におけるアクティブセットに対する推測を更新する段階
    を更に具備する
    ことを特徴とする請求項２１記載の方法。
23. 比率テストにおいてその最小値を達成するαに対応する第１のインデックスｉ_Ｂ、次の反復におけるアクティブインデックスに対する推測をＷに加える段階
    を更に具備する
    ことを特徴とする請求項２２記載の方法。
24. アクティブセットから最もネガティブな乗数を伴う制約を落とす段階
    を更に具備する
    ことを特徴とする請求項２３記載の方法。
25. 第１の反復において、もしα＜１であるならば、一致しない全ての制約を落とす段階
    を更に具備する
    ことを特徴とする請求項２４記載の方法。
26. 所望の軌道を生成する所望軌道発生器と、
    所望の軌道について線形化されたモデルを得る線形化モジュールと、
    複数の時間ステップの各々において、二次のプログラミング問題に対する問題解決として、多数の時間ステップ窓に対する所望の軌道を達成するという問題を公式化する二次のプログラミングモジュールと、
    最適の制御のプロフィールを生成するために二次のプログラミングモジュールによって確立された最適化問題を解く二次のプログラミング解決器と
    を具備し、
    二次のプログラミング解決器は、最適のアクティブセットを探索する反復アルゴリズムを使用して各時間ステップにおける二次のプログラミング問題を解き、かつ、複数の時間ステップの後の各々の時間ステップにおいて、複数の時間ステップの先行する時間ステップの最終のアクティブセットに基づいて最良のアクティブセットに対する探索を初期化する
    ことを特徴とするモデル予測制御システム。
27. アクティブセットは、最適化された問題解決において拘束している１組の制約を具備する
    ことを特徴とする請求項２６記載のシステム。
28. 二次のプログラミング解決器は、負のラグランジュ乗数を伴うたった１つの制約を落とす
    ことを特徴とする請求項２７記載のシステム。
29. 二次のプログラミング解決器は、制約のラグランジュ乗数の符号に関係なく拘束していない全ての制約を落とす
    ことを特徴とする請求項２７記載のシステム。
30. 二次のプログラミング解決器は、前記段階ｅ）の反復アルゴリズムの各々の反復においてＱＲ因数分解を実行する
    ことを特徴とする請求項２６記載のシステム。
31. 反復的なアルゴリズムの第１の反復は、完全なＱＲ因数分解を実行することを具備し、かつ、続く反復は、ＱＲ因数分解更新を実行することを具備する
    ことを特徴とする請求項３０記載のシステム。
32. 二次のプログラミング解決器は、各ＱＲ因数分解の間に退歩を検出する
    ことを特徴とする請求項３０記載のシステム。
33. 二次のプログラミング解決器は、退歩制約を落とす
    ことを特徴とする請求項３２記載のシステム。
34. 二次のプログラミング解決器は、予備処理無しに各々の反復においてアクティブ制約における線形依存を検出する
    ことを特徴とする請求項３３記載のシステム。
35. 二次のプログラミング問題を解くために最適のアクティブセットを探索する反復アルゴリズムを使用して二次のプログラミング問題を最適化するための方法であって、
    本方法は、
    初期開始点組が実行可能でないことを決定する段階と、
    実行可能な反復を取得するために問題を再公式化する段階と
    を具備する
    ことを特徴とする方法。
36. 二次のプログラミング問題を解くために最適のアクティブセットを探索する反復アルゴリズムを使用する段階
    を更に具備する
    ことを特徴とする請求項３５記載の方法。
37. 元のＱＰ問題の重み付けられた合計と実行不可能性の測定値とを最小化する段階
    を更に具備する
    ことを特徴とする請求項３５記載の方法。

Images