(57)【要約】【課題】  モデル予測制御問題を解く効率的な方法を提供すること。【解決手段】  モデル予測制御問題に基づいて、大規模疎行列を形成する。次いで、Hをまず形成せずに、Hの平方根Hrを直接形成する。次いで、大規模行列式の大規模行列の平方根(LSMルート)を、大規模疎行列をまず形成せずに、かつ複数の反復のそれぞれにおいてHrを再計算せずに、2次プログラミングソルバの複数の反復のそれぞれにおいてHrを使用して形成する。LSMルートに基づいて、大規模行列式の解を完成する。【選択図】  図１
2005038397.tif 000002
【特許請求の範囲】
【請求項１】
モデル予測制御問題を解く方法であって、
a)前記モデル予測制御問題に基づいて、大規模疎行列式LSMEを形成するステップと、
b)モデル予測制御に基づく2次プログラミング問題の解について作動装置コマンドの決定を呈示するように、制御装置の更新につき少なくとも1回、コマンドおよびフィードバックセンサを反復してサンプリングするステップと、
c)制御装置の更新につき1回のみ変化する前記LSMEの部分行列であるHの平方根である、対角に沿った各ブロックが3角であるブロック3角行列であるHrを形成するステップと、
d)各制御装置更新内における2次プログラミングソルバの複数の反復のそれぞれにおいてHrを使用して、前記LSMEの大規模疎行列の平方根であるLSMルートを形成し、前記LSMルートが、前記対角に沿った各ブロックが3角であるブロック3角行列であるステップと、
e)各反復においてLSMルートに基づいて前記LSMEの前記解を完成するステップとを含む方法。
【請求項２】
Hrが、前記複数の反復の第1において計算され、前記複数の反復のその後の反復について再計算されない、請求項1に記載の方法。
【請求項３】
前記ステップc)において、Hrが、Hをまず形成せずに形成される、請求項2に記載の方法。
【請求項４】
前記ステップc)が、QR因数分解のシーケンスで実施される、請求項3に記載の方法。
【請求項５】
H=Hr'*S*HrおよびHr=LW-Tであり、
【数１】
2005038397.tif 000034
である、請求項3に記載の方法。
【請求項６】
部分行列F、G、およびMが、QR因数分解の帰納シーケンス
【数２】
2005038397.tif 000035
によって、予測範囲の様々な将来の制御装置更新時点について生成される、請求項5に記載の方法。
【請求項７】
前記LSMルートが、LSMをまず形成せずに、前記ステップd)において形成される、請求項1に記載の方法。
【請求項８】
前記ステップd)が、通常のまたは一般化されたQR因数分解を解くことによって実施される、請求項7に記載の方法。
【請求項９】
前記ステップd)が、内点法に適応された一般化QR因数分解を解くことによって実施される、請求項7に記載の方法。
【請求項１０】
前記一般化QR因数分解が、以下の通りであり、
【数３】
2005038397.tif 000036
上式で、PインおよびSが入力であり、3角Pおよび可能であればUが出力であり、
【数４】
2005038397.tif 000037
である、請求項9に記載の方法。
【請求項１１】
前記ステップd)が、アクティブセット法に適応された一般化QR因数分解を解くことによって実施される、請求項7に記載の方法。
【請求項１２】
前記一般化QR因数分解が、
【数５】
2005038397.tif 000038
であり、上式で、PインおよびSが入力であり、PおよびSRならびに可能であればUが出力であり、
【数６】
2005038397.tif 000039
である、請求項11に記載の方法。
【請求項１３】
前記一般化QR因数分解が、通常使用されるハウスホルダ変換
【数７】
2005038397.tif 000040
のシーケンスを一般化された変換
【数８】
2005038397.tif 000041
の対応するシーケンスと置き換えることによって修正された標準的なハウスホルダ変換に基づくQRアルゴリズムにより達成され、Uは必要ではないが、所望であれば、Ui累加積である、請求項12に記載の方法。
【請求項１４】
純粋に実数または純粋に複素数を含んでいる行からなる3角行列である一般化QR因数分解の結果が、前記対角に1または-1の平方根を有する対角行列SR1/2と、純粋に実数の3角行列との積にさらに因数分解され、この対角行列の平方Srが、前記対角に1または-1を有し、SrおよびPが純粋に実数である、請求項13に記載の方法。
【請求項１５】
実数のみを使用して以下のようにLSMを因数分解するステップをさらに含み、
【数９】
2005038397.tif 000042
上式で、SおよびSRが、前記対角に沿って1または-1のみのエントリを有する一定対角行列であり、LSMルートが、各対角部分行列も3角であるブロック3角行列である、請求項11に記載の方法。
【請求項１６】
LSMEが、
【数１０】
2005038397.tif 000043
の形態を有し、上式で、zは、時間によってグループ化された各時間点の状態、制御、および状態式共役変数を含んでいるベクトルであり、mは、時間によってグループ化された不等号制約の共役変数であり、fおよびKはベクトルであり、HおよびJはバンド行列であり、選択された2次プログラムアルゴリズムに応じて、Tは対角行列またはゼロである、請求項1に記載の方法。
【請求項１７】
システムを制御する方法であって、
a)前記システムの現行条件を示す複数のセンサ信号を受信するステップと、
b)複数のコマンドを受信するステップと、
c)前記コマンドに基づいて、前記システムの所望の動的応答を決定するステップと、
d)複数の時間ステップのそれぞれにおいて、モデル予測制御の方法および前記センサ信号を使用して、2次プログラム問題に対する解として、複数時間ステップにわたるウィンドウについて所望の動的応答を達成する問題を定式化するステップと、
e)大規模行列式LSMEが、前記モデル予測問題および2次プログラミングアルゴリズムに基づいて形成される反復アルゴリズムを使用して、各時間ステップにおいて前記2次プログラミング問題を解くステップと、
f)前記反復アルゴリズムの第1反復において、Hの平方根Hrを形成するステップと、
g)前記反復アルゴリズムの前記第1反復と、前記第1反復に続く反復とを含む反復において、2次プログラミングソルバのHrを使用して、前記LSMEの大規模行列の平方根LSMルートを形成するステップと、
h)前記反復アルゴリズムの前記第1反復と、前記第1反復に続く反復とを含む反復において、LAMルートに基づいて前記LSMEの解を完成するステップとを含む方法。
【請求項１８】
プラントを制御するモデル予測制御システムであって、
該システムの現行状態を示す複数のセンサと、
コマンドに基づいて所望の動的応答を創出する所望軌跡生成装置と、
該システムの現行状態から開始して、モデル予測制御の方法を使用して2次プログラム問題に対する解として、複数の時間ステップにわたるウィンドウについて前記所望の動的応答を達成する問題を形成する2次プログラミングモジュールと、
大規模疎行列式LSMEが、前記モデル予測制御問題および2次プログラミングアルゴリズムに基づいて形成される反復アルゴリズムを使用して、各時間ステップにおいて前記2次プログラミング問題を解く2次プログラミングソルバであって、該ソルバが、前記反復アルゴリズムの第1反復において、Hの平方根Hrを形成し、該ソルバが、前記反復アルゴリズムの前記第1反復と、前記第1反復に続く複数の反復とを含む反復において、2次プログラミングソルバのHrを使用して前記LSMEの大規模疎行列の平方根LSMルートを形成し、該ソルバが、前記反復アルゴリズムの前記第1反復と、前記第1反復に続く前記複数の反復とを含む反復において、LSMルートに基づいて前記LSMEの前記解を完成するソルバとを備えるシステム。
【請求項１９】
前記ソルバが、前記第1反復においてHrを形成し、前記反復アルゴリズムの前記第1反復に続く前記複数の反復においてHrを再計算しない、請求項18に記載のシステム。
【請求項２０】
前記ソルバが、Hをまず形成せずにHrを形成する、請求項19に記載のシステム。
【請求項２１】
H=Hr'*S*HrおよびHr=LW-Tであり、
【数１１】
2005038397.tif 000044
である、請求項20に記載のシステム。
【請求項２２】
部分行列F、G、およびMが、予測範囲の様々な将来の制御装置更新時間点について、QR因数分解の帰納シーケンス
【数１２】
2005038397.tif 000045
によって生成される、請求項21に記載のシステム。
【請求項２３】
前記ソルバが、前記LSMをまず形成せずにLSMルートを形成する、請求項18に記載のシステム。
【請求項２４】
前記ソルバが、通常のまたは一般化されたQR因数分解を解くことによって前記LSMルートを形成する、請求項23に記載のシステム。
【請求項２５】
前記ソルバが、内点法について適応された一般化QR因数分解を解くことによって前記LSMルートを形成する、請求項23に記載のシステム。
【請求項２６】
前記一般化QR因数分解が、
【数１３】
2005038397.tif 000046
であり、上式で、PインおよびSが入力であり、3角Pおよび可能であればUが出力であり、
【数１４】
2005038397.tif 000047
である、請求項25に記載のシステム。
【請求項２７】
前記ソルバが、アクティブセット法について適応された一般化QR因数分解を解くことによってLSMルートを形成する、請求項23に記載のシステム。
【請求項２８】
前記一般化QR因数分解が、
【数１５】
2005038397.tif 000048
であり、上式で、PインおよびSが入力であり、PおよびSrならびに可能であればUが出力であり、
【数１６】
2005038397.tif 000049
である、請求項27に記載のシステム。
【請求項２９】
前記一般化QR因数分解が、通常使用されるハウスホルダ変換
【数１７】
2005038397.tif 000050
のシーケンスを、一般化された変換
【数１８】
2005038397.tif 000051
の対応するシーケンスと置き換えることによって修正された標準的なハウスホルダ変換に基づくQRアルゴリズムにより達成され、Uは必要ではないが、所望であればUiの累加積である、請求項28に記載のシステム。
【請求項３０】
純粋に実数または純粋に虚数を含む行からなる3角行列である一般化QR因数分解の結果が、対角が1または-1の平方根を有する対角行列SR1/2と、純粋に実数の3角行列との積にさらに因数分解され、この対角行列の平方Srが、対角上で1または-1を有し、SrおよびPが純粋に実数である、請求項29に記載のシステム。
【請求項３１】
前記ソルバが、以下のように、
【数１９】
2005038397.tif 000052
実数のみを使用してLSMを因数分解し、上式で、SおよびSRが、対角に沿って1または-1のエントリのみを有する一定対角行列であり、LSMルートが、各対角部分行列も3角行列であるブロック3角行列である、請求項27に記載のシステム。
【請求項３２】
前記LSMEが、以下の形態を有し、
【数２０】
2005038397.tif 000053
上式で、zは、時間によってグループ化された各時間点の状態、制御、および状態式共役変数を含んでいるベクトルであり、mは、時間によってグループ化された不等号制約の共役変数を含み、fおよびKはベクトルであり、HおよびJはバンド行列であり、Tは、選択された2次プログラムアルゴリズムに応じて対角行列またはゼロである、請求項18に記載のシステム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、一般に、制御システムに関し、より具体的には、モデル予測制御システムに関する。
【背景技術】
【０００２】
本発明は、米国政府契約N00421-01-2-0131の下での作業の性能について考案された。
【０００３】
モデル予測制御論理は、実際の多変数システムの制御について多くの利点を有する。多変数システムは、複数の目標(コマンドとしても知られている)を有し、それらの目標を達成するために、システムダイナミクスを変化させる複数のエフェクタを有する。多変数システムは、有意な交差結合を含むことも可能であり、主要なエフェクタが、それぞれ、複数の目標遂行することが重要である。したがって、このシステムの制御装置は、応答を分離するために、やはり多変数であるはずである。分離システムでは、目標が変化するとき、制御アルゴリズムは、他の目標を著しく変化させずに、この1つの目標を変化させるように、複数のエフェクタを設定することによって応答する。さらに、システムの交差結合は、しばしば動的であり、時間と共に変化する。この場合、いくつかの結合は永続的であり、いくつかは一過性であり、いくつかは迅速であり、いくつかは緩慢であり、または、まず一方向において結合し、次いで反転して、多方向において結合することが可能である。したがって、制御アルゴリズムの分離も、動的であるはずであり、結合の変化する性質を相殺するように時間と共に発展する。分離の他に、多変数制御論理は、システムが、安定で、安全で、耐久性であることを保証する一方で、システムダイナミクスの変化の影響および妨害の影響を低減する従来のフィードバック特性をも提供しなければならない。
【０００４】
結合多変数システムのフィードバック制御論理の設計を支持する多くの制御理論が存在し、いくつかは実施されている。この理論のほとんどは、制御されるシステムが、システムの変数またはエフェクタのいずれかに対して、本明細書では限界と呼ぶ不等号制約を有さないことを想定する。これは、理論と実施との非常に著しい不整合である。実際の実用的なシステムは、実質的で有意な限界を有する。エフェクタの範囲など、いくつかの限界は物理的であり、いくつかは、安定性、安全性、および耐久性を保証するために動作に関する。実際的な意味では、限界に到達しない場合、システムは、過度に設計されている。この不整合により、論理を多くのモードに分割することを含めて、制御法則の設計が非常に困難で、またその方法が場当たり的になる。この場合、限界または障害の条件の各順列を網羅するために、別々の制御モードが存在し、また、別々に設計された動的モード間の切替えを動的に制御するために、論理はより場当たり的になる。これらのモードを調整して、モード間に不都合な動的相互作用(2つのモード間を反復して前後にジャンプする、突然ジャンプするなど)が存在しないことを保証するために、著しく多くの場当たり的な論理が必要とされる。設計が場当たり的であることにより、制御ソフトウエアは設計が困難で、維持するのが困難になり、ソフトウエアを同じシステムに対して再使用することができる程度が大きく低減される。
【０００５】
モデル予測制御(MPC)は、多変数制御論理であり、限界を含むことが明白であり、したがって、実際のシステムとの良好な整合を提供する。MPCは、作動装置の障害など、システムの変化に実時間で応答するように構成することもできる。したがって、MPCは、限界に到達し、障害が発生した場合でも、応答を分離することができ、ならびに物理的に可能である制御アルゴリズム設計の正式な方法を多変数システムに提供する。MPC設計では、限界条件および障害に対処するための余分な制御モードの場当たり的な設計は必要ではない。これにより、制御装置ソフトウエア設計の生産性が著しく増大し、制御ソフトウエアを更新すること、そうでない場合は制御ソフトウエアを維持することがより容易になり、制御システムがより困難なタスクを達成することが可能になることが予期される。後者には、より自立的な動作、より完全に統合するシステム応答、および物理的な応力が低減された動作が含まれる。
【０００６】
モデル予測制御の実施に対する最も大きな障害は、必要とされる計算の量である。このために、化学産業において20年から30年前に始まったMPCの初期の実用化は、比較的緩慢なダイナミクスを有するシステムについてであった。化学プラント制御装置は、ほぼ5分に1度作動装置の設定を更新する。これらの初期の応用例は、当時利用可能であったコンピュータの能力を圧迫した。ガスタービンエンジンでは、更新率は、1秒当たり約40回でなければならず、約12,000倍の差であり、制御機能は、より複雑になる可能性がある。それでも、MPCは、エンジンおよび他の車両応用分野について実用的である可能性があるが、その理由は、現在、500倍から1000倍迅速なコンピュータが利用可能であるからである。それでも、システムの複雑さに応じて、MPCが明らかに必要とする計算と利用可能なコンピュータの能力との間には、10倍から30倍のギャップが依然として存在する。本発明は、モデル予測制御の計算の負担を20分の1から100分の1に低減して、ギャップをなくす計算技法を提供するので、本発明は必須である。
【０００７】
モデル予測制御
モデル予測制御理論は、化学産業において最も顕著に、実用化されてきた。これを裏付けるものとして、モデル予測制御理論は、多くの技術論文および教科書において記述されてきた。ここでそれを概述する。この概述の主な目的は、MPCの実施には、大量のオンボード計算が必要であり、また、この計算の約90%が、実時間での大規模行列の多くの解を含むことを示すことである。このセクションは、最適な式および変数がある方式で順序付けされているとき、行列式が、特定の疎でバンドの対称的な構造を有することを示す。これらの特性は、行列が、多くのゼロを含み、非ゼロのエントリが対角付近に集まり、行列の「平方根」が存在することを意味する。この特定の形態は、このセクションにおいて発展され、大規模疎行列式(LSME)の式と呼ばれる。平方根法は、計算を大きく低減するために、この構造を新規な方式で利用する。
【０００８】
モデル予測制御では、作動装置コマンドuが、最適化問題を実時間で解く結果として決定される。モデル予測制御は、あるタイプのデジタル制御装置である。デジタル制御装置は、オペレータコマンドおよび測定応答を周期的にサンプリングするコンピュータに基づく。各サンプルの結果として、デジタル制御装置は、制御装置更新と呼ばれる新しい作動装置設定を指定する。これらの設定は、次の更新まで、最もしばしば保持される。MPCは、これらの設定を作動装置軌跡の第1時間点として指定する。軌跡は、LからN+L-1まで変化する整数nで番号付けされた制御装置更新時間について、すなわちNの時間点について、u(n)を指定する。現行サンプル時間は、Lである。この作動装置軌跡は、性能指数PI
【０００９】
【数２１】
2005038397.tif 000003
【００１０】
を最小限に抑える固有の軌跡として決定される。上式で、和は、Nの未来時間点にわたって行われる。xnは、システムの動的モデルの状態ベクトルであり、unは、Nの時間点のそれぞれにおける作動装置コマンドである。行列CおよびDは、線形出力式の係数であり、この係数は、目標が存在するシステム出力を動的モデル状態および作動装置設定に関係付ける。QおよびRは、下方3角加重行列であり、ベクトルfおよびgは、システム出力の所望の軌跡と、所望のまたは公称の作動装置の値と、加重行列QならびにRとの関数である駆動項である。下方3角行列Mは、最終加重行列であり、様々な制御論理上の理由により存在する可能性がある。
【００１１】
公称システムしたがって最適化は、システムの状態動的モデルに従わなければならない。
xn+1=An・xn+Bn・un+bn
最適化の手続きは、等号制約と呼ばれる。最適化は、Nの時間点について、
限定変数yc(n)=Ccn・xn+Dcn・un+an≦Yn
に関連した物理的および動作上の限界にも従わなければならない。第1等号は、限定変数が状態動的モデルに関連付けられる方式のモデルを表す。限界の式が理論および結果的な論理の設計に明瞭に含まれることは、MPCの重要な利点である。
【００１２】
MPC設計者が、状態動的モデルと、限界モデルと、性能指数とを決定することによって、制御アルゴリズムを指定する。性能指数の設計には、出力ベクトルに含むように変数を選択することと、加重行列Q、R、およびMを選択することとが含まれる。制御ソフトウエアの設計者は、最適化問題を実時間で解くための再使用可能なソフトウエアも提供しなければならない。
【００１３】
上記の最適化問題は、2次プログラミング問題であり、条件付き数学的最適化の特別なタイプである。2次プログラミングは、汎用条件付き最適化よりはるかにより信頼性が高く、関連する制約は温和である。コンピュータ上で2次プログラミング(QP)問題を解く2つの主な方法は、アクティブセット(AS)法および内点(IP)法である。以下で記述する本発明は、どちらかに必要とされる計算量を著しく低減する。
【００１４】
MPC動作では、現行サンプル時間における状態ベクトルxLが与えられる。これは、制御装置更新時間Lまでシステム応答の測定関数である。したがって、MPC制御装置は、フィードバック制御装置である。上記の最適化の正味の出力は、u軌跡の第1時間点におけるベクトルuLであり、これは、MPC制御装置のL番目の更新の作動装置コマンドである。
【００１５】
次の作動装置コマンドuL+1を決定するために、変数が時間について変位している同じMPC最適化を実施しなければならない。制御装置の更新期間ごとに、MPC最適化を再び解かなくてはならない。次の3つのセクションは、上記のタイプの最適化を分析的な方式で最も一般的に解く方式について概述する。
【００１６】
制約をJに吸収する随伴法
最適制御を解く通常の方法では、以下のように増強性能指数を形成するために、等号制約(状態動的モデル)は、ラグランジュ乗数pnで性能指数に随伴され、不等号制約(物理的限界に関するシステム変数のモデル)は、乗数mnに随伴される。
【００１７】
【数２２】
2005038397.tif 000004
【００１８】
ラグランジュ法では、随伴制約項は、制約がないかのように増補性能指数が最適化されたときに、制約が満たされるように定式化される。随伴項も、最適ではゼロの値を有する。したがって、制約を随伴させることにより、最適値Jは変化しない。不等号制約では、随伴法は、各限界について、各時間nにおいて対応してm=0であるか、または不等号制約がその限界にあることを必要とする。
【００１９】
カルシュ-キューン-タッカー(KKT)条件
解が最適なシステム軌跡を指定する1組の代数方程式が、(KKT)条件として知られている。上記のMPC最適化では、解が、時間Lと時間N+L-1との間の描く制御装置更新について、以下のKKT条件を含むことが知られている。
【００２０】
【数２３】
2005038397.tif 000005
【００２１】
および、すべてのtおよびmは、ゼロより大きくなければならないという追加の条件。これらの式を解いた場合、解の一部であり、かつシステムダイナミクスの状態式モデルの入力として使用される作動装置コマンドu(n)は、限界を満たしながら、性能指数を最小限に抑える。
【００２２】
2次プログラミングの一般的な大規模疎行列式
Nの時間点に対するKKT条件は、以下の大規模疎行列式(LSME)に集約することができる。内点(IP)法またはアクティブセット(AS)法を使用してMPC2次プログラミング問題を解くとき、これを反復して解かなければならない。LSMEは、次の形態を有する。
【００２３】
【数２４】
2005038397.tif 000006
【００２４】
上式で、zは、時間によってグループ化された各時間点の状態x、制御u、および状態式共役変数pを含んでいるベクトルである。mは、時間によってグループ化された不等号制約の共役変数である。fおよびkはベクトルである。HおよびJは、バンド行列である。ここで
【００２５】
【数２５】
2005038397.tif 000007
【００２６】
Tは対角行列であり、行列は、時間L=0における制御装置更新について示されている。行列Hも対称的であり、これは、Hr'*Hr=Hであるように、平方根行列Hrが存在することを意味する。LSMEは、主要な2次プログラムソルバの両方に共通であるが、Tの値および右辺(RHS)ベクトルが、使用する2次プログラムソルバに依存する点が異なる。どちらの方法でも、各制御装置の更新について、LSMEを反復して解かなければならない(システムの複雑さに応じて、10回から数百回)。
【００２７】
LSMEのこの形態は、新規で特別である。これは、最適化の式(LSMEの行)およびzの変数を特別な方式で順序付けすることによって達成される。それにより、結果的なLSME行列は、疎、バンド、および対称的になる。
【００２８】
内点法
2次プログラムを解く内点法では、Tは対角である。各要素は、各時間点について、対応するスラック変数tと共役変数mとの積の平方根である。Jは、疎バンド行列であり、すべての限界の式を含む。内点法の最も一般的な族では、Tおよび右辺ベクトルは、反復ごとに2回変化する。右辺ベクトルは、以前の反復の解から決定される。HおよびJは、同じ制御装置更新期間の反復では一定のままである。最適な作動装置設定の解は、問題の複雑さに応じて、10回から数百回の反復を必要とする。さらなる説明については、具体的にはスラック変数が排除されている場合である、参考文献1を参照されたい。
【００２９】
アクティブセット法
アクティブセット法は、どの限界がアクティブであるかを探索する。不等号制約のいくつかのみJc*z<=Ymaxが、アクティブである。アクティブセットqpソルバの各反復は、制約のどれがアクティブであるかを「推測」する。「アクティブな」ベクトルが、アクティブな制約の指数を含むとする。アクティブなとき、不等式は等式になり、Jc(アクティブ,;)*z=Ymax(アクティブ)であり、システムは、アクティブな限界について実行される。LSMEと比較して、J=Jc(アクティブ,:)、K=Ymax(アクティブ)、およびT=0である。所与の制御装置更新期間では、「アクティブ」のいくつかの反復が、適切なアクティブセットを探索する際に評価される。探索戦略は、LSMEの以前の反復の解の結果に基づいて、「アクティブ」の新しい反復を決定する。Hは、所与の制御装置更新の反復について、一定のままである。
【００３０】
この段階で重要な観測は、内点ソルバまたはアクティブセットソルバを使用するかにかかわらず、制御装置が更新されるたびに、MPC2次プログラムを解くために、LSMEのある形態を多数回解かなければならないということである。実際の制御装置では、Nおよび他のファクタに応じて、制御装置の更新につき、LSMEを10回から100回解く。
【発明の開示】
【発明が解決しようとする課題】
【００３１】
ここで記述する新規なアルゴリズムは、大規模行列式を因数分解する平方根法を提供することによって、この従来の技術を改善する。LSMEなどの行列式を解くことには、まず、行列を上方3角行列と、1つまたは複数の対角行列と、下方3角行列(またはその反対)との積に因数分解して、2つが3角行列を含み、残りが対角行列を含むより簡単な行列式を解くことが含まれる。
【課題を解決するための手段】
【００３２】
以下でより詳細に記述する他の新規な特徴の中でも、本発明は、大規模疎行列を形成することさえ必要とせずに、LSMルートを直接形成する新規な方法を使用する。本発明は、また、Hをまず形成することを必要とせずに、制御装置の更新につきHrを一度形成する新規な方法を使用して、2次プログラム解法を通してHrを使用する。
【００３３】
本発明のこれらおよび他の新規な特徴により、モデル予測制御の計算負担が、10分の1から30分の1に低減される。このような計算負担の低減は、ガスタービンエンジンなど、はるかにより迅速なダイナミクスを有するより複雑なシステムについて、実用的なモデル予測制御を使用する。
【００３４】
本発明の他の利点は、添付の図面と関連して考慮するとき、以下の詳細な記述を参照することによって理解することができる。
【発明を実施するための最良の形態】
【００３５】
本発明のアルゴリズムは、大規模疎行列式を因数分解する平方根法を提供する。LSMEなどの行列式を解くことには、まず、行列を上方3角行列と、1つまたは複数の対角行列と、下方3角行列(またはこの反対)との積に因数分解して、2つが3角行列を含み、残りが対角行列を含むより簡単な行列式を解くことが含まれる。予測されるMPC問題のサイズでは、因数分解は、組み合わされた対角行列と両方の3角行列とを解くより、1次から2次程度より複雑である。本発明は、この因数分解をはるかにより迅速に実施するアルゴリズムである。
【００３６】
記述の最高レベルでは、アルゴリズムは、以下によって優れた計算効率を達成する。
1.式および変数を順序付けして、対称的なバンドLSMEを得る。
2.制御装置が更新されるたびに、解かれるこのLSMEを多数回解く必要性を実現する。
3.LSME行列LSMが対称的なので、LSM=LSMc'*LSMcであり、かつLSMcが下方3角行列であるように、複素平方根行列LSMcを有することを実現する。
4.LSM=LSMルート'*S*LSMルートであり、Sが対角に沿って1または-1のエントリのみを有する一定対角であり、かつLSMルートが下方3角であるように、平方根行列のこの概念を一般化して、実数のみを使用する(計算効率のために)。
5.LSMを形成することさえ必要とせずに、LSMルートを直接形成する新規な方法を使用する。
6.所与の制御装置の更新時に、行列Hが、2次問題解法を通じて一定のままであり、かつ平方根Hrを有することを実現する。
7.Hを形成することを必要とせずに、Hrを直接形成する新規な方法を使用する。これは、より小さい周知のQR因数分解のシーケンスを含む。
8.それぞれがHrを再度使用する2つの一般化QR因数分解の一方を解くことによって、2次問題ソルバの各反復についてLSMルートを形成する新規な方法を使用する。
9.一般化QR因数分解を実施する新規なアルゴリズム。
10.LSMをステップ4の形態に因数分解した後、LSMEの解を完成するために、標準的な方法を使用する。
【００３７】
これらのステップは、LSMEの効率的な因数分解について2つの重要な新規な手段を実現する。第1は、部分行列Hの効率的な因数分解であり、これは、制御装置の更新につき1度実施するだけでよい(ステップ1、6、および7)。第2は、LSME全体をより迅速に因数分解するために、Hの因数分解を再使用することであり(ステップ8および9)、これは、QP法および問題のサイズに応じて、制御装置の更新につき10回から200回(QPの反復につき1回または2回)実施する必要がある。
【００３８】
制御装置の更新につき1回:Hを因数分解する
HおよびJの部分行列:A、B、C、D、Q1/2、R1/2、およびCcならびにDcを入力する。これらは、この制御装置の更新のMPC問題を記述する。Hを以下のように因数分解する。
【００３９】
ステップI:標準的なQR因数分解アルゴリズムを使用して、Nの別々でより小さいQR因数分解を計算する。
【００４０】
【数２６】
2005038397.tif 000008
【００４１】
n=Nから開始して、N-1、N-2、・・・、1まで時間を後退して作業する。左辺行列は、因数分解への入力である。右辺行列は、結果である。標準的な方法では、直交行列オルトを生成するかは随意選択である。ここでは必要ではない。行列A、B、C、D、Q1/2、およびR1/2は、MNのように、事前に指定されている。次いで、第1QR因数分解は、MN-1を決定し、これは、時間N-1の第2因数分解に入力される、などである。したがって、時間を後退してすべてのNを機能的に実施することができる。非制約MPC問題(JおよびKが0に等しい場合)を解く際のQR因数分解のこの帰納シーケンスおよびその値は、制御理論の教科書において既知である。Hを因数分解する際の使用(ステップ2)は、新規である。制御装置の更新につきHを1度だけ因数分解をすることを可能にする方式で因数分解されたHを再使用する(ステップ3、4、および5)。
【００４２】
ステップ2:下方3角行列であるLおよびブロック対角である行列Wを以下のように形成する。やはり、行列Sの対角要素を形成して、記憶する。対角要素のエントリは1であるが、エントリが-1であるベクトルzの共役変数pに対応する列を除く。
【００４３】
【数２７】
2005038397.tif 000009
【００４４】
効率的な行列計算の標準的な方法は、L、S、およびWの同じ情報をよりコンパクトな方式で記憶する。ここでは、理解を補助するために、疎行列の形態を使用する。
【００４５】
このアルゴリズムでは、Hを決定する必要はない。しかし、Hが所望される場合、
H=Hr'*S*Hr,左式でHr=LW-T
で計算することができる。
【００４６】
計算を保存するために、Hrは、ここでも明らかには形成されない。代わりに、L、S、およびWを別々に使用する。
【００４７】
各2次プログラムソルバ反復:LSMEを因数分解して解く
QPソルバ探索法から解く必要がある次のLSMEを記述するパラメータを入力する:T、K、およびf。一般に、2次プログラム解は、解の構造化された反復探索を含む。LSMEを指定するために使用される探索法およびパラメータは、ASおよびIPに関して異なる。したがって、平方根アルゴリズムの残りのステップは、IP法およびAS法について異なる。
【００４８】
ステップ3(IP):次のセクションにおいて記述する一般化QRアルゴリズムを使用して、以下の一般化QR因数分解を実施する。
【００４９】
【数２８】
2005038397.tif 000010
【００５０】
右側の角括弧の行列は、入力である。他の行列は、出力である。Uは、ここでは計算する必要はない。Pは、下方3角行列である。Tが対角で、JおよびLがバンドなので、Pもバンドであり、計算を低減するために、特別なQR因数分解を使用して、バンド構造を利用することができる。Uが正規直交ではないので、これは、新規な一般化因数分解である。代わりに、Uは、より一般的な特性を有する。
【００５１】
【数２９】
2005038397.tif 000011
【００５２】
ここで、LSMEの行列(LSMと呼ぶ)を形成する必要はない。しかし、所望であれば、以下のように計算することができる。
【００５３】
【数３０】
2005038397.tif 000012
【００５４】
ステップ4(IP):中間の式
【００５５】
【数３１】
2005038397.tif 000013
【００５６】
を以下のサブステップを含むブロック後退代入により、w1およびw2について解く。
4.1.対角行列の標準的な方法を使用して、w2についてTTw2=Kを解く。
4.2.対角行列の標準的な方法を使用して、q1についてT・q1=w2を解く。
4.3.WおよびJの構造を利用する方法を使用して、q2=W・f-JTq1を解く。
4.3.バンド行列の標準的な後退代入法を使用して、w1についてPTw1=q2を解く。
【００５７】
ステップ5(IP):最後の式
【００５８】
【数３２】
2005038397.tif 000014
【００５９】
を、サブステップを使用するブロック前進代入法により、mzおよびmについて解く。
5.1.Sの構造を利用する方法を使用して、q3=S・w1を解く。
5.2.バンド行列の標準的な前進代入法を使用して、P・mz=q3を解く。
5.3.TおよびJの構造を利用する方法を使用して、q4についてq4=T-TJ・mzを解く。
5.3.対角行列の標準的な方法を使用して、mについてT・m=w2-q4を解く。
5.4変数pが必要な場合:Wの構造を利用するように修正された標準的な後退代入法を使用して、z=WTmzを解く(pは、異なるzの要素のみである)。
【００６０】
ここで、LSMEを解く。変数mおよびz(またはmz)は、IP法の探索法に対する出力である。
【００６１】
ステップ3(AS):次のセクションで記述する一般化QRアルゴリズムを使用して、以下の一般化QR因数分解を実施する。
【００６２】
【数３３】
2005038397.tif 000015
【００６３】
右辺の角括弧の行列は、入力である。他の行列は、出力である。ここで、Uを計算する必要はない。Pは、上方3角である。残念ながら、LTのバンド構造は、転置されたときに失われている。したがって、QR因数分解は、計算を保存するためにバンド構造を利用することができず、Pはバンドではない。Uが正規直交ではないので、これは、新規な一般化因数分解である。代わりに、Uは、より一般的な特性を有する。
UTS・U=S
【００６４】
L-TJTは、対角より下にある非ゼロのエントリのバンドを除いて、上方3角行列である。一般化QR因数分解は、このバンドより下の下方3角部分がすでにゼロであるということを利用するはずである
【００６５】
行列LSMを決定する必要はない。しかし、以下のように計算することができる。
【００６６】
【数３４】
2005038397.tif 000016
【００６７】
行列Tは、AS法ではゼロである。
【００６８】
ステップ4(AS):中間の式
【００６９】
【数３５】
2005038397.tif 000017
【００７０】
を以下のステップを使用して、ブロック後退代入によりw1およびw2について解く。
4.1.標準的な後退代入法を使用して、w2についてLT=w2=W・fを解く。
4.2.Sの構造を利用する方法を使用して、q1=S・w2を解く。
4.2.標準的な前進代入法を使用して、q2についてL・q2=q1を解く。
4.3.Jのバンド構造を利用する方法を使用して、q3=(K+J・q2)を解く。
4.4標準的な.前進代入法を使用して、w1についてPTw1=q3を解く。
【００７１】
ステップ5(AS):最後の式
【００７２】
【数３６】
2005038397.tif 000018
【００７３】
を、サブステップを使用して、ブロック前進代入法によりmおよびmzについて解く。
5.1.Srの構造を利用する方法を使用して、q4=-SR・w1を解く。
5.2.標準的な後退代入法を使用して、mについてP・m=q4を解く。
5.3.Jの構造を利用する方法を使用して、q5=JTmを解く。
5.4.標準的な後退代入法を使用して、q6について、LTq6=q5を解く。
6.Sの構造を利用する方法を使用して、q7=S(w2-q6)を解く。
6.標準的な前進代入法を使用して、mzについてL・mz=q7を解く。
7.変数pが必要である場合:Wの構造を利用するように修正された標準的な後退代入法を使用して、z=WTmzを解く。
【００７４】
一般化QR因数分解アルゴリズム
新規な一般化QR因数分解が、ステップ3(ASまたはIP)において必要である。形態P=U*Pインの従来のQR因数分解を決定するために、周知のアルゴリズムが利用可能である。Pインは、入力行列であり、Uは、UTU=Iの特性を有し、Pは、上方3角または下方3角である。記述する方法は、ハウスホルダ変換のシーケンスを使用する既知のQR法の族の新規な拡張である。
【００７５】
平方根法のステップ3の因数分解は、Sr1/2P=S1/2U*Pインを必要とする。PインおよびSは入力であり、SrおよびSは対角エントリが1または-1である対角行列であり、Pは下方3角または上方3角であり、Srは出力であり、UはUTS*U=Sの特性を有する。Pインが実数である場合、Pも実数であることに留意されたい。Uは、複素数成分を含む可能性がある。この因数分解は、以下の変換がハウスホルダ変換の代わりに使用される場合、ハウスホルダ変換のシーケンスを使用する周知のQR因数分解法と同じ方式で実施することができる。
【００７６】
【数３７】
2005038397.tif 000019
【００７７】
上式で、eiは、ゼロの列ベクトルであるが、i番目の位置にある1のエントリを除く。xiおよびeiは、従来のQR因数分解と同じ手続きに従って確定される。行列およびベクトルを使用した上記の記述は、単に説明を目的としている。計算線形代数の当業者なら、これを実施するアルゴリズムは、Uおよびvを実際には形成しないことによって計算を保存することができることを認識している。そのような当業者なら、変換のシーケンスを構成して、入力行列に基づいて各xを確定し、変換を適用して下方3角または上方3角の出力行列Pを生成する方式をも認識している。出力行列Srは、対角行列であり、各要素は、対応するswとSの要素との積である。さらなる説明上の補助として、Sが恒等行列に等しい場合、すべてのswは1であり、新しい変換は、従来のハウスホルダ変換になる。ハウスホルダ変換の代入以外に、ステップ3(ASおよびIP)において因数分解を生成するために、従来のQR因数分解アルゴリズムを使用することができ、また、ステップ3IPにおいて、随意選択である構造を利用するために使用することができる。ステップ3IPでは、Pインの構造は、恒等行列に等しいSrとなり、したがってSrを省略することができることに留意されたい。
【００７８】
平方根因数分解の理論的根拠
一般化QR因数分解の理論
一般化QR因数分解において確定されたUは、以下の特性を有する。
【００７９】
【数３８】
2005038397.tif 000020
【００８０】
Uは、以下の特性も有する。
【００８１】
【数３９】
2005038397.tif 000021
【００８２】
上式で、xiは、vの確定に使用された。
【００８３】
したがって、この変換は、σに設定された1つを除いて、列行列のすべての要素をゼロにする。行列を3角形態に変換するためにこの後者の特性を使用することは、周知である。これが実施されたとき、PインからPへの変換のシーケンスが存在する。
【００８４】
【数４０】
2005038397.tif 000022
【００８５】
シーケンスのUおよびeについて確定された構造により、これは、以下の形態に再構成することができる。
【００８６】
【数４１】
2005038397.tif 000023
【００８７】
SrおよびSwは、Sr-1/2S1/2=Sw1/2として関係付けられる。SwをSrで置き換えることにより、以下を得る。
【００８８】
【数４２】
2005038397.tif 000024
【００８９】
上記で、最後の式は、その左側の先行する各辺にそれぞれの転置を乗算した結果である。最後の式は、これが所望の変換であることを示す。
【００９０】
Hを因数分解する理論
Hの因数分解の理論は、下式
【００９１】
【数４３】
2005038397.tif 000025
【００９２】
を得るように、式1の両辺に各時間点のそれぞれの転置を乗算することによって形成された恒等式に基づく。
【００９３】
これらの恒等式は、制御理論の教科書において既知である。HおよびLSMを因数分解するためにそれを使用することが新規である。Hの因数分解は、式2および3においてL、S、およびWの定義を使用して、W-1LTSLW-Tの表現を記号により展開して、次いでこれらの恒等式を使用することによって実証することができる。
【００９４】
IP法についてLSMEを因数分解する理論
式5の因数分解は、
【００９５】
【数４４】
2005038397.tif 000026
【００９６】
を意味する。これを再構成して、必要な関係を得る。
W-1PTS・P・W-T-JTT-tT-TJ=W-1LTSL・W-T=H
【００９７】
式8を式7に代入して、S-1=Sを使用する場合、以下を得る。
【００９８】
【数４５】
2005038397.tif 000027
【００９９】
W-1によって列の第1ブロックをスケーリングして、mzに代入すると、以下を得る。
【０１００】
【数４６】
2005038397.tif 000028
【０１０１】
行列の乗算を記号により実施すると、次のようになる。
【０１０２】
【数４７】
2005038397.tif 000029
【０１０３】
これは、IPの平方根アルゴリズムが、式5の一般化QR因数分解を実施するアルゴリズムが機能することを想定して、LSMEを解くことを示す。
【０１０４】
AS法のLSMEを因数分解する理論
式9の因数分解は、以下を意味する。
【０１０５】
【数４８】
2005038397.tif 000030
【０１０６】
式12を式11に代入して、S-1=Sを使用する場合、以下を得る。
【０１０７】
【数４９】
2005038397.tif 000031
【０１０８】
W-1によって行の第2ブロックをスケーリングして、mzに代入すると、以下を得る。
【０１０９】
【数５０】
2005038397.tif 000032
【０１１０】
行列の乗算を記号により実施すると、以下のようになる。
【０１１１】
【数５１】
2005038397.tif 000033
【０１１２】
これは、式9の一般化QR因数分解を実施するアルゴリズムが機能することを想定して、ASの平方根アルゴリズムが、LSMEを解くことを示す。
【０１１３】
平方根法の計算効率
この因数分解の利点は、以下の通りである。
【０１１４】
LSME行列は、上方3角行列と、下方3角行列と、対角行列との積に因数分解される。
【０１１５】
次いで、ブロック後退代入工程でLSMEを解くことができる。ブロック上方3角行列には、ブロック下方3角行列の前進代入が続く。
【０１１６】
LSMおよびHを決定する必要はない。
【０１１７】
L、M、およびPに必要なコンピュータワード長は、計算精度の所与のレベルについて、H、およびLSMに必要とされる精度のわずかに半分である。倍精度の代わりに単精度コンピュータ演算を使用する能力により、計算速度が約2倍になる。
【０１１８】
Hの因数分解は、制御装置の更新につき1回行うだけでよい。QPソルバの各反復は、この因数分解を再使用する。
【０１１９】
QP反復ごとの因数分解の完成は、バンド構造を利用する。IP法では完全に使用し、AS法では部分的に使用する。
【０１２０】
LおよびPの残りの構造を利用する特別なQRおよび3角行列のソルバで、さらなる計算効率を達成することができる。
【０１２１】
Hを因数分解するこのアルゴリズムは、以下の計算上の利点を有する。
【０１２２】
M行列およびF行列は、下方3角であり、これにより、具体的にはNが大きいとき、M行列およびF行列に関する行列演算の計算が大きく低減される。
【０１２３】
L行列およびS行列は、疎バンド行列であり、やはりこれにより、計算が低減される。
【０１２４】
制御理論の「平方根」法の多くを適用することができる。これらのいくつかは、計算をさらに低減するために、行列A、B、C、D、Q、およびR内の構造を利用する。
【０１２５】
Hを因数分解するために必要な計算量は、Nに関して線形にのみ増大する。
【０１２６】
HおよびJの反復構造を利用する。
【０１２７】
図1は、上記で記述したアルゴリズムを実施するための、本発明によるモデル予測制御システム10の概略図である。システム10は、センサ16からのフィードバックとコマンドとに基づいて、作動装置14でプラント12を制御する。システム10は、コマンドに基づいて所望の軌跡を生成する所望軌跡生成装置20を含む。システム10は、上記で記述した方式でKKTの様々な行列およびベクトルを決定するためのモジュール22を含めて、2次プログラミングモジュール21を含む。KKT条件は、センサ16からのフィードバックに基づいたシステム10の現行状態の関数でもある。2次プログラミングモジュール21は、上記で記述したHrを決定するためのモジュール24をさらに含む。修正2次プログラミングソルバ26が、プラント12およびシステム10の状態が所望通りになるように、上記で記述したアルゴリズムに従って、作動装置14に対する作動装置コマンドを生成する。
【０１２８】
図2は、標準的な方式の2次プログラミング(QP)ソルバがモデル予測制御装置(MPC)に埋め込まれる標準的な方式と、新規な平方根法が標準的なQPソルバに埋め込まれる方式とを示す。図2のアルゴリズムは、制御装置が更新されるたびに実行される。図のステップ50は、モデル予測制御を2次プログラミング問題に変換し、ステップ52は、QP解から、現行作動装置コマンドを含めて、次のN制御装置更新について最適な軌跡を抽出する。すべての他のステップが、QPソルバに関連付けられる。一般的なQPソルバは、制約を満たす最適解の構造化された反復探索からなる。これは、ステップ56において、従来の反復から解を分析して、LSMEのバージョンを設定することと(ここで提示した構造化形態では一般的ではないが)、LSMEを解くことと、ステップ62において、結果を分析して、解が必要な精度に達したかを決定することとからなる。図2に示すように、平方根法は、探索が開始される前に、ステップ54においてHの平方根であるHrを事前に計算するための新しい論理を追加する。また、ステップ58においてLSMルートを形成することを含めて、LSMEの解を、ステップ60の計算がより効率的な新規なアルゴリズムと置き換える。
【０１２９】
平方根法は、図2の上部ブロックとQPソルバとのインタフェースを変化させる。ステップ50において、制御下にあるシステムのコマンドを着け取り、それらのコマンドに対する所望の応答を定式化し、システムの現行状態を示すセンサ信号を受信して、KKT条件に従ってQP問題を定式化する。KKT条件の式は、一般に、3つの行列式にまとめられる。1つは、性能指数に関連し、1つは、システムダイナミクスモデルに関連し、1つは、限界条件に関連する。最適軌跡の以前の解は、KKT条件において様々な行列を形成する際に使用することが可能である。標準的な方法では、QPは、次いで、各反復ごとに、異なる構造を有するLSMEのバージョンを形成する。対照的に、平方根法は、KKT条件の様々な変数、式、および行列が、特定の順序および構造のLSMEにまとめられるという点で、ステップ56においてQPのその部分の修正もする。
【０１３０】
特許の法令および法体系の規定によれば、上記で記述した例示的な構成は、本発明の好ましい実施形態を表すと見なされる。しかし、本発明は、その精神または範囲から逸脱せずに、具体的に示し記述したものとは異なる方式で実施することができることに留意されたい。
【図面の簡単な説明】
【０１３１】
【図１】本発明によるモデル予測制御システムの概略図である。
【図２】図1のシステムによって使用される本発明による方法のフローチャートである。
【符号の説明】
【０１３２】
MPC  モデル予測制御
LSME  大規模疎行列式
QP  2次プログラミング
IP  内点
AS  アクティブセット
【図１】
2005038397.tif 000054
【図２】
2005038397.tif 000055