(57)【要約】  非線形、時変、及びパラメータ変動動的プロセスあるいはシステムを支配できる微分または差分方程式を推定する方法及びシステム。これらの方程式を推定する前記方法及びシステムは、観測した出力と、必要に応じて、方程式の入力データとの推定値に基づいてもよい。前記方法及びシステムは、非線形、時変、及びパラメータ変動差分方程式を用いて表せる任意のシステムあるいはプロセスとともに利用できる。また、システム制御、故障検出、状態推定及び予測、並びにシステムまたは方法における変化へのその適応において用いられるシステムあるいは方法の詳細な挙動を記述する際に、差分方程式の自動抽出にも利用できる。【選択図】図５
2013504133.tif 000002
【特許請求の範囲】
【請求項１】
非線形、線形パラメータ変動、及び時変モデルを同定するシステムであって、
線形部分空間アルゴリズムを含み、該線形部分空間アルゴリズムは、第１の出力信号と第１の入力信号の少なくとも一方を入力し、前記第１の出力信号及び前記第１の入力信号は、該第１の入力信号とスケジューリング関数の間のクロネッカー積、状態関数と前記スケジューリング関数の間のクロネッカー積、及び実質的に同等な信号のうちの少なくとも１つと併用される、前記線形部分空間アルゴリズムの第２の出力信号及び拡張入力信号の少なくとも一方を含む、システム。
【請求項２】
パラメータ推定値及び状態推定値を精緻化する前記線形部分空間アルゴリズムの少なくとも１回の反復をさらに含む、請求項１に記載のシステム。
【請求項３】
１回の反復の収束を判定する尤度関数の評価をさらに含む、請求項２に記載のシステム。
【請求項４】
最適な状態次数を、尤度計算及びＡＩＣ計算の少なくとも一方を用いて最適な状態次数を選択する、請求項１に記載のシステム。
【請求項５】
前記システムにおける少なくとも一つの故障あるいは変化を判定する尤度比テストをさらに含む、請求項１に記載のシステム。
【請求項６】
初期化時に、状態関数と、スケジューリング関数と状態関数の間のクロネッカー積との少なくとも一方が存在しない、請求項１に記載のシステム。
【請求項７】
不安定な計算をもたらすデータ利用時における、該データを無視する異常値打ち切りをさらに含む、請求項２に記載のシステム。
【請求項８】
ＡＲＸモデル化計算の監視により前記異常値打ち切りを実現する、請求項７に記載のシステム。
【請求項９】
計算を安定化させる実験計画の利用をさらに含む、請求項１に記載のシステム。
【請求項１０】
非線形、時変、及びパラメータ変動処理を支配する方程式のセットを推定するシステムであって、
第１の入力と、
第２の入力と、
フィードバックボックスと、
時間遅延ボックスと、
を備え、
前記第１の入力及び前記第２の入力は、前記フィードバックボックスを通って前記時間遅延ボックスに至り、出力を生成するシステム。
【請求項１１】
前記システムは、連続する出力からの変動が所定の閾値に達するまで反復する、請求項１１に記載のシステム。
【請求項１２】
前記システムの出力は、前記システムの反復の間、前記第１の入力及び前記第２の入力のうち少なくとも一方として用いられる、請求項１１に記載のシステム。
【請求項１３】
前記システムは、連続する出力からの前記変化が所定の閾値に達するまで反復する、請求項１１に記載のシステム。
【請求項１４】
前記時間遅延ボックスは、所定のサンプル持続時間である、請求項１１に記載のシステム。
【請求項１５】
前記時間遅延ボックスは、複数の線形時不変係数の少なくとも１つを有する、請求項１４に記載のシステム。
【請求項１６】
前記フィードバックボックスは、複数の基底関数、複数のスケジューリング関数、及び複数の状態関数のそれぞれのクロネッカー積のうち少なくとも１つを有する、請求項１１に記載のシステム。
【請求項１７】
前記第１の入力は、複数のスケジューリングパラメータのうちの少なくとも１つである、請求項１１に記載のシステム。
【請求項１８】
第２の入力は、複数の基底関数のうちの少なくとも１つである、請求項１１に記載のシステム。
【発明の詳細な説明】
【背景技術】
【０００１】
測定された出力データ及び入力される可能性のあるデータに基づき、非線形かつ時変な動的プロセスあるいはシステムをモデル化する技術は、新たな技術分野である。この技術は、その応用や理論の分野によるが、統計学では時系列分析、機械工学ではシステム同定、心理学では縦断的解析、財務分析では予測とそれぞれ呼ばれることがある。
【０００２】
従来、部分空間システム同定法が導入され、フィードバックに関するシステムの最適な方法、双線形システムを含む非線形システムの方法の研究、及び線形パラメータ変動（ＬＰＶ）システムなどにおいて著しく発展し改良がなされてきた。部分空間方法は、収束しない反復非線形パラメータ最適化を回避し、高次・大規模システムに有用な、数値的に安定な方法を用いる。
【０００３】
時変非線形システムの分野においては、必ずしも望ましい結果が得られるとは限らない状況で研究が続けられてきた。非線形システムのモデル化に線形部分空間法の直接的拡張を用いているという点で、この研究は現在の最先端技術の典型である。この手法は、過去と未来を、過去の入力及び出力の非線形関数の線形結合として表現する。この手法の結果の１つとして、過去と未来の次元が、測定され利用される過去の入力、出力、状態、及び遅延において指数関数的に増加する。これらの変数の一部のみを用いる場合、過去の次元は１０４以上、あるいは１０６以上に達することもある。典型的な産業プロセスでは、過去の次元は１０９あるいは１０１２を優に超える可能性がある。このような極端な数字からは、非効率な開発または結果しか得られない。
【０００４】
他の技術では、モデルにおける非線形項を推定するための反復部分空間手法を用いるため、必要とされる計算量は極めて少ない。この手法は発見的アルゴリズムを含み、ランダムスケジューリング関数、すなわちホワイトノイズ特性を有するＬＰＶシステムの場合における高精度なモデル同定に用いられてきた。しかしながら、大抵のＬＰＶシステムにおいては、スケジューリング関数が通常特定のアプリケーションにより決定され、特徴がしばしば極めて非ランダムであるという問題がある。非ランダムなスケジューリング関数の場合の精度を高める試みがなされたいくつかの変形例においては、結果的にモデル化の精度を実質的に高めることはできなかった。
【０００５】
さらに一般的な状況において、非線形システムの同定は、一般的非線形正準変量分析（ＣＶＡ）法として知られる一般的な問題がある。この問題は、簡単な非線形差分方程式によって表されるカオス的非線形システムであるローレンツアトラクタにより例証された。それゆえ、過去と未来の非線形関数は、プロセスの状態を表すように決定され、プロセスの状態は、同様にシステムの非線形状態方程式を表すために用いられる。この方法では、過去と未来の必要とされる非線形関数の数が、周知の通り指数関数的に増加するため、実行可能な計算の実現が大きな難題の１つとなっている。一般的非線形システムにも当てはまる、システム同定の問題の解決策を見出す際にも、しばしばこのような難題に直面する。
【０００６】
それゆえ、後述する例示的な実施形態においては、観測の「大規模サンプル」が得られる場合に著しい改善がなされ、また最適な結果を生み出すことのできる方法及びシステムを表現できる。更には、この方法は「場当たり的な」方法ではなく、最適な統計的方法である。
【発明の概要】
【０００７】
例示的な一実施形態は、非線形、時変、及びパラメータ変動動的プロセスを利用する方法を説明する。該方法は、時変要素を含むシステムの縮小モデルを生成するために用いられてよい。前記方法は、状態空間差分方程式を展開するステップと、出力及び拡張入力について、差分方程式を線形時不変システムとして表現するステップと、状態方程式の係数を推定するステップと、を備えてもよい。
【０００８】
他の例示的な実施形態において、非線形、時変、及びパラメータ変動プロセスを支配する方程式のセットを推定するシステムを説明してもよい。該システムは、第１の入力と、第２の入力と、フィードバックボックスと、時間遅延ボックスと、を含んでもよい。また、前記システムにおいて、前記第１の入力と前記第２の入力は、前記フィードバックボックスを通って前記時間遅延ボックスに至り、出力を生成する。
【発明を実施するための形態】
【０００９】
本発明の態様を、以下の説明及び本発明の具体的な実施形態に関する関連図面に開示する。代替の実施形態は、本発明の精神及び請求項の範囲から逸脱せずに考案され得ることは当業者に明らかであろう。また、本発明の例示的な実施形態における周知の要素は、本発明の関連する詳細を不明瞭にすることのないよう、詳細に説明されないか、または省略される。
【００１０】
本明細書で用いる「例示的」という用語は、「例、事例、または実例として」を意味する。ここに説明されているいずれの実施形態もこれに限定されず、単なる例示にすぎない。なお、ここで説明されている実施形態は、必ずしも他の実施形態よりも望ましい、または有利と解釈されない。また、「本発明の実施形態」、「実施形態」、または「本発明」という用語は、発明の全ての実施形態が検討された特徴、利点または動作モードを含むことを必要としない。
【００１１】
さらに、本発明の多くの実施形態は、たとえばコンピューティングデバイスの要素によって実行されるべき一連の行為に関して説明する。ここで説明する様々な行為は、特定の回路（たとえば、特定用途向け集積回路（ＡＳＩＣ））によって、１つまたは複数のプロセッサによって実行されるプログラム命令によって、あるいはその両方の組合せによって、実行できることを認識されよう。加えて、ここで説明するこれら一連の行為は、実行すると、ここで説明する機能を関連するプロセッサに実行させるコンピュータ命令の対応するセットを記憶した任意の形態のコンピュータ可読記憶媒体内で全体として実施されるものとみなすことができる。したがって、本発明の様々な態様は、すべてが請求する主題の範囲内に入るように熟慮された、いくつかの異なる形態で実施することが可能である。加えて、ここで説明する実施形態の各々について、そのような実施形態の対応する形態は、ここでは、たとえば、説明する行為を実行する「ように構成された論理」として説明することがある。
【図面の簡単な説明】
【００１２】
（図面の簡単な説明に対応する記載なし）
例示的な図１〜図６を広く参照すると、時変、パラメータ変動、及び非線形差分方程式の経験的モデル化のための方法及びシステムについて説明できる。該方法及びシステムは、様々な結果を提供するように実施及び利用することができ、効率的に実施できる。
【００１３】
例示的な図１に、例示的な実施形態に係る時変、パラメータ変動、及び非線形差分方程式の経験的モデル化のための方法論のフローチャートが示される。ここで、１０２において、時変及びパラメータ変動差分方程式、ならびに、必要に応じて、非線形状態空間差分方程式のセットを用いてもよい。１０４において、これらの方程式は、選択した１組の基底関数について展開してもよい。例えば、非線形入出力方程式は、xt及びutの多項式に展開してもよい。次に、１０６において、それぞれの差分方程式は、出力yt及び拡張された入力ut（入力utに加え、入力utと、スケジューリング関数ptと、状態xtとにおける例えば多項式の基底関数を含む）に関して、例えば、線形時不変として表現されてもよい。
【００１４】
例示的な図２に、差分方程式の線形パラメータ変動システムを利用できる例示的なフローチャートが示される。本実施形態では、２０２において、下記に示す式１及び式２に示すような線形パラメータ変動状態空間方程式のセットを用いてよい。
【数１】
2013504133.tif 000003
【数２】
2013504133.tif 000004
【００１５】
そして、２０４において、例えば次式３及び４に示すような、スケジューリング関数pt、状態xt、及び入力utにおける多項式に関して、状態空間差分方程式を展開してもよい。
【数３】
2013504133.tif 000005
【数４】
2013504133.tif 000006
【００１６】
2013504133.tif 000007
【数５】
2013504133.tif 000008
【数６】
2013504133.tif 000009
【００１７】
2013504133.tif 000010
【００１８】
例示的な図３に、反復部分空間同定のために実施される反復アルゴリズムのフローチャートが示される。本実施形態において、３０２に示すように、非線形差分方程式は、付加的基底関数に展開され、また拡張入力utを用いて線形時不変形式として表現できる。いくつかの例では、これは、出力、状態、及びスケジューリング関数に関する非線形基底関数を含む。
【００１９】
2013504133.tif 000011
【００２０】
2013504133.tif 000012
【００２１】
３０８において、収束を確認してもよい。ここで、対数尤度関数の値の変化と、反復k-1及び反復kにおける状態次数を比較してもよい。いくつかの例では、状態次数は同一であり、対数尤度関数の変化は選択された閾値εよりも小さい。この閾値は、多くの場合、１よりも小さく、例えば０．０１であり、その後に反復は終了あるいは停止されてもよい。その他の場合においては、その値は選択された閾値εより大きい場合、上記ステップ３０６に戻り、反復k+1を行ってもよい。反復k+1を行った後、収束を再び確認してもよい。
【００２２】
別の実施形態においては、別の手法を用いて、例えば、線形時不変システムのために開発された部分空間方法を用いるシステムで自己相関された誤差及びフィードバックの場合において、未知のパラメータの最適あるいは所望の推定値を直接的かつ簡単に求めてもよい。これは、状態推定値について反復の要望につながる可能性のある、異なる形式の問題を表現することによって成し得る。しかしながら、反復回数は極めて少なく、またシステム及びその開発をさらに簡素化することにより、確率的ノイズが除去されてもよい。
【００２３】
例えば、システム行列が、下記式に示すような形式のスケジューリングパラメータρt=(ρt(1) ρt(2) … ρt(s))Tのベクトルの時変関数である線形システムについて考える。
【数７】
2013504133.tif 000013
【数８】
2013504133.tif 000014
【００２４】
スケジューリングパラメータのアフィン依存性に関しては、状態空間行列は、下記式９〜１２のような形式を有してもよい。
【数９】
2013504133.tif 000015
【数１０】
2013504133.tif 000016
【数１１】
2013504133.tif 000017
【数１２】
2013504133.tif 000018
【００２５】
上記式において、左辺の行列は、アフィン形式と呼ばれることもあるスケジューリングパラメータρtによって特定される右辺の線形結合(1 ρt(1) ρt(2) … ρt(s))Tとして表現されることに注目してもよい。さらに説明すると、この表記法は、システム行列A=[A0 A1 … As]=[A0 Aρ]に用いられ、B、C、Dについても同様に用いられる。
【００２６】
さらに他の例示的な実施形態において、ＬＰＶシステムのクラスのためのシステム同定法は、多くの適用可能性と経済的価値がある。そのようなシステムの例としては、これらに限定されないが、例えば、航空機や船舶のような、空気力学的あるいは流体力学的な車両、自動車エンジン力学、タービンエンジン力学、及び、攪拌槽型反応器や蒸留塔などの化学プロセスが挙げられる。特徴の１つとしては、ある任意の操作点ρtにおいて、システム力学は線形システムと称されることがある。スケジューリングパラメータρtは、操作点変数の複合非線形関数であってもよく、その例としては、これに限定されないが、速度、圧力、温度、燃料流量などが挙げられ、これらは、未知の可能性のある定数行列A、B、C、Dにおけるシステム力学を特徴付ける、既知あるいは正確に測定された変数であってもよい。そのような任意の操作点変数を知ることによって、ρtを計算可能あるいは決定可能であるものとしてもよい。例えば、自動車エンジンのＬＰＶモデルは、ベクトルρtの要素を様々な操作点変数の非常に複雑な非線形関数として明示するＬＰＶ状態空間方程式を含んでもよい。本明細書に記載の例示的な実施形態では、システム同定計算を行う際に、スケジューリングパラメータρtを利用できさえすればよい。これにより、リアルタイム利用、あるいはリアルタイム制御やフィルタリングのような用途に対する要求条件が緩和される。
【００２７】
以上の説明を簡単にまとめると、下記式に示すように、スケジューリングパラメータρtを入力ut及び状態xtに関連付けることにより、ＬＰＶ方程式を時不変形式で表すことできる。
【数１３】
2013504133.tif 000019
【数１４】
2013504133.tif 000020
【００２８】
2013504133.tif 000021
【数１５】
2013504133.tif 000022
【数１６】
2013504133.tif 000023
【００２９】
2013504133.tif 000024
【００３０】
また、ここでは、上記の式を、下記式１７及び１８のように記述してもよい。
【数１７】
2013504133.tif 000025
【数１８】
2013504133.tif 000026
【００３１】
2013504133.tif 000027
【００３２】
さらに他の例示的な実施形態では、ＬＰＶシステムは、既知のパラメータ変動関数ρtを含んでもよい非線形内部フィードバックによる線形時不変システムとして表現してもよい（Vincent Verdult, Nonlinear System Identification: Section 2.1: A State‐Space Approach、2002, Ph. D. Thesis, University Of Twente, The Netherlandsを参照。この全内容は参照により本明細書に組み込まれている。）。この例示的な実施形態において、階数riのシステム行列Piは、式１９に示されるように、1≦i≦sとなる各iについて、特異値分解を用いて因数分解されてもよい。
【数１９】
2013504133.tif 000028
【００３３】
さらに、式２０〜式２３に示されるように、その数量は次のように定義できる。
【数２０】
2013504133.tif 000029
【数２１】
2013504133.tif 000030
【数２２】
2013504133.tif 000031
【数２３】
2013504133.tif 000032
【００３４】
2013504133.tif 000033
【数２４】
2013504133.tif 000034
【数２５】
2013504133.tif 000035
【数２６】
2013504133.tif 000036
【数２７】
2013504133.tif 000037
【００３５】
上記式に関して例示的な図４を参照し、図４及び式２６に示すフィードバック構成のように、任意の時刻tにおけるフィードバックftの出力ztへの影響がないと仮定できるとすると、線形分数変換（ＬＦＴ）記述におけるパラメータ依存性は存在しない（K. Zhou, J. Doyle, and K. Glover (1996), Robust and Optimal Control, Prentice-Hall、Inc., Section 10.2を参照。この全内容は参照により本明細書に組み込まれている。）。図４に示したように、ボックス４０４及びボックス４０２の２つのボックスを有するシステム４００を備えてよい。ボックス４０４は、メモリレス非線形フィードバックシステムであってよく、ボックス４０２は線形時不変システムであってよい。したがって、ある例示的な実施形態では、パラメータ依存性はアフィンとなり得る。
【００３６】
2013504133.tif 000038
【００３７】
また、下記に示すように、これは、スケジューリングパラメータベクトルρにおいて状態方程式が線形に成り得る式１５及び１６に示すようなＬＰＶ形式に等しい。
【００３８】
2013504133.tif 000039
【数２８】
2013504133.tif 000040
【数２９】
2013504133.tif 000041
【００３９】
続いて、上記の因数を定義するために式１９を用いる場合、式２８及び式２９は、式１５及び式１６と同様であってよい。
【００４０】
2013504133.tif 000042
【数３０】
2013504133.tif 000043
【数３１】
2013504133.tif 000044
【００４１】
2013504133.tif 000045
【００４２】
2013504133.tif 000046
【００４３】
別の例示的な実施形態において、既知のスケジューリング関数とＬＴＩサブシステムの状態に関する非線形フィードバックシステムのＬＴＩサブシステムの実施に続き、ＬＴＩシステム行列及び状態ベクトルを決定してもよい。本実施形態は、ＬＴＩシステムを記述するＬＴＩ状態空間行列及びＬＴＩシステム状態の両者について反復決定法を用いてよい。
【００４４】
2013504133.tif 000047
【００４５】
2013504133.tif 000048
【００４６】
2013504133.tif 000049
【００４７】
2013504133.tif 000050
【００４８】
2013504133.tif 000051
【００４９】
さらなる例示的な実施形態において、追加入力の影響を、反復アルゴリズムによってトレースしてもよい。さらに、以下の例示的なステップは、例示的な図５に示す表や、例示的な図６におけるフローチャートに反映されてもよい。
【００５０】
2013504133.tif 000052
【００５１】
例示的なステップ６０４において、訂正された未来を計算してもよい。未来出力に対する未来入力の影響は、ＡＲＸモデルを用いて決定し、出力から差し引いてもよい。この影響は、時刻tの未来において入力がないと仮定すると、時刻tの状態から求められる未来出力の計算であってもよい。
【００５２】
例示的なステップ６０６において、正準変量分析（ＣＶＡ）を実行あるいは計算してもよい。ここで、過去と訂正された未来とにおけるＣＶＡを計算してもよい。さらに、過去と訂正された未来とにおける共分散行列を計算してもよい。これは、ＡＲＸモデルを求めるための、過去の共分散の次数からなる過去・未来の結合共分散行列におけるＳＶＤに類似すると言えよう。このステップの結果、「メモリ」と呼ばれるシステムの状態の推定値が得られる。
【００５３】
例示的なステップ６０８において、状態方程式を用いて回帰を実行してもよい。ステップ６０６の「メモリ」はあたかもデータのように用いることができ、ノイズ過程の状態空間行列及び共分散行列の結果としての推定値が得られる。これらの推定値は、ＡＲＸやＳＳモデルのパラメータ値について事前仮定を行わない状態空間モデルの漸近的なＭＬ推定値であってもよい。
【００５４】
2013504133.tif 000053
【００５５】
再びステップ６０２を参照すると、ＡＲＸモデルフィッティングにおいて、拡張出力の次元は、dimuからdimua＝dimu+dimp(dimx+dimu)へと増加する。ＡＲＸモデルをフィットすると、同定されたＡＲＸ次数lagpは、非線形入力項に起因し、また出力及び拡張入力変数に見られる統計的に有意な動的特性に依存して、かなり高次となり得る。この計算は、lagpが考えられる最大のＡＲＸ次数となる場合のdimua*lagpの次元であるデータ共分散行列に関するＳＶＤ計算を含んでもよい。ＳＶＤに利用できる計算は、60*(dimua*lagp)3の次数からなるため、(dimp(dimx+dimu)/dimu)3に比例して増加する。
【００５６】
2013504133.tif 000054
【００５７】
2013504133.tif 000055
【００５８】
別の例示的な実施形態においては、ＬＰＶの場合の反復アルゴリズムの収束について説明できる。本実施形態では、非線形差分方程式の他の形式にほぼ同様の手法を用いてよい。
【００５９】
2013504133.tif 000056
【００６０】
さらに、ここで説明する方法及びシステムは、ＥＭアルゴリズムに照らして説明してもよい。それは、両者間にある平行処理が存在する可能性があるためである。ＬＰＶアルゴリズムの収束を説明するため、以下、ＧＮと記し、その全内容が参照により本明細書に組み込まれる、Ｇｉｂｓｏｎ及びＮｉｎｎｅｓｓ（２００５）における開発（S. Gibson及びB. Ninness、「Robust maximum-likelihood estimation of multivariable dynamic systems」、Automatica、第４１巻、Ｎｏ．５、１６６７〜１６８２ページ、２００５）について、ＬＰＶアルゴリズムに適した様々な変形例を用いて説明する。Ｇｉｂｓｏｎ及びＮｉｎｎｅｓｓ（２００５）から引用されたすべての式番号は、例えば、（２３ＧＮ）のように、ＧＮ論文中に含まれる番号の次に記されるＧＮを含む。
【００６１】
この例示的な実施形態において、ＬＰＶアルゴリズムを求めるＧＮの説明において、次のような置き換えが可能である。すなわち、ＬＴＩ状態方程式をＬＰＶ状態方程式に置き換え、また、欠測データについて、カルマンスムーザによる状態推定値を、反復アルゴリズムにおける「メモリ」ベクトルmtに置き換える。この結果は有意となり得る。何故ならば、ここに述べる反復アルゴリズムによれば、ＧＮの場合と同様に、線形システムに対して、大規模サンプルにおいて１ステップでグローバルなＭＬパラメータ推定値を得ることができるためである。一方、線形システムに対して、ＭＬ解を求めるためにＥＭアルゴリズムを複数回反復してもよい。
【００６２】
2013504133.tif 000057
【００６３】
2013504133.tif 000058
【００６４】
したがって、代わりにＬＰＶアルゴリズムを用いることにより、高速収束を得てもよい。よって、本実施形態においては、上記「欠測データ」をカルマンスムーザの推定値とする代わりに、部分空間手法によってＣＶＡ状態推定値mtまたは「メモリ」を欠測データとしてもよい。メモリmtは、図５に明記された入出力ベクトルを用いた反復アルゴリズムの例示的なステップ６０６において求められる訂正された未来及び過去の正準変量分析によって得られた状態ベクトルの推定値である。これは、最後に推定されたパラメータ値におけるカルマンスムーザ状態推定値というよりは、むしろ反復kにおける出入力データに関連するグローバルなＭＬパラメータ推定値におけるカルマンフィルタ状態推定値に類似してもよい。差異があるとすれば、以下のような理由からである。すなわち、ＣＶＡ方法は、メモリmtの推定値が出入力データのみに依存し得るように過去を必要条件として訂正された未来に関して尤度関数を表現し、またこの出入力データの分布は、状態の平滑化された推定値よりも、むしろ反復kにおいて用いられる出入力データに関連するＭＬ推定値に依存しているからである。ξtT=[xt+1TytT]及びztT=[xtTutT]による（２６ＧＮ）の実際の条件付き尤度関数pθ(ξt/zt)は、すべてのＥＭ計算に含まれ得るものである。これは、ＧＮのＬｅｍｍａ３．３のように、状態空間方程式のパラメータを推定するＣＶＡアルゴリズムの例示的なステップ６０８に関する同様の尤度関数であってもよい。相違点としては、ＣＶＡアルゴリズムの例示的なステップ６０８においては、反復kにおける出入力データに関連する真にグローバルなＭＬ推定値についての期待値であってよいのに対し、ＧＮ推定値は、前回の反復で得られたパラメータ値についての期待値である点である。
【００６５】
2013504133.tif 000059
【００６６】
2013504133.tif 000060
【００６７】
2013504133.tif 000061
【００６８】
2013504133.tif 000062
【００６９】
さらに別の例示的な実施形態では、ＬＰＶシステム同定アルゴリズムが最大尤度解の近くまで、幾何学的に収束することを例証してもよい。ここで、線形システムの結果を展開してもよい。同様の手法はＬＰＶシステムにも有効であるが、下記に示す式は時間に依存し、非常に複雑になる可能性がある。
【００７０】
2013504133.tif 000063
【数３２】
2013504133.tif 000064
【数３３】
2013504133.tif 000065
【００７１】
2013504133.tif 000066
【数３４】
2013504133.tif 000067
【００７２】
次に、xtに関する式３４の右辺を再帰的に代入することにより、下記式３５が得られる。
【数３５】
2013504133.tif 000068
【００７３】
2013504133.tif 000069
【００７４】
2013504133.tif 000070
【００７５】
2013504133.tif 000071
【数３６】
2013504133.tif 000072
【００７６】
2013504133.tif 000073
【００７７】
他の例示的な実施例では、反復線形部分空間計算の収束は、ＬＰＶ差分方程式の安定性、より具体的には、ここに述べるＬＰＶ線形部分空間システム同定の安定性の影響を受けることがある。状態遷移行列の固有値のすべてが安定的であれば、時不変線形状態空間差分方程式のセットは安定的となり、安定的となるのはその場合に限る。そのため、固有値は、例えば１未満である。ＬＰＶの場合はさらに複雑であるが、本発明の例示的な実施形態の目的として、サンプル時間ごとの拡張率または縮小率が、式９から得られたＬＰＶ状態遷移行列の各固有値により状態ベクトルxtのそれぞれの固有ベクトル成分に対して与えられる。ここでは、下記式３７のように表される。
【数３７】
2013504133.tif 000074
【００７８】
式３７において、行列A=(A0 Aρ)=(A0 A1 … As)を未知定数行列と見なすことができる。したがって、遷移行列A(ρt)は、0≦i≦sに対する行列Aiの線形結合[1; ρt]であってもよいことが明らかである。よって、0≦i≦sに対する行列Aiを任意に選択することで、特定のサンプル時間tにおいて不安定な固有値を生成可能なρtの可能値が示される場合がある。このことが散発的にのみ及び／またはごく限られた連続した観測時間だけ起きると仮定する、あるいはそれが起こる場合、いかなる問題や誤差も発生しない可能性がある。
【００７９】
2013504133.tif 000075
【００８０】
さらに他の例示的な実施形態では、不安定性は、予測システム応答が急速に大きくなり得る期間を生じさせることがある。例えば、固有値がすべて１未満で有界である場合、予測システム応答は有界であり得る。これに対して、推定誤差が大きいと、例えば、１つの固有値は、ある期間において２に等しくなり、さらにその期間において、予測応答によってほぼすべてのサンプルが２倍となり得る。したがって、０３０＝（１０３）１０≒（２１０）１０＝２１００とすると、無意味な結果をもたらし得る計算において、わずか１００サンプルで３０桁の精度を損なう可能性がある。
【００８１】
したがって、上記によれば、ＬＰＶ遷移行列が、例えば、延長した時間に対して不安定であるような期間に関連し得る数値的精度が著しく損なわれる状況も有り得る。また、アルゴリズムは、大規模サンプルの最適解において安定及び収束するため、アルゴリズムの初期化及びサンプルサイズが問題となる場合がある。無限の精度で計算が可能であると仮定すると、不良条件の問題は避けられるが、小数第１５位または３０位の精度では、例えば航空機翼のフラッタのある実際のデータセットは、さらなる検討が有益であろう。
【００８２】
2013504133.tif 000076
【００８３】
他の例においては、例えば、スケジューリングパラメータ値がA(ρt)において重大な不安定性をもたらす期間のように、ft|qtのある成分、つまり、未来入力qtの効果に対して補正された未来出力ftが極めて大きくなることがある。したがって、共分散行列ft|qtの計算において、そのような大きな数値は、重大な不良条件及び数値的精度における損失をもたらすことがある。
【００８４】
不安定性が発生する時間を反復アルゴリズム及び他の部分空間方法を用いて任意に編集できる例においては、不安定領域の異常値の編集を行ってもよい。例えば、訂正された未来を計算できるＡＲＸ共分散計算において、重大な不安定性を含む期間を判定し、計算から除外してもよい。
【００８５】
実験計画を扱う例においては、操作点変数の軌道が特定できれば、いくつかの方法でスケジューリングパラメータをスケジューリングし、システム同定を向上させることができる。計算上の不安定性を回避できる初期領域を選択し、ＬＰＶパラメータ(A, B, C, D)の十分正確な初期推定値を求めてもよい。さらに、これは、計算上の不安定性をもたらすことのない十分な精度でアルゴリズムの状態推定値を初期化するために、他の領域において利用されてもよい。
【００８６】
2013504133.tif 000077
【００８７】
さらに、ここで説明した方法及びシステムの多数の例示的な潜在的応用例においては、産業上の問題点について今のところ実現不可能な既存の方法と比較して、提案されたアルゴリズムをより好ましく実施することが期待できる。また、多くの場合、可能な限り短時間で、低コスト、及び低資源で、作動空間の特定されたグローバルな領域に対して所望の忠実度の結果を得られる実験を設計することが望まれる。反復アルゴリズムの線形部分空間法はプロシージャーに基づいた最大尤度であるため、ＬＴＩシステム同定に関して設計開発してもよい。また、それは、動的周波数応答関数のような数量についての信頼帯を含む、推定外乱モデルにより確率論的モデルを同定するため、外乱プロセスについての少ない事前情報を用いて、必要なサンプルサイズ及びシステムの入力励起を開発できる。
【００８８】
さらに他の例示的な実施形態において、ここで説明したＬＰＶ方法及びシステムは、非線形システムに拡張してもよい。例えば、いくつかの複雑な非線形システムを、ここで説明したＬＰＶ部分空間システム同定法を適用するのに十分な略ＬＰＶ形式で表現できることを示してもよい。
【００８９】
一例では、一般的非線形、時変、及びパラメータ変動動的システムは、式３８及び式３９に示すような非線形状態方程式のシステムにより記述できる。
【数３８】
2013504133.tif 000078
【数３９】
2013504133.tif 000079
【００９０】
式３８及び式３９において、xtは状態ベクトル、utは入力ベクトル、ytは出力ベクトル、νtはホワイトノイズ測定ベクトルであってよい。ある例示的な実施形態では、「パラメータ変動」システムに対応するため、時変パラメータであり得る「スケジューリング」変数ρtは、システムの現在の操作点を表してもよい。関数f(・)及びh(・)のごく一般的なクラスは、必ずしも連続する必要のない加法的ボレル関数により表してもよい。
【００９１】
簡単な方法では、Ｒｕｇｈ、Ｓｅｃｔｉｏｎ  ６．３（W. J. Rugh、Nonlinear System Theory: The Volterra/ Wiener Approach. Baltimore、MD: Johns Hopkins大学出版、１９８１を参照。この全内容は参照により本明細書により組み込まれている。）、ここで、ρt及びνtは省略してもよい、に記載のテーラー展開を採用する関数により、下式４０及び４１が与えられる。
【００９２】
【数４０】
2013504133.tif 000080
【数４１】
2013504133.tif 000081
【００９３】
2013504133.tif 000082
【００９４】
このように、式４０及び式４１は、非線形関数f(・)及びh(・)の多項式展開であってもよい。なお、非線形方程式は、様々な目的のために容易に計算可能な積の和のみを含む近似多項式のような、比較的簡単な形式の非線形関数を含んでもよい。しかしながら、自己相関された誤差の存在下における係数の経験的推定については、たとえ部分空間法を用いたとしても、y、u、及びxが低次元であると、問題が困難となり得る。部分空間法では、利用できる「過去」の次元とともに行列次元が指数関数的に増大することがある。これは、式４０の左辺のxtについて状態方程式４０の右辺のxt[k]へ代入を繰り返すことにより式４０が展開される際に起こり得る。したがって、tをt-1で置き換えた式４０の右辺全体は、推定ＡＲＸモデル指数として過去に通常通り選択された次数である冪指数lagpまで累乗してよい。比較的低次の過去、及び低次元のxt及びutについて、加法的項数は指数関数的に増加してもよい．
【００９５】
しかしながら、さらなる例示的な実施形態及び後述のＲｕｇｈにおいては、式３９及び式４０は、カーレマン双線形化により、下記式４２に示す状態変数における双線形ベクトル微分方程式へ変換してもよい。
【数４２】
2013504133.tif 000083
また、入力指数及び積変数は下記式４３のように表される。
【数４３】
2013504133.tif 000084
【００９６】
さらに他の実施形態では、状態・アフィン形式を表す式４０は、下記式４４及び４５のように書き換えてもよい。
【数４４】
2013504133.tif 000085
【数４５】
2013504133.tif 000086
【００９７】
2013504133.tif 000087
【数４６】
2013504133.tif 000088
【数４７】
2013504133.tif 000089
【００９８】
式４６及び式４７は、明示的なＬＰＶ形式である。したがって、十分に平滑化され、変数の近傍に明確に定義されたべき級数展開が存在するような関数f(xt, uu, ρt, νt)及びh(xt, uu, ρt, νt)に関して、プロセスのＬＰＶ近似値が存在してもよい。
【００９９】
2013504133.tif 000090
【０１００】
前述の説明及び添付図面では、本発明に係る原理、好適な実施形態及び動作モードを示したものである。しかしながら、本発明は、上記の特定な実施形態によって範囲を限定されるものではない。上述の実施形態に加えられる様々な変形は、当業者には明らかとなろう。
【０１０１】
したがって、上述の実施形態は、限定的ではなく例示的なものであると考えるべきであり、以下に述べる特許請求の範囲に定義されるように、本発明の範囲を逸脱しない限り、当業者によってこれらの実施形態は変形なされ得ると認識すべきである。
【図１】
2013504133.tif 000091
【図２】
2013504133.tif 000092
【図３】
2013504133.tif 000093
【図４】
2013504133.tif 000094
【図５】
2013504133.tif 000095
【図６】
2013504133.tif 000096
【国際調査報告】
2013504133.tif 000097
2013504133.tif 000098