(57)【要約】【課題】システムのモデル構築にあたって、閉ループデータを使用可能であり、数値的に安定なモデル構築方法を提供する。【解決手段】対象システムを線形システムと仮定して、その対象システムの入出力データから仮定的にモデル化した状態空間モデルのパラメータを求めて前記対象システムのシステム同定を行う、システムモデル構築方法において、部分空間法を用いて仮定した状態空間モデルの時刻についての多項式を未来部分へ拡張した係数行列を特異値分解した行列に対して、Ｅｉｇｅｎｓｙｓｔｅｍ  Ｒｅａｌｉｚａｔｉｏｎ  Ａｌｇｏｒｉｔｈｍ（ＥＲＡ）処理を行うことで、前記状態空間モデルのパラメータを算出して前記対象システムをモデル化するシステムモデル構築方法を用いる。【選択図】なし【特許請求の範囲】
【請求項１】
対象システムを線形システムと仮定して、その対象システムの入力データ及び出力データから、状態空間モデルのパラメータを求めて前記対象システムのシステム同定を行う、システムモデル構築方法であって、
部分空間法を用いて仮定した状態空間モデルの時刻についての多項式を未来部分へ拡張した係数行列を特異値分解した行列に対して、Ｅｉｇｅｎｓｙｓｔｅｍ  Ｒｅａｌｉｚａｔｉｏｎ  Ａｌｇｏｒｉｔｈｍ処理を行うことで前記状態空間モデルのパラメータを算出して前記対象システムをモデル化する、システムモデル構築方法。
【請求項２】
下記のステップ（Ａ）乃至（Ｈ）の手順により、
対象システムを線形システムと仮定してモデル化した状態空間モデルの多項式のパラメータを求め、前記対象システムのシステム同定化を行う、請求項１に記載のシステム同定方法。
（Ａ）同定を行う上記プロセスプラントの入力データ及び出力データを収集してそのプロセスプラントを線形なシステムと仮定する。
（Ｂ）ｙ（ｔ）∈Ｒnyはシステムの出力、ｕ（ｔ）∈Ｒnuはシステムの入力、ｎｙはｙの次数、ｎｕはｕの次数であり、ｘ（ｔ）∈Ｒnは次数ｎの状態ベクトルであり、ｅ（ｔ）∈Ｒnyは平均値がゼロで共分散行列が未知あるイノベーション過程であり、行列Ａ∈Ｒn×n、行列Ｂ∈Ｒn×nu、行列Ｃ∈Ｒny×n、行列Ｋ∈Ｒn×nyである、下記式（２）の時刻ｔに対する状態空間モデルを構築する。
【数１】
2007260504.tif 000057
（Ｃ）上記（２）式を、下記（４）乃至（６）式により変形して下記（７）式とする。
【数２】
2007260504.tif 000058
（Ｄ）上記の入力ｕ（ｔ）及び出力ｙ（ｔ）とともに上記状態ベクトルｘ（ｔ）の近似値を求めるためのカルマンディレクターマルコフパラメータを最小二乗法により求める。
（Ｅ）上記カルマンプレディクターマルコフパラメータを用いて、状態空間モデルを未来部分（ｔ＋ｆ−１）へ拡張する。
（Ｆ）上記カルマンプレディクターマルコフパラメータを要素として有する下記式（２８）のハンケル行列を定義し、ｋ＝０の際の前記ハンケル行列を、上記サンプリングデータに対する相関関数を用いて最小二乗法により推定する。
（Ｇ）推定された規格化係数を含む行列に特異値分解を施す正準相関解析を行う。（なお、ｐはｔ以下の正数を示す。）
【数３】
2007260504.tif 000059
（Ｈ）上記の特異値分解を施した係数行列に、正準相関解析による正規化のための規格化行列を含むＥｉｇｅｎｓｙｓｔｅｍ  Ｒｅａｌｉｚａｔｉｏｎ  Ａｌｇｏｒｉｔｈｍ処理を適用して、上記の（４）式及び（５）式の左辺で表されるパラメータ及び上記Ｃを算出する。
（Ｉ）上記（Ｈ）で算出したパラメータ及び上記Ｃから、上記（４）乃至（７）式により、上記のパラメータＡ，Ｂ，Ｋを算出する。
【請求項３】
請求項１又は２に記載のシステムモデル構築方法により得た上記状態空間モデル及び上記パラメータを用いて、入力データに対する推定出力データを求めるとともに、上記対象システムにおいて前記入力データに対する実測出力データを測定し、前記推定出力データと前記実測出力データとを比較して、上記状態空間モデル及び上記パラメータの精度を確認する、システムモデル構築方法。
【請求項４】
請求項１乃至３のいずれかに記載のシステムモデル構築方法により求めた上記空間状態モデル及び上記パラメータを用いて、上記対象システムへの入力データに対する推定出力データを求める、システム出力推定方法。
【請求項５】
上記対象システムで実際に得られた実プロセス出力データと、希望する希望出力データとを比較し、その誤差が小さくなるように実プロセスへの入力値を計算するシステムにおいて、
請求項４に記載のシステム出力推定方法により推定した推定出力データが、前記希望出力データに近くなるように、上記対象システムに与える入力データを計算して、その入力データを前記対象システムへ入力する、システム運用方法。
【請求項６】
上記対象システムは化学プラントにおいて１つ以上の製造材料から１つ以上の生成物を生じるシステムであり、
上記入力データが化学プラントに導入する製造材料の供給量、成分比、温度、圧力のうちの少なくとも１つの数値であり、上記出力データが生成物の生産量、成分比、温度、圧力のうちの少なくとも１つの数値であり、
上記対象システムにおける必要な生成物の製造計画を運用する、請求項５に記載のシステム運用方法。
【請求項７】
上記対象システムの上記入力データ及び出力データを入力する入力手段と、
請求項１又は３に記載のモデル構築方法により、上記入力データ及び出力データからモデル化した上記状態空間モデルのパラメータＡ，Ｂ，Ｃ，Ｋを算出してシステム同定化を行うシステム同定化手段と、
算出された上記のパラメータＡ，Ｂ，Ｃ，Ｋからなる状態空間モデルを出力する出力手段とを有するシステムモデル化装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、入力と出力とを有するシステムをモデル化する方法、及び、そのモデルを用いて入力に対する出力を推定する方法に関し、特に、化学プラントのシステムをモデル化する方法に関する。
【背景技術】
【０００２】
従来、化学プラントにおいて、生産量を最適な量に調整するためには、反応塔の動きや特徴を把握するために、操作可能なバルブを動かして、反応塔の挙動、すなわち温度変化や製品の組成変化などを把握するプラントテストを行う。通常、バルブは調節器によって自動的に制御されているが、プラントテストを行う際には、その調節器による制御を切り離して、自動的に制御される制御モードから、人間が操作するマニュアルモードにして、バルブを直接操作してデータを採取する。
【０００３】
しかし、このようなマニュアルモードにすると、プラントを不安定にさせたり、製品の組成が悪くなったりする運用面での問題が発生しやすくなり、現実にスタッフが操作するにしても、テストは困難となることがあった。このため、出力に応じて入力を調製する閉ループとなる制御モードのままで調製を行えるようにする必要があった。
【０００４】
これらの問題を避けるために、実際にプラントに実装されている制御を切り離して行うテストをしなくても、シミュレーションによって適切な入力値を得られるようにするために、入出力データに基づいてシステムの伝達関数をモデル化した、その関数のパラメータを推定し、システム同定するシステムモデルの構築が行われている。このモデル化した伝達関数を用いて入力値に対する出力値の応答を調べることにより、プラントに実装されている制御を切り離して行うプラントテストを行わずに済む。また、得られた伝達関数のパラメータから、全体としてはブラックボックス状態であるシステムの内容を解析することもできる。
【０００５】
ただし、このような同定を行おうとするシステムは単純な入力から出力が得られるものだけではなく、化学プラントにおいて出力である生産量や成分等に応じて随時、操作端であるバルブを調節器を用いて操作することで出力を一定に収束させる場合や、生成物の一部が原料に還流される場合など、出力が入力にフィードバックされるリサイクルプロセスを有するものもある。このため、単純にモデル化できるものではなく、種種の方法が検討されている。
【０００６】
古典法としては、構築した線形多項式のモデルを観測データに当てはめたときに生じる予測誤差の分散を小さくする最尤法を用いて線形多項式のパラメータを決定する予測誤差法（ｐｒｅｄｉｃｔｉｏｎ  ｅｒｒｏｒ  ｍｅｔｈｏｄ：ＰＥＭ）がある。この手法では、閉ループデータを扱うことができ、精度が高いが、収束計算が必要であるために、数値計算的には楽ではなく、入力と出力とが複数ある多変数を扱おうとすると煩雑になってしまう。
【０００７】
また別の古典法として、構築した線形多項式のモデルのパラメータを、最小二乗法を用いて算出する自己回帰法（ａｕｔｏｒｅｇｒｅｓｓｉｖｅ  ｍｏｄｅｌ  ｗｉｔｈ  ｅｘｏｇｅｎｏｕｓ  ｉｎｐｕｔ：ＡＲＸ）がある。この手法は閉ループデータを用いることができ、収束計算が不要であるので計算しやすいが、入力と出力とが複数ある多変数を扱おうとすると煩雑になり、また、上記ＰＥＭより予測精度が落ちてしまう。
【０００８】
これらの煩雑さを解決するために、データを行列化した行列の部分空間を計算することで、収束計算を行わずに、入力と出力とが複数ある多変数を扱うことができ、必要な計算量を減らした部分空間法（ｓｕｂｓｐａｃｅ  ｉｄｅｎｔｉｆｉｃａｔｉｏｎ  ｍｅｔｈｏｄ：ＳＩＭ）が開発された。しかし、この手法は閉ループデータに対して用いると、閉ループを含めたシステムの入出力を計算してしまい、システムだけのモデル化ができないという欠点を有している。
【０００９】
またこれらとは別に、状態変数と行列とを用いた状態空間モデルを構築して、入力と出力とが複数ある多変数を扱うことができるＯｂｓｅｒｖｅｒ／Ｋａｌｍａｎ  ｆｉｌｔｅｒ  Ｉｄｅｎｔｉｆｉｃａｔｉｏｎ（ＯＫＩＤ）という手法が開発された。これは、次数が高いケースに限り閉ループデータを用いることができ、多変数を容易に扱うことができるが、この状態空間モデルの精度は簡易的であり、ＰＥＭに匹敵する精度を求めるには多くのパラメータを用いる必要あるため、使用しにくい。
【００１０】
これらを解決するために、閉ループデータを使用可能であるように拡張した部分空間法として、非特許文献１に記載のＳｕｂｓｐａｃｅ  ｉｄｅｎｔｉｆｉｃａｔｉｏｎ  ｍｅｔｈｏｄ（ＳＳＡＲＸ）が開発されている。この手法は具体的には以下のようなステップによるものである。
【００１１】
（ａ）対象システムが線形なシステムであると仮定する。
（ｂ）出力ｙ、入力ｕ、状態変数ｘ、イノベーション過程ｅ、行列Ａ、Ｂ，Ｃ，Ｋを用いて、前記対象システムの入出力を下記式（２）の状態空間モデルを仮に構築する。
【数４】
2007260504.tif 000002
（ｃ）前記（２）式をカルマンプレディクター形式に変形する。
（ｄ）変形した式中のカルマンプレディクターマルコフパラメータを最小二乗法により算出する。
（ｅ）部分空間法の準備として状態空間モデルを未来部分へ拡張する。
（ｆ）前記状態空間モデルの係数を、サンプリングデータから最小二乗法を用いて求める。
（ｇ）前記（ｆ）の係数行列に特異値分解を施す。
（ｈ）特異値分解した結果から、時刻ｔにおける状態変数ｘを割り当てる。
（ｉ）時刻（ｔ＋１）における状態変数ｘを、前記（ｃ）乃至（ｇ）と同様に求める。
（ｊ）前記（ｂ）の（２）式の第二式から、パラメータＣを最小二乗法で求める。
（ｋ）前記（ｂ）の（２）式の第一式から、パラメータＡ，Ｂ，Ｋを最小二乗法で求める。
【００１２】
【非特許文献１】１３ｔｈ  ＩＦＡＣ  Ｓｙｍｐｏｓｉｕｍ  ｏｎ  Ｓｙｓｔｅｍ  Ｉｄｅｎｔｉｆｉｃａｔｉｏｎ  Ｐ１６２５−１６３０
【発明の開示】
【発明が解決しようとする課題】
【００１３】
しかしながら、上記ＳＳＡＲＸの方法では、係数行列を上記（ｇ）で特異値分解した後、上記（ｊ）で逆行列を用いた最小二乗法により得たパラメータＣを用いて生じる誤差ｅを求め、さらに最小二乗法を用いて残りのパラメータを求めているが、この、最小二乗法では逆行列計算をしている。しかし、逆行列計算では、ゼロに割り算や発散の可能性を有するために、数値が不安定になる問題があった。
【００１４】
そこでこの発明は、対象システムのモデル構築にあたって、閉ループデータを使用可能であり、より数値的に安定なモデル構築方法を提供することを目的とする。
【課題を解決するための手段】
【００１５】
この発明は、対象システムを線形システムと仮定して、その対象システムの入出力データから仮定的にモデル化した状態空間モデルのパラメータを求めて前記対象システムのシステム同定を行う、システムモデル構築方法において、
部分空間法を用いて仮定した状態空間モデルの時刻についての多項式を未来部分へ拡張した係数行列を特異値分解した行列に対して、Ｅｉｇｅｎｓｙｓｔｅｍ  Ｒｅａｌｉｚａｔｉｏｎ  Ａｌｇｏｒｉｔｈｍ（ＥＲＡ）処理を行うことで、前記状態空間モデルのパラメータを算出して前記対象システムをモデル化するシステムモデル構築方法により、上記の課題を解決したのである。
【００１６】
すなわち、従来のＳＳＡＲＸでは特異値分解した後に逆行列を使用する最小二乗法を繰り返していたが、この発明では、逆行列を用いずに行列計算するＥＲＡ処理により直接にパラメータを算出することで、数値的に安定してモデル構築を行うことができる。また、単純な部分空間法と違い、未来部分へ拡張した係数行列を用いるため、出力値が入力値にフィードバックされる閉ループデータを使用可能となっている。
【発明の効果】
【００１７】
この発明にかかるモデル構築方法を用いることにより、収束計算の必要なく、多変数が関与し、閉ループであるシステムのモデル構築を、安定的に行うことができる。このように出来たモデルにより、従来のＳＳＡＲＸモデルと同等以上の精度を得ることができる。
【発明を実施するための最良の形態】
【００１８】
以下、この発明について詳細に説明する。この発明は、対象システムを線形システムと仮定して、その対象システムの入力データ及び出力データ（以下、まとめて「入出力データ」と表記する。）から仮定的にモデル化した状態空間モデルのパラメータを求めて前記対象システムのシステム同定を行う、システムモデル構築方法である。
【００１９】
この発明でモデル化する上記対象システムとは、入力に対して出力を有するシステムであり、その入力と出力とをいずれも数値化できるものである。例えば、質量、エネルギー量、体積、温度、圧力、などである。この、入出力を数値化した入出力サンプリングデータを複数得て、それらのデータが線形を示す線形システムと仮定できるシステムに対して、この発明を用いる。この発明にかかるシステムでは、出力が入力にフィードバックされ、出力が未来の出力に影響を及ぼす閉ループデータであっても、モデル化することができる。
【００２０】
この発明では、上記の入出力サンプリングデータから上記対象システムを線形なシステムであると仮定して、状態空間モデルを作成する。この状態空間モデルは、時刻ｔに対する状態変数ｘ、入力ｕ、出力ｙ、及び誤差ノイズｅの、行列との積からなる多項式であり、時刻（ｔ＋１）における状態変数ｘを表現することで、フィードバックが可能な形式となっている。このパラメータである行列を一部まとめると、項数を減らし、解析しやすい形式に変形することができる。この形式は、カルマンフィルターを用いた予測形式であり、カルマンプレディクター形式と呼ばれるものである。
【００２１】
まとめたパラメータのうち、出力ｙに対する状態変数ｘの係数行列と、時刻（ｔ＋１）の状態変数に対する時刻ｔの状態変数ｘ及び入力ｕの係数行列とからなるマルコフパラメータを、入力ｕと出力ｙとの行列ｚから、最小二乗法で算出する。なお、ここまでの手順は、いわゆるＡＲＸモデル、ＯＫＩＤモデルと同一の手順である。
【００２２】
次に、上記の状態空間モデルを、時刻ｔから、未来時刻（ｔ＋ｆ−１）にまで拡張する。この未来時刻に関する係数行列を含む多項式について、その係数行列を、上記入出力サンプリングデータについて最小二乗法することにより求め、この係数行列を特異値分解する。なお、ここまでの手順はＳＳＡＲＸと同じである。
【００２３】
上記の特異値分解とは行列計算の一手法であり、ｏ×ｑの行列Ｘを、ｏ次直交行列Ｕと、ｑ次直交行列Ｖと、ｏ×ｑ次行列であって対角要素にのみ値を有する行列ＳとからなるＸ＝ＵＳＶTと分解する手法である。
【００２４】
さらに、上記の特異値分解を行った係数行列に対して、Ｅｉｇｅｎｓｙｓｔｅｍ  Ｒｅａｌｉｚａｔｉｏｎ  Ａｌｇｏｒｉｔｈｍ処理（以下、「ＥＲＡ処理」と略す。）を行うことで前記状態空間モデルのパラメータを算出することで、先に立てた状態空間モデルの多項式の係数行列を求めて前記対象システムをモデル化することができる。ここで前記ＥＲＡ処理とは、特異値分解を行った固有行列に相当する行列を求める方法である。具体例を挙げると、上記のマルコフパラメータに相当する係数行列、Ａ，Ｂ，Ｃについて、行列の要素がＣＡKＢであるハンケル行列：ＨKを定義し、このハンケル行列の固有行列を求める。固有行列であれば何乗しても元の行列と変わらないことから、まずマルコフパラメータのうちのＡを求めることができる。また残りの行列ＢとＣとは、固有行列であるＡを切り分けた残りとして算出することができる。また別の方法として、正方行列でなければならない固有行列を求めなくても、特異行列を求めることからＡ，Ｂ，Ｃを算出する特異値分解法（ＳＶＤ）が挙げられる。以下、特異値分解法を用いたＥＲＡ処理を前提に記述する。なお、ここで記述したＡ，Ｂ，Ｃは後述する式（１）のマルコフパラメータである。
【００２５】
このＥＲＡ処理は、従来から知られたものであるが、行列に対して特異値分解を行うためには、上記のＸ＝ＵＳＶTの形式に分解するために、上記入出力サンプリングデータの行列である規格化係数Ｗ1とＷ2とを、行列の乗算の先頭と末尾とに有するように変形し、この規格化係数を正準相関解析（ｃａｎｏｎｉｃａｌ  ｃｏｒｒｅｌａｔｉｏｎ  ａｎａｌｙｓｉｓ：以下、ＣＣＡと略す。）重みとして有する係数行列に対して上記ＥＲＡ処理を行うことにより、逆行列を用いた最小二乗法を使用することなく、パラメータを対応させて求めることができる。
【００２６】
上記のＣＣＡ重みを用いることで、対象となる行列を、計算上左と右とから別個の規格化された行列とにバランスよく切り分けると、左側に切り分けられた部分が、状態を観測できるかという部分に相当し、右側に切り分けられた部分が、状態を制御することができるかという部分に相当するため、プロセスの状態を観測できかつ制御できるようにすることができるという利点がある。逆に、左右の一方にのみ切り分けが偏ると、観測又は制御の一方だけが優位にできるということになってしまう。
【００２７】
これらのパラメータにより、上記対象システムを上記状態空間モデルでモデル化した多項式が完成し、この多項式に入力データｕを導入することで、それに対する推定の出力データｙを求めることができる。上記対象システムを運用するにあたって、この発明にかかるシステムモデル構築方法で同定したモデルを用いて、このように入力データに対する上記対象システムの推定出力データを推定することを、複数の入力データにおいて行うと、要求する出力データが得られる入力データを調べることができる。すなわち、好適な出力データが得られる入力データを上記対象システムに導入することができるようになり、効率のよい上記対象システムの運用が可能となる。
【００２８】
また、推定された推定出力データと、同様の入力データを実際に上記対象システムに対して実行させた実測出力データとを比較して、得られた上記の状態空間モデルとパラメータとの精度を確認してもよい。精度とはすなわち適合性であって、入力データに対する推定出力データが、実際の出力データにどれほど適合するかを示す。その精度が不十分であると判断される場合には、さらに入出力サンプリングデータを測定し直して、又は追加して、上記の手順によるモデル化を行うことで、さらに精度を向上させるとよい。
【００２９】
このようなシステムモデル構築方法を適用する上記対象システムとしては、例えば、化学プラントにおいて、複数種の製造材料から複数種の生成物を生じる多次元システムが挙げられる。このようなシステムは解析すべき変数が多くなるために、従来の方法ではモデル化が困難であるが、この発明にかかる方法を用いることで安定した数値でのモデル化を行うことができる。
【００３０】
このような化学プラントのシステムにおいて、モデル化する入力データとしては、原材料の供給量、成分比、温度、圧力などの、数値化できる値のうち、出力への関与が無視できない要素を選択するとよい。ここで原材料とは単純な製造原料だけではなく、反応温度等を調整するために熱源として用いる水蒸気などの、間接的に用いる材料も含んでよい。上記対象システムに関与する要素であれば考慮されうるためである。一方、出力データとしては、生成物の生産量、成分比、温度、圧力などが挙げられ、これらの中で入力による関与が無視できない要素を選択するとよい。ここで生成物とは、求める製品だけではなく、副生物や、分離して廃棄する廃液や廃棄物等も含む。
【００３１】
この発明にかかるシステムモデル構築方法を利用する方法としては、この発明にかかるシステムモデル構築方法により、上記入出力データからモデル化した上記状態空間モデルのパラメータＡ，Ｂ，Ｃ，Ｋを算出してシステム同定化を行うシステム同定化手段と、上記対象システムの上記入出力データを入力する入力手段と、算出された上記のパラメータＡ，Ｂ，Ｃ，Ｋからなる状態空間モデルを表示する出力手段とを有し、上記入出力データを入力することで、モデル化を行う一体のシステム同定化装置として利用する方法が挙げられる。また、このようにモデル化を行うシステム同定化装置は、さらに、追加で入力データを入力することで、モデル化した状態空間モデルを用いて、追加で入力された入力データに対応する推定出力データを算出できるものであると、より実際の上記対象システムに利用しやすくなり、好ましい。
【００３２】
ここで上記入力手段としては、キーボード等の入力デバイスによる他、上記対象システム内の製造原料の物性値を測定するセンサーや、生成物の物性値を測定するセンサーを用いて、直接に数値を入力するものでもよい。また、上記出力手段としては、ディスプレイ、プリンタなどの一般的な出力デバイスが挙げられる。また、上記のシステム同定化手段は、この発明にかかるモデル構築方法を実現するプログラムを実行するコンピュータである。また、上記物性値として用いることができる上記入力データとしては、上記対象システムで用いる製造原料の供給量、成分比、温度、圧力のうちの、少なくとも一つの数値が挙げられる。さらに、上記物性値として用いることのできる上記出力データとしては、上記対象システムで得られる生成物の生成量、成分比、温度、圧力のうちの、少なくとも一つの数値が挙げられる。
【００３３】
また、より具体的にこの発明にかかるシステムモデル構築方法を利用する方法としては、上記対象システムの上記入出力データ、及び目的とする出力データである希望出力データを入力する入力手段と、この発明にかかるモデル構築方法により、上記入出力データからモデル化した上記状態空間モデルのパラメータＡ，Ｂ，Ｃ，Ｋを算出してシステム同定化を行うシステム同定化手段と、算出された上記のパラメータＡ，Ｂ，Ｃ，Ｋからなる状態空間モデルを用いて、新たな入力データに対する推定出力データを算出する出力算出手段と、上記希望出力データと上記推定出力データとを比較し、その結果により前記新たな入力データの適合性を判断する適合性判断手段と、上記適合性が基準条件を満たした前記新たな入力データを、上記対象システムにおいて実行する実行手段とを有するシステム運用装置として利用する方法が挙げられる。
【００３４】
この場合の、上記入力手段としては、上記で列挙した入力デバイスやセンサーに加えて、上記対象システムの運用上必要な、目的とする希望出力データを、それまでの上記対象システムの運用実績、規模などの過去の情報から自動的に算出するものも含む。また、上記出力手段としては、上記の出力デバイスの他に、直接に上記対象システム内の製造原料供給量や熱量、圧力等を調製するバルブや弁等を調製する調整機も含む。これにより、上記の基準条件を満たした入力データを、実際の上記対象システムに実行させることができる。
【００３５】
また、上記出力算出手段は上記状態空間モデルによる行列計算可能なプログラムを実行するコンピュータである。また、上記適合性判断手段は、上記推定出力データと、上記希望出力データとを比較し、上記対象システムの運用上、その差が許容できるものであるか否かを、上記基準条件を満たすか否かにより判断するプログラムを実行するコンピュータである。ここで、上記基準条件は、上記入力手段から入力されるものであり、上記対象システムの管理者が直接入力するものであってもよいし、それまでの上記対象システムの運用実績から算出されるものでもよい。
【００３６】
さらに、上記システム同定化手段、上記出力算出手段、上記適合性判断手段は、それぞれの役割を担ったプログラムが、単一のコンピュータ上で実行されるものでも良い。これらの役割を担い、かつ希望出力データの近似値を実際の出力データとして得ることができるように、上記入力手段を調製するシステム制御プログラムとすると、この発明にかかるシステムモデル構築方法を有用に上記対象システムに適用することができる。
【００３７】
このようなシステム運用装置は、上記の化学プラントのシステムにおいて使用することができ、上記と同様に、モデル化する入力データ及び出力データとしては、製造材料及び姿勢物の供給量、成分比、温度、圧力などの、数値化できる値のうち、入力では出力への関与が無視できない要素を、出力では入力からの関与が無視できない要素を選択するとよい。
【００３８】
以下、この発明でモデルを構築するための具体的なパラメータの同定手順を説明する。
まず、下記式（１）のようなイノベーションモデル形式で表される線形時不変な状態空間モデルを想定する。なお、イノベーションモデルとは、推測できない外乱データを含むモデルをいう。
【００３９】
【数５】
2007260504.tif 000003
【００４０】
ここで、ｙ（ｔ）∈Ｒnyはシステムの出力、ｕ（ｔ）∈Ｒnuはシステムの入力、ｎｙはｙの次数、ｎｕはｕの次数であり、ｘ（ｔ）∈Ｒnは次数ｎのベクトルである状態変数であり、ｅ（ｔ）∈Ｒnyは平均値がゼロで共分散行列が未知あるイノベーション過程であり、それぞれの行列は、行列Ａ∈Ｒn×n、行列Ｂ∈Ｒn×nu、行列Ｃ∈Ｒny×n、行列Ｋ∈Ｒn×nyである。前記イノベーション過程ｅ（ｔ）は、時々刻々システム内に入ってくる外乱データであり、共分散行列が未知であるため、入出力データから推定されるものである。ここで、離散時間システムに対象を絞ることで、同時刻における入力ｕから出力ｙに直結する直達項Ｄをゼロとおき、下記式（２）のモデルにより議論を進める。
【００４１】
【数６】
2007260504.tif 000004
【００４２】
上記式（２）を、下記式（３）のようなプレディクター形式に整理し直す。
なお、プレディクター形式とは、入力ｕ（ｔ）と出力ｙ（ｔ）とをまとめて新たな入力ｚ（ｔ）として扱うオブザーバー形式の中でも、特に、後述するカルマンゲインＫを用いて表現したカルマンプレディクター形式のことをいう。すなわち、以下の解析では、システムの入出力を一つにまとめたものとみなす観測器として解析する。なお、この式（３）は、下記の式（４）乃至（６）を前提とする。
【００４３】
【数７】
2007260504.tif 000005
【数８】
2007260504.tif 000006
【数９】
2007260504.tif 000007
【数１０】
2007260504.tif 000008
【００４４】
ここで、上記式（４）が含むパラメータＫは、上記式（２）におけるイノベーション過程ｅの係数行列であるが、これは、十分な時間経過がされると収束し定数となる定常カルマンゲインと呼ばれる行列であるので、（Ａ−）（上記式（４）の左辺を示す。本文中の表記上、行列記号の上にバーを有する記号を以下同様に表す。）は漸近的に安定であるとみなすことができる。従って、上記式（３）から状態ベクトルｘは、不安定なイノベーション過程を排除し、それより以前の時点における状態変数ｘと入力ｕ、出力ｙとで表すことが出来、具体的には、下記式（７）のように書くことができる。ここでｐは正数である。
【００４５】
【数１１】
2007260504.tif 000009
【００４６】
ここで、（Ａ−）が漸近的に安定であるので、十分大きなある正数ｐに対して、上記式（７）の右辺第一項は（Ａ−）p＝０と近似することが出来る。従って右辺第一項を０と近似することができ、状態ベクトルｘ（ｔ）は、過去の入力ｕ（ｔ）、出力ｙ（ｔ）、及びカルマンプレディクターマルコフパラメータ（Ａ−）及び（Ｂ−）から、下記式（８）により推定することが出来る。
【００４７】
【数１２】
2007260504.tif 000010
【００４８】
ここで、（ｘ∧）（上記式（８）の左辺を示す。本文中の表記上、行列記号の上に∧を有する記号を以下同様に示す。）は、状態変数ｘの推定値を示す。次に、上記式（２）を、上記式（８）で推定された状態変数で置き換えることで、出力ｙを下記式（９）のように近似することが出来る。
【００４９】
【数１３】
2007260504.tif 000011
【００５０】
入出力の時系列データに対して、上記式（９）は、下記式（１０）のようにまとめることができる。なお、下記式（１０）中の記号は下記式（１１）乃至（１４）の通りであり、式中Ｎはデータの数であり、正数である。
【００５１】
【数１４】
2007260504.tif 000012
【数１５】
2007260504.tif 000013
【数１６】
2007260504.tif 000014
【数１７】
2007260504.tif 000015
【数１８】
2007260504.tif 000016
【００５２】
また、上記式（１０）はｙtについての近似式であるが、ｕtとｅtに関しても上記（１３）と同様の行列を定義することとする。上記式（１０）における、カルマンプレディクターマルコフパラメータＣ（Ｌ−）zは、下記式（１５）で示される最小二乗解を有する。この計算手順は、上記式（１０）に対して、両辺にＺｐの逆行列を右からかけることで求まる。この捜査は最小二乗法と等価である。また、ｅtとＺpとの積は、ｅtがイノベーション、すなわち白色の性質であるために、積の期待値が０となる。
なお、下記式（１５）中のＲは、上記入出力データのサンプリングデータに対する相関行列であり、下記式（１６）及び（１７）で表される。
【００５３】
【数１９】
2007260504.tif 000017
【数２０】
2007260504.tif 000018
【数２１】
2007260504.tif 000019
【００５４】
この式（１５）により推定される、カルマンプレディクターマルコフパラメータ要素は、上記式（５）及び上記式（１１）より、下記式（１８）のように表すことが出来る。なお、上記式（１５）は、いわゆるＡＲＸモデルに対応するものである。
【００５５】
【数２２】
2007260504.tif 000020
【００５６】
ここで、従来のＯＫＩＤアルゴリズムの場合は、上記式（１５）のマルコフパラメータから、システムのマルコフパラメータを導き出し、そのマルコフパラメータからハンケル行列を作成し、そのハンケル行列にＥＲＡ処理を施して最終的にシステムの行列Ａ，Ｂ，Ｃ，Ｋを算出する。しかし、この発明で用いるシステムモデル構築方法では、上記のように、入出力をｚ（ｔ）としてまとめたオブザーバであるカルマンプレディクターマルコフパラメータを、未来部分へ状態空間モデルを拡張する部分空間法を用いて、ハンケル行列を求めてＥＲＡ処理する。外見上の手順は従来のＯＫＩＤよりも煩雑であるが、直接に行列Ａ，Ｂ，Ｃ，Ｋを求めることができ、精度が高いため、結果として数値解析的に安定となる。
【００５７】
このように推定した上記式（３）の状態空間モデルについて、時間を、基本となる時刻ｔから未来時刻（ｔ＋ｆ−１）まで拡張して、部分空間法によるカルマンプレディクターマルコフパラメータの解析を行う。その手順は以下の通りである。まず、上記式（３）を下記式（１９）のように拡張する。それぞれの記号の内容は下記式（２０）乃至（２５）のようになる。
【数２３】
2007260504.tif 000021
【数２４】
2007260504.tif 000022
【数２５】
2007260504.tif 000023
【数２６】
2007260504.tif 000024
【数２７】
2007260504.tif 000025
【数２８】
2007260504.tif 000026
【数２９】
2007260504.tif 000027
【００５８】
上記式（１５）において推定されたカルマンプレディクターマルコフパラメータは、上記式（２２）の行列（Φ−）を構成する要素に利用する。また、状態関数ｘからなる状態ベクトルＸ（ｔ）は、可制御行列（Ｌ−）z及び、過去データＺpから、下記式（２６）のように推定される。
【００５９】
【数３０】
2007260504.tif 000028
【００６０】
この上記式（２６）を用いると、上記式（１９）は下記式（２７）のように書くことが出来る。この式は上記式（１０）を、未来時刻（ｔ＋ｆ−１）まで拡張したものに相当する。
【００６１】
【数３１】
2007260504.tif 000029
【００６２】
次に、カルマンプレディクターマルコフパラメータの解析のために、カルマンプレディクターマルコフパラメータからなるハンケル行列を定義して、その推定を行う。このハンケル行列は下記式（２８）で定義される、ｆｎy×ｐ（ｎu＋ｎy）の行列である。
【００６３】
【数３２】
2007260504.tif 000030
【００６４】
上記の式（２８）において、ｋ＝０とおくと、上記式（２８）は、下記式（２９）となる。
【００６５】
【数３３】
2007260504.tif 000031
【００６６】
次に、上記式（２７）において、未来の入力に寄与する部分を、未来の出力に関する部分から取り除く。すなわち、右辺のＺf-1の項を左辺へ移項する。この出力に関する部分を下記式（３０）の通りに（Ｙ〜）ｆと定義する。
【００６７】
【数３４】
2007260504.tif 000032
【００６８】
また一方で、上記サンプリングデータに対する相関関数である上記式（１７）を用いて最小二乗法により、（Ｈ−）0は下記式（３１）のように推定することができる。ただし、（Ｒ∧）fpは下記式（３２）で定義される相関関数である。
【００６９】
【数３５】
2007260504.tif 000033
【数３６】
2007260504.tif 000034
【００７０】
また、未来時刻に関する時間要素ｆにおけるサンプリングデータに対する相関行列として、（Ｒ∧）ffを下記式（３３）のように定義すると、同様に、サンプリングデータに対する相関行列である上記式（１７）と上記式（２７）とを用いて、（Ｈ−）0を下記式（３４）のように変形することができる。ただし、規格化係数Ｗ1、Ｗ2は下記式（３５）であり、Ｇは下記式（３６）である。
【数３７】
2007260504.tif 000035
【００７１】
【数３８】
2007260504.tif 000036
【数３９】
2007260504.tif 000037
【数４０】
2007260504.tif 000038
【００７２】
ここで、上記式（３６）は信号（Ｙ〜）fとＺpとの正準相関行列（ｃａｎｏｎｉｃａｌ  ｃｏｒｒｅｌａｔｉｏｎ  ｍａｔｒｉｘ）となっている。このＧに対して特異値分解を施す正準相関解析（ｃａｎｏｎｉｃａｌ  ｃｏｒｒｅｌａｔｉｏｎ  ａｎａｌｙｓｉｓ、以下、「ＣＣＡ」と略す。）を行う。システムの次元をｎと仮定して特異値分解を実施すると、下記式（３７）のようになる。
【００７３】
【数４１】
2007260504.tif 000039
【００７４】
ここで、ＵとＶとは、それぞれ正規直行行列であり、Ｓは対角成分が特異値からなる対角行列で、その特異値は対角成分に沿ってその値が順に小さくなっていくものである。また、ＵnとＶnは、それぞれＵ及びＶの最初のｎ列からなる行列である。さらに、ＳnはＳのｎ個の特異値からなる対角行列である。これにより上記式（３４）から、（Ｈ−∧）0は下記式（３８）のように書き表すことができる。
【数４２】
2007260504.tif 000040
【００７５】
このように特異値分解を施した係数行列に対して、正規化するための規格化係数であるＷ1とＷ2とを、一次元において正規化する際の分母に相当する逆行列として含む、すなわちＣＣＡ重みを含んだ形式で、以下のＥＲＡ処理を行う。
【００７６】
まず、行列（Ｈ−）0をムーアペンローズ形疑似逆行列Ｈ†を用いて記述できるように構築すると下記式（３９）のようになる。
【００７７】
【数４３】
2007260504.tif 000041
【００７８】
ここで、ＵnTＵn＝Ｉn、及びＶnTＶn＝Ｉnという関係を用いて、上記式（３８）に上記式（３９）の関係を適用させると、下記式（４０）のようになる。
【００７９】
【数４４】
2007260504.tif 000042
【００８０】
次に、上記式（３９）において行列（Ｈ−）0を記述するために用いる、（Γ−）fと（Ｌ−）zの選び方の一つに、並行実現の考え方から、（Γ−）f＝Ｗ1-1ＵnＳn1/2並びに、（Ｌ−）z＝Ｓn1/2ＶnTＷ2のように決める方法が挙げられる。また、下記式（４１）の関係が成立する。
【００８１】
【数４５】
2007260504.tif 000043
【００８２】
さらに、Ｉiをｉ次の単位行列とし、Ｏiをｉ次の零行列とする。ｎuは入力変数の数、ｎyは出力変数の数であることから、次の下記式（４２）及び（４３）のような行列を定義する。
【００８３】
【数４６】
2007260504.tif 000044
【数４７】
2007260504.tif 000045
【００８４】
これらの定義を用いて、ＣＣＡ重み付きのＥＲＡ処理を行うと、その過程及び結果は、下記式（４４）のようになる。
【００８５】
【数４８】
2007260504.tif 000046
【００８６】
なお、上記式の三行目から四行目への計算は、上記式（４１）を適用している。上記式（４４）から、カルマンプレディクターマルコフパラメータのうち、（Ａ−）、（Ｂ−）、Ｃは下記式（４５）乃至（４７）のように計算される。
【００８７】
【数４９】
2007260504.tif 000047
【数５０】
2007260504.tif 000048
【数５１】
2007260504.tif 000049
【００８８】
ただし、（Ａ−）はハンケル行列である（Ｈ−）1を含んでおり、別途この（Ｈ−）1を求める必要がある。ハンケル行列の定義式である上記式（２８）より、（Ｈ−）1は次のようになる。
【数５２】
2007260504.tif 000050
【００８９】
この上記式（４８）の二行目は、行列（Ｈ−）0の可制御行列内の２番目からｐ番目までの列ブロックを選び、（Ｈ−）１の可制御行列内の１番目からｐ−１番目までの列ブロックに割り当て、最後のｐブロックは（Ａ−）pを０と近似できることから、零行列を割り当てることで求められる。
【００９０】
従って、最初に定義した上記式（１）の状態空間モデルのカルマンプレディクターマルコフパラメータは、定義した上記式（４）乃至（６）から、下記式（４９）乃至（５２）のようにして導くことができる。なお、表現法はＭＡＴＬＡＢ（登録商標）の記述法を用いている。
【００９１】
【数５３】
2007260504.tif 000051
【数５４】
2007260504.tif 000052
【数５５】
2007260504.tif 000053
【数５６】
2007260504.tif 000054
【００９２】
このようにしてカルマンプレディクターマルコフパラメータを求めることで、この発明にかかるシステムモデル構築方法を実現することができる。この発明にかかる方法によることで、従来のＳＳＡＲＸに比べて最小二乗法の使用回数を減らし、数値的により安定したモデルの同定と、それによる出力データの推定が可能となる。
【実施例１】
【００９３】
以下、この発明を用いた実施例について具体例を挙げて説明する。
閉ループデータの適用検証用システムとして、Ｖｅｒｈａｅｇｅｎ，１９９３で取り上げられた例題を用いる。このシステムは、上記式（２）の形式でモデル化される、下記式（５３）のパラメータを有する５次のシステムである。
【００９４】
【数５７】
2007260504.tif 000055
【００９５】
このパラメータから、状態空間モデルで記述したコントローラは上記式（２）におけるパラメータが次の式（５４）の通りに表すことができる。なお、コントローラとは、出力を入力にフィードバックし、出力を望み通りに調節するものであり、後述する図３中のＫ（ｑ）に相当する。
【００９６】
【数５８】
2007260504.tif 000056
【００９７】
以上を決定し、単位円のプロットを行った。その結果を図１に示す。なお、図１中の曲線が単位円であり、横軸は実数軸、縦軸は虚数軸を表している。この対象プロセスに、モンテカルロ法を用いて、参照信号ｒは分散１のガウス性白色信号で、イノベーション過程ｅは分散１／９のガウス性白色信号として作成した１２００点からなる閉ループデータを、１００集合作成して、１００個のデータセットをそれぞれ作る。
【００９８】
このデータセットを、実施例としてこの発明にかかるシステムモデル構築方法によりモデル化したモデル（図１中「ＳＯＫＩＤ」と表記する。）に導入し、その固有値を計算し、図１にプロットした。なお、このモデルについて、過去のラグｐ及び未来のラグｆの次数を決定する方法としては、ＦＯＥ（Ｆｉｎａｌ  ｏｕｔｐｕｔ  ｅｒｒｏｒ  ｃｒｉｔｅｒｉｏｎ）（Ｚｈｕ，２００１）により決定を行った。なお、ＦＯＥとは、実測値と同定されたモデルの推定出力値とを比較して、その誤差の二乗平均値がもっとも小さくなる値を探索する方法をいう。ただし、上記式（５３）の次数ｎは５であるが、ｎはｐ又はｆのいずれか小さい方の値より大きくすることはできないという条件を満たす必要がある。１００個のデータセットに対して計算した結果、過去のラグｐ及び未来のラグｆをそれぞれ２０とした。
【００９９】
同様に、比較例として、従来のＳＳＡＲＸ及びＯＫＩＤによりモデル化したモデルに、上記と同様のデータセットを導入し、その固有値を計算して図１にプロットした。ただしＯＫＩＤについては部分空間法ではないのでＦ＝１とした。また、行列計算ソフトＭａｔｌａｂ（登録商標）に搭載されている、システム同定ツールボックスの部分空間法関数であるＮ４ＳＩＤによりモデル化したモデルにおいて、同様のデータセットを導入し、その固有値を計算して図１にプロットした。
【０１００】
その結果、この発明にかかるシステムモデル構築方法を用いたＳＯＫＩＤの固有値のばらつきは、ＳＳＡＲＸとほぼ同程度であり、誤差がＳＳＡＲＸと同程度まで少ないことが示されている。また、ＯＫＩＤはこれらと比較するとばらつきが大きくて、誤差が多いことがわかる、また、Ｎ４ＳＩＤは、ばらつきは小さいが、偏りが生じており、閉ループの同定がうまくいっていないことがわかる。したがってこの発明にかかるシステムモデル構築方法によるモデル化は、従来の方法と比較して、同等以上の精度でのシミュレーションが可能であることがわかった。
【実施例２】
【０１０１】
図２に記載のプラントについて、実データを用いてシステムのモデル化を行った。まずこのプラントについて説明する。このプラントは、製造材料ガスＦ（ＦｒｅｓｈＦｅｅｄ）の量をバルブ１１で調製し、このバルブが、製造材料ガスが熱交換器１２及び熱交換器１３を通過した後に設けた温度コントローラからの指令により調製され、適切な温度に温められるようにして、反応器Ｒ１に導入されるようになっている。反応機Ｒ１の底から反応した高温のガスは、途中で分離させ、一部は熱交換器１３で反応器Ｒ１に導入する前のガスと熱交換して冷やした後に、反応器Ｒ２に導入し、残りは反応器Ｒ１から反応器Ｒ２へ直接導入している。ここで、反応器Ｒ１から反応器Ｒ２へ直接導入するラインのバルブ１４は、熱交換器１３を通過したガスと通過してないガスとの混合ガスの温度を温度コントローラＴＣ２により制御され、反応器Ｒ２で望ましい温度になるようにバルブ１４が調製されている。さらに、反応器Ｒ２から出た製品ガスＰ（Ｐｒｏｄｕｃｔ）が有する熱量は、熱交換器１２において、前記製造材料ガスに受け渡されている。
【０１０２】
なお、反応器Ｒ１及びＲ２からの排出温度は、温度センサーＴＩ１及びＴＩ２で計測しており、それぞれへの入力を外乱変数ＤＶ１、ＤＶ２とする。また、バルブ１４への操作信号を外乱変数ＤＶ３とする。
【０１０３】
また、それぞれの温度コントローラＴＣ１及びＴＣ２は、温度の指示値を測定し，オペレータが与えた設定値と指示値との偏差が減少するように、バルブ１１及び１４へ操作すべき信号値を計算し、それぞれのバルブへ信号線を介して出力操作する。具体的には、ＰＩＤコントローラを用いる。
【０１０４】
ここで、出力ｙ（ｔ）を温度コントローラＴＣ１が検知するガスの温度の指示値とし、入力ｕ（ｔ）をバルブ１１で調製するガスの供給量として対象システムを定義する。このシステムの動作は、図３のように表すことができる。熱交換器１２及び熱交換器１３がシステムＧ（ｑ）であり、これに外乱変数であるＤＶ１乃至ＤＶ３をまとめた外乱ｖ（ｔ）が関与して出力ｙ（ｔ）がなされる。またこの出力ｙ（ｔ）がフィードバックされ、設定値ｒ（ｔ）と出力ｙ（ｔ）との偏差が減少するようにコントローラＫ（ｑ）が入力ｕ（ｔ）を制御する。
【０１０５】
そのときの入力ｕ（ｔ）と出力ｙ（ｔ）との対応は図６のようになっている。横軸は１分周期で測定したデータ点数を示し、縦軸は、ｙ（ｔ）については温度（℃）、ｕ（ｔ）についてはコントローラ出力（％）となる。また一方で、コントローラの制御をＯＦＦにしてプラントテストを行ったデータ、すなわち、開ループで得たデータを図４に示す。この、プロセスの入力ｕ（ｔ）を上げ下げして、プロセスの応答を取ることで得られたプラントテストデータから、得られた応答を基に解析した結果を、周波数領域でのステップ応答モデルとして図５の太線に示す。これをプロセスの真の値として扱い、以下の手順でモデル化した値と比較する。
【０１０６】
次に、図２のプロセスで、コントローラＫ（ｑ）を用いた通常運転、すなわち閉ループでの運用データを図６に示す。
【０１０７】
一方で、図６に示す閉ループデータを用いてプロセスの同定を行い、周波数領域でのステップ応答モデルを求めた。実施例１と同様にＦＯＥ規範を用いて、ＦＯＥが最小になるようにふさわしい次数ｎを９として、この発明にかかるシステムモデル構築方法（ＳＯＫＩＤ）で推定した出力値を図５中実線で示す。なお、データのサンプリング周期は１分間で２５００点であり、後述する比較例も同様である。
【０１０８】
また、比較例として、実施例１と同様に、ＳＳＡＲＸ、ＯＫＩＤ、Ｍａｔｌａｂ−Ｎ４ＳＩＤ（ＣＶＡ重み）でモデル化したモデルについて、上記のＳＯＫＩＤについての実施例と同様にＦＯＥ規範を用いて、ＦＯＥが最小となる次数を選択しモデル化して推定した、周波数に対する出力値を求めた。その結果を図５に示す。なお、ダッシュドット線がＳＳＡＲＸ、点線はＯＫＩＤ、破線がＭａｔｌａｂ−Ｎ４ＳＩＤである。
【０１０９】
さらに、図５において実線で表される真のプロセスモデルに対する、各々のモデル化による推定値の偏差の絶対値を計算したものを図７に示す。同様に、実線はこの発明にかかるシステムモデル構築方法（ＳＯＫＩＤ）であり、ダッシュドット線はＳＳＡＲＸ、点線はＯＫＩＤ、破線がＭａｔｌａｂ−Ｎ４ＳＩＤである。この絶対値について、周波数に対する二乗誤差平均を計算すると、ＳＯＫＩＤが０．２２、ＳＳＡＲＸが０．３４、ＯＫＩＤは０．２９、Ｎ４ＳＩＤは０．３１となった。すなわち、このシステムのモデル化にあたっては、この発明にかかるシステムモデル構築方法によるモデル化を行って推定した場合の推定値が、最も真のプロセスの値に近い値が得られた。
【図面の簡単な説明】
【０１１０】
【図１】実施例における真の値と推定値の固有値分布を示すグラフ
【図２】実施例２でモデル化するプラントの概略図
【図３】実施例２における出力からのフィードバックの仕組みを示す閉ループの概略図
【図４】実施例２におけるコントローラの制御をＯＦＦにした開ループでの入力ｕと出力ｙのサンプリングデータを示すグラフ
【図５】実施例２におけるゲインの周波数応答の真の値と推定値を示すグラフ
【図６】実施例２で用いるコントローラの制御をＯＮにした閉ループデータでの外乱変数を含むサンプリングデータを示すグラフ
【図７】図５における真の値との偏差の絶対値を示すグラフ
【符号の説明】
【０１１１】
ｒ（ｔ）  設定値
ｕ（ｔ）  入力
ｖ（ｔ）  外乱
ｙ（ｔ）  出力
ＤＶ１，ＤＶ２，ＤＶ３  外乱変数
Ｇ（ｑ）  システム
Ｋ（ｑ）  コントローラ
Ｆ  製造材料ガス
Ｐ  製品ガス
Ｒ１，Ｒ２  反応器
ＴＣ１，ＴＣ２  温度コントローラ
ＴＩ１，ＴＩ２  温度センサー
１１，１４  バルブ
１２，１３  熱交換器
【図１】
2007260504.tif 000060
【図２】
2007260504.tif 000061
【図３】
2007260504.tif 000062
【図４】
2007260504.tif 000063
【図５】
2007260504.tif 000064
【図６】
2007260504.tif 000065
【図７】
2007260504.tif 000066