JP4815391B2

JP4815391B2 - モデルパラメータ推定演算装置及び方法、モデルパラメータ推定演算処理プログラム並びにそれを記録した記録媒体

Info

Publication number: JP4815391B2
Application number: JP2007129381A
Authority: JP
Inventors: 信行友近
Original assignee: Kobe Steel Ltd
Current assignee: Kobe Steel Ltd
Priority date: 2007-05-15
Filing date: 2007-05-15
Publication date: 2011-11-16
Anticipated expiration: 2027-05-15
Also published as: JP2008287344A

Description

本発明は、動的モデルを有するシステムに対する理論的・経験的知識を活用するシステム同定法を用いて、当該動的モデルのモデルパラメータを推定演算するモデルパラメータの推定演算装置及び方法、モデルパラメータ推定演算プログラム並びに当該プログラムを記録したコンピュータにより読取可能な記録媒体に関する。

あるシステムの制御系を設計したり、その挙動をシミュレートしたりするような場合には、システムの挙動を模擬できるモデルが必要となる。例えば、モデルベースのコントローラを設計する場合にはモデルが不可欠であるし、ＰＩＤ（Proportional Integral Derivative）補償器を用いる場合でも事前に制御系の挙動をシミュレートする際にモデルが必要となる。

一般に、モデルを構築する方法は、２つ（ないしは３つ）に大別される。

１つの方法は、いわゆる第一原理に基づいたモデリング（ホワイトボックスモデリングとか、物理モデリングと呼ばれる。）であり、システムを支配する物理化学的な自然法則（運動方程式、電磁界方程式、物質収支、熱収支、化学反応方程式など）に基づいてモデリングを行う方法である。大規模なシステムの場合は、サブシステムに分割してそれぞれモデリングを行い、それらを統合してシステム全体のモデルを構築するアプローチである。

もう１つの方法は、ブラックボックスモデリングと呼ばれるもので、システムに関する物理化学的な法則や先験的知識を利用しないアプローチである。システムをブラックボックスと見なして、観測された（有限個の）入出力データからモデルのパラメータを推定する方法であり、制御工学の分野ではシステム同定という。システムが比較的大規模あるいは複雑であり、物理化学法則のみによってその動特性が記述できないような場合に有効である。

さらには、（モデリング手法を３つに大別する場合には）３つ目の分類として、グレーボックスモデリングが存在する。これは、ホワイトボックスモデリングとブラックボックスモデリングの中間に位置づけられるもので、例えば、入出力データに加えて、物理化学的な自然法則を部分的に利用するようなモデリング手法である。

さて、モデリングの対象となるシステムには、多リンク構造のロボット（マニピュレータ）のように基本的な挙動が運動方程式で表現できるものもあるが、熱や化学反応を伴うプラント（例えば、化学プラントや燃焼炉など）のように、その動特性が自然法則から導き難いものもある。この場合、入出力データに基づいてモデルを構築するブラックボックスモデリングに頼らざるを得ない。

しかしながら、ブラックボックスモデリングにおいては、例えば、入力信号に対するＰＥ性（Persistently Exciting:持続励振性）などのように満たすべき条件が存在する。ところが、現実のプロセスにおいては、システム同定用入出力データを採取する実験機会が限られていたり、印加できる同定入力の波形や大きさが限られている場合も多い。また、入出力データには、少なからず外乱やノイズも存在し、場合によっては、無視できないほど大きな外乱が存在する場合もある。例えば、ごみ焼却炉のように、不特定で多種多様のものを燃焼させるようなプロセスでモデルを構築する場合、ごみ供給機の速度入力に対して、実際に炉内に投入されるごみの種類や量はばらついてしまう。従って、ごみ供給機の速度を入力とし、炉内の温度を出力とするようなモデルを構築する場合、入出力データには大きな外乱が存在することになる。しかも、データ採取の機会も限られているため、外乱の影響を打ち消すのに十分な量のデータを取得することもできない。

このように、入出力データに無視できない外乱やノイズが存在し、印加できる同定入力の波形や大きさが制約を受け、入出力データを採取する実験機会も限られているようなプロセスでは、高精度にモデリングするのに十分な（質的・量的に十分な）データを揃えることができないため、モデリング作業は試行錯誤の繰り返しに頼らざるを得ず、モデルを構築するまでに莫大な時間を要していた。例えば、データの取得から、データ区間の選定、サンプリングタイムの調整、データの前処理（例えば、異常値の除去、フィルタリング等を含む。）、モデル次数の選定などを何度も何度も繰り返す必要があった。また、このように入出力データの質と量が限られている中で、高精度で信頼性のあるモデルを得ることは非常に困難であった。

特開２００４−２７２９１６号公報。足立修一著，「ＭＡＴＬＡＢによる制御のためのシステム同定」，東京電機大学出版局，ｐｐ．１〜５，１９９６年１１月３０日発行。片山徹著，「システム同定入門（システム制御情報ライブラリー９）」，朝倉書店，ｐｐ．４〜７，１９９４年５月２５日発行。 Herbert J. A. F. Tulleken, "Grey-box Modelling and Identification Using Physical Knowledge and Bayesian Techniques", Automatica, Vol. 29, No. 2, pp.285-308, 1993. 片山徹著，「システム同定−部分空間法からのアプローチ−」，朝倉書店，ｐｐ．１０１〜１０７，１８５〜２４０，足立修一著，「ＭＡＴＬＡＢによる制御のための上級システム同定」，東京電機大学出版局，ｐｐ．１７９〜１８１，２００４年３月１０日。 B. D. O. Anderson et al., 荒木光彦訳，「制御装置の低次元化―その考え方と手法−」，システム／制御／情報，システム制御情報学会，Ｖｏｌ．３４，Ｎｏ．９，ｐｐ．４９１〜４９９，平成２年９月。大日方五郎，「平衡実現、ハンケルノルム及び連分数展開を用いた制御系の低次元化」，システム／制御／情報，システム制御情報学会，Ｖｏｌ．３４，Ｎｏ．９，ｐｐ．５０８〜５１４，平成２年９月。

以上の問題点を解決するために、先験的知識を補うことによって、モデルを高精度化しようとしたものに以下の先行技術がある。

例えば、非特許文献１においては、モデルの安定性と定常利得の符号に関する先験的知識を制約条件として考慮したモデリング技術が開示されている。しかしながら、プロセスに対して考慮すべき先見的情報は、定常利得の符号や安定性だけではない。定常利得や応答速度のおよその値（取り得る値の範囲）、任意の入力に対する応答波形の概略は、理論的・経験的に既知の場合も多いが、本手法では、これらの先験的知識を考慮に入れることができない。従って、本手法では、定常利得や応答速度、任意の入力に対する応答波形に関する先験的知識に合致するモデルが得られるとは限らず、モデルの精度や信頼性が不十分であるという問題点があった。

また、特許文献１においては、まず、既存のモデリング手法を使って、入出力データから初期モデルを求めた後、この初期モデルのモデルパラメータとの誤差（重みつき誤差）が小さくなるように、かつ、定常利得や安定性、システムの応答に関する制約条件を満たすようにモデルパラメータを変更し、この変更されたモデルパラメータを新たなモデルとして求める手法が開示されている。しかしながら、変更後のモデルパラメータが、初期モデルのモデルパラメータと値が近いからといって、モデルの動特性が近いとは限らない。例えば、次式の１入力１出力で次数１の状態空間モデルについて考える。

ここで、出力データｙ（ｋ）、入力データｕ（ｋ）、状態パラメータｘ（ｋ）、モデルパラメータａ_０，ｂ_０，ｃ_０，ｄ_０（それぞれスカラーである。）に対して、初期モデルが（ａ_０，ｂ_０，ｃ_０，ｄ_０）＝（０．１，１．０，１．０，０．０）だったとする。この初期モデルに対して、ｂ_０，ｃ_０のパラメータが異なる２つのモデル、モデル１（ａ_０，ｂ_０，ｃ_０，ｄ_０）＝（０．１，０．９，０．９，０．０）、モデル２（ａ_０，ｂ_０，ｃ_０，ｄ_０）＝（０．１，０．２，５．０，０．０）を定義する。初期モデルとパラメータ空間上で値が近いのは、モデル１の方であるが、初期モデルと動特性が近いのはモデル２の方である。この場合、初期モデルとモデル２は全く同じモデルを表している。

そもそも、状態空間表現の場合、モデルのパラメータは一意に決まらない（無数の値をとり得る）ので、モデルのパラメータの近さを評価関数とすることは意味がない。初期モデルを求めるときのみ入出力データを用いており、先見的情報に基づく制約条件をできるだけ満たすようにモデルパラメータを最適化計算する際には、初期モデルのみを規範としており、入出力データへの適合性を全く考慮していない。入出力データへの適合性は初期モデルを通して間接的に考慮されているだけであり、入出力データがモデリングに反映されているか否かは、初期モデルの精度に依存している。

本発明の目的は以上の問題点を解決し、入出力データに基づいて離散時間状態空間モデルのモデルパラメータを推定演算するときに、質や量の不十分な入出力データからでも、精度と信頼性の高いモデルを構築することができるモデルパラメータ推定演算装置及び方法、モデルパラメータ推定演算処理プログラム並びにそれを記録した記録媒体を提供することにある。

第１の発明に係るモデルパラメータ推定演算装置は、入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算装置において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件と、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件とを有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つに関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算手段を備えたことを特徴とする。

第２の発明に係るモデルパラメータ推定演算装置は、入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算装置において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件を有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つ、及び
（ｃ）マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方
に関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算手段を備えたことを特徴とする。

上記モデルパラメータ推定演算装置において、上記演算手段における制約条件はさらに、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件を有することを特徴とする。

上記モデルパラメータ推定演算装置において、上記評価関数は、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差に関する重み付きユークリッドノルムと、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差に関する重み付きユークリッドノルムと
の和で表されたことを特徴とする。

さらに、上記モデルパラメータ推定演算装置において、上記規範となるマルコフパラメータベクトルは、上記入力系列に対する応答波形の上下限制約条件の中間値に基づいて設定されたことを特徴とする。

またさらに、上記モデルパラメータ推定演算装置において、上記制約条件はステップ応答波形に関する制約条件であって、上記制約条件がインパルス応答系列であるマルコフパラメータで構成されたことを特徴とする。

また、上記モデルパラメータ推定演算装置において、上記演算手段は、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算したときに、演算後の離散時間状態空間モデルが最適化計算の際の制約条件を満たさなかった場合に、演算後の離散時間状態空間モデルから構成されるマルコフパラメータを上記規範として更新して、上記評価関数を最小化するようにマルコフパラメータベクトルを演算する制約条件付き最適化計算を繰り返すことを特徴とする。

第３の発明に係るモデルパラメータ推定演算方法は、入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算方法において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件と、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件とを有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つに関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算ステップを含むことを特徴とする。

第４の発明に係るモデルパラメータ推定演算方法は、入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算方法において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件を有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つ、及び
（ｃ）マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方
に関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算ステップを含むことを特徴とする。

上記モデルパラメータ推定演算方法において、上記演算ステップはさらに、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件を有することを特徴とする。

上記モデルパラメータ推定演算方法において、上記評価関数は、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差に関する重み付きユークリッドノルムと、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差に関する重み付きユークリッドノルムと
の和で表されたことを特徴とする。

さらに、上記モデルパラメータ推定演算方法において、上記入力系列に対する応答波形の上下限制約条件の中間値に基づいて設定されたことを特徴とする。

またさらに、上記モデルパラメータ推定演算方法において、上記制約条件はステップ応答波形に関する制約条件であって、上記制約条件がインパルス応答系列であるマルコフパラメータで構成されたことを特徴とする。

また、上記モデルパラメータ推定演算方法において、上記演算ステップは、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算したときに、演算後の離散時間状態空間モデルが最適化計算の際の制約条件を満たさなかった場合に、演算後の離散時間状態空間モデルから構成されるマルコフパラメータを上記規範として更新して、上記評価関数を最小化するようにマルコフパラメータベクトルを演算する制約条件付き最適化計算を繰り返すことを含むことを特徴とする。

第５の発明に係るモデルパラメータ推定演算処理プログラムは、上記モデルパラメータ推定演算方法の演算ステップを含むことを特徴とする。

第６の発明に係るコンピュータで読取可能な記録媒体は、上記モデルパラメータ推定演算処理プログラムを記録したことを特徴とする。

従って、本発明に係るモデルパラメータ推定演算装置及び方法によれば、以下の特有の効果を有する。
（１）任意の入力に対する応答波形に関して制約条件を考慮できるので、幅広い先験的情報を活用してモデリングすることができる。その結果、質が悪く量の少ない入出力データからでも、精度と信頼性の高いモデルを求めることができる。
（２）精度と信頼性の高いモデルを使って制御系を設計したり、シミュレーションを実施したりできるので、制御系やシミュレーションの精度・信頼性も向上する。
（３）モデリングが困難で試行錯誤を特に要する多入力多出力系にも適用できる。
（４）マルコフパラメータ、又はその差分に対し、上下限制約を考慮しているので、あるいは、評価関数によってマルコフパラメータ又はその差分の絶対値が小さくなるように考慮しているので、マルコフパラメータから構成される応答波形が滑らかになるため、マルコフパラメータからモデルパラメータ（後述するＡ、Ｂ、Ｃ行列）を求める際に、制約条件が守られやすくなる。その結果、モデリングのパラメータ調整に要する手間が少なくなる。
（５）規範となるマルコフパラメータを上下限制約の中間に設定したり、任意に波形を設定することによって、規範モデル（後述する（２０）式）が適切に求まらないシステムに対しても、応答波形の制約条件を満たすモデルが求まりやすくなる。精度の悪い規範モデルに変に引っ張られることが少ない。
（６）得られたモデルは必然的に先験的情報に合致したものとなる。従来法では、入出力データからモデルを求めてから、そのモデルが先験的情報に合致するか検証し、合致しない場合には、再度、モデリング工程を（最初から、あるいは、途中から）やり直す必要があった。特に、外乱やノイズの大きなプロセス、印加できる入力波形が限られたプロセス、採取できるデータ数が少ないプロセスにおいては、先験的情報に合致したモデルを構築することが困難であり、入出力データの採取・データ区間の選定・サンプリングタイムの調整・データの前処理などを試行錯誤しながら繰り返す必要があり、モデリング工程に多大な時間と労力を要した。先験的情報を予め考慮に入れることによって、モデリング工程における試行錯誤・やり直しの回数を大幅に減らすことが可能となり、短期間で精度と信頼性の高いモデルを構築できるようになる。その結果、制御系構築やシミュレーションによる検討のリードタイムが短縮化され、早い段階からそれらの効果を享受することができる。

以下、本発明に係る実施形態について図面を参照して説明する。なお、以下の各実施形態において、同様の構成要素については同一の符号を付している。

図１は本発明の一実施形態に係るモデルパラメータ推定演算装置１０の構成を示すブロック図である。本実施形態のモデルパラメータ推定演算装置１０は、図１に示すように、ディジタル計算機を含むように構成され、例えば図２のモデルパラメータ推定演算処理を実行することにより、モデルパラメータ推定演算結果を表示して出力する。モデルパラメータ推定演算装置１０は、入出力データに基づいて任意の入力数と出力数を有するシステムのモデルパラメータを推定して演算するモデルパラメータ推定演算装置であり、
（ａ）当該モデルパラメータ推定演算装置１０の動作及び処理を演算及び制御するコンピュータのＣＰＵ（中央演算処理装置）２０と、
（ｂ）オペレーションプログラムなどの基本プログラム及びそれを実行するために必要なデータを格納するＲＯＭ（読み出し専用メモリ）２１と、
（ｃ）ＣＰＵ２０のワーキングメモリとして動作し、図２のモデルパラメータ推定演算処理において必要なパラメータやデータを一時的に格納するＲＡＭ（ランダムアクセスメモリ）２２と、
（ｄ）例えばハードディスクメモリで構成され、入力パラメータのデータやシミュレーション結果のデータなどのデータを格納するデータメモリ２３と、
（ｅ）例えばハードディスクメモリで構成され、ＣＤ−ＲＯＭドライブ装置４５を用いて読みこんだ図２のモデルパラメータ推定演算処理プログラムを格納するプログラムメモリ２４と、
（ｆ）所定のデータや指示コマンドを入力するためのキーボード４１に接続され、キーボード４１から入力されたデータや指示コマンドを受信して所定の信号変換などのインターフェース処理を行ってＣＰＵ２０に伝送するキーボードインターフェース３１と、
（ｇ）ＣＲＴディスプレイ４３上で指示コマンドを入力するためのマウス４２に接続され、マウス４２から入力されたデータや指示コマンドを受信して所定の信号変換などのインターフェース処理を行ってＣＰＵ２０に伝送するマウスインターフェース３２と、
（ｈ）ＣＰＵ２０によって処理されたデータや設定指示画面などを表示するＣＲＴディスプレイ４３に接続され、表示すべき画像データをＣＲＴディスプレイ４３用の画像信号に変換してＣＲＴディスプレイ４３に出力して表示するディスプレイインターフェース３３と、
（ｉ）ＣＰＵ２０によって処理されたデータ及び所定のモデルパラメータ推定演算結果などを印字するプリンタ４４に接続され、印字すべき印字データの所定の信号変換などを行ってプリンタ４４に出力して印字するプリンタインターフェース３４と、
（ｊ）図２のモデルパラメータ推定演算処理プログラムが記憶されたＣＤ−ＲＯＭ４５ａからモデルパラメータ推定演算処理プログラムのプログラムデータを読み出すＣＤ−ＲＯＭドライブ装置４５に接続され、読み出されたモデルパラメータ推定演算処理プログラムのプログラムデータを所定の信号変換などを行ってプログラムメモリ２４に転送するドライブ装置インターフェース３５と、
（ｋ）入力信号発生器５１からのアナログ入力信号をＡ／Ｄ変換するＡ／Ｄ変換器２５ａと、上記アナログ入力信号に応答してモデリング対象システム５０から出力されるアナログ出力信号を検出して出力する出力信号検出器５２からのアナログ出力信号をＡ／Ｄ変換するＡ／Ｄ変換器２５ｂとを含む入力インターフェース２５とを備える。
ここで、これらの回路２０−２５及び３１−３４はバス３０を介して接続される。なお、Ａ／Ｄ変換器２５ａ，２５ｂによりＡ／Ｄ変換された入力データ及び出力データはバス３０を介してデータメモリ２３に格納される。

次いで、本実施形態に係るモデルパラメータ推定演算装置１０において用いるモデルパラメータ推定演算処理の推定演算方法について以下に詳述する。

この推定演算方法の要旨は、以下の処理を実行することを特徴としている。すなわち、任意の入力に対する応答波形に関する制約条件を、インパルス応答系列であるマルコフ（Marcov）パラメータの関数で構成し、少なくともマルコフパラメータの大きさ（又は、マルコフパラメータの差分の大きさ）に関する制約条件を構成し、上記制約条件の下で、最小化すべき評価関数をＪとし、マルコフパラメータＧ（ｋ）（ｋ＝１，２，…，Ｎ；Ｎは入出力データ長である。）の各要素から構成されるマルコフパラメータベクトルをθとし、規範とするマルコフパラメータ

の各要素から構成されるマルコフパラメータベクトルをθ_ｒｅｆとし、出力データｙ（ｋ）から構成される出力データベクトルをＹ_ｄａｔａとし、入力データｕ（ｋ）から構成される行列をΨとしたときに、上記各評価関数及び上記各ベクトルは次式で表される。

ここで、Ｐ＜Ｎであり、ｋ＞Ｐにおいて、Ｇ（ｋ）≒Ｏ（零行列）である。また、Ｇｒｓ（ｋ）は行列Ｇ（ｋ）の行展開であり次式のように表される。

に対して次式で表される。

ここで、

は

の行展開である。また、次式のように表される。

ここで、重み行列Ｗ_１，Ｗ_２においてＷ_１≠Ｏ（零行列）又はＷ_２≠Ｏ（零行列）である。すなわち、重み行列Ｗ_１，Ｗ_２のうち少なくともどちらかは零行列ではなく、一方は零行列であってもよい。本発明の本実施形態に係るモデルパラメータ推定演算処理は、上記制約条件の下で、評価関数Ｊを最小化するマルコフパラメータベクトルθを算出し、マルコフパラメータベクトルθから離散時間状態空間モデルのパラメータを求めることを特徴としている。

次いで、具体的なモデルパラメータ推定演算方法について以下に詳述する。

例えば実際のプラント装置から取得した入出力データ（ｌ入力数ｍ出力；ここで、ｌ，ｍはそれぞれ自然数である。）に対し、非特許文献１に記載されているように、フィルタリング、異常値の除去、デシメーション、データの切り出し等によって前処理された入力データ系列（入力データ行列）Ｕ_ｄａｔａ及び出力データ系列（出力データ行列）Ｙｄａｔａを次式で表す。

ここで、

はｌ次の実数ベクトルを表す。なお、

はｌ次の実数ベクトルの集合を示す。また、

はサンプル時点ｋにおけるｍ次の出力データベクトルを表す。なお、ベクトルや行列の右上の添え字Ｔは、ベクトルや行列の転置を表すものとする。

また、モデリング対象となるシステムについて、ある既知の次式の入力系列Ｕ_ｒを当該モデルのシステムに印加したときの出力系列に関し、その上限値Ｙ_ｍａｘ及び下限値Ｙ_ｍｉｎが先験的情報として分かっていたとする。

ここで、

はサンプル時点ｋにおけるｌ次の既知入力ベクトルを表す。上記上限値Ｙ_ｍａｘは次式で表される。

ここで、

はサンプル時点ｋにおけるｍ次の出力ベクトルの上限値を表す。また、上記下限値Ｙ_ｍｉｎは次式で表される。

ここで、

はサンプル時点ｋにおけるｍ次の出力ベクトルの下限値を表す。

さて、モデリング対象となるシステムのモデル式を次式で表す。

ここで、

はサンプル時点ｋにおけるｎ次の状態ベクトルである。また、

はサンプル時点ｋにおける出力ベクトル（モデルによる計算値）であり、

は、それぞれｎ行ｎ列、ｎ行ｌ列、ｍ行ｎ列の定数行列である。なお、

はｍ行ｎ列の実数行列の集合を表す。（７）式のシステムのマルコフパラメータ（インパルス応答系列）

を用いると、よく知られているようにサンプル時点ｋにおける出力ベクトルは次式で表される。

ここで、マルコフパラメータ（行列）Ｇ（ｋ）の各要素を次式とする。

また、マルコフパラメータ（行列）Ｇ（ｋ）を行展開した縦ベクトルを

と定義し、入力ベクトルｕ（ｋ）の転置ベクトルを出力数（ｍ個）だけブロック対角に並べた行列を

と定義すると、（９）式を次式に書き換えることができる。

さらに、

と定義すると、（２）式の入力系列Ｕ_ｄａｔａが加えられたときのｋ＝１〜Ｎにおける出力ベクトルのモデル計算値Ｙ_ｍは次式で表される。

さて、求めるモデルが安定、すなわち、Ａのすべての固有値が負であれば、（８）式のマルコフパラメータは十分大きなサンプル時点において０（零行列）に収束する。そこで、サンプル時点Ｐを越える時点においては、（８）式のマルコフパラメータはＯ（零行列）とおいて差し支えないものとする。すなわち、ｋ＞Ｐにおいて、マルコフパラメータ（行列）Ｇ（ｋ）＝Ｏ（零行列）が成り立つとする。一般に（２）式と（３）式の入出力データを採取するときには、モデリング対象となるシステムの収束時間よりも長く入出力データを採取するので、入出力データ長ＮはＰよりも大きくなる（Ｎ＞Ｐ）。以上を考慮すると、（１５’）〜（１７’）式は、それぞれ次式で表される。

さて、（４）式の任意の入力系列

に対する応答系列は同様に

先験的情報（ある入力系列が加わったときの出力応答の上下限）を満たすためには（２０）式の

が上下限に入ればよいので、先験的情報を満たすための制約条件式は次式で表される。

一方、規範（参照）となるマルコフパラメータを

とし、（１０）〜（１１）式と同様に行展開した縦ベクトル

を用いて、

と定義する。

ここで、

は、例えば、入出力データ系列から部分空間法などの既存のモデリング手法を用いて求めたモデル

に対して、

として求める。
ここで、制約条件（２０）式の下で、評価関数

又は

を最小化するマルコフパラメータθを求める。

ここで、ｗ_１は、規範となるマルコフパラメータと同定するモデルのマルコフパラメータの誤差を小さくするための重み係数であり、ｗ_２は、入出力データへのフィッティング度合いを調整する重み係数であり、それぞれ０以上の実数である（０でもよいが、ｗ_１とｗ_２がともに０であることはない）。また、Ｗ_１は、規範となるマルコフパラメータと同定するモデルのマルコフパラメータの誤差を小さくするための重み行列であり、それぞれの要素に異なった重みをつけることができる。Ｗ_２は、入出力データへのフィッティング度合いを調整する重み行列であり、それぞれの要素に異なった重みをつけることができる。重み行列Ｗ_１とＷ_２は、通常、対角行列にすることが多く、その対角要素はそれぞれ０以上の実数である。また、重み行列Ｗ_１とＷ_２は零行列でもよいが、重み行列Ｗ_１とＷ_２がともに零行列となることはない。

なお、

はユークリッドノルムを表し、任意の縦ベクトルｑに対し、

である。

また、制約条件（２０）式の下で、（２５）式の評価関数Ｊ、又は（２５’）式の評価関数Ｊを最小化するマルコフパラメータθを求める問題は、二次計画問題として知られ、市販のパッケージソフト等を用いて解くことができる。

求まったマルコフパラメータベクトル

に対して、各要素Ｇ_ｒｓ（ｋ）からマルコフパラメータを構成し（（１１）式から（１０）式を求める。）、ホー・カルマン（Ho-Kalman）の実現問題（例えば、非特許文献３及び４参照。）に従って、（７）式のＡ、Ｂ、Ｃ行列を求める。このとき、モデルの次数を低次に設定してしまうと、（２０）式のような満たすべき制約条件を満たさなくなってしまう可能性がある。そこで、この段階では、高次のモデルを求めておく。

そして、（７）式のＡ、Ｂ、Ｃ行列から、平衡実現に基づくモデル低次元化手法（例えば、非特許文献６及び７参照。）により、低次のモデル式に近似する。

さて、上記では制約条件式として（２０）式を考慮したが、任意の応答波形に関する制約条件式だけでは、マルコフパラメータの各要素ｇ_ｉｊ（ｋ）（ｉ＝１〜ｍ；ｊ＝１〜ｌ）の大きさが非現実的に大きく求まってしまったり（マルコフパラメータの絶対値が大きい＝ステップ応答が急激に変化する）、その差分

（ここで、

）が大きくなりステップ応答波形が尖ったり（ギザギザ状になったり）する可能性がある（図４参照。）。

例えば、簡単のために１入力１出力系においてステップ応答波形に関して制約条件を設定した場合を考える。制約条件（２０）式の下で、評価関数Ｊ（（２５）式又は（２５’）式で表される。）を最小化して求まったマルコフパラメータθが図９のようであったとする。

図９から明らかなように、ステップ応答は上下限制約Ｙ_ｍａｘとＹ_ｍｉｎの間に存在しているが、マルコフパラメータの大きさに制約を設けないと、図９（ａ）のように、ステップ応答波形が急激に変化したり、マルコフパラメータの差分に制約を設けないと、図９（ａ）のようにステップ応答波形がギザギザ状になったりする可能性があることがわかる。その結果、最終的に求まるモデル（（７）式のＡ、Ｂ、Ｃ行列）が不必要に高次なものになるおそれがあるし、逆に低次のモデルでフィッティングしようとすると応答波形に関する制約条件を満たさなくなってしまうおそれがある。

そこで、ステップ応答波形を急激に変化させないようにするためには、（２０）式の制約条件とは別に、マルコフパラメータの各要素ｇ_ｉｊ（ｋ）に対して制約条件を次式のごとく設定すればよい。

この制約条件は、インパルス応答系列

の上下限、言い換えると、次式のステップ応答系列の差分の上下限について、先験的情報がある場合に有効である。

ここで、

ただし、

である。

例えば、すべての入力からすべての出力へのステップ応答波形の利得がすべて正で、単調増加である（逆応答や振動がない）ことが先験的情報として分かっている場合には、（２６）式において、マルコフパラメータθ_ｍｉｎ＝Ｏ（零ベクトル）とすればよい。

また、ステップ応答波形を尖らさない（ギザギザ状にしない）ためには、（２０）式の制約条件とは別に、マルコフパラメータの差分

（ただし、

）に対して制約条件を次式のごとく設定すればよい。

ここで、Ｉ_ｍｌは（ｍ×ｌ）行（ｍ×ｌ）列の単位行列であり、Δθ_ｍａｘ，Δθ_ｍｉｎはマルコフパラメータの差分の上下限制約ベクトルである。

さらに、（２０）式の制約条件において、ステップ応答以外の応答波形（例えばインパルス応答など）に関する制約条件を定めたような場合で、ステップ応答系列

に関しても制約条件を設定しておく場合には、（２０）式の制約条件とは別に、制約条件式として次式を考慮すればよい。

ここで、Θ_ｍａｘ，Θ_ｍｉｎはステップ応答系列

に関する上下限制約ベクトルである。

なお、（２０）式、（２６）式、（２７）式、（２９）式の制約条件は任意に組み合わせ可能であり、例えば、すべての制約条件を考慮する場合には、次式のようにすればよい。

（３１）式の制約条件の下で（２５）式又は（２５’）式を最小化する二次計画問題を解けばよい。

複数の入力パターンに対する応答波形に関して制約条件をつけることも可能であり、ｉ番目の入力パターンによって構成される（１６）式の行列を

とし、応答波形の上下限を定めた（５）式、（６）式のベクトルをそれぞれＹ_{ｉ，ｍａｘ}，Ｙ_{ｉ，ｍｉｎ}とすると、Ｌ個の入力パターンに関する応答波形の制約条件式は次式で表される。

図２は図１のモデルパラメータ推定演算装置１０によって実行されるモデルパラメータ推定演算処理を示すフローチャートである。また、図３は図２のモデルパラメータ推定演算処理を実行するときの一実施例であって、ごみ焼却炉における蒸気弁解度（入力波形データ（ａ））と蒸気流量（出力波形データ（ｂ））とを示す波形図である。以下、図２のモデルパラメータ推定演算処理のフローチャートを参照して、本実施形態に係るモデルパラメータ推定演算処理について以下に説明する。ここでは、図３に示すような入出力データ系列（ごみ焼却炉における蒸気弁開度と蒸気流量）を用いたケースを例に説明する。

［ステップＳ１］
モデリング対象システムに対して入力信号発生器５１からアナログ入力信号を印加し、これに応答してモデリング対象システム５０から出力されるアナログ出力信号を出力信号検出器５２により検出する。アナログ入力信号は入力インターフェース２５内のＡ／Ｄ変換器２５ａに入力されてＡ／Ｄ変換された後、データメモリ２３に格納される。また、アナログ出力信号は入力インターフェース２５内のＡ／Ｄ変換器２５ｂに入力されてＡ／Ｄ変換された後、データメモリ２３に格納される。

［ステップＳ１Ａ］
次いで、先験的情報に基づき、ステップ応答領域の上下限制約を設定する。ここでは、図４に示すように、ステップ応答の上限と下限の波形をそれぞれ設定した。図４（ａ）は実施例において、ステップ応答の上限波形及び下限波形、最適化計算によって求めたマルコフ列積算値（数６６や数７３のように、マルコフパラメータを積算した値である。）及び規範モデルのマルコフ列積算値を示す波形図であり、図４（ｂ）は実現問題によるモデル（４次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。求めるモデルのステップ応答波形は、この上限と下限に挟まれる領域に存在することになる。なお、上下限制約の設定にあたっては、図１のマウス４２又はキーボード４１（入力装置）を用いて入力する。ここで、例えばキーボード４１を用いて折線の座標を入力してもよい。先験的情報を記憶している装置から、自動的にステップ応答領域のデータ（上下限制約データ）を入力してデータメモリ２３に格納してもよい。

［ステップＳ２］
ステップ応答波形の制約条件だけでは、ステップ応答波形がギザギザ状になることも考えられる。特に、入出力データに大きな外乱や雑音が乗っている場合には、その可能性が大きくなってギザギザ状になると考えられる（図４参照）。ステップ応答波形がギザギザ状になると、ホー・カルマンの実現問題で求めた（７）式のシステムが、あるいは、それを低次元化したシステムが、ステップ応答の上下限値から外れてしまうおそれがある。

そこで、ステップ応答波形（＝マルコフパラメータの積算値）が滑らかになるように、マルコフパラメータ（＝インパルス応答系列）（あるいはマルコフパラメータの差分（

；ただし

））にも上下限制約を設定する。

つまり、（２０）式の制約条件式と（２７）式の制約条件式を考慮する。ここでは、マルコフパラメータの差分の絶対値が０．００２以下になるように、かつ、サンプル時点ｋが大きくなるほど（Ｐに近づくほど）、マルコフパラメータの差分の絶対値の上限値が小さくなるようにした（図１０）。図１０は本発明に係る変形例においてマルコフパラメータの差分値の絶対値を時間経過につれて小さくするように設定した一例を示す図である。図１０に示すように、マルコフパラメータの差分の絶対値を小さくすることによって、サンプル時点ｋのマルコフパラメータがサンプル時点ｋ−１に対して急激に変化しないようになる。言い換えれば、マルコフパラメータの積算値であるステップ応答波形の傾きが急激に変化しなくなる。ステップ応答波形の傾きが急激に変化しないということは、ステップ応答波形に尖った部分やギザギザ部分が存在しにくくなり、ステップ応答波形を滑らかにする効果がある。

また、本モデリング対象は、安定なシステムであるため、ステップ応答波形は発散せずにある値に収束する。たとえ、ステップ応答波形に振動成分や逆応答があったとしても、その波形の変化は徐々に緩やかになりながらある値に収束すると考えられる（図１１）。図１１（ａ）はステップ入力の入出力波形データの一例を示し、図１１（ｂ）は上記ステップ入力に対する出力波形データを示す波形図である。

図１１から明らかなように、ステップ応答波形の差分（すなわちマルコフパラメータ）は、巨視的に見れば徐々に小さくなることが考えられ、さらに、その差分（すなわちマルコフパラメータの差分）も巨視的に見れば徐々に小さくなると考えられる。

そこで、このような先験的情報を活用して、マルコフパラメータの差分の絶対値が（巨視的に）徐々に小さくなるように制約条件を定めることによって、最適化計算によって求まるマルコフパラメータの信頼度を上げる効果がある。また、マルコフパラメータからホー・カルマンの手法によって（７）式のＡ、Ｂ、Ｃ行列を求める際や、平衡実現に基づくモデル低次元化手法により、低次のモデル式に近似する際にも、元のマルコフパラメータ（二次計画問題を解いて求めたマルコフパラメータ）ｇ_ｉｊ（ｋ）の持っていた情報の欠落量を少なくし、応答波形に関する制約条件式（（２０）式）を満たしやすくする効果がある。もし、二次計画問題を解いた結果、サンプル時点ｋが大きくなってからマルコフパラメータ（又はその差分）が大きくなるような（ステップ入力に対して時間が経過してから大きく応答波形が変化するような）解が求まったとすると、ホー・カルマンの手法によって（７）式のＡ、Ｂ、Ｃ行列を求める際に不安定なモデルが求まったり、応答波形に関する制約条件式（（２０）式）から逸脱するようなモデルが求まったりする可能性が高くなる。これは、安定かつ低次のシステムにおいて、時間が経ってからステップ応答が急激に変化することは考えにくいためである。

なお、マルコフパラメータもしくはその差分の上下限の設定にあたっては、図１のように、マウスで波形４２を入力し、又はキーボード４１で折線の座標を入力するなど入力装置を用いて入力してもよいし、先験的情報を記憶している記憶装置から、自動的にステップ応答領域のデータを入力してデータメモリ２３に格納してもよい。

［ステップＳ３］
既存のモデリング手法によって規範となるモデルを同定する。ここでは、部分空間法（例えば、特許文献４及び５参照。）によって求めた２次のモデルを規範モデルとした。

［ステップＳ４］
次いで、（２２）〜（２４）式に従って、規範モデルのマルコフパラメータを構成する。規範モデルのマルコフパラメータの積算値（（２２）式のθ_ｒｅｆの積算値）

を図５に示す。なお、本実施例は１入力１出力系なので当該システムの関数は次式で表される。

図５（ａ）は実施例において、ステップ応答の上限波形及び下限波形、最適化計算によって求めたマルコフ列積算値（マルコフパラメータの積算値）、及び規範モデルのマルコフ列積算値を示す波形図であり、図５（ｂ）は図２のモデルパラメータ推定演算処理によって得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。なお、規範モデルマルコフパラメータの積算値は、ステップ応答波形に対する上限値と下限値（（５）式のＹ_ｍａｘ（ｋ）及び（６）式のＹ_ｍｉｎ（ｋ）に対応する。）の制約を満たしていないことに注意する。

［ステップＳ５］
（２５）式の評価関数の重み係数（ｗ_１，ｗ_２）を設定する。規範モデルのマルコフパラメータへのフィッティングを考慮しない場合には、ｗ_１＝０としてもよい。ここでは、重み係数ｗ_１＝ｗ_２＝１とした。なお、（２５’）式の評価関数においては、重み行列Ｗ_１、Ｗ_２をともに単位行列とした場合に相当する。

［ステップＳ６］
ステップＳ１Ａ及びＳ２の制約条件の下、（２５）式の評価関数を最小化するマルコフパラメータを二次計画問題を解くことによって求める。求めたマルコフパラメータの積算値（（１５）式のθの積算値

に対応する。）を図５に示す。図５から明らかなように、求めたマルコフパラメータの積算値は、ステップ応答波形に対する上限値と下限値（（５）式のＹ_ｍａｘ（ｋ）、（６）式のＹ_ｍｉｎ（ｋ）に相当）の制約を満たしており、意図通りにマルコフパラメータが求まっていることが確認できる。

［ステップＳ７］
ホー・カルマンの手法で実現問題を解くことによって、（７）式のＡ、Ｂ、Ｃ行列を求める。ここでは１０次のモデルを求めた。求めたＡ、Ｂ、Ｃ行列を用いてステップ応答波形を描画したものを図５（ｂ）に示す。求めたＡ、Ｂ、Ｃ行列は、ステップ応答波形に関する制約条件を満たしていることが確認できる。

［ステップＳ８−Ｓ１１］
平衡実現に基づき、モデルを低次元化する。求まったモデルが制約条件を満たすなど許容範囲であることを確認し、同定されたモデルパラメータをデータメモリ２３に格納し、ＣＲＴディスプレイ４３に表示して出力し、モデリングを終了する。なお、本実施例では、１０次のモデルを２次のモデルに低次元化した。ステップＳ９でＮＯのときは、規範マルコフ列を所定の方法で更新した後ステップＳ５に戻り、上記の処理を繰り返す。

以上のモデルパラメータ推定演算処理によって、求まったモデルを図５（ｂ）に示す。図５（ｂ）から明らかなように、最終的に得られたモデル（低次元化されたモデル）のステップ応答波形は予め定めた上下限制約の間に入っており、意図通りにモデルが求まったことを確認できる。

なお、ステップＳ２の制約条件を加えない場合には、マルコフパラメータから安定してモデルパラメータが求まらなかった。すなわち、マルコフパラメータからホー・カルマンの実現問題によって１０次のモデルを求めたが不安定なモデルとなってしまった。試行錯誤でパラメータを調整した結果、ホー・カルマンの方法で低次（４次）のモデルを求め、２次のモデルに低次元化した（図４（ｂ）参照。）。ステップ応答が上下限制約に入ったものの、ステップＳ２の制約条件を付け加えた場合に比べて、モデリングに時間を要することとなった。

また、ステップＳ２において、マルコフパラメータそのものに制約条件をつけてもよい。つまり、インパルス応答に制約をつけたものになる。ここでは、マルコフパラメータの絶対値が０．０１以下になるように、かつ、サンプル時点ｋが大きくなるほど（Ｐに近づくほど）、マルコフパラメータの絶対値の上限値が小さくなるようにした。以上によって求めたモデルを図６に示す。図６（ａ）は実施例においてステップ応答の上限波形及び下限波形、マルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、図６（ｂ）は図２のモデルパラメータ推定演算処理においてマルコフパラメータに制約条件を付与したときに得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。図６（ｂ）から明らかなように、求まったモデルは、ステップ応答に関する上下限制約を満たしていることがわかる。

図７（ａ）は実施例において、ステップ応答の上限波形及び下限波形、マルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、図７（ｂ）は図２のモデルパラメータ推定演算処理において規範モデルのマルコフパラメータを求めるときにステップ応答の上下限制約の中間値を用いたときに得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。ここでは、中間値として平均値を用いた。

さらに、ステップＳ３〜Ｓ４において、規範となるマルコフパラメータを求める際、ステップ応答の上下限制約の中間値を用いてもよい（図７（ａ））。つまり、ステップ応答の上下限制約の中間値を差分したものを規範マルコフパラメータとする。求まったモデルは、図７（ｂ）から明らかなように、ステップ応答に関する上下限制約を満たしていることが確認できる。あるいは、直接、規範となるマルコフパラメータを指定してもよい。例えば、ステップ応答の上下限制約の領域内で、先験的情報に基づき適当な波形を入力あるいは手書きし、その差分を規範マルコフパラメータとして採用してもよい。

さらに、ステップＳ８において、求まったモデルが制約条件を満たさない場合、求まったモデルでマルコフパラメータを構成し、これを（２２）式の最適化計算の際の規範マルコフパラメータθ_ｒｅｆとして更新してもよい。求まったモデルを規範モデルとして更新することによって、徐々にステップ応答制約に入るようにモデルが改善することが期待できる。

さらに、ステップ８で平衡実現に基づきモデルを低次元化する際、周波数重み付きの平衡実現（例えば、非特許文献６及び７参照。）を用いてもよい。外乱や雑音の影響を受けがちな高周波成分を無視し、応答波形の特徴的な低周波成分を重視して低次元化できるため、低次元化後のモデルが応答波形制約から外れる可能性が小さくなると考えられる。

本実施例では、ステップ応答波形に関する制約条件を扱ったが、任意の入力に対する出力応答波形に関して、制約条件を与えてもよい。例えば、正弦波入力に関する応答波形制約を考慮することによって、任意の周波数における利得や位相に関する制約を考慮することができる。例えば、角周波数（角速度）ω［ｒａｄ／ｓ］における利得の上限値と下限値がそれぞれＧ_ｍａｘ、Ｇ_ｍｉｎであり、位相遅れの上限値と下限値がそれぞれφ_ｍａｘ、φ_ｍｉｎであることが先験的に分かっているような場合には、図１２に示すように周期２π／ωで振幅１の正弦波入力に対し、振幅がＧ_ｍｉｎからＧ_ｍａｘの間で、かつ、正弦波入力に対して位相遅れがφ_ｍｉｎとφ_ｍａｘの間となるような正弦波形がすべて包含できるように上限値と下限値（（５）式のＹ_ｍａｘ（ｋ）、（６）式のＹ_ｍｉｎ（ｋ）に相当）を定めればよい。ここで、図１２は本発明に係る変形例において正弦波入力に関する応答波形制約を考慮したときの出力波形データ（ａ）及び入出力波形データ（ｂ）を示す波形図である。

また、（３１）式や（３２）式で示したように複数の制約条件を組み合わせることができる。複数の先験的情報を活用することによって、より高精度で信頼性の高いモデルを求めることができる。

さらに、マルコフパラメータに関して、その大きさに対する制約条件式（（２６）式）、その差分の大きさに対する制約条件式（（２７）式）、その積和の大きさに対する制約条件式（（２９）式）を考慮する代わりに、マルコフパラメータの大きさ、差分の大きさ、積和の大きさに関して、評価関数で考慮することによって、それぞれの大きさを抑制してもよい。つまり、（２５）式や（２５’）式の評価関数Ｊの代わりに、次式を用いてもよい。

なお、重み係数ｗ_３〜ｗ_５はそれぞれ０以上の実数であり、重み行列Ｗ_３〜Ｗ_５はそれぞれ対角行列で、その対角要素はそれぞれ０以上の実数である。もちろん、次式のように組み合わせてもよい。

（２１）式の制約条件の下で、（３３）〜（４０）式のいずれかの評価関数を最小化する問題（二次計画問題）を解いてマルコフパラメータを求めればよい。（２６）式、（２７）式、（２９）式のように制約条件として考慮した場合には、不適切な制約条件を設定した場合には、すべての制約を満たす解が存在しなくなる可能性もあるが、（３３）〜（４０）式のように評価関数として考慮することによって、解が求まる可能性を高くすることができる。また、マルコフパラメータの大きさ、差分の大きさ、積和の大きさに対して、それぞれできるだけ大きさを抑制したいが、明確な制約条件を定め難いときには、評価関数として考慮することが有効である。

さらには、マルコフパラメータの大きさ、差分の大きさ、積和の大きさに対して、あるものは制約条件として、また、あるものは評価関数として考慮するようにしてもよいし、制約条件と評価関数の両方で考慮するようにしてもよい。制約条件と評価関数の両方で考慮することによって、マルコフパラメータの大きさ、差分の大きさ、積和の大きさに関して、制約条件に入ることを保証しながら、かつ、それらの大きさができるだけ小さくなるような解を求めることができる。

＜変形例＞
上記の実施形態では、求めたマルコフパラメータからホー・カルマンの手法を用いて（７）式のモデルパラメータＡ、Ｂ、Ｃ行列を求めたが、マルコフパラメータをＡ、Ｂ、Ｃ行列に変換することなく、マルコフパラメータをそのまま用いて、制御系を設計したり、シミュレーションに用いたりすることができる。例えば、任意の入力系列

に対する応答をシミュレーションする場合は、出力系列を

とし、

とすると、出力のシミュレーションの結果は次式で表される。

図８にシミュレーション結果を示す。図７の事例で求めたマルコフパラメータθを用いて図８（ａ）のような入力波形に対する応答を求めた結果を図８（ｂ）に示す。ここで、マルコフパラメータからＡ，Ｂ，Ｃ行列を求めることなく、シミュレーションできることが分かる。

また、マルコフパラメータを用いて制御系を設計する場合には、例えば、モデル予測制御手法を用いればよい。制御対象となるシステムのインパルス応答系列

あるいは、ステップ応答系列

があれば制御系を構成できるので、ホー・カルマンの手法を用いて（７）式のモデルパラメータＡ、Ｂ、Ｃ行列を求める必要がない。

以上によれば、ホー・カルマンの手法を用いて（７）式のモデルパラメータＡ、Ｂ、Ｃ行列を求める際、あるいは、平衡実現によってモデルを低次元化する際に生じる誤差（情報の欠落）がない（応答波形に関する制約条件を逸脱することもない）ので、それだけ高精度なシミュレーションや制御系の設計ができる可能性がある。また、それだけモデルを求める際に必要な労力や計算量も少なくてすむ。

以上の実施形態又は実施例においては、図２のモデルパラメータ推定演算処理プログラムデータをＣＤ−ＲＯＭ４５ａに格納して実行するときにプログラムメモリ２４にロードして実行しているが、本発明はこれに限らず、ＣＤ−Ｒ、ＣＤ−ＲＷ、ＤＶＤ、ＭＯなどの光ディスク又は光磁気ディスクの記録媒体、もしくは、フロッピーディスクなどの磁気ディスクの記録媒体など種々の記録媒体に格納してもよい。これらの記録媒体は，コンピュータで読み取り可能な記録媒体である。また、図２のモデルパラメータ推定演算処理のデータを予めプログラムメモリ２４に格納して当該処理を実行してもよい。

以上詳述したように、本発明に係るモデルパラメータ推定演算装置及び方法によれば、以下の特有の効果を有する。
（１）任意の入力に対する応答波形に関して制約条件を考慮できるので、幅広い先験的情報を活用してモデリングすることができる。その結果、質が悪く量の少ない入出力データからでも、精度と信頼性の高いモデルを求めることができる。
（２）精度と信頼性の高いモデルを使って制御系を設計したり、シミュレーションを実施したりできるので、制御系やシミュレーションの精度・信頼性も向上する。
（３）モデリングが困難で試行錯誤を特に要する多入力多出力系にも適用できる。
（４）マルコフパラメータ、又はその差分に対し、上下限制約を考慮しているので、あるいは、評価関数によってマルコフパラメータ又はその差分の絶対値が小さくなるように考慮しているので、マルコフパラメータから構成される応答波形が滑らかになるため、マルコフパラメータからモデルパラメータ（Ａ、Ｂ、Ｃ行列）を求める際に、制約条件が守られやすくなる。その結果、モデリングのパラメータ調整に要する手間が少なくなる。
（５）規範となるマルコフパラメータを上下限制約の中間に設定したり、任意に波形を設定することによって、規範モデル（（２０）式）が適切に求まらないシステムに対しても、応答波形の制約条件を満たすモデルが求まりやすくなる。精度の悪い規範モデルに変に引っ張られることが少ない。
（６）得られたモデルは必然的に先験的情報に合致したものとなる。従来法では、入出力データからモデルを求めてから、そのモデルが先験的情報に合致するか検証し、合致しない場合には、再度、モデリング工程を（最初から、あるいは、途中から）やり直す必要があった。特に、外乱やノイズの大きなプロセス、印加できる入力波形が限られたプロセス、採取できるデータ数が少ないプロセスにおいては、先験的情報に合致したモデルを構築することが困難であり、入出力データの採取・データ区間の選定・サンプリングタイムの調整・データの前処理などを試行錯誤しながら繰り返す必要があり、モデリング工程に多大な時間と労力を要した。先験的情報を予め考慮に入れることによって、モデリング工程における試行錯誤・やり直しの回数を大幅に減らすことが可能となり、短期間で精度と信頼性の高いモデルを構築できるようになる。その結果、制御系構築やシミュレーションによる検討のリードタイムが短縮化され、早い段階からそれらの効果を享受することができる。

本発明の一実施形態に係るモデルパラメータ推定演算装置１０の構成を示すブロック図である。図１のモデルパラメータ推定演算装置１０によって実行されるモデルパラメータ推定演算処理を示すフローチャートである。図２のモデルパラメータ推定演算処理を実行するときの一実施例であって、ごみ焼却炉における蒸気弁解度（入力波形データ（ａ））と蒸気流量（出力波形データ（ｂ））とを示す波形図である。（ａ）は実施例において、ステップ応答の上限波形及び下限波形、最適化計算によって求めたマルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、（ｂ）は実現問題によるモデル（４次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。（ａ）は実施例において、ステップ応答の上限波形及び下限波形、最適化計算によって求めたマルコフ列積算値（マルコフパラメータの積算値）、及び規範モデルのマルコフ列積算値を示す波形図であり、（ｂ）は図２のモデルパラメータ推定演算処理によって得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。（ａ）は実施例においてステップ応答の上限波形及び下限波形、マルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、（ｂ）は図２のモデルパラメータ推定演算処理においてマルコフパラメータに制約条件を付与したときに得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。（ａ）は実施例において、ステップ応答の上限波形及び下限波形、マルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、（ｂ）は図２のモデルパラメータ推定演算処理において規範モデルのマルコフパラメータを求めるときにステップ応答の上下限制約の中間値を用いたときに得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。シミュレーション結果であって、（ａ）は入力波形の波形図であり、（ｂ）はマルコフパラメータを用いて図８（ａ）のような入力波形に対する結果を示す波形図である。簡単のために１入力１出力系においてステップ応答波形に関して制約条件を設定した場合において、制約条件（２０）式の下で、評価関数Ｊ（（２５）式又は（２５’）式で表される。）を最小化して求まったマルコフパラメータθを示す図である。本発明に係る変形例においてマルコフパラメータの差分値の絶対値を時間経過につれて小さくするように設定した一例を示す図である。（ａ）はステップ入力の入出力波形データの一例を示し、（ｂ）は上記ステップ入力に対する出力波形データを示す波形図である。本発明に係る変形例において正弦波入力に関する応答波形制約を考慮したときの出力波形データ（ａ）及び入出力波形データ（ｂ）を示す波形図である。

符号の説明

１０…モデルパラメータ推定演算装置、
２０…ＣＰＵ、
２１…ＲＯＭ、
２２…ＲＡＭ、
２３…データメモリ、
２４…プログラムメモリ、
２５…入力インターフェース、
２５ａ，２５ｂ…Ａ／Ｄ変換器、
３０…バス、
３１…キーボードインターフェース、
３２…マウスインターフェース、
３３…ディスプレイインターフェース、
３４…プリンタインターフェース、
３５…ドライブ装置インターフェース、
４１…キーボード、
４２…マウス、
４３…ＣＲＴディスプレイ、
４４…プリンタ、
４５…ＣＤ−ＲＯＭドライブ装置、
４５ａ…ＣＤ−ＲＯＭ、
５０…モデリング対象システム、
５１…入力信号発生器、
５２…出力信号検出器。

Claims

入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算装置において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件と、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件とを有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つに関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算手段を備えたことを特徴とするモデルパラメータ推定演算装置。
入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算装置において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件を有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つ、及び
（ｃ）マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方
に関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算手段を備えたことを特徴とするモデルパラメータ推定演算装置。
上記演算手段における制約条件はさらに、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件を有することを特徴とする請求項２記載のモデルパラメータ推定演算装置。
上記評価関数は、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差に関する重み付きユークリッドノルムと、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差に関する重み付きユークリッドノルムと
の和で表されたことを特徴とする請求項１乃至３のうちのいずれか１つに記載のモデルパラメータ推定演算装置。
上記規範となるマルコフパラメータベクトルは、上記入力系列に対する応答波形の上下限制約条件の中間値に基づいて設定されたことを特徴とする請求項１乃至４のうちのいずれか１つに記載のパラメータ推定装置。
上記制約条件はステップ応答波形に関する制約条件であって、上記制約条件がインパルス応答系列であるマルコフパラメータで構成されたことを特徴とする請求項１乃至５記載のうちのいずれか１つに記載のモデルパラメータ推定演算装置。
上記演算手段は、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算したときに、演算後の離散時間状態空間モデルが最適化計算の際の制約条件を満たさなかった場合に、演算後の離散時間状態空間モデルから構成されるマルコフパラメータを上記規範として更新して、上記評価関数を最小化するようにマルコフパラメータベクトルを演算する制約条件付き最適化計算を繰り返すことを特徴とする請求項１乃至６のうちのいずれか１つに記載のモデルパラメータ推定演算装置。
入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算方法において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件と、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件とを有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つに関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算ステップを含むことを特徴とするモデルパラメータ推定演算方法。
入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算方法において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件を有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つ、及び
（ｃ）マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方
に関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算ステップを含むことを特徴とするモデルパラメータ推定演算方法。
上記演算ステップにおける制約条件はさらに、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件を有することを特徴とする請求項９記載のモデルパラメータ推定演算方法。
上記評価関数は、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差に関する重み付きユークリッドノルムと、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差に関する重み付きユークリッドノルムと
の和で表されたことを特徴とする請求項８乃至１０のうちのいずれか１つに記載のモデルパラメータ推定演算方法。
上記規範となるマルコフパラメータベクトルは、上記入力系列に対する応答波形の上下限制約条件の中間値に基づいて設定されたことを特徴とする請求項８乃至１１のうちのいずれか１つに記載のパラメータ推定方法。
上記制約条件はステップ応答波形に関する制約条件であって、上記制約条件がインパルス応答系列であるマルコフパラメータで構成されたことを特徴とする請求項８乃至１２記載のうちのいずれか１つに記載のモデルパラメータ推定演算方法。
上記演算ステップは、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算したときに、演算後の離散時間状態空間モデルが最適化計算の際の制約条件を満たさなかった場合に、演算後の離散時間状態空間モデルから構成されるマルコフパラメータを上記規範として更新して、上記評価関数を最小化するようにマルコフパラメータベクトルを演算する制約条件付き最適化計算を繰り返すことを含むことを特徴とする請求項８乃至１３のうちのいずれか１つに記載のモデルパラメータ推定演算方法。
請求項８乃至１４のうちのいずれか１つに記載のモデルパラメータ推定演算方法の演算ステップを含むことを特徴とするモデルパラメータ推定演算処理プログラム。
請求項１５記載のモデルパラメータ推定演算処理プログラムを記録したことを特徴とするコンピュータにより読取可能な記録媒体。