(57)【要約】【課題】離散時間状態空間モデルのシステムにおいて、質や量の不十分な入出力波形データからでも、精度と信頼性の高いモデルパラメータを推定演算してモデルを構築する。【解決手段】モデルパラメータ推定演算装置１０において、任意の入力データに対する応答波形である出力データに関する制約条件をインパルス応答系列であるマルコフパラメータの関数で構成する。モデルパラメータ推定演算装置１０は、入出力データに基づいて、制約条件の下で、マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、出力データを含む出力データベクトルと、入力データを含む出力データベクトルに離散時間状態空間モデルのモデル関数を乗算してなる乗算ベクトルとの誤差とに関する評価関数を最小化するようにマルコフパラメータベクトルを演算し、マルコフパラメータベクトルに基づいて離散時間状態空間モデルのモデルパラメータを演算する。【選択図】図２
2008287344.tif 000002
【特許請求の範囲】
【請求項１】
入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算装置において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件と、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件とを有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つに関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算手段を備えたことを特徴とするモデルパラメータ推定演算装置。
【請求項２】
入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算装置において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件を有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つ、及び
（ｃ）マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方
に関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算手段を備えたことを特徴とするモデルパラメータ推定演算装置。
【請求項３】
上記演算手段における制約条件はさらに、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件を有することを特徴とする請求項２記載のモデルパラメータ推定演算装置。
【請求項４】
上記評価関数は、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差に関する重み付きユークリッドノルムと、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差に関する重み付きユークリッドノルムと
の和で表されたことを特徴とする請求項１乃至３のうちのいずれか１つに記載のモデルパラメータ推定演算装置。
【請求項５】
上記規範となるマルコフパラメータベクトルは、予め別の方法で構築されたモデル式に基づいて演算されたことを特徴とする請求項１乃至４のうちのいずれか１つに記載のモデルパラメータ推定演算装置。
【請求項６】
上記規範となるマルコフパラメータベクトルは、上記入力系列に対する応答波形の上下限制約条件の中間値に基づいて設定されたことを特徴とする請求項１乃至４のうちのいずれか１つに記載のパラメータ推定装置。
【請求項７】
上記制約条件はステップ応答波形に関する制約条件であって、上記制約条件がインパルス応答系列であるマルコフパラメータで構成されたことを特徴とする請求項１乃至６記載のうちのいずれか１つに記載のモデルパラメータ推定演算装置。
【請求項８】
上記演算手段は、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算したときに、演算後の離散時間状態空間モデルが最適化計算の際の制約条件を満たさなかった場合に、演算後の離散時間状態空間モデルから構成されるマルコフパラメータを上記規範として更新して、上記評価関数を最小化するようにマルコフパラメータベクトルを演算する制約条件付き最適化計算を繰り返すことを特徴とする請求項１乃至７のうちのいずれか１つに記載のモデルパラメータ推定演算装置。
【請求項９】
上記演算手段は、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算することに代えて、上記マルコフパラメータベクトルをそのまま所定の処理に用いることを特徴とする請求項１乃至８記載のモデルパラメータ推定演算装置。
【請求項１０】
入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算方法において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件と、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件とを有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つに関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算ステップを含むことを特徴とするモデルパラメータ推定演算方法。
【請求項１１】
入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算方法において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件を有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つ、及び
（ｃ）マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方
に関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算ステップを含むことを特徴とするモデルパラメータ推定演算方法。
【請求項１２】
上記演算ステップにおける制約条件はさらに、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件を有することを特徴とする請求項１１記載のモデルパラメータ推定演算方法。
【請求項１３】
上記評価関数は、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差に関する重み付きユークリッドノルムと、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差に関する重み付きユークリッドノルムと
の和で表されたことを特徴とする請求項１０乃至１２のうちのいずれか１つに記載のモデルパラメータ推定演算方法。
【請求項１４】
上記規範となるマルコフパラメータベクトルは、予め別の方法で構築されたモデル式に基づいて演算されたことを特徴とする請求項１０乃至１３のうちのいずれか１つに記載のモデルパラメータ推定演算方法。
【請求項１５】
上記規範となるマルコフパラメータベクトルは、上記入力系列に対する応答波形の上下限制約条件の中間値に基づいて設定されたことを特徴とする請求項１０乃至１３のうちのいずれか１つに記載のパラメータ推定方法。
【請求項１６】
上記制約条件はステップ応答波形に関する制約条件であって、上記制約条件がインパルス応答系列であるマルコフパラメータで構成されたことを特徴とする請求項１０乃至１５記載のうちのいずれか１つに記載のモデルパラメータ推定演算方法。
【請求項１７】
上記演算ステップは、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算したときに、演算後の離散時間状態空間モデルが最適化計算の際の制約条件を満たさなかった場合に、演算後の離散時間状態空間モデルから構成されるマルコフパラメータを上記規範として更新して、上記評価関数を最小化するようにマルコフパラメータベクトルを演算する制約条件付き最適化計算を繰り返すことを含むことを特徴とする請求項１０乃至１６のうちのいずれか１つに記載のモデルパラメータ推定演算方法。
【請求項１８】
上記演算ステップは、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算することに代えて、上記マルコフパラメータベクトルをそのまま所定の処理に用いることを含むことを特徴とする請求項１０乃至１７記載のモデルパラメータ推定演算方法。
【請求項１９】
請求項１０乃至１８のうちのいずれか１つに記載のモデルパラメータ推定演算方法の演算ステップを含むことを特徴とするモデルパラメータ推定演算処理プログラム。
【請求項２０】
請求項１９記載のモデルパラメータ推定演算処理プログラムを記録したことを特徴とするコンピュータにより読取可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、動的モデルを有するシステムに対する理論的・経験的知識を活用するシステム同定法を用いて、当該動的モデルのモデルパラメータを推定演算するモデルパラメータの推定演算装置及び方法、モデルパラメータ推定演算プログラム並びに当該プログラムを記録したコンピュータにより読取可能な記録媒体に関する。
【背景技術】
【０００２】
あるシステムの制御系を設計したり、その挙動をシミュレートしたりするような場合には、システムの挙動を模擬できるモデルが必要となる。例えば、モデルベースのコントローラを設計する場合にはモデルが不可欠であるし、ＰＩＤ（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., 荒木光彦訳，「制御装置の低次元化―その考え方と手法−」，システム／制御／情報，システム制御情報学会，Ｖｏｌ．３４，Ｎｏ．９，ｐｐ．４９１〜４９９，平成２年９月。
【非特許文献７】大日方五郎，「平衡実現、ハンケルノルム及び連分数展開を用いた制御系の低次元化」，システム／制御／情報，システム制御情報学会，Ｖｏｌ．３４，Ｎｏ．９，ｐｐ．５０８〜５１４，平成２年９月。
【発明の開示】
【発明が解決しようとする課題】
【００１１】
以上の問題点を解決するために、先験的知識を補うことによって、モデルを高精度化しようとしたものに以下の先行技術がある。
【００１２】
例えば、非特許文献１においては、モデルの安定性と定常利得の符号に関する先験的知識を制約条件として考慮したモデリング技術が開示されている。しかしながら、プロセスに対して考慮すべき先見的情報は、定常利得の符号や安定性だけではない。定常利得や応答速度のおよその値（取り得る値の範囲）、任意の入力に対する応答波形の概略は、理論的・経験的に既知の場合も多いが、本手法では、これらの先験的知識を考慮に入れることができない。従って、本手法では、定常利得や応答速度、任意の入力に対する応答波形に関する先験的知識に合致するモデルが得られるとは限らず、モデルの精度や信頼性が不十分であるという問題点があった。
【００１３】
また、特許文献１においては、まず、既存のモデリング手法を使って、入出力データから初期モデルを求めた後、この初期モデルのモデルパラメータとの誤差（重みつき誤差）が小さくなるように、かつ、定常利得や安定性、システムの応答に関する制約条件を満たすようにモデルパラメータを変更し、この変更されたモデルパラメータを新たなモデルとして求める手法が開示されている。しかしながら、変更後のモデルパラメータが、初期モデルのモデルパラメータと値が近いからといって、モデルの動特性が近いとは限らない。例えば、次式の１入力１出力で次数１の状態空間モデルについて考える。
【００１４】
【数１】
2008287344.tif 000003
【００１５】
ここで、出力データｙ（ｋ）、入力データｕ（ｋ）、状態パラメータｘ（ｋ）、モデルパラメータａ０，ｂ０，ｃ０，ｄ０（それぞれスカラーである。）に対して、初期モデルが（ａ０，ｂ０，ｃ０，ｄ０）＝（０．１，１．０，１．０，０．０）だったとする。この初期モデルに対して、ｂ０，ｃ０のパラメータが異なる２つのモデル、モデル１（ａ０，ｂ０，ｃ０，ｄ０）＝（０．１，０．９，０．９，０．０）、モデル２（ａ０，ｂ０，ｃ０，ｄ０）＝（０．１，０．２，５．０，０．０）を定義する。初期モデルとパラメータ空間上で値が近いのは、モデル１の方であるが、初期モデルと動特性が近いのはモデル２の方である。この場合、初期モデルとモデル２は全く同じモデルを表している。
【００１６】
そもそも、状態空間表現の場合、モデルのパラメータは一意に決まらない（無数の値をとり得る）ので、モデルのパラメータの近さを評価関数とすることは意味がない。初期モデルを求めるときのみ入出力データを用いており、先見的情報に基づく制約条件をできるだけ満たすようにモデルパラメータを最適化計算する際には、初期モデルのみを規範としており、入出力データへの適合性を全く考慮していない。入出力データへの適合性は初期モデルを通して間接的に考慮されているだけであり、入出力データがモデリングに反映されているか否かは、初期モデルの精度に依存している。
【００１７】
本発明の目的は以上の問題点を解決し、入出力データに基づいて離散時間状態空間モデルのモデルパラメータを推定演算するときに、質や量の不十分な入出力データからでも、精度と信頼性の高いモデルを構築することができるモデルパラメータ推定演算装置及び方法、モデルパラメータ推定演算処理プログラム並びにそれを記録した記録媒体を提供することにある。
【課題を解決するための手段】
【００１８】
第１の発明に係るモデルパラメータ推定演算装置は、入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算装置において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件と、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件とを有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つに関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算手段を備えたことを特徴とする。
【００１９】
第２の発明に係るモデルパラメータ推定演算装置は、入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算装置において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件を有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つ、及び
（ｃ）マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方
に関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算手段を備えたことを特徴とする。
【００２０】
上記モデルパラメータ推定演算装置において、上記演算手段における制約条件はさらに、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件を有することを特徴とする。
【００２１】
上記モデルパラメータ推定演算装置において、上記評価関数は、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差に関する重み付きユークリッドノルムと、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差に関する重み付きユークリッドノルムと
の和で表されたことを特徴とする。
【００２２】
また、上記モデルパラメータ推定演算装置において、上記規範となるマルコフパラメータベクトルは、予め別の方法で構築されたモデル式に基づいて演算されたことを特徴とする。
【００２３】
さらに、上記モデルパラメータ推定演算装置において、上記規範となるマルコフパラメータベクトルは、上記入力系列に対する応答波形の上下限制約条件の中間値に基づいて設定されたことを特徴とする。
【００２４】
またさらに、上記モデルパラメータ推定演算装置において、上記制約条件はステップ応答波形に関する制約条件であって、上記制約条件がインパルス応答系列であるマルコフパラメータで構成されたことを特徴とする。
【００２５】
また、上記モデルパラメータ推定演算装置において、上記演算手段は、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算したときに、演算後の離散時間状態空間モデルが最適化計算の際の制約条件を満たさなかった場合に、演算後の離散時間状態空間モデルから構成されるマルコフパラメータを上記規範として更新して、上記評価関数を最小化するようにマルコフパラメータベクトルを演算する制約条件付き最適化計算を繰り返すことを特徴とする。
【００２６】
さらに、上記モデルパラメータ推定演算装置において、上記演算手段は、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算することに代えて、上記マルコフパラメータベクトルをそのまま所定の処理に用いることを特徴とする。
【００２７】
第３の発明に係るモデルパラメータ推定演算方法は、入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算方法において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件と、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件とを有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つに関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算ステップを含むことを特徴とする。
【００２８】
第４の発明に係るモデルパラメータ推定演算方法は、入出力データに基づいて、離散時間状態空間モデルのモデルパラメータを推定演算するモデルパラメータ推定演算方法において、
インパルス応答系列であるマルコフパラメータの関数で構成された任意の入力系列に対する応答波形に関する制約条件を有し、
上記入出力データに基づいて、上記制約条件の下で、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差と、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差と
のうち、少なくともいずれか１つ、及び
（ｃ）マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方
に関する評価関数を最小化するように上記離散時間状態空間モデルのモデルパラメータを演算する演算ステップを含むことを特徴とする。
【００２９】
上記モデルパラメータ推定演算方法において、上記演算ステップはさらに、
マルコフパラメータの大きさと、マルコフパラメータの差分との少なくとも一方に関する制約条件を有することを特徴とする。
【００３０】
上記モデルパラメータ推定演算方法において、上記評価関数は、
（ａ）マルコフパラメータベクトルと規範となるマルコフパラメータベクトルとの誤差に関する重み付きユークリッドノルムと、
（ｂ）上記出力データベクトルと、マルコフパラメータと上記入力データにより計算される出力予測値（推定値）との誤差に関する重み付きユークリッドノルムと
の和で表されたことを特徴とする。
【００３１】
また、上記モデルパラメータ推定演算方法において、上記規範となるマルコフパラメータベクトルは、予め別の方法で構築されたモデル式に基づいて演算されたことを特徴とする。
【００３２】
さらに、上記モデルパラメータ推定演算方法において、上記入力系列に対する応答波形の上下限制約条件の中間値に基づいて設定されたことを特徴とする。
【００３３】
またさらに、上記モデルパラメータ推定演算方法において、上記制約条件はステップ応答波形に関する制約条件であって、上記制約条件がインパルス応答系列であるマルコフパラメータで構成されたことを特徴とする。
【００３４】
また、上記モデルパラメータ推定演算方法において、上記演算ステップは、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算したときに、演算後の離散時間状態空間モデルが最適化計算の際の制約条件を満たさなかった場合に、演算後の離散時間状態空間モデルから構成されるマルコフパラメータを上記規範として更新して、上記評価関数を最小化するようにマルコフパラメータベクトルを演算する制約条件付き最適化計算を繰り返すことを含むことを特徴とする。
【００３５】
さらに、上記モデルパラメータ推定演算方法において、上記演算ステップは、上記マルコフパラメータベクトルに基づいて上記離散時間状態空間モデルのパラメータを演算することに代えて、上記マルコフパラメータベクトルをそのまま所定の処理に用いることを含むことを特徴とする。
【００３６】
第５の発明に係るモデルパラメータ推定演算処理プログラムは、上記モデルパラメータ推定演算方法の演算ステップを含むことを特徴とする。
【００３７】
第６の発明に係るコンピュータで読取可能な記録媒体は、上記モデルパラメータ推定演算処理プログラムを記録したことを特徴とする。
【発明の効果】
【００３８】
従って、本発明に係るモデルパラメータ推定演算装置及び方法によれば、以下の特有の効果を有する。
（１）任意の入力に対する応答波形に関して制約条件を考慮できるので、幅広い先験的情報を活用してモデリングすることができる。その結果、質が悪く量の少ない入出力データからでも、精度と信頼性の高いモデルを求めることができる。
（２）精度と信頼性の高いモデルを使って制御系を設計したり、シミュレーションを実施したりできるので、制御系やシミュレーションの精度・信頼性も向上する。
（３）モデリングが困難で試行錯誤を特に要する多入力多出力系にも適用できる。
（４）マルコフパラメータ、又はその差分に対し、上下限制約を考慮しているので、あるいは、評価関数によってマルコフパラメータ又はその差分の絶対値が小さくなるように考慮しているので、マルコフパラメータから構成される応答波形が滑らかになるため、マルコフパラメータからモデルパラメータ（後述するＡ、Ｂ、Ｃ行列）を求める際に、制約条件が守られやすくなる。その結果、モデリングのパラメータ調整に要する手間が少なくなる。
（５）規範となるマルコフパラメータを上下限制約の中間に設定したり、任意に波形を設定することによって、規範モデル（後述する（２０）式）が適切に求まらないシステムに対しても、応答波形の制約条件を満たすモデルが求まりやすくなる。精度の悪い規範モデルに変に引っ張られることが少ない。
（６）得られたモデルは必然的に先験的情報に合致したものとなる。従来法では、入出力データからモデルを求めてから、そのモデルが先験的情報に合致するか検証し、合致しない場合には、再度、モデリング工程を（最初から、あるいは、途中から）やり直す必要があった。特に、外乱やノイズの大きなプロセス、印加できる入力波形が限られたプロセス、採取できるデータ数が少ないプロセスにおいては、先験的情報に合致したモデルを構築することが困難であり、入出力データの採取・データ区間の選定・サンプリングタイムの調整・データの前処理などを試行錯誤しながら繰り返す必要があり、モデリング工程に多大な時間と労力を要した。先験的情報を予め考慮に入れることによって、モデリング工程における試行錯誤・やり直しの回数を大幅に減らすことが可能となり、短期間で精度と信頼性の高いモデルを構築できるようになる。その結果、制御系構築やシミュレーションによる検討のリードタイムが短縮化され、早い段階からそれらの効果を享受することができる。
【発明を実施するための最良の形態】
【００３９】
以下、本発明に係る実施形態について図面を参照して説明する。なお、以下の各実施形態において、同様の構成要素については同一の符号を付している。
【００４０】
図１は本発明の一実施形態に係るモデルパラメータ推定演算装置１０の構成を示すブロック図である。本実施形態のモデルパラメータ推定演算装置１０は、図１に示すように、ディジタル計算機を含むように構成され、例えば図２のモデルパラメータ推定演算処理を実行することにより、モデルパラメータ推定演算結果を表示して出力する。モデルパラメータ推定演算装置１０は、入出力データに基づいて任意の入力数と出力数を有するシステムのモデルパラメータを推定して演算するモデルパラメータ推定演算装置であり、
（ａ）当該モデルパラメータ推定演算装置１０の動作及び処理を演算及び制御するコンピュータのＣＰＵ（中央演算処理装置）２０と、
（ｂ）オペレーションプログラムなどの基本プログラム及びそれを実行するために必要なデータを格納するＲＯＭ（読み出し専用メモリ）２１と、
（ｃ）ＣＰＵ２０のワーキングメモリとして動作し、図２のモデルパラメータ推定演算処理において必要なパラメータやデータを一時的に格納するＲＡＭ（ランダムアクセスメモリ）２２と、
（ｄ）例えばハードディスクメモリで構成され、入力パラメータのデータやシミュレーション結果のデータなどのデータを格納するデータメモリ２３と、
（ｅ）例えばハードディスクメモリで構成され、ＣＤ−ＲＯＭドライブ装置４５を用いて読みこんだ図２のモデルパラメータ推定演算処理プログラムを格納するプログラムメモリ２４と、
（ｆ）所定のデータや指示コマンドを入力するためのキーボード４１に接続され、キーボード４１から入力されたデータや指示コマンドを受信して所定の信号変換などのインターフェース処理を行ってＣＰＵ２０に伝送するキーボードインターフェース３１と、
（ｇ）ＣＲＴディスプレイ４３上で指示コマンドを入力するためのマウス４２に接続され、マウス４２から入力されたデータや指示コマンドを受信して所定の信号変換などのインターフェース処理を行ってＣＰＵ２０に伝送するマウスインターフェース３２と、
（ｈ）ＣＰＵ２０によって処理されたデータや設定指示画面などを表示するＣＲＴディスプレイ４３に接続され、表示すべき画像データをＣＲＴディスプレイ４３用の画像信号に変換してＣＲＴディスプレイ４３に出力して表示するディスプレイインターフェース３３と、
（ｉ）ＣＰＵ２０によって処理されたデータ及び所定のモデルパラメータ推定演算結果などを印字するプリンタ４４に接続され、印字すべき印字データの所定の信号変換などを行ってプリンタ４４に出力して印字するプリンタインターフェース３４と、
（ｊ）図２のモデルパラメータ推定演算処理プログラムが記憶されたＣＤ−ＲＯＭ４５ａからモデルパラメータ推定演算処理プログラムのプログラムデータを読み出すＣＤ−ＲＯＭドライブ装置４５に接続され、読み出されたモデルパラメータ推定演算処理プログラムのプログラムデータを所定の信号変換などを行ってプログラムメモリ２４に転送するドライブ装置インターフェース３５と、
（ｋ）入力信号発生器５１からのアナログ入力信号をＡ／Ｄ変換するＡ／Ｄ変換器２５ａと、上記アナログ入力信号に応答してモデリング対象システム５０から出力されるアナログ出力信号を検出して出力する出力信号検出器５２からのアナログ出力信号をＡ／Ｄ変換するＡ／Ｄ変換器２５ｂとを含む入力インターフェース２５とを備える。
ここで、これらの回路２０−２５及び３１−３４はバス３０を介して接続される。なお、Ａ／Ｄ変換器２５ａ，２５ｂによりＡ／Ｄ変換された入力データ及び出力データはバス３０を介してデータメモリ２３に格納される。
【００４１】
次いで、本実施形態に係るモデルパラメータ推定演算装置１０において用いるモデルパラメータ推定演算処理の推定演算方法について以下に詳述する。
【００４２】
この推定演算方法の要旨は、以下の処理を実行することを特徴としている。すなわち、任意の入力に対する応答波形に関する制約条件を、インパルス応答系列であるマルコフ（Marcov）パラメータの関数で構成し、少なくともマルコフパラメータの大きさ（又は、マルコフパラメータの差分の大きさ）に関する制約条件を構成し、上記制約条件の下で、最小化すべき評価関数をＪとし、マルコフパラメータＧ（ｋ）（ｋ＝１，２，…，Ｎ；Ｎは入出力データ長である。）の各要素から構成されるマルコフパラメータベクトルをθとし、規範とするマルコフパラメータ
【数２】
2008287344.tif 000004
の各要素から構成されるマルコフパラメータベクトルをθｒｅｆとし、出力データｙ（ｋ）から構成される出力データベクトルをＹｄａｔａとし、入力データｕ（ｋ）から構成される行列をΨとしたときに、上記各評価関数及び上記各ベクトルは次式で表される。
【００４３】
【数３】
2008287344.tif 000005
【数４】
2008287344.tif 000006
【００４４】
ここで、Ｐ＜Ｎであり、ｋ＞Ｐにおいて、Ｇ（ｋ）≒Ｏ（零行列）である。また、Ｇｒｓ（ｋ）は行列Ｇ（ｋ）の行展開であり次式のように表される。
【００４５】
【数５】
2008287344.tif 000007
に対して次式で表される。
【００４６】
【数６】
2008287344.tif 000008
【数７】
2008287344.tif 000009
【００４７】
ここで、
【数８】
2008287344.tif 000010
は
【数９】
2008287344.tif 000011
の行展開である。また、次式のように表される。
【００４８】
【数１０】
2008287344.tif 000012
【００４９】
【数１１】
2008287344.tif 000013
【００５０】
ここで、重み行列Ｗ１，Ｗ２においてＷ１≠Ｏ（零行列）又はＷ２≠Ｏ（零行列）である。すなわち、重み行列Ｗ１，Ｗ２のうち少なくともどちらかは零行列ではなく、一方は零行列であってもよい。本発明の本実施形態に係るモデルパラメータ推定演算処理は、上記制約条件の下で、評価関数Ｊを最小化するマルコフパラメータベクトルθを算出し、マルコフパラメータベクトルθから離散時間状態空間モデルのパラメータを求めることを特徴としている。
【００５１】
次いで、具体的なモデルパラメータ推定演算方法について以下に詳述する。
【００５２】
例えば実際のプラント装置から取得した入出力データ（ｌ入力数ｍ出力；ここで、ｌ，ｍはそれぞれ自然数である。）に対し、非特許文献１に記載されているように、フィルタリング、異常値の除去、デシメーション、データの切り出し等によって前処理された入力データ系列（入力データ行列）Ｕｄａｔａ及び出力データ系列（出力データ行列）Ｙｄａｔａを次式で表す。
【００５３】
【数１２】
2008287344.tif 000014
【数１３】
2008287344.tif 000015
【００５４】
ここで、
【数１４】
2008287344.tif 000016
はｌ次の実数ベクトルを表す。なお、
【数１５】
2008287344.tif 000017
はｌ次の実数ベクトルの集合を示す。また、
【数１６】
2008287344.tif 000018
はサンプル時点ｋにおけるｍ次の出力データベクトルを表す。なお、ベクトルや行列の右上の添え字Ｔは、ベクトルや行列の転置を表すものとする。
【００５５】
また、モデリング対象となるシステムについて、ある既知の次式の入力系列Ｕｒを当該モデルのシステムに印加したときの出力系列に関し、その上限値Ｙｍａｘ及び下限値Ｙｍｉｎが先験的情報として分かっていたとする。
【００５６】
【数１７】
2008287344.tif 000019
【００５７】
ここで、
【数１８】
2008287344.tif 000020
はサンプル時点ｋにおけるｌ次の既知入力ベクトルを表す。上記上限値Ｙｍａｘは次式で表される。
【００５８】
【数１９】
2008287344.tif 000021
【００５９】
ここで、
【数２０】
2008287344.tif 000022
はサンプル時点ｋにおけるｍ次の出力ベクトルの上限値を表す。また、上記下限値Ｙｍｉｎは次式で表される。
【００６０】
【数２１】
2008287344.tif 000023
【００６１】
ここで、
【数２２】
2008287344.tif 000024
はサンプル時点ｋにおけるｍ次の出力ベクトルの下限値を表す。
【００６２】
さて、モデリング対象となるシステムのモデル式を次式で表す。
【００６３】
【数２３】
2008287344.tif 000025
【００６４】
ここで、
【数２４】
2008287344.tif 000026
はサンプル時点ｋにおけるｎ次の状態ベクトルである。また、
【００６５】
【数２５】
2008287344.tif 000027
はサンプル時点ｋにおける出力ベクトル（モデルによる計算値）であり、
【００６６】
【数２６】
2008287344.tif 000028
【００６７】
【数２７】
2008287344.tif 000029
【００６８】
【数２８】
2008287344.tif 000030
は、それぞれｎ行ｎ列、ｎ行ｌ列、ｍ行ｎ列の定数行列である。なお、
【数２９】
2008287344.tif 000031
はｍ行ｎ列の実数行列の集合を表す。（７）式のシステムのマルコフパラメータ（インパルス応答系列）
【数３０】
2008287344.tif 000032
を用いると、よく知られているようにサンプル時点ｋにおける出力ベクトルは次式で表される。
【００６９】
【数３１】
2008287344.tif 000033
【００７０】
ここで、マルコフパラメータ（行列）Ｇ（ｋ）の各要素を次式とする。
【００７１】
【数３２】
2008287344.tif 000034
【００７２】
また、マルコフパラメータ（行列）Ｇ（ｋ）を行展開した縦ベクトルを
【数３３】
2008287344.tif 000035
と定義し、入力ベクトルｕ（ｋ）の転置ベクトルを出力数（ｍ個）だけブロック対角に並べた行列を
【数３４】
2008287344.tif 000036
と定義すると、（９）式を次式に書き換えることができる。
【００７３】
【数３５】
2008287344.tif 000037
【００７４】
さらに、
【数３６】
2008287344.tif 000038
【数３７】
2008287344.tif 000039
【数３８】
2008287344.tif 000040
と定義すると、（２）式の入力系列Ｕｄａｔａが加えられたときのｋ＝１〜Ｎにおける出力ベクトルのモデル計算値Ｙｍは次式で表される。
【００７５】
【数３９】
2008287344.tif 000041
【００７６】
さて、求めるモデルが安定、すなわち、Ａのすべての固有値が負であれば、（８）式のマルコフパラメータは十分大きなサンプル時点において０（零行列）に収束する。そこで、サンプル時点Ｐを越える時点においては、（８）式のマルコフパラメータはＯ（零行列）とおいて差し支えないものとする。すなわち、ｋ＞Ｐにおいて、マルコフパラメータ（行列）Ｇ（ｋ）＝Ｏ（零行列）が成り立つとする。一般に（２）式と（３）式の入出力データを採取するときには、モデリング対象となるシステムの収束時間よりも長く入出力データを採取するので、入出力データ長ＮはＰよりも大きくなる（Ｎ＞Ｐ）。以上を考慮すると、（１５’）〜（１７’）式は、それぞれ次式で表される。
【００７７】
【数４０】
2008287344.tif 000042
【００７８】
【数４１】
2008287344.tif 000043
【００７９】
【数４２】
2008287344.tif 000044
【００８０】
さて、（４）式の任意の入力系列
【数４３】
2008287344.tif 000045
に対する応答系列は同様に
【００８１】
【数４４】
2008287344.tif 000046
【００８２】
【数４５】
2008287344.tif 000047
【００８３】
【数４６】
2008287344.tif 000048
【００８４】
【数４７】
2008287344.tif 000049
【００８５】
先験的情報（ある入力系列が加わったときの出力応答の上下限）を満たすためには（２０）式の
【数４８】
2008287344.tif 000050
が上下限に入ればよいので、先験的情報を満たすための制約条件式は次式で表される。
【００８６】
【数４９】
2008287344.tif 000051
【００８７】
一方、規範（参照）となるマルコフパラメータを
2008287344.tif 000052
とし、（１０）〜（１１）式と同様に行展開した縦ベクトル
【数５０】
2008287344.tif 000053
を用いて、
【数５１】
2008287344.tif 000054
と定義する。
【００８８】
ここで、
【数５２】
2008287344.tif 000055
は、例えば、入出力データ系列から部分空間法などの既存のモデリング手法を用いて求めたモデル
【数５３】
2008287344.tif 000056
に対して、
【数５４】
2008287344.tif 000057
として求める。
ここで、制約条件（２０）式の下で、評価関数
【数５５】
2008287344.tif 000058
又は
【数５６】
2008287344.tif 000059
を最小化するマルコフパラメータθを求める。
【００８９】
ここで、ｗ１は、規範となるマルコフパラメータと同定するモデルのマルコフパラメータの誤差を小さくするための重み係数であり、ｗ２は、入出力データへのフィッティング度合いを調整する重み係数であり、それぞれ０以上の実数である（０でもよいが、ｗ１とｗ２がともに０であることはない）。また、Ｗ１は、規範となるマルコフパラメータと同定するモデルのマルコフパラメータの誤差を小さくするための重み行列であり、それぞれの要素に異なった重みをつけることができる。Ｗ２は、入出力データへのフィッティング度合いを調整する重み行列であり、それぞれの要素に異なった重みをつけることができる。重み行列Ｗ１とＷ２は、通常、対角行列にすることが多く、その対角要素はそれぞれ０以上の実数である。また、重み行列Ｗ１とＷ２は零行列でもよいが、重み行列Ｗ１とＷ２がともに零行列となることはない。
【００９０】
なお、
【数５７】
2008287344.tif 000060
はユークリッドノルムを表し、任意の縦ベクトルｑに対し、
【数５８】
2008287344.tif 000061
である。
【００９１】
また、制約条件（２０）式の下で、（２５）式の評価関数Ｊ、又は（２５’）式の評価関数Ｊを最小化するマルコフパラメータθを求める問題は、二次計画問題として知られ、市販のパッケージソフト等を用いて解くことができる。
【００９２】
求まったマルコフパラメータベクトル
【数５９】
2008287344.tif 000062
に対して、各要素Ｇｒｓ（ｋ）からマルコフパラメータを構成し（（１１）式から（１０）式を求める。）、ホー・カルマン（Ho-Kalman）の実現問題（例えば、非特許文献３及び４参照。）に従って、（７）式のＡ、Ｂ、Ｃ行列を求める。このとき、モデルの次数を低次に設定してしまうと、（２０）式のような満たすべき制約条件を満たさなくなってしまう可能性がある。そこで、この段階では、高次のモデルを求めておく。
【００９３】
そして、（７）式のＡ、Ｂ、Ｃ行列から、平衡実現に基づくモデル低次元化手法（例えば、非特許文献６及び７参照。）により、低次のモデル式に近似する。
【００９４】
さて、上記では制約条件式として（２０）式を考慮したが、任意の応答波形に関する制約条件式だけでは、マルコフパラメータの各要素ｇｉｊ（ｋ）（ｉ＝１〜ｍ；ｊ＝１〜ｌ）の大きさが非現実的に大きく求まってしまったり（マルコフパラメータの絶対値が大きい＝ステップ応答が急激に変化する）、その差分
【数６０】
2008287344.tif 000063
（ここで、
【数６１】
2008287344.tif 000064
）が大きくなりステップ応答波形が尖ったり（ギザギザ状になったり）する可能性がある（図４参照。）。
【００９５】
例えば、簡単のために１入力１出力系においてステップ応答波形に関して制約条件を設定した場合を考える。制約条件（２０）式の下で、評価関数Ｊ（（２５）式又は（２５’）式で表される。）を最小化して求まったマルコフパラメータθが図９のようであったとする。
【００９６】
図９から明らかなように、ステップ応答は上下限制約ＹｍａｘとＹｍｉｎの間に存在しているが、マルコフパラメータの大きさに制約を設けないと、図９（ａ）のように、ステップ応答波形が急激に変化したり、マルコフパラメータの差分に制約を設けないと、図９（ａ）のようにステップ応答波形がギザギザ状になったりする可能性があることがわかる。その結果、最終的に求まるモデル（（７）式のＡ、Ｂ、Ｃ行列）が不必要に高次なものになるおそれがあるし、逆に低次のモデルでフィッティングしようとすると応答波形に関する制約条件を満たさなくなってしまうおそれがある。
【００９７】
そこで、ステップ応答波形を急激に変化させないようにするためには、（２０）式の制約条件とは別に、マルコフパラメータの各要素ｇｉｊ（ｋ）に対して制約条件を次式のごとく設定すればよい。
【００９８】
【数６２】
2008287344.tif 000065
【００９９】
この制約条件は、インパルス応答系列
【数６３】
2008287344.tif 000066
の上下限、言い換えると、次式のステップ応答系列の差分の上下限について、先験的情報がある場合に有効である。
【０１００】
【数６４】
2008287344.tif 000067
ここで、
【数６５】
2008287344.tif 000068
【数６６】
2008287344.tif 000069
ただし、
【数６７】
2008287344.tif 000070
である。
【０１０１】
例えば、すべての入力からすべての出力へのステップ応答波形の利得がすべて正で、単調増加である（逆応答や振動がない）ことが先験的情報として分かっている場合には、（２６）式において、マルコフパラメータθｍｉｎ＝Ｏ（零ベクトル）とすればよい。
【０１０２】
また、ステップ応答波形を尖らさない（ギザギザ状にしない）ためには、（２０）式の制約条件とは別に、マルコフパラメータの差分
【数６８】
2008287344.tif 000071
（ただし、
【数６９】
2008287344.tif 000072
）に対して制約条件を次式のごとく設定すればよい。
【０１０３】
【数７０】
2008287344.tif 000073
【０１０４】
【数７１】
2008287344.tif 000074
【０１０５】
ここで、Ｉｍｌは（ｍ×ｌ）行（ｍ×ｌ）列の単位行列であり、Δθｍａｘ，Δθｍｉｎはマルコフパラメータの差分の上下限制約ベクトルである。
【０１０６】
さらに、（２０）式の制約条件において、ステップ応答以外の応答波形（例えばインパルス応答など）に関する制約条件を定めたような場合で、ステップ応答系列
【数７２】
2008287344.tif 000075
【数７３】
2008287344.tif 000076
に関しても制約条件を設定しておく場合には、（２０）式の制約条件とは別に、制約条件式として次式を考慮すればよい。
【０１０７】
【数７４】
2008287344.tif 000077
【０１０８】
【数７５】
2008287344.tif 000078
【０１０９】
ここで、Θｍａｘ，Θｍｉｎはステップ応答系列
【数７６】
2008287344.tif 000079
に関する上下限制約ベクトルである。
【０１１０】
なお、（２０）式、（２６）式、（２７）式、（２９）式の制約条件は任意に組み合わせ可能であり、例えば、すべての制約条件を考慮する場合には、次式のようにすればよい。
【０１１１】
【数７７】
2008287344.tif 000080
【０１１２】
（３１）式の制約条件の下で（２５）式又は（２５’）式を最小化する二次計画問題を解けばよい。
【０１１３】
複数の入力パターンに対する応答波形に関して制約条件をつけることも可能であり、ｉ番目の入力パターンによって構成される（１６）式の行列を
【数７８】
2008287344.tif 000081
とし、応答波形の上下限を定めた（５）式、（６）式のベクトルをそれぞれＹｉ，ｍａｘ，Ｙｉ，ｍｉｎとすると、Ｌ個の入力パターンに関する応答波形の制約条件式は次式で表される。
【０１１４】
【数７９】
2008287344.tif 000082
【０１１５】
図２は図１のモデルパラメータ推定演算装置１０によって実行されるモデルパラメータ推定演算処理を示すフローチャートである。また、図３は図２のモデルパラメータ推定演算処理を実行するときの一実施例であって、ごみ焼却炉における蒸気弁解度（入力波形データ（ａ））と蒸気流量（出力波形データ（ｂ））とを示す波形図である。以下、図２のモデルパラメータ推定演算処理のフローチャートを参照して、本実施形態に係るモデルパラメータ推定演算処理について以下に説明する。ここでは、図３に示すような入出力データ系列（ごみ焼却炉における蒸気弁開度と蒸気流量）を用いたケースを例に説明する。
【０１１６】
［ステップＳ１］
モデリング対象システムに対して入力信号発生器５１からアナログ入力信号を印加し、これに応答してモデリング対象システム５０から出力されるアナログ出力信号を出力信号検出器５２により検出する。アナログ入力信号は入力インターフェース２５内のＡ／Ｄ変換器２５ａに入力されてＡ／Ｄ変換された後、データメモリ２３に格納される。また、アナログ出力信号は入力インターフェース２５内のＡ／Ｄ変換器２５ｂに入力されてＡ／Ｄ変換された後、データメモリ２３に格納される。
【０１１７】
［ステップＳ１Ａ］
次いで、先験的情報に基づき、ステップ応答領域の上下限制約を設定する。ここでは、図４に示すように、ステップ応答の上限と下限の波形をそれぞれ設定した。図４（ａ）は実施例において、ステップ応答の上限波形及び下限波形、最適化計算によって求めたマルコフ列積算値（数６６や数７３のように、マルコフパラメータを積算した値である。）及び規範モデルのマルコフ列積算値を示す波形図であり、図４（ｂ）は実現問題によるモデル（４次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。求めるモデルのステップ応答波形は、この上限と下限に挟まれる領域に存在することになる。なお、上下限制約の設定にあたっては、図１のマウス４２又はキーボード４１（入力装置）を用いて入力する。ここで、例えばキーボード４１を用いて折線の座標を入力してもよい。先験的情報を記憶している装置から、自動的にステップ応答領域のデータ（上下限制約データ）を入力してデータメモリ２３に格納してもよい。
【０１１８】
［ステップＳ２］
ステップ応答波形の制約条件だけでは、ステップ応答波形がギザギザ状になることも考えられる。特に、入出力データに大きな外乱や雑音が乗っている場合には、その可能性が大きくなってギザギザ状になると考えられる（図４参照）。ステップ応答波形がギザギザ状になると、ホー・カルマンの実現問題で求めた（７）式のシステムが、あるいは、それを低次元化したシステムが、ステップ応答の上下限値から外れてしまうおそれがある。
【０１１９】
そこで、ステップ応答波形（＝マルコフパラメータの積算値）が滑らかになるように、マルコフパラメータ（＝インパルス応答系列）（あるいはマルコフパラメータの差分（
【数８０】
2008287344.tif 000083
；ただし
【数８１】
2008287344.tif 000084
））にも上下限制約を設定する。
【０１２０】
つまり、（２０）式の制約条件式と（２７）式の制約条件式を考慮する。ここでは、マルコフパラメータの差分の絶対値が０．００２以下になるように、かつ、サンプル時点ｋが大きくなるほど（Ｐに近づくほど）、マルコフパラメータの差分の絶対値の上限値が小さくなるようにした（図１０）。図１０は本発明に係る変形例においてマルコフパラメータの差分値の絶対値を時間経過につれて小さくするように設定した一例を示す図である。図１０に示すように、マルコフパラメータの差分の絶対値を小さくすることによって、サンプル時点ｋのマルコフパラメータがサンプル時点ｋ−１に対して急激に変化しないようになる。言い換えれば、マルコフパラメータの積算値であるステップ応答波形の傾きが急激に変化しなくなる。ステップ応答波形の傾きが急激に変化しないということは、ステップ応答波形に尖った部分やギザギザ部分が存在しにくくなり、ステップ応答波形を滑らかにする効果がある。
【０１２１】
また、本モデリング対象は、安定なシステムであるため、ステップ応答波形は発散せずにある値に収束する。たとえ、ステップ応答波形に振動成分や逆応答があったとしても、その波形の変化は徐々に緩やかになりながらある値に収束すると考えられる（図１１）。図１１（ａ）はステップ入力の入出力波形データの一例を示し、図１１（ｂ）は上記ステップ入力に対する出力波形データを示す波形図である。
【０１２２】
図１１から明らかなように、ステップ応答波形の差分（すなわちマルコフパラメータ）は、巨視的に見れば徐々に小さくなることが考えられ、さらに、その差分（すなわちマルコフパラメータの差分）も巨視的に見れば徐々に小さくなると考えられる。
【０１２３】
そこで、このような先験的情報を活用して、マルコフパラメータの差分の絶対値が（巨視的に）徐々に小さくなるように制約条件を定めることによって、最適化計算によって求まるマルコフパラメータの信頼度を上げる効果がある。また、マルコフパラメータからホー・カルマンの手法によって（７）式のＡ、Ｂ、Ｃ行列を求める際や、平衡実現に基づくモデル低次元化手法により、低次のモデル式に近似する際にも、元のマルコフパラメータ（二次計画問題を解いて求めたマルコフパラメータ）ｇｉｊ（ｋ）の持っていた情報の欠落量を少なくし、応答波形に関する制約条件式（（２０）式）を満たしやすくする効果がある。もし、二次計画問題を解いた結果、サンプル時点ｋが大きくなってからマルコフパラメータ（又はその差分）が大きくなるような（ステップ入力に対して時間が経過してから大きく応答波形が変化するような）解が求まったとすると、ホー・カルマンの手法によって（７）式のＡ、Ｂ、Ｃ行列を求める際に不安定なモデルが求まったり、応答波形に関する制約条件式（（２０）式）から逸脱するようなモデルが求まったりする可能性が高くなる。これは、安定かつ低次のシステムにおいて、時間が経ってからステップ応答が急激に変化することは考えにくいためである。
【０１２４】
なお、マルコフパラメータもしくはその差分の上下限の設定にあたっては、図１のように、マウスで波形４２を入力し、又はキーボード４１で折線の座標を入力するなど入力装置を用いて入力してもよいし、先験的情報を記憶している記憶装置から、自動的にステップ応答領域のデータを入力してデータメモリ２３に格納してもよい。
【０１２５】
［ステップＳ３］
既存のモデリング手法によって規範となるモデルを同定する。ここでは、部分空間法（例えば、特許文献４及び５参照。）によって求めた２次のモデルを規範モデルとした。
【０１２６】
［ステップＳ４］
次いで、（２２）〜（２４）式に従って、規範モデルのマルコフパラメータを構成する。規範モデルのマルコフパラメータの積算値（（２２）式のθｒｅｆの積算値）
【数８２】
2008287344.tif 000085
を図５に示す。なお、本実施例は１入力１出力系なので当該システムの関数は次式で表される。
【０１２７】
【数８３】
2008287344.tif 000086
【０１２８】
図５（ａ）は実施例において、ステップ応答の上限波形及び下限波形、最適化計算によって求めたマルコフ列積算値（マルコフパラメータの積算値）、及び規範モデルのマルコフ列積算値を示す波形図であり、図５（ｂ）は図２のモデルパラメータ推定演算処理によって得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。なお、規範モデルマルコフパラメータの積算値は、ステップ応答波形に対する上限値と下限値（（５）式のＹｍａｘ（ｋ）及び（６）式のＹｍｉｎ（ｋ）に対応する。）の制約を満たしていないことに注意する。
【０１２９】
［ステップＳ５］
（２５）式の評価関数の重み係数（ｗ１，ｗ２）を設定する。規範モデルのマルコフパラメータへのフィッティングを考慮しない場合には、ｗ１＝０としてもよい。ここでは、重み係数ｗ１＝ｗ２＝１とした。なお、（２５’）式の評価関数においては、重み行列Ｗ１、Ｗ２をともに単位行列とした場合に相当する。
【０１３０】
［ステップＳ６］
ステップＳ１Ａ及びＳ２の制約条件の下、（２５）式の評価関数を最小化するマルコフパラメータを二次計画問題を解くことによって求める。求めたマルコフパラメータの積算値（（１５）式のθの積算値
【数８４】
2008287344.tif 000087
に対応する。）を図５に示す。図５から明らかなように、求めたマルコフパラメータの積算値は、ステップ応答波形に対する上限値と下限値（（５）式のＹｍａｘ（ｋ）、（６）式のＹｍｉｎ（ｋ）に相当）の制約を満たしており、意図通りにマルコフパラメータが求まっていることが確認できる。
【０１３１】
［ステップＳ７］
ホー・カルマンの手法で実現問題を解くことによって、（７）式のＡ、Ｂ、Ｃ行列を求める。ここでは１０次のモデルを求めた。求めたＡ、Ｂ、Ｃ行列を用いてステップ応答波形を描画したものを図５（ｂ）に示す。求めたＡ、Ｂ、Ｃ行列は、ステップ応答波形に関する制約条件を満たしていることが確認できる。
【０１３２】
［ステップＳ８−Ｓ１１］
平衡実現に基づき、モデルを低次元化する。求まったモデルが制約条件を満たすなど許容範囲であることを確認し、同定されたモデルパラメータをデータメモリ２３に格納し、ＣＲＴディスプレイ４３に表示して出力し、モデリングを終了する。なお、本実施例では、１０次のモデルを２次のモデルに低次元化した。ステップＳ９でＮＯのときは、規範マルコフ列を所定の方法で更新した後ステップＳ５に戻り、上記の処理を繰り返す。
【０１３３】
以上のモデルパラメータ推定演算処理によって、求まったモデルを図５（ｂ）に示す。図５（ｂ）から明らかなように、最終的に得られたモデル（低次元化されたモデル）のステップ応答波形は予め定めた上下限制約の間に入っており、意図通りにモデルが求まったことを確認できる。
【０１３４】
なお、ステップＳ２の制約条件を加えない場合には、マルコフパラメータから安定してモデルパラメータが求まらなかった。すなわち、マルコフパラメータからホー・カルマンの実現問題によって１０次のモデルを求めたが不安定なモデルとなってしまった。試行錯誤でパラメータを調整した結果、ホー・カルマンの方法で低次（４次）のモデルを求め、２次のモデルに低次元化した（図４（ｂ）参照。）。ステップ応答が上下限制約に入ったものの、ステップＳ２の制約条件を付け加えた場合に比べて、モデリングに時間を要することとなった。
【０１３５】
また、ステップＳ２において、マルコフパラメータそのものに制約条件をつけてもよい。つまり、インパルス応答に制約をつけたものになる。ここでは、マルコフパラメータの絶対値が０．０１以下になるように、かつ、サンプル時点ｋが大きくなるほど（Ｐに近づくほど）、マルコフパラメータの絶対値の上限値が小さくなるようにした。以上によって求めたモデルを図６に示す。図６（ａ）は実施例においてステップ応答の上限波形及び下限波形、マルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、図６（ｂ）は図２のモデルパラメータ推定演算処理においてマルコフパラメータに制約条件を付与したときに得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。図６（ｂ）から明らかなように、求まったモデルは、ステップ応答に関する上下限制約を満たしていることがわかる。
【０１３６】
図７（ａ）は実施例において、ステップ応答の上限波形及び下限波形、マルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、図７（ｂ）は図２のモデルパラメータ推定演算処理において規範モデルのマルコフパラメータを求めるときにステップ応答の上下限制約の中間値を用いたときに得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。ここでは、中間値として平均値を用いた。
【０１３７】
さらに、ステップＳ３〜Ｓ４において、規範となるマルコフパラメータを求める際、ステップ応答の上下限制約の中間値を用いてもよい（図７（ａ））。つまり、ステップ応答の上下限制約の中間値を差分したものを規範マルコフパラメータとする。求まったモデルは、図７（ｂ）から明らかなように、ステップ応答に関する上下限制約を満たしていることが確認できる。あるいは、直接、規範となるマルコフパラメータを指定してもよい。例えば、ステップ応答の上下限制約の領域内で、先験的情報に基づき適当な波形を入力あるいは手書きし、その差分を規範マルコフパラメータとして採用してもよい。
【０１３８】
さらに、ステップＳ８において、求まったモデルが制約条件を満たさない場合、求まったモデルでマルコフパラメータを構成し、これを（２２）式の最適化計算の際の規範マルコフパラメータθｒｅｆとして更新してもよい。求まったモデルを規範モデルとして更新することによって、徐々にステップ応答制約に入るようにモデルが改善することが期待できる。
【０１３９】
さらに、ステップ８で平衡実現に基づきモデルを低次元化する際、周波数重み付きの平衡実現（例えば、非特許文献６及び７参照。）を用いてもよい。外乱や雑音の影響を受けがちな高周波成分を無視し、応答波形の特徴的な低周波成分を重視して低次元化できるため、低次元化後のモデルが応答波形制約から外れる可能性が小さくなると考えられる。
【０１４０】
本実施例では、ステップ応答波形に関する制約条件を扱ったが、任意の入力に対する出力応答波形に関して、制約条件を与えてもよい。例えば、正弦波入力に関する応答波形制約を考慮することによって、任意の周波数における利得や位相に関する制約を考慮することができる。例えば、角周波数（角速度）ω［ｒａｄ／ｓ］における利得の上限値と下限値がそれぞれＧｍａｘ、Ｇｍｉｎであり、位相遅れの上限値と下限値がそれぞれφｍａｘ、φｍｉｎであることが先験的に分かっているような場合には、図１２に示すように周期２π／ωで振幅１の正弦波入力に対し、振幅がＧｍｉｎからＧｍａｘの間で、かつ、正弦波入力に対して位相遅れがφｍｉｎとφｍａｘの間となるような正弦波形がすべて包含できるように上限値と下限値（（５）式のＹｍａｘ（ｋ）、（６）式のＹｍｉｎ（ｋ）に相当）を定めればよい。ここで、図１２は本発明に係る変形例において正弦波入力に関する応答波形制約を考慮したときの出力波形データ（ａ）及び入出力波形データ（ｂ）を示す波形図である。
【０１４１】
また、（３１）式や（３２）式で示したように複数の制約条件を組み合わせることができる。複数の先験的情報を活用することによって、より高精度で信頼性の高いモデルを求めることができる。
【０１４２】
さらに、マルコフパラメータに関して、その大きさに対する制約条件式（（２６）式）、その差分の大きさに対する制約条件式（（２７）式）、その積和の大きさに対する制約条件式（（２９）式）を考慮する代わりに、マルコフパラメータの大きさ、差分の大きさ、積和の大きさに関して、評価関数で考慮することによって、それぞれの大きさを抑制してもよい。つまり、（２５）式や（２５’）式の評価関数Ｊの代わりに、次式を用いてもよい。
【０１４３】
【数８５】
2008287344.tif 000088
【数８６】
2008287344.tif 000089
【数８７】
2008287344.tif 000090
【０１４４】
【数８８】
2008287344.tif 000091
【数８９】
2008287344.tif 000092
【数９０】
2008287344.tif 000093
【０１４５】
なお、重み係数ｗ３〜ｗ５はそれぞれ０以上の実数であり、重み行列Ｗ３〜Ｗ５はそれぞれ対角行列で、その対角要素はそれぞれ０以上の実数である。もちろん、次式のように組み合わせてもよい。
【０１４６】
【数９１】
2008287344.tif 000094
【数９２】
2008287344.tif 000095
【０１４７】
（２１）式の制約条件の下で、（３３）〜（４０）式のいずれかの評価関数を最小化する問題（二次計画問題）を解いてマルコフパラメータを求めればよい。（２６）式、（２７）式、（２９）式のように制約条件として考慮した場合には、不適切な制約条件を設定した場合には、すべての制約を満たす解が存在しなくなる可能性もあるが、（３３）〜（４０）式のように評価関数として考慮することによって、解が求まる可能性を高くすることができる。また、マルコフパラメータの大きさ、差分の大きさ、積和の大きさに対して、それぞれできるだけ大きさを抑制したいが、明確な制約条件を定め難いときには、評価関数として考慮することが有効である。
【０１４８】
さらには、マルコフパラメータの大きさ、差分の大きさ、積和の大きさに対して、あるものは制約条件として、また、あるものは評価関数として考慮するようにしてもよいし、制約条件と評価関数の両方で考慮するようにしてもよい。制約条件と評価関数の両方で考慮することによって、マルコフパラメータの大きさ、差分の大きさ、積和の大きさに関して、制約条件に入ることを保証しながら、かつ、それらの大きさができるだけ小さくなるような解を求めることができる。
【０１４９】
＜変形例＞
上記の実施形態では、求めたマルコフパラメータからホー・カルマンの手法を用いて（７）式のモデルパラメータＡ、Ｂ、Ｃ行列を求めたが、マルコフパラメータをＡ、Ｂ、Ｃ行列に変換することなく、マルコフパラメータをそのまま用いて、制御系を設計したり、シミュレーションに用いたりすることができる。例えば、任意の入力系列
【数９３】
2008287344.tif 000096
に対する応答をシミュレーションする場合は、出力系列を
【数９４】
2008287344.tif 000097
とし、
【数９５】
2008287344.tif 000098
【数９６】
2008287344.tif 000099
とすると、出力のシミュレーションの結果は次式で表される。
【０１５０】
【数９７】
2008287344.tif 000100
【０１５１】
図８にシミュレーション結果を示す。図７の事例で求めたマルコフパラメータθを用いて図８（ａ）のような入力波形に対する応答を求めた結果を図８（ｂ）に示す。ここで、マルコフパラメータからＡ，Ｂ，Ｃ行列を求めることなく、シミュレーションできることが分かる。
【０１５２】
また、マルコフパラメータを用いて制御系を設計する場合には、例えば、モデル予測制御手法を用いればよい。制御対象となるシステムのインパルス応答系列
【数９８】
2008287344.tif 000101
あるいは、ステップ応答系列
【数９９】
2008287344.tif 000102
【数１００】
2008287344.tif 000103
があれば制御系を構成できるので、ホー・カルマンの手法を用いて（７）式のモデルパラメータＡ、Ｂ、Ｃ行列を求める必要がない。
【０１５３】
以上によれば、ホー・カルマンの手法を用いて（７）式のモデルパラメータＡ、Ｂ、Ｃ行列を求める際、あるいは、平衡実現によってモデルを低次元化する際に生じる誤差（情報の欠落）がない（応答波形に関する制約条件を逸脱することもない）ので、それだけ高精度なシミュレーションや制御系の設計ができる可能性がある。また、それだけモデルを求める際に必要な労力や計算量も少なくてすむ。
【０１５４】
以上の実施形態又は実施例においては、図２のモデルパラメータ推定演算処理プログラムデータをＣＤ−ＲＯＭ４５ａに格納して実行するときにプログラムメモリ２４にロードして実行しているが、本発明はこれに限らず、ＣＤ−Ｒ、ＣＤ−ＲＷ、ＤＶＤ、ＭＯなどの光ディスク又は光磁気ディスクの記録媒体、もしくは、フロッピーディスクなどの磁気ディスクの記録媒体など種々の記録媒体に格納してもよい。これらの記録媒体は，コンピュータで読み取り可能な記録媒体である。また、図２のモデルパラメータ推定演算処理のデータを予めプログラムメモリ２４に格納して当該処理を実行してもよい。
【産業上の利用可能性】
【０１５５】
以上詳述したように、本発明に係るモデルパラメータ推定演算装置及び方法によれば、以下の特有の効果を有する。
（１）任意の入力に対する応答波形に関して制約条件を考慮できるので、幅広い先験的情報を活用してモデリングすることができる。その結果、質が悪く量の少ない入出力データからでも、精度と信頼性の高いモデルを求めることができる。
（２）精度と信頼性の高いモデルを使って制御系を設計したり、シミュレーションを実施したりできるので、制御系やシミュレーションの精度・信頼性も向上する。
（３）モデリングが困難で試行錯誤を特に要する多入力多出力系にも適用できる。
（４）マルコフパラメータ、又はその差分に対し、上下限制約を考慮しているので、あるいは、評価関数によってマルコフパラメータ又はその差分の絶対値が小さくなるように考慮しているので、マルコフパラメータから構成される応答波形が滑らかになるため、マルコフパラメータからモデルパラメータ（Ａ、Ｂ、Ｃ行列）を求める際に、制約条件が守られやすくなる。その結果、モデリングのパラメータ調整に要する手間が少なくなる。
（５）規範となるマルコフパラメータを上下限制約の中間に設定したり、任意に波形を設定することによって、規範モデル（（２０）式）が適切に求まらないシステムに対しても、応答波形の制約条件を満たすモデルが求まりやすくなる。精度の悪い規範モデルに変に引っ張られることが少ない。
（６）得られたモデルは必然的に先験的情報に合致したものとなる。従来法では、入出力データからモデルを求めてから、そのモデルが先験的情報に合致するか検証し、合致しない場合には、再度、モデリング工程を（最初から、あるいは、途中から）やり直す必要があった。特に、外乱やノイズの大きなプロセス、印加できる入力波形が限られたプロセス、採取できるデータ数が少ないプロセスにおいては、先験的情報に合致したモデルを構築することが困難であり、入出力データの採取・データ区間の選定・サンプリングタイムの調整・データの前処理などを試行錯誤しながら繰り返す必要があり、モデリング工程に多大な時間と労力を要した。先験的情報を予め考慮に入れることによって、モデリング工程における試行錯誤・やり直しの回数を大幅に減らすことが可能となり、短期間で精度と信頼性の高いモデルを構築できるようになる。その結果、制御系構築やシミュレーションによる検討のリードタイムが短縮化され、早い段階からそれらの効果を享受することができる。
【図面の簡単な説明】
【０１５６】
【図１】本発明の一実施形態に係るモデルパラメータ推定演算装置１０の構成を示すブロック図である。
【図２】図１のモデルパラメータ推定演算装置１０によって実行されるモデルパラメータ推定演算処理を示すフローチャートである。
【図３】図２のモデルパラメータ推定演算処理を実行するときの一実施例であって、ごみ焼却炉における蒸気弁解度（入力波形データ（ａ））と蒸気流量（出力波形データ（ｂ））とを示す波形図である。
【図４】（ａ）は実施例において、ステップ応答の上限波形及び下限波形、最適化計算によって求めたマルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、（ｂ）は実現問題によるモデル（４次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。
【図５】（ａ）は実施例において、ステップ応答の上限波形及び下限波形、最適化計算によって求めたマルコフ列積算値（マルコフパラメータの積算値）、及び規範モデルのマルコフ列積算値を示す波形図であり、（ｂ）は図２のモデルパラメータ推定演算処理によって得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。
【図６】（ａ）は実施例においてステップ応答の上限波形及び下限波形、マルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、（ｂ）は図２のモデルパラメータ推定演算処理においてマルコフパラメータに制約条件を付与したときに得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。
【図７】（ａ）は実施例において、ステップ応答の上限波形及び下限波形、マルコフ列積算値及び規範モデルのマルコフ列積算値を示す波形図であり、（ｂ）は図２のモデルパラメータ推定演算処理において規範モデルのマルコフパラメータを求めるときにステップ応答の上下限制約の中間値を用いたときに得られた、実現問題によるモデル（１０次）による推定波形及び低次元化モデル（２次）による推定波形を示す波形図である。
【図８】シミュレーション結果であって、（ａ）は入力波形の波形図であり、（ｂ）はマルコフパラメータを用いて図８（ａ）のような入力波形に対する結果を示す波形図である。
【図９】簡単のために１入力１出力系においてステップ応答波形に関して制約条件を設定した場合において、制約条件（２０）式の下で、評価関数Ｊ（（２５）式又は（２５’）式で表される。）を最小化して求まったマルコフパラメータθを示す図である。
【図１０】本発明に係る変形例においてマルコフパラメータの差分値の絶対値を時間経過につれて小さくするように設定した一例を示す図である。
【図１１】（ａ）はステップ入力の入出力波形データの一例を示し、（ｂ）は上記ステップ入力に対する出力波形データを示す波形図である。
【図１２】本発明に係る変形例において正弦波入力に関する応答波形制約を考慮したときの出力波形データ（ａ）及び入出力波形データ（ｂ）を示す波形図である。
【符号の説明】
【０１５７】
１０…モデルパラメータ推定演算装置、
２０…ＣＰＵ、
２１…ＲＯＭ、
２２…ＲＡＭ、
２３…データメモリ、
２４…プログラムメモリ、
２５…入力インターフェース、
２５ａ，２５ｂ…Ａ／Ｄ変換器、
３０…バス、
３１…キーボードインターフェース、
３２…マウスインターフェース、
３３…ディスプレイインターフェース、
３４…プリンタインターフェース、
３５…ドライブ装置インターフェース、
４１…キーボード、
４２…マウス、
４３…ＣＲＴディスプレイ、
４４…プリンタ、
４５…ＣＤ−ＲＯＭドライブ装置、
４５ａ…ＣＤ−ＲＯＭ、
５０…モデリング対象システム、
５１…入力信号発生器、
５２…出力信号検出器。
【図１】
2008287344.tif 000104
【図２】
2008287344.tif 000105
【図３】
2008287344.tif 000106
【図４】
2008287344.tif 000107
【図５】
2008287344.tif 000108
【図６】
2008287344.tif 000109
【図７】
2008287344.tif 000110
【図８】
2008287344.tif 000111
【図９】
2008287344.tif 000112
【図１０】
2008287344.tif 000113
【図１１】
2008287344.tif 000114
【図１２】
2008287344.tif 000115