"装入された鉄鉱石とコークスを含む原材料から、複数の化学反応からなる還元反応により溶銑を製造する高炉の操業において、前記原材料が固相、液相及び気相の部分からなる装入物の状態を含む高炉の内部状態の指標値であって、高炉の操業状態を表す特徴量である高炉操業特徴量の時間変化を、メッシュで複数の領域に分割した高炉内部それぞれの領域での反応を記述するモデルを用いたシミュレーションにより演算して、高炉の操業状態を制御する高炉操業状態制御システムであって、高炉に設けられた複数のセンサ及び高炉を制御する制御手段と、前記センサ及び前記制御手段から取得した、高炉の操業状態に関するデータである高炉操業情報と、高炉内部の前記固体相、液体相及び気体相の各部分について、それぞれの時間的分解能で、またはそれぞれの時間的分解能および空間的分解能で前記シミュレーションを行い、前記高炉操業特徴量を演算する高炉操業特徴量演算装置と、前記高炉操業情報及び前記高炉操業特徴量に基づいて、高炉の熱バランスとマスバランスを含む操業プロセスに関する高炉の操業状態の安定性を評価するための定量的な指標である操業安定性評価情報を演算して取得し、前記操業安定性評価情報より判別される高炉の操業状態から、高炉の操業を安定させる前記制御手段の操作量を決定する操業制御装置と、を備え、前記センサ及び前記制御手段と、前記高炉操業特徴量演算装置と、前記操業制御装置との相互連携により、動的システムである高炉の安定性評価に基づき高炉の操業状態を制御することを特徴とする、高炉操業状態制御システム。前記操業制御装置は、操業動作点における高炉の操業安定性評価モデルを取得し、前記操業安定性評価モデルより高炉の操業状態の安定性を表す操業安定性評価情報を取得する安定性評価部と、前記操業安定性評価情報に基づいて、高炉の操業状態が安定状態あるいは不安定状態のいずれにあるかを判別し、高炉の操業状態を安定させる前記制御手段の操作量を決定する制御方法決定部と、を備えることを特徴とする、請求項1に記載の高炉操業状態制御システム。前記安定性評価部は、前記高炉操業情報及び前記高炉操業特徴量に基づいて、高炉の動的システム特性を表す状態空間表現モデルの状態空間行列を同定し、当該状態空間行列を前記操業安定性評価モデルとして用い、制御理論に基づき定量的な指標である操業安定性評価情報を演算して取得することを特徴とする、請求項2に記載の高炉操業状態制御システム。前記安定性評価部は、前記状態空間表現モデルの状態空間行列のランクを操業安定性評価情報として演算して取得し、当該ランクが可制御であることを示す所定値であるか否かに基づいて、高炉の操業状態の安定性を評価することを特徴とする、請求項3に記載の高炉操業状態制御システム。前記安定性評価部は、前記状態空間表現モデルの状態空間行列の固有値を操業安定性評価情報として演算して取得し、当該固有値に基づいて、高炉の操業状態の安定性を評価することを特徴とする、請求項3に記載の高炉操業状態制御システム。前記安定性評価部は、同定された状態空間行列で表される状態空間表現モデルの周波数特性に基づいて、ナイキストの安定判別法を用いて高炉の操業状態の安定性を評価することを特徴とする、請求項3に記載の高炉操業状態制御システム。前記安定性評価部は、同定された状態空間行列で表される状態空間表現モデルに関し、リアプノフの安定判別法を用いて高炉の操業状態の安定性を評価することを特徴とする、請求項3に記載の高炉操業状態制御システム。高炉操業特徴量演算装置は、境界条件に外乱を加えた情報を、前記高炉操業特徴量を算出するための高炉の非定常物理モデルへの入力情報とし、当該入力情報に基づき前記非定常物理モデルを用いて算出された前記高炉操業特徴量から取得される外乱の影響を出力情報として取得して、前記入力情報および前記出力情報を前記操業制御装置へ出力し、前記操業制御装置の前記安定性評価部は、前記高炉操業情報、前記高炉操業特徴量、前記入力情報および前記出力情報に基づいて、4SID法を用いて、入出力信号のダイナミクスを近似する状態空間表現モデルを導出し、当該状態空間表現モデルを前記操業安定性評価モデルとして用い、制御理論に基づき定量的な指標である操業安定性評価情報を演算して取得することを特徴とする、請求項2に記載の高炉操業状態制御システム。装入された鉄鉱石とコークスを含む原材料から、複数の化学反応からなる還元反応により溶銑を製造する高炉の操業において、前記原材料が固相、液相及び気相の部分からなる装入物の状態を含む高炉の内部状態の指標値であって、高炉の操業状態を表す特徴量である高炉操業特徴量の時間変化を、メッシュで複数の領域に分割した高炉内部それぞれの領域での反応を記述するモデルを用いたシミュレーションにより演算して、高炉の操業状態を制御する操業制御装置であって、高炉に設けられた複数のセンサ及び高炉を制御する制御手段から取得した、高炉の操業状態に関するデータである高炉操業情報と、高炉操業特徴量演算装置により、高炉内部の前記固体相、液体相及び気体相の各部分について、それぞれの時間的分解能で、またはそれぞれの時間的分解能および空間的分解能で前記シミュレーションを行って演算された前記高炉操業特徴量とに基づいて、高炉の熱バランスとマスバランスを含む操業プロセスに関する高炉の操業状態の安定性を評価するための定量的な指標である操業安定性評価情報を演算して取得し、前記操業安定性評価情報より判別される高炉の操業状態から、高炉の操業を安定させる前記制御手段の操作量を決定することを特徴とする、操業制御装置。装入された鉄鉱石とコークスを含む原材料から、複数の化学反応からなる還元反応により溶銑を製造する高炉の操業において、前記原材料が固相、液相及び気相の部分からなる装入物の状態を含む高炉の内部状態の指標値であって、高炉の操業状態を表す特徴量である高炉操業特徴量の時間変化を、メッシュで複数の領域に分割した高炉内部それぞれの領域での反応を記述するモデルを用いたシミュレーションにより演算して、高炉の操業状態を制御する高炉操業状態制御システムによる高炉操業状態制御方法であって、高炉に設けられた複数のセンサ及び高炉を制御する制御手段から、高炉の操業状態に関するデータである高炉操業情報を取得するステップと、高炉操業特徴量演算装置によって、高炉内部の前記固体相、液体相及び気体相の各部分について、それぞれの時間的分解能で、またはそれぞれの時間的分解能および空間的分解能で前記シミュレーションを行い、前記高炉操業特徴量を演算するステップと、前記高炉操業情報及び前記高炉操業特徴量に基づいて、操業制御装置によって、高炉の熱バランスとマスバランスを含む操業プロセスに関する高炉の操業状態の安定性を評価するための定量的な指標である操業安定性評価情報を演算して取得し、前記操業安定性評価情報より判別される高炉の操業状態から、高炉の操業を安定させる前記制御手段の操作量を決定するステップと、を含み、前記センサ及び前記制御手段と、前記高炉操業特徴量演算装置と、前記操業制御装置とを相互連携させて、上記ステップを制御タイミングにおいて行うことにより、動的システムである高炉の安定性評価に基づき高炉の操業状態を制御することを特徴とする、高炉操業状態制御方法。装入された鉄鉱石とコークスを含む原材料から、複数の化学反応からなる還元反応により溶銑を製造する高炉の操業において、前記原材料が固相、液相及び気相の部分からなる装入物の状態を含む高炉の内部状態の指標値であって、高炉の操業状態を表す特徴量である高炉操業特徴量の時間変化を、メッシュで複数の領域に分割した高炉内部それぞれの領域での反応を記述するモデルを用いたシミュレーションにより演算して、高炉の操業状態を制御する操業制御装置の機能をコンピュータに実現させるためのプログラムであって、高炉に設けられた複数のセンサ及び高炉を制御する制御手段から取得した、高炉の操業状態に関するデータである高炉操業情報と、高炉操業特徴量演算装置により、高炉内部の前記固体相、液体相及び気体相の各部分について、それぞれの時間的分解能で、またはそれぞれの時間的分解能および空間的分解能で前記シミュレーションを行って演算された前記高炉操業特徴量とに基づいて、高炉の熱バランスとマスバランスを含む操業プロセスに関する高炉の操業状態の安定性を評価するための定量的な指標である操業安定性評価情報を演算して取得し、前記操業安定性評価情報より判別される高炉の操業状態から、高炉の操業を安定させる前記制御手段の操作量を決定することを特徴とする、操業制御装置として機能させるためのプログラム。ここで、高炉の内部の状態を精密に予測するためにはシミュレーションの時間刻みを細かく設定する必要がある一方、高炉の内部で進行する現象を予測するためにはシミュレーションを行う期間例えば、何時間後の状態までシミュレートするか等をある程度大きく設定することが求められる。従って、上記特許文献1~3に記載されているような従来のシミュレーションでは、ある期間のシミュレーションを行うために多大なリソース及び多大な計算時間が必要となり、高炉の内部状態をリアルタイムに予測することができないという問題があった。このため、従来は、高炉を線形定常モデルとして取り扱い、演算負荷を軽減させていたが、操業目標や操業動作点の変更や外乱に対する高炉の動的挙動を計算することはできなかった。したがって、高炉の、グローバルな操業目標や操業動作点の変更や外乱に対する大局的な操業状態を定量的に評価することもなされていなかった。そこで、本発明は、上記問題に鑑みてなされたものであり、本発明の目的とするところは、高炉の内部状態を示す操業の特徴量をシミュレーションにより求めて、高炉の操業状態の安定性を定量的かつ従来よりも短時間に評価することが可能な、高炉操業状態制御システム、操業制御装置、高炉操業状態制御方法およびプログラムを提供することにある。上記課題を解決するために、本発明のある観点によれば、装入された鉄鉱石とコークスを含む原材料から、複数の化学反応からなる還元反応により溶銑を製造する高炉の操業において、原材料が固相、液相及び気相の部分からなる装入物の状態を含む高炉の内部状態の指標値であって、高炉の操業状態を表す特徴量である高炉操業特徴量の時間変化を、メッシュで複数の領域に分割した高炉内部それぞれの領域での反応を記述するモデルを用いたシミュレーションにより演算して、高炉の操業状態を制御する高炉操業状態制御システムが提供される。高炉操業状態制御システムは、高炉に設けられた複数のセンサ及び高炉を制御する制御手段と、センサ及び制御手段から取得した、高炉の操業状態に関するデータである高炉操業情報と、高炉内部の固体相、液体相及び気体相の各部分について、それぞれの時間的分解能で、またはそれぞれの時間的分解能および空間的分解能でシミュレーションを行い、高炉操業特徴量を演算する高炉操業特徴量演算装置と、高炉操業情報及び高炉操業特徴量に基づいて、高炉の熱バランスとマスバランスを含む操業プロセスに関する高炉の操業状態の安定性を評価するための定量的な指標である操業安定性評価情報を演算して取得し、操業安定性評価情報より判別される高炉の操業状態から、高炉の操業を安定させる制御手段の操作量を決定する操業制御装置と、を備え、センサ及び制御手段と、高炉操業特徴量演算装置と、操業制御装置との相互連携により、動的システムである高炉の安定性評価に基づき高炉の操業状態を制御することを特徴とする。ここで、操業制御装置は、操業動作点における高炉の大局的なシステム特性を表す操業安定性評価モデルを取得し、操業安定性評価モデルより高炉の大局的な操業状態の安定性を表す操業安定性評価情報を取得する安定性評価部と、操業安定性評価情報に基づいて、高炉の操業状態が安定状態あるいは不安定状態のいずれにあるかを判別し、高炉の操業状態を安定させる制御手段の操作量を決定する制御方法決定部と、を備えてもよい。また、安定性評価部は、高炉操業情報及び高炉操業特徴量に基づいて、高炉の動的システム特性を表す状態空間表現モデルの状態空間行列を同定し、当該状態空間行列を操業安定性評価モデルとして用い、制御理論に基づき定量的な指標である操業安定性評価情報を演算して取得してもよい。安定性評価部は、状態空間表現モデルの状態空間行列のランクを操業安定性評価情報として演算して取得し、当該ランクが可制御であることを示す所定値であるか否かに基づいて、高炉の操業状態の安定性を評価してもよい。あるいは、安定性評価部は、状態空間表現モデルの状態空間行列の固有値を操業安定性評価情報として演算して取得し、当該固有値に基づいて、高炉の操業状態の安定性を評価してもよい。または、安定性評価部は、同定された状態空間行列で表される状態空間表現モデルの周波数特性に基づいて、ナイキストの安定判別法を用いて高炉の操業状態の安定性を評価してもよい。さらに、安定性評価部は、同定された状態空間行列で表される状態空間表現モデルに関し、リアプノフの安定判別法を用いて高炉の操業状態の安定性を評価してもよい。高炉操業特徴量演算装置は、境界条件に外乱を加えた情報を、高炉操業特徴量を算出するための高炉の非定常物理モデルへの入力情報とし、当該入力情報に基づき非定常物理モデルを用いて算出された高炉操業特徴量から取得される外乱の影響を出力情報として取得して、入力情報および出力情報を操業制御装置へ出力してもよい。このとき、操業制御装置の安定性評価部は、高炉操業情報、高炉操業特徴量、入力情報および出力情報に基づいて、4SID法を用いて、入出力信号のダイナミクスを近似する状態空間表現モデルを導出し、当該状態空間表現モデルを操業安定性評価モデルとして用い、制御理論に基づき定量的な指標である操業安定性評価情報を演算して取得してもよい。また、上記課題を解決するために、本発明の別の観点によれば、装入された鉄鉱石とコークスを含む原材料から、複数の化学反応からなる還元反応により溶銑を製造する高炉の操業において、原材料が固相、液相及び気相の部分からなる装入物の状態を含む高炉の内部状態の指標値であって、高炉の操業状態を表す特徴量である高炉操業特徴量の時間変化を、メッシュで複数の領域に分割した高炉内部それぞれの領域での反応を記述するモデルを用いたシミュレーションにより演算して、高炉の操業状態を制御する操業制御装置が提供される。操業制御装置は、高炉に設けられた複数のセンサ及び高炉を制御する制御手段から取得した、高炉の操業状態に関するデータである高炉操業情報と、高炉操業特徴量演算装置により高炉内部の固体相、液体相及び気体相の各部分について、それぞれの時間的分解能で、またはそれぞれの時間的分解能および空間的分解能でシミュレーションを行って演算された高炉操業特徴量とに基づいて、高炉の熱バランスとマスバランスを含む操業プロセスに関する高炉の操業状態の安定性を評価するための定量的な指標である操業安定性評価情報を演算して取得し、操業安定性評価情報より判別される高炉の操業状態から、高炉の操業を安定させる制御手段の操作量を決定することを特徴とする。さらに、上記課題を解決するために、本発明の別の観点によれば、装入された鉄鉱石とコークスを含む原材料から、複数の化学反応からなる還元反応により溶銑を製造する高炉の操業において、原材料が固相、液相及び気相の部分からなる装入物の状態を含む高炉の内部状態の指標値であって、高炉の操業状態を表す特徴量である高炉操業特徴量の時間変化を、メッシュで複数の領域に分割した高炉内部それぞれの領域での反応を記述するモデルを用いたシミュレーションにより演算して、高炉の操業状態を制御する高炉操業状態制御システムによる高炉操業状態制御方法が提供される。高炉操業状態制御方法は、高炉に設けられた複数のセンサ及び高炉を制御する制御手段から、高炉の操業状態に関するデータである高炉操業情報を取得するステップと、高炉操業特徴量演算装置によって、高炉内部の固体相、液体相及び気体相の各部分について、それぞれの時間的分解能で、またはそれぞれの時間的分解能および空間的分解能でシミュレーションを行い、高炉操業特徴量を演算するステップと、高炉操業情報及び高炉操業特徴量に基づいて、操業制御装置によって、高炉の熱バランスとマスバランスを含む操業プロセスに関する高炉の操業状態の安定性を評価するための定量的な指標である操業安定性評価情報を演算して取得し、操業安定性評価情報より判別される高炉の操業状態から、高炉の操業を安定させる制御手段の操作量を決定するステップと、を含み、センサ及び制御手段と、高炉操業特徴量演算装置と、操業制御装置とを相互連携させて、上記ステップを制御タイミングにおいて行うことにより、動的システムである高炉の安定性評価に基づき高炉の操業状態を制御することを特徴とする。また、上記課題を解決するために、本発明の別の観点によれば、装入された鉄鉱石とコークスを含む原材料から、複数の化学反応からなる還元反応により溶銑を製造する高炉の操業において、原材料が固相、液相及び気相の部分からなる装入物の状態を含む高炉の内部状態の指標値であって、高炉の操業状態を表す特徴量である高炉操業特徴量の時間変化を、メッシュで複数の領域に分割した高炉内部それぞれの領域での反応を記述するモデルを用いたシミュレーションにより演算して、高炉の操業状態を制御する操業制御装置の機能をコンピュータに実現させるためのプログラムであって、高炉に設けられた複数のセンサ及び高炉を制御する制御手段から取得した、高炉の操業状態に関するデータである高炉操業情報と、高炉操業特徴量演算装置により、高炉内部の固体相、液体相及び気体相の各部分について、それぞれの時間的分解能で、またはそれぞれの時間的分解能および空間的分解能でシミュレーションを行って演算された高炉操業特徴量とに基づいて、高炉の熱バランスとマスバランスを含む操業プロセスに関する高炉の操業状態の安定性を評価するための定量的な指標である操業安定性評価情報を演算して取得し、操業安定性評価情報より判別される高炉の操業状態から、高炉の操業を安定させる制御手段の操作量を決定することを特徴とする、操業制御装置として機能させるためのプログラムが提供される。以上説明したように本発明によれば、高炉の内部状態を示す操業の特徴量を、高炉内部の固体相、液体相及び気体相の各部分について、それぞれの時間的分解能で、またはそれぞれの時間的分解能および空間的分解能でシミュレーションにより求めることにより、高炉の操業状態の安定性を定量的かつ従来よりも短時間に評価することが可能な、高炉操業状態制御システム、操業制御装置、高炉操業状態制御方法およびプログラムを提供することができる。以下に添付図面を参照しながら、本発明の好適な実施の形態について詳細に説明する。なお、本明細書及び図面において、実質的に同一の機能構成を有する構成要素については、同一の符号を付することにより重複説明を省略する。なお、説明は、以下の順序で行うものとする。1高炉という動的システムについて1-1高炉内で進行している反応の概略1-2システムとしての取り扱い2第1の実施形態2-1高炉操業状態制御システムの概略について2-2高炉操業状態制御装置の全体構成について2-3高炉操業状態制御方法の全体的な流れについて2-5高炉操業特徴量の演算について2-5-1高炉操業特徴量演算部の構成について2-5-2高炉操業特徴量の演算方法について2-6高炉操業の安定性評価について2-6-1高炉操業状態制御システムの構成概念2-6-2操業制御部の構成について2-6-3高炉操業の安定性評価3本発明の実施形態に係る高炉操業状態制御装置のハードウェア構成について4実施例高炉という動的システムについてまず、図1を参照しながら、本発明の実施形態で着目する高炉という動的システムについて、簡単に説明する。図1は、本発明の実施形態で着目する高炉について説明するための説明図である。<高炉内で進行している反応の概略>図1に示したように、高炉は、円筒の徳利形状を有する竪型炉の一種であり、炉の頭頂部以下、炉頂部とも称する。から供給される原料と、炉の下方に設けられた羽口から供給される熱風により生成される還元性ガスとが反応する反応装置として機能する。炉頂部から供給される原料としては、主に、鉄鉱石や焼結鉱以下では、特に断りのない限り、鉄鉱石及び焼結鉱をまとめて鉄鉱石と称することとする。などの鉄酸化物、コークス、石灰石等がある。鉄鉱石は、高炉における反応で生成される銑鉄の鉄源となるものであり、コークスは、鉄鉱石の還元剤及び原料を溶解するための熱源として機能するだけでなく、高炉内の通気性を保持する役割を有している。また、石灰石は、鉄鉱石の脈石成分と反応して低溶融点を持ち流動性のよいスラグを生成するために添加される媒溶剤として機能する。高炉の内部では、図1に示したように、鉄鉱石及び石灰石からなる層と、コークスからなる層とが交互に積層されている。これらの原料は、図1に示したような積層状態を維持しつつ、炉の下方へと移動していく。また、図1に示した羽口からは、熱風及びコークスの補完還元剤として機能する微粉炭とが供給され、かかる熱風により微粉炭やコークスがガス化して、一酸化炭素や水素等からなる高温の還元性ガスが生成される。この高温の還元性ガスは、炉内を移動する上昇気流となって炉頂部へと吹き昇っていく。この還元性ガスにより炉内の鉄鉱石は還元されていき間接還元、更に、還元性ガスが有する熱によって固体から液体へと変化する。液体となった鉄分は、コークス層内を滴下しながらコークスの炭素によって更に還元され直接還元、炭素を5程度含む溶銑となる。図1に示した融着帯では、半溶融状態にある鉄分の間に固体コークスがスリット状に存在している部分であり、主にこの融着帯において、上述のような鉄分の相変化が生じている。このように、高炉という反応装置では、固体、液体、気体が共存して反応がダイナミック動的に進行している。また、高炉内で進行している反応には、後述するように吸熱反応と発熱反応とが存在しているため、この相反する2種類の反応間での熱バランス及びマスバランスを調節しながら、銑鉄を製造する操業が行われている。<システムとしての取り扱い>以上説明したように、高炉の内部では、固体、液体、気体の3つの相状態が共存しており、互いに関連しあいながら複数の化学反応が進行しており、高炉の内部で生じている現象は、非定常的な現象非定常現象であるといえる。また、高炉の内部で起こっている非定常現象は、秒単位で進行するものから、分単位、時間単位、日単位、週単位で進行するものまで、各種のものが存在している。従って、以下で説明する本発明の実施形態では、高炉という反応装置を、一つの動的なシステムとして取り扱うこととする。この際、高炉という動的システムは、固体状態にある物質を原料とした反応のサブシステム、液体状態にある物質を原料とした反応のサブシステム、気体状態にある物質を原料とした反応のサブシステムという、複数3つのサブシステムからなるシステムであると考えることができる。ここで、各サブシステムでは、高炉内の様々な部位で上述のような様々な時間単位換言すれば、時間分解能の反応が異なるダイナミクスで進行している。従って、高炉という動的システムは、より詳細には、異なる分解能を有する複数のサブシステムが相互に作用しあう異種間相互作用系動的システムであるといえる。また、異なる観点から考えると、高炉という異種間相互作用系多分解能動的システムは、複数の化学反応が進行しているシステムであって、かかる化学反応の一つ一つを、高炉というシステム全体を構成するサブシステムと捉えることもできる。更に異なる観点から考えると、高炉という異種間相互作用系多分解能動的システムで生じている様々な現象は、以下で説明するように、複数の偏微分方程式に基づく分布定数系システムであって、かかる偏微分方程式の一つ一つを、高炉というシステム全体を構成するサブシステムと捉えることも可能である。第1の実施形態<高炉操業状態制御システムの概略について>続いて、図2を参照しながら、本発明の第1の実施形態に係る高炉操業状態制御システムについて説明する。図2は、本実施形態に係る高炉操業状態制御システムの一例を示した説明図である。本実施形態に係る高炉操業状態制御システム1は、例えば図2に示したように、制御対象であり、先に簡単に説明したような各種反応が内部で進行している高炉3と、高炉3及び高炉3の近傍に設置されている各種のセンサ5と、高炉3を制御するための各種の制御手段7と、高炉操業状態制御装置10と、を含む。各種のセンサ5は、高炉3に関する各種のデータを計測する。このようなセンサ5として、例えば、高炉に装入される各種原料を秤量する秤量センサ、装入された原料の高さを測定する測定装置、高炉に供給される熱風の送風量、送風温度、送風湿分等を計測する各種センサ、炉頂部の気体の圧力を計測する圧力センサ、高炉の炉壁や炉床等の各部位の温度の分布を測定する、通常は多数の温度計、冷却水の給水圧力、給排水流量、排水温度等を計測する各種センサ等を挙げることができる。また、高炉3及び高炉3の近傍に設置されるセンサ5は、上記例に限定されるものではなく、上記のもの以外にも、出銑温度を測定するための各種温度計や、銑鉄中Si量、銑鉄中S量、スラグ組成、炉頂ガス成分等を特定するための各種分析装置等を挙げることができる。これらセンサ5によって測定された各種データは、後述する高炉操業状態制御装置10に出力され、高炉操業状態制御装置10における各種処理に利用される。制御手段7は、高炉3における製銑プロセスを制御するために用いられる各種の設備や装置である。このような制御手段7として、例えば、原料装入用ベルトコンベアの制御装置、原料装入用ホッパーや原料装入用シュートの制御装置、熱風を高炉に供給するための熱風炉の制御装置、冷却水の供給制御装置等を挙げることができる。また、本実施形態に係る制御手段7は、上記例に限定されるわけではなく、高炉3における製銑プロセスを制御するために用いられるあらゆる設備や装置が、本実施形態に係る制御手段7に相当する。これら制御手段7は、後述する高炉操業状態制御装置10から出力される制御信号に基づいて動作し、高炉3の操業制御運転制御を行う。高炉操業状態制御装置10は、高炉3及び高炉3の近傍に設けられた各種センサ5から取得した各種データを取得し、取得したこれらのデータを利用して、高炉の操業状態に関する特徴量を算出する。その後、高炉操業状態制御装置10は、取得したデータや算出した特徴量に基づいて、高炉操業状態の安定性に関する評価値を算出したり、高炉の制御方法を決定したりする。高炉操業状態制御装置10は、高炉の制御方法を決定した場合には、決定した制御方法に対応する制御信号を生成して、各種制御手段7に出力する。ここで、高炉操業状態の安定性とは、高炉の内部全体に亘る熱的バランスやマスバランスを含む各化学的、および物理的なプロセスのバランスに関する安定性を意味する。本実施形態に係る高炉操業状態制御装置10の詳細については、以下で改めて詳細に説明する。以上、図2を参照しながら、本実施形態に係る高炉操業状態制御システム1の概略について説明した。<高炉操業状態制御装置の全体構成について>次に、図3を参照しながら、本実施形態に係る高炉操業状態制御装置の全体構成について説明する。図3は、本実施形態に係る高炉操業状態制御装置の全体構成を示したブロック図である。本実施形態に係る高炉操業状態制御装置10は、以下で説明するように、高炉3の操業状態を表す特徴量の演算を行う高炉操業特徴量演算機能と、高炉3の操業状態の安定性を判別し、高炉3が安定して操業するように制御する高炉操業制御機能と、を有している。すなわち、本実施形態に係る高炉操業状態制御装置10は、高炉操業特徴量の演算を行う高炉操業特徴量演算装置の一例であり、高炉3の操業状態の安定性を判別し、高炉3が安定して操業するように制御する操業制御装置の一例でもある。本実施形態に係る高炉操業状態制御装置10は、例えば図3に示したように、高炉操業情報取得部101と、高炉操業特徴量演算部103と、操業制御部105と、表示制御部107と、通信制御部109と、記憶部111と、を主に備える。高炉操業情報取得部101は、例えば、CPUCentralProcessingUnit、ROMReadOnlyMemory、RAMRandomAccessMemory等により実現される。高炉操業情報取得部101は、後述する通信制御部109を介して、高炉3及び高炉3の近傍に設置されている各種センサ5から、当該センサ5が計測した各種測定データを取得する。また、高炉操業情報取得部101は、高炉操業の制御を行っている各種装置例えば、制御手段7等や制御コンピュータ等から、高炉の設計値や現在設定されている操業条件等の各種情報を取得することもできる。高炉操業情報取得部101が取得した各種データや各種情報は、高炉3の操業状態に関する様々な情報を含むものである。高炉操業情報取得部101は、取得したこれらのデータや情報を、高炉操業情報として後述する高炉操業特徴量演算部103及び操業制御部105に出力する。なお、高炉操業情報取得部101は、取得した高炉操業情報を、当該情報を取得した日時に関する時刻情報等と関連付け、履歴情報として後述する記憶部111に格納してもよい。高炉操業特徴量演算部103は、例えば、CPU、ROM、RAM等により実現される。高炉操業特徴量演算部103は、高炉操業情報取得部101が取得した高炉操業情報に基づいて、高炉の操業状態を表す特徴量の演算を行う。具体的には、高炉操業特徴量演算部103は、各種のセンサ等から入力される高炉操業情報及び後述する記憶部111等に格納されている各種文献値やデータベース等を利用して、着目する化学反応に関する反応速度を算出する。その後、高炉操業特徴量演算部103は、算出した反応速度を利用して、固体相、液体相、気体相それぞれの温度、移動速度、成分密度、成分体積、圧力等といった、各相の状態等を示すところの高炉の操業状態の指標であり、当該操業状態を特徴づける特徴量を算出する。すなわち、本実施形態に係る高炉操業特徴量演算部103は、高炉の操業状態をシミュレートするシミュレータとして機能する処理部であるともいえる。ここで、上述のような高炉操業特徴量の演算を行うためには、以下で説明するように、時間に依存する項を含む偏微分方程式を考慮する必要がある。かかる場合において、本実施形態に係る高炉操業特徴量演算部103では、固体相、液体相、気体相それぞれに関する特徴量を演算する際に、各相における時間刻みを同一幅とせず各相における時間刻みを互いに独立な値として取り扱う。このとき各相の時間刻みは、相内の現象が時間変化する速さ換言すれば、時定数に基づき、所望の演算精度が得られるように設定すればよい。高炉操業特徴量演算部103は、高炉操業特徴量の演算が終了すると、得られた高炉操業特徴量を、後述する操業制御部105に出力する。また、高炉操業特徴量演算部103は、得られた高炉操業特徴量を後述する表示制御部107に出力して、高炉操業状態制御装置10等が備える表示部に演算結果を表示させてもよい。更に、高炉操業特徴量演算部103は、得られた高炉操業特徴量に、当該特徴量を算出した日時等に関する時刻情報を関連付けて、履歴情報として記憶部111等に記録してもよい。なお、本実施形態に係る高炉操業特徴量演算部103については、以下で改めて詳細に説明する。操業制御部105は、例えば、CPU、ROM、RAM等により実現される。操業制御部105は、各種センサ5により計測された各種測定データや高炉3の操業条件等の各種情報等からなる高炉操業情報と、高炉操業特徴量演算部103により演算された高炉操業特徴量とに基づき、高炉3の操業状態の安定性を判別する。この際、操業制御部105は、適宜記憶部111に記憶された情報を参照して情報を取得することができ、処理を実行する過程において必要な情報を格納することもできる。本実施形態における操業制御部105は、高炉3の操業状態の安定性を、制御理論に基づき操業安定性評価情報を用いて定量的に評価する。操業安定性評価情報とは、例えば、高炉3のダイナミクスを表現する状態空間表現モデルから取得できる、高炉3の大局的な非定常システム特性を表す情報操業安定性評価モデルから、システムの安定性判別に関する制御理論に基づき取得される定量的な情報である。操業制御部105は、従来の高炉3の制御における、例えば、溶銑の温度を制御目標値に一定に制御するために装入物分布条件及び羽口への装入条件等の操作量を操作するといった制御を行うのではなく、高炉3という動的システム全体の安定性を定量的に評価し、当該システムが安定して操業するように操作することにより、所望の状態を実現する。操業制御部105は、高炉3の制御方法を決定すると、後述する通信制御部109を介して、各種制御手段7を操作する操作量が出力される。操作量に基づき各種制御手段7が操作されることにより、高炉3における各種プロセスを制御するための設備や装置が動作され、高炉3の操業状態を所望の状態に近づける。また、操業制御部105により取得された操業安定性評価モデルや操業安定性評価情報、安定性評価に基づき決定された高炉3の制御方法等は、後述する表示制御部107へ出力してもよい。これにより、これらの情報を高炉操業状態制御装置10の利用者が把握することが可能となる。なお、本実施形態に係る操業制御部105については、以下で改めて詳細に説明する。表示制御部107は、例えば、CPU、ROM、RAM等により実現される。表示制御部107は、高炉操業特徴量演算部103による演算結果、操業制御部105による安定性評価値の演算結果や制御方法の決定結果換言すれば、制御するための操作量の決定結果等を、高炉操業状態制御装置10が備えるディスプレイ等の出力装置に表示する際の表示制御を行う。また、表示制御部107は、後述する通信制御部109を介して、高炉操業状態制御装置10の外部に設けられた各種装置に設けられた出力装置に各種結果を表示する際の表示制御を行うこともできる。これにより、高炉操業状態制御装置10の利用者は、高炉の状態や高炉の制御等に関する操作情報を、その場で把握することが可能となる。通信制御部109は、例えば、CPU、ROM、RAM、通信装置等により実現される。通信制御部109は、本実施形態に係る高炉操業状態制御装置10と、外部に設けられた各種装置との間で行われる通信の制御を行う。ここで、高炉操業状態制御装置10の外部に設けられた各種装置として、各種センサ5や各種制御手段7を挙げることができる。また、通信制御部109は、これらの装置以外にも、高炉操業状態制御装置10が各種ネットワークを介して接続可能なサーバやコンピュータ等の装置や、高炉操業状態制御装置10に各種端子により直接接続された装置との間で行われる通信を制御可能である。記憶部111は、本実施形態に係る高炉操業状態制御装置10が備えるストレージ装置の一例である。記憶部111には、本実施形態に係る高炉操業状態制御装置10が、何らかの処理を行う際に保存する必要が生じた様々なパラメータや処理の途中経過等、または、各種のデータベースやプログラム等が、適宜格納されている。この記憶部111は、高炉操業情報取得部101、高炉操業特徴量演算部103、操業制御部105、表示制御部107、通信制御部109等が、自由に読み書きを行うことが可能である。以上、本実施形態に係る高炉操業状態制御装置10の機能の一例を示した。上記の各構成要素は、汎用的な部材や回路を用いて構成されていてもよいし、各構成要素の機能に特化したハードウェアにより構成されていてもよい。また、各構成要素の機能を、CPU等が全て行ってもよい。従って、本実施形態を実施する時々の技術レベルに応じて、適宜、利用する構成を変更することが可能である。なお、上述のような本実施形態に係る高炉操業状態制御装置の各機能を実現するためのコンピュータプログラムを作製し、パーソナルコンピュータ等に実装することが可能である。また、このようなコンピュータプログラムが格納された、コンピュータで読み取り可能な記録媒体も提供することができる。記録媒体は、例えば、磁気ディスク、光ディスク、光磁気ディスク、フラッシュメモリなどである。また、上記のコンピュータプログラムは、記録媒体を用いずに、例えばネットワークを介して配信してもよい。<高炉操業状態制御方法の全体的な流れについて>続いて、図4を参照しながら、本実施形態に係る高炉操業状態制御方法の全体的な流れについて説明する。図4は、本実施形態に係る高炉操業状態制御方法の全体的な流れの一例を示した流れ図である。本実施形態に係る高炉操業状態制御装置10の高炉操業情報取得部101は、高炉3及び高炉3の近傍に設けられた各種センサ5や、高炉操業の制御を行っている各種装置から測定データ、高炉の各部の寸法値や形状、高炉に設定されている操業条件等を含む高炉操業情報を取得するステップS11。高炉操業情報取得部101は、取得した高炉操業情報を、高炉操業特徴量演算部103及び操業制御部105に出力する。次に、高炉操業特徴量演算部103は、取得した高炉操業情報等を利用して、高炉の操業状態を特徴づける高炉操業特徴量の演算を実施するステップS13。高炉操業特徴量演算部103は、得られた高炉操業特徴量を、操業制御部105に出力する。続いて、操業制御部105は、高炉操業特徴量演算部103により演算された高炉操業特徴量や、高炉操業情報等を利用して、高炉操業の安定性評価を実施するステップS15。その後、操業制御部105は、高炉操業特徴量、高炉操業情報、安定性評価の結果等を利用して、高炉の制御方法を決定するステップS17。本実施形態に係る高炉操業状態制御装置10は、このような流れで処理を行うことで、高炉の操業状態を制御することができる。なお、高炉操業特徴量の演算処理、高炉操業の安定性評価処理、高炉の制御方法の決定処理の詳細な流れについては、以下で改めて詳細に説明する。また、図4においては、高炉操業の安定性評価処理及び高炉の制御方法の決定処理の双方が実施される場合の一例について図示しているが、安定性評価処理と操業方法の決定処理のいずれか一方が実施されてもよい。また、高炉操業の安定性評価処理及び高炉の制御方法の決定処理は、並行して実施されてもよい。以上、本実施形態に係る高炉操業状態制御方法の全体的な流れについて、簡単に説明した。<高炉操業特徴量の演算について>[高炉操業特徴量演算部の構成について]次に、図5~図10を参照しながら、本実施形態に係る高炉操業状態制御装置10が備える高炉操業特徴量演算部103の構成について、詳細に説明する。図5は、本実施形態に係る高炉操業特徴量演算部の構成を示したブロック図である。本実施形態に係る高炉操業特徴量演算部103は、図5に示したように、演算制御部121と、パラメータ設定部123と、反応速度算出部125と、特徴量算出部127と、を更に備える。演算制御部121は、例えば、CPU、ROM、RAM等により実現される。演算制御部121は、本実施形態に係る高炉操業特徴量演算部103で実施される高炉操業特徴量の演算処理を統括する制御部である。すなわち、演算制御部121は、後述するパラメータ設定部123、反応速度算出部125及び特徴量算出部127と連携しながら、これらの処理部における各種処理の実施状況を統括する。演算制御部121は、高炉操業情報取得部101から高炉操業情報が出力されると、まず、後述するパラメータ設定部123に対して、高炉操業特徴量の演算処理で用いられる各種パラメータの設定を要請する。パラメータ設定部123により各種パラメータの設定値が通知されると、演算制御部121は、設定されたパラメータにより規定される演算上の時間に関して、着目する化学反応式それぞれの反応速度の算出を後述する反応速度算出部125に要請する。反応速度算出部125により着目する化学反応式における反応速度が算出されると、演算制御部121は、着目している時間における高炉操業特徴量の算出を、後述する特徴量算出部127に要請する。ここで、本実施形態に係る高炉操業特徴量演算部103では、固体相、液体相、気体相の3つの相状態それぞれについて、高炉操業特徴量の演算を行う。以下で詳述するように、各相における高炉操業特徴量の演算を行うためには、時間に依存する項を含む偏微分方程式を考慮しなければならない。そこで、演算制御部121は、パラメータ設定部123により設定された時間刻みの設定値に基づいて、反応速度や特徴量の算出を行う時間を更新していく。以下では、図6を参照しながら、演算制御部121による時間更新処理について、具体的に説明する。図6は、本実施形態に係る演算制御部121における時間更新方法を説明するための説明図である。本実施形態に係る高炉操業特徴量演算処理では、固体相、液体相、気体相それぞれの高炉操業特徴量を算出する際に、各相における時間刻み換言すれば、時間を更新するまでの間隔が互いに独立に設定される。以下の説明では、固体相の時間刻み、液体相の時間刻み、気体相の時間刻みを、それぞれΔtS、ΔtL、ΔtGと表すこととし、ΔtL=2ΔtS、ΔtG=2ΔtL=4ΔtSの関係が成立している場合を示す一例である。なお、各相の時間刻みは、相内の現象が時間変化する速さ換言すれば、時定数に基づき、所望の演算精度が得られるように設定すればよい。いま、時間t=0において、固体相、液体相、気体相の演算ここでいう演算は、反応速度の算出処理及び高炉操業特徴量の算出処理の双方を含む。を開始したとする。かかる演算により、それぞれ固体相_1とラベリングされた状態、液体相_1とラベリングされた状態、気体相_1とラベリングされた状態に関する反応速度や特徴量が算出されることとなる。時間t=0における演算が終了すると、演算制御部121は、演算時間tを更新する処理を実施する。図32の例の場合、一番短い時間刻みはΔtSであるため、演算制御部121は、時間tを0ΔtS=ΔtSに更新して、後述する反応速度算出部125及び特徴量算出部127に処理開始を要請する。かかる場合において、液体相及び気体相については時間を更新するタイミングが到来していないため、反応速度算出部125及び特徴量算出部127は、固体相に関しては固体相_2とラベリングされた状態の演算のみを開始する。また、各演算処理において液体相や気体相の状態を参照する必要がある場合には、液体相_1とラベリングされた状態、及び、気体相_1とラベリングされた状態の値を用いて処理が実施されることとなる。固体相_2とラベリングされた状態の演算が終了すると、演算制御部121は、演算時間tをΔtSΔtS=2ΔtS=ΔtLに更新して、反応速度算出部125及び特徴量算出部127に処理開始を要請する。かかる場合において、気体相については時間を更新するタイミングが到来していないため、反応速度算出部125及び特徴量算出部127は、固体相に関しては固体相_3とラベリングされた状態、液体相に関しては液体相_2とラベリングされた状態の演算を開始するが、気体相に関する新たな演算は開始しない。また、各演算処理において気体相の状態を参照する必要がある場合には、気体相_1とラベリングされた状態の値を用いて処理が実施されることとなる。以下、同様の処理が継続して実施され、時間t=4ΔtS=2ΔtL=ΔtGとなった時点で、初めて気体相_2とラベリングされた状態の演算が開始されることとなる。このように、本実施形態に係る高炉操業特徴量の演算処理では、演算上の時間は固体相、液体相、気体相の3相に共通して流れている換言すれば、同一の時間軸を用いて各相の演算が定義されているものの、各相における新たな状態の演算を開始するタイミングは、個別の時間刻みに応じて規定されている。その結果、本実施形態に係る高炉操業特徴量の演算処理の流れでは、時間に関して多重のループが存在することとなる。このように、演算制御部121は、後述するパラメータ設定部123が設定した時間刻みに基づいて時間更新処理を行いながら、高炉操業特徴量の演算処理を統括する。また、演算制御部121は、パラメータ設定部123、反応速度算出部125及び特徴量算出部127から出力された各種データを、データ出力源以外の処理部に通達する仲介を行う。これにより、パラメータ設定部123、反応速度算出部125及び特徴量算出部127は、それぞれが出力した各種データを参照しつつ、個々の処理を実施することが可能となる。演算制御部121は、設定された演算期間シミュレーション時間内での演算処理が終了すると、得られた演算結果すなわち、高炉操業特徴量に関する情報を、操業制御部105に出力する。これにより、本実施形態に係る操業制御部105は、高炉操業特徴量を利用して、高炉の操業安定性の評価や、高炉操業の制御等を実施することが可能となる。また、演算制御部121は、得られた高炉操業特徴量に関するデータを、当該データを生成した日時等の時刻情報と関連付けて、履歴情報として記憶部111等に格納してもよい。パラメータ設定部123は、例えば、CPU、ROM、RAM等により実現される。パラメータ設定部123は、高炉操業特徴量演算部103で行われる高炉操業特徴量の演算処理で用いられる各種のパラメータを設定する。この際、パラメータ設定部123は、記憶部111等に格納されている各種のデータベースや文献値等に関する情報や、高炉操業状態制御装置10の使用者例えば高炉操業オペレータや操業管理者が、当該高炉操業状態制御装置10に付属するキーボード、タッチパネル、各種操作ボタン等の入力装置を介して入力した各種の設定情報や、各種センサ5や各種制御手段7等から取得した各種データや、高炉操業状態制御装置10に接続可能な外部装置から取得した情報等を適宜利用することができる。パラメータ設定部123が設定する主なパラメータとしては、例えば、着目する高炉内の成分とその種類の個数、着目する化学反応式、固体相、液体相、気体相の各相の時間刻み幅、演算単位となる高炉内の空間の区分であるメッシュの大きさの設定値メッシュの設定値等を挙げることができる。これらのパラメータのうち、各相の時間刻み幅、メッシュの設定値、着目する化学反応式及び成分について、以下で簡単に説明する。これらのパラメータは、例えばオペレータ等の入力に基づき設定される。本実施形態に係る高炉操業特徴量の演算処理では、後述するように、高炉の内部状態を、固体相液体相気体相の3つの相状態に分類して、各相状態について高炉操業特徴量を算出する。この際、本実施形態に係るパラメータ設定部123では、気体相の高炉操業特徴量の時間刻み幅を、液体相及び固体相の高炉操業特徴量の時間刻み幅よりも長く設定する。また、本実施形態に係る高炉操業特徴量の演算処理において、パラメータ設定部123は、図7に示したように高炉内部を複数のメッシュに区分する。本実施形態に係る高炉操業特徴量演算部103は、かかるメッシュを演算単位として、固体相液体相気体相のモデルいわゆる物理モデルを考慮する。なお、図7に示したメッシュの区分からも明らかなように、本実施形態に係る高炉操業特徴量の演算処理においては、高炉の高さ方向と高炉の半径方向の2つの方向を考慮した2次元物理モデルを例示するが、かかる2次元物理モデルに対して更に高炉の炉周方向の角度方向を追加した3次元メッシュを演算単位とする3次元物理モデルを考慮することも可能である。ここでメッシュの個数すなわち、高炉内部をどのように区切るかについては、高炉操業特徴量の演算にかけることができるパーソナルコンピュータの演算能力やメモリ容量等のリソース及び演算時間、並びに、演算対象とする高炉の大きさ等に応じて、操業制御の精度等を考慮して適宜設定すればよい。また、図7に示したように、各メッシュの形状は矩形とすることが好ましいが、かかる形状に限定されるわけではない。更に、高炉内部を矩形のメッシュで区切る場合、高炉の炉壁を含むメッシュでは、メッシュの形状と当該メッシュに属する高炉の形状とが表記上は一致しないこととなる。しかしながら、メッシュの面積又は体積を、当該メッシュに属する実際の面積又は体積と設定することで、見かけ上の形状の不一致を解消することが可能である。また、パラメータ設定部123は、高炉操業特徴量演算部103で実施される演算処理で着目する化学反応式を設定する。炉内で進行している化学反応を表現する化学反応式はいくつかの組み合わせを仮定することができるが、以下においては、その一つの組み合わせの例として、図8に示した10個の化学反応式に着目する。すなわち、本実施形態に係るパラメータ設定部123は、図8に示した10個の化学反応式に着目する。図8に示したように、着目する化学反応式は、原料である鉄鉱石がFeOまで還元される一連の反応反応式1,2,9,10、固体状態にあるFeOが溶解する反応反応式3、溶融FeOがFeに還元される直接還元反応反応式4、原料である鉄鉱石を還元する還元性ガスが生成される反応反応式5,6,7、水性ガスの変成反応反応式8に大別される。図8に例示したように、高炉内で生じている化学反応には、吸熱反応も存在すれば、発熱反応も存在する。このような化学反応のうち、吸熱反応が多く進行するようになれば、高炉内の温度炉内温度が低下することとなり、出銑量の低下や、製造される銑鉄の品質の低下等を招くこととなる。従って、高炉を安定的に操業させるためには、これら化学反応の熱バランスを制御することが重要となる。パラメータ設定部123は、着目すべき化学反応式を設定すると、次に、着目する化学成分を選択する。以下で説明する具体例では、図8に示した多数の化学成分のうち、Fe2O3S、Fe3O4S、FeOS、FeOL、FeL、CS、COGCO2G、H2G、H2OGについて着目する。なお、各化学式の末尾に記載されたSは、化学式で表される成分が固体状態Solidにあることを示しており、末尾に記載されたLは、化学式で表される成分が液体状態Liquidにあることを示しており、末尾に記載されたGは、化学式で表される成分が気体状態Gasにあることを示している。なお、上記化学反応式及び成分はあくまでも一例であって、図8に示した10種類以外の反応に着目してもよいし、図8に示した10種類の反応の中から更に選択すべき反応式を抽出してもよい。また、N2やSiなどの他の成分について着目してもよい。また、パラメータ設定部123は、以上説明したようなパラメータの他にも、高炉操業特徴量を演算する際に設定される各種のパラメータを設定することが可能である。以上、パラメータ設定部123が設定するパラメータの一例について、簡単に説明した。パラメータ設定部123は、このようなパラメータを設定すると、設定したパラメータを表す情報を、演算制御部121に出力する。また、パラメータ設定部123は、設定したパラメータを表す情報を、後述する反応速度算出部125及び特徴量算出部127に直接出力してもよい。また、パラメータ設定部123は、設定したパラメータを表す情報を、パラメータを設定した日時等の時刻情報と関連付けて、履歴情報として記憶部111等に格納してもよい。反応速度算出部125は、例えば、CPU、ROM、RAM等により実現される。反応速度算出部125は、パラメータ設定部123が設定したパラメータを利用して、着目している化学反応それぞれの反応速度を算出する。反応速度算出部125は、反応速度を算出する際に、記憶部111等に格納されている各種のデータベースや文献値等に関する情報を適宜利用することが可能である。また、反応速度算出部125は、各種ネットワークに接続された情報検索サーバ等で管理されている様々な情報を、反応速度の算出に利用することも可能である。本実施形態に係る反応速度算出部125は、原料であるFe2O3が還元性ガスであるH2又はCOによってFe2O3→Fe3O4→FeO→Feと還元される一連の反応速度を算出する際に、図9A~図9Cに例示した未反応核界面モデルを利用する。未反応核界面モデルの考え方や当該モデルにおける仮定事項、及び、未反応核界面モデルに基づく一連の還元反応速度式は、例えば、非特許文献2及び非特許文献3等に開示されている。以下では、まず、図9Aに示した未反応2界面モデルFe2O3→Fe3O4→FeOを一例として、未反応核界面モデルの考え方及び仮定事項、並びに、還元反応速度式について説明する。原料である鉄鉱石のガス還元反応は、気体と固体との反応である気固反応に分類され、固体粒子径はほぼ一定に保たれ、粒子内部に未反応核未反応の部分が存在し、未反応核の外側に還元反応生成物層が形成される反応であると考えられている。かかる還元反応を考慮するために、未反応核界面モデルというモデルが考えられている。この未反応核界面モデルは、粒子内の界面において反応が起こり、その反応界面が内部に向かって移動していくというモデルである。ここで、未反応固体物質の空隙率が非常に小さい場合には、反応ガスが粒子内部まで浸透して行かず、反応が外表面近傍でまず起こり、その反応が次第に粒子内部に向かって移動していくという、トポケミカルtopochemical反応となる。この際の反応素過程は、1粒子の周囲に存在するガス境膜を通って還元性ガスが粒子表面へ移動するガス境膜内拡散過程、2還元反応生成物層内に存在する気孔を介して粒子表面から反応界面へ還元性ガスが移動する気孔内拡散過程、3反応界面における化学反応過程、の3過程が直列に並んで進行すると考えることができる。ここで、鉄鉱石の還元反応のような気固反応では、固体粒子内における反応成分の濃度分布を考えた場合、濃度分布の時間的変化速度は小さく、比較的短い時間内では、あたかも定常的な濃度分布が成立しているとみなすことができる。このような擬定常状態近似が成立している場合には、上記3つの反応素過程の速度が等しいとして取り扱うことが可能となる。従って、本実施形態に係る反応速度算出部125では、このような擬定常状態近似が成立しているとして、鉄鉱石の還元に関する反応速度を算出する。鉄鉱石の還元に関する反応速度を算出する際に、反応速度算出部125は、まず、未反応核界面モデルの界面における還元反応速度を算出する。本実施形態に係る未反応核2界面モデルでは、図9Aに示したように、未反応核であるFe2O3層の周囲に第1の還元反応生成物層であるFe3O4層が存在し、このFe3O4層の周囲に第2の還元反応生成物層であるFeO層が存在し、FeO層の周囲にガス境膜が存在するモデルを仮定する。Fe2O3層、Fe3O4層及びFeO層は、図9Aに示したように同心球として存在するものと仮定し、FeO層の半径すなわち、鉄鉱石粒子の半径をr0と表し、Fe2O3層の半径をr1と表し、Fe3O4層とFeO層との界面の半径をr2と表すこととする。また、以下で説明する未反応核2界面モデルでは、各パラメータについて、以下のような表記を行うこととする。ここで、第1界面とは、反応式1又は反応式9で表される化学反応が進行している界面であり、第2界面とは、反応式2又は反応式10で表される化学反応が進行している界面である。反応式1及び反応式2に着目した場合の未反応核2界面モデル[H2]:還元性ガスH2のモル分率[H2O]:H2Oガスのモル分率kc,1:反応式1の還元反応の速度定数[cm/s]kc,2:反応式2の還元反応の速度定数[cm/s]ke,1:反応式1の還元反応の平衡定数ke,2:反応式2の還元反応の平衡定数反応式9及び反応式10に着目した場合の未反応核2界面モデル[CO]:還元性ガスCOのモル分率[CO2]:CO2ガスのモル分率kc,9:反応式9の還元反応の速度定数[cm/s]kc,10:反応式10の還元反応の速度定数[cm/s]ke,9:反応式9の還元反応の平衡定数ke,10:反応式10の還元反応の平衡定数共通するパラメータy1:第1界面に関する界面半径比y2:第2界面に関する界面半径比Yb:粒子周囲のガスのモル分率Y0:鉄鉱石粒子表面におけるガスのモル分率Y1:第1界面でのガスのモル分率Y2:第2界面でのガスのモル分率Ye,1:第1界面でのガスの平衡モル分率Ye,2:第2界面でのガスの平衡モル分率VC,1:第1還元反応生成物層の界面での還元反応速度[mol/s]VC,2:第2還元反応生成物層の界面での還元反応速度[mol/s]VD,1:第1還元反応生成物層での還元反応速度[mol/s]VD,2:第2還元反応生成物層での還元反応速度[mol/s]kf:ガス境膜内の物質移動係数[m/s]TG:還元性ガスの温度[K]PG:還元性ガスの圧力[Pa]RG:気体定数[J/molK]ここで、Rf[s/m3]は、ガス境膜内の物質移動抵抗であり、対応する反応速度をVfmol/sと表すこととする。RC,1[s/m3]は、第1還元反応生成物層の界面での化学反応抵抗であり、RC,2[s/m3]は、第2還元反応生成物層の界面での化学反応抵抗である。また、RD,1[s/m3]は、第1還元反応生成物層における気孔内ガスの有効拡散抵抗であり、RD,2[s/m3]は、第2還元反応生成物層における気孔内ガスの有効拡散抵抗である。ガス境膜内の物質移動抵抗Rfの抵抗をFと表すこととすると、抵抗Fは、ガス境膜内の物質移動係数kfに関連する値となる。また、第1還元反応生成物層の界面及び第2還元反応生成物層での界面での化学反応抵抗RC,1,RC,2に基づく抵抗をそれぞれ新たにA1、A2と表記すると、かかる化学反応抵抗A1、A2は、還元反応の速度定数kc,1,kc,2と、各反応式の化学反応の平衡定数K1,K2とに関連する値となる。更に、第1還元反応生成物層及び第2還元反応生成物層における気孔内ガスの有効拡散抵抗RD,1、RD,2に基づく有効拡散抵抗をそれぞれ新たにB1,B2と表記すると、かかる有効拡散抵抗B1,B2は、第1還元反応生成物層及び第2還元反応生成物層の気孔内ガスの有効拡散係数De,1,De,2に関連する値となる。本実施形態に係る反応速度算出部125は、演算パラメータr0、r1、r2、De,1、De,2、kf、PG、RG、TG、[H2]、[H2O]、εに基づいて、各抵抗における具体的な抵抗と、駆動力となる平衡モル分率と、を算出して、各抵抗における還元反応速度を算出する。ここで、反応速度算出部125は、上記演算パラメータのうち、高炉操業情報に含まれているものについては、高炉操業情報に記載されている値を、演算処理の初期値とする。また、反応速度算出部125は、上記演算パラメータのうち高炉操業情報に含まれていないものについては、各種のデータベースや文献値等に関する情報や、高炉操業状態制御装置10の使用者がキーボード、タッチパネル、各種操作ボタン等の入力装置を介して入力した各種の設定情報や、高炉操業状態制御装置10に接続可能な外部装置から取得した情報に基づいて決定した値を、演算処理の初期値とする。また、時間刻みの更新が行われた場合には、更新前の演算により算出された値が転用可能な演算パラメータについては、算出値を利用して時間刻みの更新後の演算パラメータが決定される。なお、反応速度算出部125において利用される具体的な演算式については、反応工学における公知の方法により定式化することが可能であるが、以下では、反応式1及び反応式2における反応速度を算出する際に用いられる式について、適宜例を挙げて示すものとする。反応速度算出部125は、設定された演算パラメータを利用して、速度定数kc,1,kc,2、及び、平衡定数ke,1,ke,2を算出する。速度定数kc,1,kc,2、及び、平衡定数ke,1,ke,2は、還元性ガスの温度TGをパラメータとして含む指数を用いて指数関数で表されることが各種の非特許文献例えば、非特許文献2,3,5等に開示されており、反応速度算出部125は、これらの文献等に記載されている式を用いて、速度定数kc,1,kc,2、及び、平衡定数ke,1,ke,2を算出可能である。また、反応速度算出部125は、設定された演算パラメータを利用して、以下の式101及び式102により、界面半径比y1,y2を算出する。反応速度算出部125は、これらの値を算出すると、演算パラメータ及び算出した値を利用して、以下の式103~式108により、総括物質移動係数W2を算出する。他方、反応速度算出部125は、演算パラメータ及び算出した値を利用して、反応の駆動力となる平衡モル分率Ye,1H,Ye,2Hを算出する。この平衡モル分率は、着目している反応に関与する還元性ガスの平衡モル分率が、平衡状態にある全ガス成分のモル分率の和に占める割合として定義することができる。反応速度算出部125は、かかる平衡モル分率を、着目している反応に関与するガス成分本例では、H2ガス及びH2Oガスのモル分率[H2],[H2O]と、平衡定数ke,1,ke,2を用いて算出する。反応速度算出部125は、以上説明したような値を算出すると、以下の式109及び式110を利用して、H2ガスを還元性ガスとして用いた場合における、第1界面での還元反応速度VC,1Hと、第2界面での還元反応速度VC,2Hとを算出する。反応速度算出部125は、各界面での還元反応速度を算出すると、以下の式111及び式112に基づいて、反応式1における反応速度R1及び反応式2における反応速度R2を算出する。なお、以下の式111におけるV1,H2は、VC,1Hに対応しており、以下の式112におけるV2,H2は、VC,2Hに対応している。また、COガスを還元性ガスとして利用する反応式9及び反応式10における反応速度についても、同様にして算出することが可能である。すなわち、反応速度算出部125は、反応式1及び反応式2における場合と同様にして、設定された演算パラメータを利用して、速度定数kc,9,kc,10、及び、平衡定数ke,9,ke,10を算出する。反応式9及び反応式10についても、速度定数kc,1,kc,2、及び、平衡定数ke,9,ke,10は、還元性ガスの温度TGをパラメータとして含む指数を用いて指数関数で表されることが各種の非特許文献例えば、非特許文献4,5等に開示されており、反応速度算出部125は、これらの文献等に記載されている式を用いて、速度定数kc,9,kc,10、及び、平衡定数ke,9,ke,10を算出可能である。かかる値を算出することで、反応速度算出部125は、反応式1及び反応式2の場合と同様にして、以下の式113~式116により、各抵抗を算出することができる。その結果、反応速度算出部125は、上記式107及び式108を利用して、総括物質移動係数W2を算出することができる。他方、反応速度算出部125は、反応式1及び反応式2における場合と同様にして、演算パラメータ及び算出した値を利用して、反応の駆動力となる平衡モル分率Ye,1C,Ye,2Cを算出する。すなわち、反応速度算出部125は、平衡モル分率Ye,1C,Ye,2Cを、着目している反応に関与するガス成分本例では、COガス及びCO2ガスのモル分率[CO],[CO2]と、平衡定数ke,9,ke,10を用いて算出する。反応速度算出部125は、以上説明したような値を算出すると、以下の式117及び式118を利用して、COガスを還元性ガスとして用いた場合における、第1界面での還元反応速度VC,1Cと、第2界面での還元反応速度VC,2Cとを算出する。反応速度算出部125は、各界面での還元反応速度を算出すると、以下の式119及び式120に基づいて、反応式9における反応速度R9及び反応式10における反応速度R10を算出する。なお、以下の式119におけるV1,COは、VC,1Cに対応しており、以下の式120におけるV2,COは、VC,2Cに対応している。以上のようにして、本実施形態に係る反応速度算出部125は、原料である鉄鉱石の還元反応の反応速度を算出することができる。また、図8に示した他の気固反応の反応速度についても、適切な界面数を考慮した未反応核界面モデルを考慮し、速度定数kcや平衡定数keに関する文献値等を適宜利用することで、同様にして反応速度を算出することが可能である。また、文献に記載されている反応速度の算出式そのものを、利用することも可能である。以上、図9Aを参照しながら、未反応核2界面モデルを利用した反応速度の算出方法の一例について、具体的に説明した。なお、本実施形態に係る反応速度算出部125は、先だって述べたように、上記未反応核2界面モデル以外のモデルを利用して反応速度を算出することもできる。以下では、図9Bを参照しながら、未反応核1界面モデルを用いた場合の反応速度の算出方法について、簡単に説明する。なお、以下で説明する未反応核1界面モデルは、図9Bに示したように、考慮する界面数を1としたモデルであり、以下の反応式11及び反応式12で表される化学反応Fe2O3の還元反応を考慮するものである。Fe2O33H2→2Fe3H2O反応式11Fe2O33CO→2Fe3CO反応式12かかる未反応核1界面モデルを利用する場合、反応速度算出部125は、以下の式121に基づいて、還元反応速度VC,1H,VC,1Cを算出することができる。なお、以下の式121では、VC,1H,VC,1CをまとめてVC,1として1つの式で表記している。ここで、上記式121において、パラメータA1は、以下の式122で表される化学反応抵抗であり、パラメータB1は、以下の式123で表される有効拡散抵抗であり、パラメータFは、以下の式124で表される物質移動抵抗である。また、反応速度算出部125は、算出したこれらの抵抗を利用して、以下の式125で表される総括物質移動係数W1を算出することができる。なお、以下の式122におけるパラメータK1は、上記反応式11又は反応式12の平衡定数である。反応速度算出部125は、上記式120に基づいて還元反応速度を算出すると、先だって説明した未反応核2界面モデルと同様にして、着目している化学反応反応式11、反応式12の反応速度を算出することができる。続いて、図9Cを参照しながら、未反応核3界面モデルを用いた場合の反応速度の算出方法について、簡単に説明する。なお、以下で説明する未反応核3界面モデルは、図9Cに示したように、考慮する界面数を3としたモデルであり、以下の反応式13~反応式18で表される一連の化学反応を考慮するものである。ここで、以下の反応式13~反応式15は、還元性ガスとしてH2ガスが関与する場合の化学反応式であり、以下の反応式16~反応式18は、還元性ガスとしてCOガスが関与する場合の化学反応式である。かかる未反応核3界面モデルを利用する場合、反応速度算出部125は、以下の式126~式128に基づいて、還元反応速度VC,1H,VC,2H,VC,3H,VC,1C,VC,2C,VC,3Cを算出することができる。なお、以下の式126~式128では、VC,nH,VC,nCをまとめてVC,nとして1つの式で表記している。ここで、上記式126は、反応式13及び反応式16に基づく第1界面における還元反応速度であり、上記式127は、反応式14及び反応式17に基づく第2界面における還元反応速度であり、上記式128は、反応式15及び反応式18に基づく第3界面における還元反応速度である。また、上記式126~式128において、パラメータAnは、以下の式129で表される化学反応抵抗であり、パラメータBnは、以下の式130で表される有効拡散抵抗であり、パラメータFは、以下の式131で表される物質移動抵抗である。なお、以下の式129におけるパラメータKnは、各界面n=1~3にて生じている化学反応の平衡定数である。また、以下の式129及び式130におけるパラメータRnは、各界面における還元率であり、以下の式132で表されるものである。なお、以下の式132においてR4=0である。反応速度算出部125は、算出したこれらの抵抗を利用して、以下の式133で表される総括物質移動係数W3を算出することができる。反応速度算出部125は、上記式126~式128に基づいて還元反応速度を算出すると、先だって説明した未反応核2界面モデルと同様にして、着目している化学反応の反応速度を算出することができる。以上、図9A~図9Cを参照しながら、未反応核界面モデルを利用した酸化鉄の還元反応の反応速度の算出処理について、詳細に説明した。なお、図8に示した反応式1,2,9,10以外の化学反応の反応速度は、各種非特許文献等に基づいて公知の方法により算出することが可能である。例えば、図8の反応式3に示した、固体相FeOの液体相FeOへの相変化反応の反応速度R3は、例えば、溶融温度における平衡条件式により与えることができる。また、図8の反応式4に示した溶融FeOのCによる直接還元反応の反応速度R4は、例えば、非特許文献6に基づいて算出することができる。また、図8の反応式5に示したレースウエイにおけるO2によるCの燃焼反応の反応速度R5は、本実施形態に係る高炉操業特徴量の演算処理におけるパラメータとして与えることができる。更に、図8の反応式6に示したCO2によるカーボンソリューション反応の反応速度R6は、例えば、非特許文献7及び非特許文献8に基づいて算出することができる。また、図8の反応式7に示したH2Oによるカーボンソリューション反応の反応速度R7は、例えば、非特許文献7及び非特許文献8に基づいて算出することができる。また、図8の反応式8に示した水性ガス変性反応シフト反応の反応速度R8は、反応式1の反応速度R1、反応式2の反応速度R2及び反応式7の反応速度R7を利用して、R8=1/3R1R2-R7として算出することが可能である。反応速度算出部125は、以上説明したようにして着目する化学反応式における反応速度を算出すると、算出した反応速度を表す情報を、演算制御部121に出力する。なお、反応速度算出部125は、算出した反応速度を表す情報を、特徴量算出部127に直接出力してもよい。また、反応速度算出部125は、算出した反応速度を表す情報を、反応速度を算出した日時等の時刻情報に関連付けて、履歴情報として記憶部111等に格納してもよい。特徴量算出部127は、例えば、CPU、ROM、RAM等により実現される。特徴量算出部127は、高炉操業情報、パラメータ設定部123により設定された各種パラメータ、及び、反応速度算出部125により算出された反応速度に基づいて、高炉の操業状態を表す特徴量である高炉操業特徴量を算出する。この際、特徴量算出部127は、パラメータ設定部123により設定されたメッシュを演算単位として、固体相液体相気体相の各相状態モデルを利用し、高炉操業特徴量の算出をオンラインで行う。ここで、固体相モデル液体相モデル気体相モデルの各物理モデルは、物質収支式、熱収支式、移動速度決定式の3種類の偏微分方程式から構成される。以下、各相の物理モデルについて、具体的な式を示しながら説明する。○固体相モデルまず、固体相モデルについて説明する。固体相モデルは、以下の式141で表される物質収支式と、式142で表される熱収支式と、式143及び式144で表される移動速度決定式からなる。ここで、下記式143及び式144で表される移動速度決定式は、いわゆるkinematicmodelにおける定式を利用している。ここで、上記式141において、方程式中の各パラメータの意味は、以下の通りである。αi:固体成分iの体積分率vS:固体相流れの速度ベクトル[m/s]αij:固体成分iの化学反応式jにおける係数mS,i:固体成分iの単位モルあたりの質量[kg/mol]di:固体成分iの密度[kg/m3]Rj:化学反応式jの反応速度t:時間[s]ε:空隙率また、上記式142において、方程式中の各パラメータの意味は、以下の通りである。なお、式141と共通するパラメータについては、記載を省略する。TS:固体相温度[K]TL:液体相温度[K]TG:気体相温度[K]ρS:固体相密度[kg/m3]=ΣdiαiCS:固体相比熱[J/kgK]ηS:固体相への熱分配率kS:固体相の熱拡散率[W/mKm]UGS:気体相から固体相への熱伝達率[W/m2Km2]ULS:液体相から固体相への熱伝達率[W/m2Km2]また、上記式143及び式144において、方程式中の各パラメータの意味は、以下の通りである。vS,x:固体相移動速度のx成分[m/s]vS,y:固体相移動速度のy成分[m/s]BS:比例定数α:固体相の体積分率の総和=Σαi○液体相モデル次に、液体相モデルについて説明する。液体相モデルは、以下の式145で表される物質収支式と、式146で表される熱収支式と、式147及び式148で表される移動速度決定式からなる。ここで、下記式147及び式148で表される移動速度決定式は、いわゆるkinematicmodelにおける定式を利用している。ここで、上記式145において、方程式中の各パラメータの意味は、以下の通りである。γi:液体成分iの密度[kg/m3]vL:液体相流れの速度ベクトル[m/s]bij:液体成分iの化学反応式jにおける係数mL,i:液体成分iの単位モルあたりの質量[kg/mol]di:固体成分iの密度[kg/m3]Rj:化学反応式jの反応速度t:時間[s]また、上記式146において、各パラメータの意味は、以下の通りである。TS:固体相温度[K]TL:液体相温度[K]TG:気体相温度[K]γL:液体相密度[kg/m3]=ΣγiCL:液体相比熱[J/kgK]ηL:液体相への熱分配率kL:液体相の熱拡散率[W/mKm]USL:固体相から液体相への熱伝達率[W/m2Km2]UGL:気体相から液体相への熱伝達率[W/m2Km2]また、上記式147及び式148において、方程式中の各パラメータの意味は、以下の通りである。vL,x:液体相移動速度のx成分[m/s]vL,y:液体相移動速度のy成分[m/s]BL:比例定数γ:液体相の体積分率の総和=Σγi○気体相モデル次に、気体相モデルについて説明する。気体相モデルは、以下の式149で表される物質収支式と、式150で表される熱収支式と、式151~式153で表される移動速度決定式からなる。ここで、下記式151~式153で表される移動速度決定式は、いわゆるErgun式を利用している。ここで、上記式149において、方程式中の各パラメータの意味は、以下の通りである。βi:気体成分iの密度[kg/m3]vG:気体相流れの速度ベクトル[m/s]kG:気体拡散係数[m2/s]aij:気体成分iの化学反応式jにおける係数mG,i:気体成分iの単位モルあたりの質量[kg/mol]di:固体成分iの密度[kg/m3]Rj:化学反応式jの反応速度t:時間[s]ε:空隙率また、上記式150において、方程式中の各パラメータの意味は、以下の通りである。TS:固体相温度[K]TL:液体相温度[K]TG:気体相温度[K]βG:気体相密度[kg/m3]=ΣβiCG:気体相比熱[J/kgK]ηG:気体相への熱分配率kG:気体相の熱拡散率[W/mKm]USG:固体相から気体相への熱伝達率[W/m2Km2]ULG:液体相から気体相への熱伝達率[W/m2Km2]また、上記式151及び式153において、方程式中の各パラメータの意味は、以下の通りである。P:気体圧力[Pa]vρ,G:気体相の質量速度ベクトル[kg/m2s]vρ,G,x:気体相の質量速度のx成分[kg/m2s]vρ,G,y:気体相の質量速度のy成分[kg/m2s]g:重力加速度[m/s2]μG:気体の粘度[kgm/s]φG:粒子の形状係数ρG:気体の密度[kg/m3]DP:粒子径m上記モデルを構成する方程式の形から明らかなように、各モデルを構成する方程式には、時間に依存する項時間tで変数を偏微分した項が含まれており、非線形の偏微分方程式になっていることがわかる。従って、本実施形態に係る特徴量算出部127は、各物理モデル固体相モデル液体相モデル気体相モデルにおいて、物質収支式、熱収支式、移動速度決定式の3種類の方程式を連立させて、高炉操業特徴量を算出する。特徴量算出部127が算出する高炉操業特徴量として、例えば以下のような特徴量を挙げることができる。固体相に関する高炉操業特徴量固体相速度、各固体成分の体積分率、固体相温度など液体相に関する高炉操業特徴量液体相速度、各液体成分の密度、液体相温度など気体相に関する高炉操業特徴量気体相速度、気体相圧力、各気体成分の密度、気体相温度など○各モデルの取り扱いここで、実際に高炉内で生じている現象について着目する。実際の高炉操業において、羽口から供給される熱風のガスによって生じるガス流れは、秒のオーダー数秒程度で炉頂部まで到達する。一方、銑鉄の原料となる鉄鉱石等の固体物は、炉内滞留時間が10時間程であることが知られている。このような動特性換言すれば、現象の時定数の違いを考慮すると、高炉全体の現象を考慮するにあたって支配的な時定数は、固体相の動特性で決まると考えることができる。そこで、本実施形態に係る特徴量算出部127は、固体相モデル及び液体相モデルについては、上記方程式をそのまま利用して、時間tに関する偏微分を取り扱う非定常モデルとして取り扱う一方で、動特性が秒オーダー程度と考えられる気体相モデルについては、以下のような取り扱いを行う。ここで、非定常モデルとは、当該モデルで記述する現象が時間変化する項や時間発展する項を含んでいる数式モデルをいう。すなわち、上述のような動特性に関する考察から、気体相の物質収支の動特性は、本実施形態に係る高炉操業特徴量演算部103において着目する化学反応式の反応速度が律速になると考えられる。そこで、本実施形態に係る特徴量算出部127は、式149に示した気体相モデルの物質収支式において、時間tの偏微分に関する項≒0とみなして、定常モデルとして物質収支式を取り扱う。同様に、気体相の熱容量は固体相の熱容量と比べて十分に小さいと考えることができるため、式150に示した気体相モデルの熱収支式において、時間tの偏微分に関する項≒0とみなして、定常モデルとして熱収支式を取り扱う。これにより、演算負荷の大きな気体相モデルを定常モデルとして取り扱うことが可能となり、演算負荷の低減、及び、更なる計算時間の高速化を図ることができる。ここで、定常モデルとは、当該モデルで記述する現象が時間変化する項や時間発展する項を含まない数式モデルをいう。○方程式の解法について本実施形態に係る特徴量算出部127が実施する上記式141~式153の解法は、特に限定されるものではないが、例えば、以下のような方法を用いることが可能である。例えば、式141~式153に示した連立方程式の解を求める際に、逆行列による代数的計算手法を採用するAX=bを満たす変数行列Xを、X=bA-1として算出することができる。また、気体相モデルにおける式151~式153移動速度決定式を用いて、気体相の流速分布を算出する際には、いわゆる流れ関数法を採用することで、式151~式153で定義されるErgun式を、連立方程式として取り扱い、かつ、収束計算を行うことで方程式を満たす流れ関数を算出することが可能となる。その結果、気体相の流速分布を算出する際に、流れのポテンシャルΦを容易に算出することが可能となる。また、本実施形態に係る特徴量算出部127では、気体相における高炉操業特徴量を算出する際に、各種の物性値データベースから取得した物性値気体相物性値を利用する。この際、特徴量算出部127は、着目した気体相物性値の温度特性を、物性値データベースから取得した値を非線形関数で同定してもよい。本実施形態に係る特徴量算出部127は、以上のようにして、固体相液体相気体相それぞれにおける高炉操業特徴量を算出すると、算出した高炉操業特徴量を表す情報を、演算制御部121に出力する。また、特徴量算出部127は、算出した高炉操業特徴量に、当該高炉操業特徴量を算出した日時等に関する時刻情報を関連付けて、記憶部111に履歴情報として格納してもよい。以上、本実施形態に係る高炉操業特徴量演算部の構成について、詳細に説明した。[高炉操業特徴量の演算方法について]続いて、図10~図13を参照しながら、本実施形態に係る高炉操業特徴量演算部103により実施される高炉操業特徴量の演算方法の流れについて、詳細に説明する。図10~図13は、本実施形態に係る高炉操業特徴量の演算方法の流れの一例を示した流れ図である。○全体的な流れについてまず、図10を参照しながら、本実施形態に係る高炉操業特徴量の演算方法の全体的な流れの一例について説明する。なお、先だって説明したように、本実施形態に係る高炉操業特徴量演算部103では、固体相に関する演算の時間刻み、液体相に関する演算の時間刻み、及び、気体相に関する演算の時間刻みを互いに独立に設定し、各相の時間刻みを多重にループさせることで、計算時間の高速化を図っている。以下で説明する流れ図では、液体相に関する演算の時間刻みと、気体相に関する演算の時間刻みとを同一の値に設定し、固体相に関する演算の時間刻みと、液体相及び気体相に関する演算の時間刻みとが二重にループしている場合の流れを図示している。高炉操業特徴量演算部103の演算制御部121は、高炉操業情報取得部101から出力された高炉操業情報を取得すると、まず、パラメータ設定部123に各種パラメータの設定を要請する。要請を受けたパラメータ設定部123は、高炉内部のメッシュへの区分や、演算に用いられる状態変数等、処理に利用する各種の状態変数やパラメータについて、値を初期化するステップS101。その後、パラメータ設定部123は、固体相液体相気体相それぞれの時間刻みと、演算シミュレーションの終了時間を設定するステップS103。本説明の場合では、パラメータ設定部123は、液体相及び気体相に関する時間刻みΔt1、固体相に関する時間刻みΔt2、及び、演算の終了時間t_endをそれぞれ設定する。この際、時間刻みΔt1及びΔt2の間には、Δt1>Δt2の関係が成立している。続いて、パラメータ設定部123は、設定した時間刻みΔt1及びΔt2を利用して、時間更新の二重ループ回数nを設定するステップS105。具体的には、パラメータ設定部123は、intΔt1/Δt2を算出し、得られた整数値を二重ループ回数nの値として設定する。その後、パラメータ設定部123は、高炉操業特徴量の演算処理における時間tを初期化して、t=0とするステップS107。パラメータ設定部123は、これらの処理が終了すると、設定したパラメータを表す情報を、演算制御部121に出力する。演算制御部121は、パラメータ設定部123から出力された時間に関するパラメータ時間t、時間刻みΔt1及びΔt2、二重ループ回数、終了時間に基づいて演算上の時間更新を管理しつつ、パラメータ設定部123、反応速度算出部125及び特徴量算出部127における処理の制御を実施する。まず、演算制御部121は、反応速度算出部125に対して、着目している化学反応式の反応速度の算出を要請する。反応速度算出部125は、演算制御部121から指定された演算時間t=0と、パラメータ設定部123により設定された各種のパラメータや演算初期値に基づいて、着目している化学反応式の反応速度を算出するために利用される各種物理量と反応速度とを算出するステップS109。反応速度算出部125は、算出した演算時刻tにおける反応速度を表す情報を、演算制御部121に出力する。次に、演算制御部121は、高炉操業特徴量を算出する際に利用する二重ループカウンタのカウント数を表すパラメータiを初期化i=1と設定するステップS111。その後、演算制御部121は、特徴量算出部127に対して、パラメータ設定部123により設定された各種のパラメータ等と、反応速度算出部125により算出された演算時間tにおける反応速度とを通知して、高炉操業特徴量の算出を要請する。特徴量算出部127は、演算制御部121からの要請を受けて、まず、固体相モデルを利用して固体相に関する演算処理を実施するステップS113。特徴量算出部127は、固体相に関する演算処理が終了すると、固体相に関する演算結果を、演算制御部121に出力する。演算制御部121は、特徴量算出部127から固体相に関する演算結果が出力されると、二重ループカウンタのカウント値iが二重ループ回数n以下であるか否かすなわち、i≦nであるか否かを判断するステップS115。演算制御部121は、i≦nが成立している場合には、二重ループカウンタのカウント値iを1増加させた後ステップS117、演算時間tを更新して、tΔt2とするステップS119。その後、演算制御部121は、新たな演算時間t=tΔt2を特徴量算出部127に通知して、固体相の演算処理を特徴量算出部127に要請する。すると、特徴量算出部127は、再びステップS113に戻って、演算時間t=tΔt2における固体相の演算処理を開始する。他方、演算制御部121は、i≦nが成立していない場合すなわち、i>nであった場合には、特徴量算出部127にi≦nが成立しなくなった演算時間tを通知して、液体相に関する演算処理を開始するように要請する。特徴量算出部127は、通知された演算時間tと、パラメータ設定部123により設定された各種のパラメータ等と、反応速度算出部125により算出された演算時間tにおける反応速度とを利用して、液体相モデルに基づいて、液体相に関する高炉操業特徴量の演算処理を実施するステップS121。その後、特徴量算出部127は、液体相に関する高炉操業特徴量の演算結果を演算制御部121に出力するとともに、通知された演算時間tと、パラメータ設定部123により設定された各種のパラメータ等と、反応速度算出部125により算出された演算時間tにおける反応速度とを利用して、気体相モデルに基づいて、気体相に関する高炉操業特徴量の演算処理を実施するステップS123。特徴量算出部127は、気体相に関する高炉操業特徴量の演算処理が終了すると、得られた気体相に関する高炉操業特徴量の演算結果を、演算制御部121に出力する。演算制御部121は、特徴量算出部127から気体相に関する高炉操業特徴量の演算結果が出力されると、それまでに得られた固体相液体相気体相のそれぞれに関する高炉操業特徴量を、グラフとしてプロットしたり、高炉操業特徴量演算部103の外部例えば、操業制御部105や表示制御部107に出力したりするステップS125。その後、演算制御部121は、演算時間tが終了時間になったか否かすなわち、t≦t_endが成立するか否かを判定するステップS127。t≦t_endが成立している場合には、演算制御部121は、その時点での演算時間tを反応速度算出部125に通知して、演算時間tにおける反応速度の算出を要請する。反応速度算出部125は、ステップS109に戻って、演算制御部121から通知された演算時間tにおける反応速度を算出する。他方、t≦t_endが成立していない場合すなわち、t>t_endが成立している場合には、所定の演算時間が経過したと判断して、高炉操業特徴量の算出を終了する。○固体相の演算処理について次に、図11を参照しながら、図10におけるステップS113で実施される固体相の演算処理の流れについて、簡単に説明する。特徴量算出部127は、演算制御部121から固体相の演算処理の実施を要請されると、まず、固体相モデルの物質収支式例えば、上記式141と、移動速度決定式例えば、上記式143及び式144とを主に利用して、固体相速度式141における変数vSを算出するステップS131。続いて、特徴量算出部127は、固体相モデルの物質収支式例えば、上記式141と、移動速度決定式例えば、上記式143及び式144とを主に利用して、各固体成分の体積分率式141等における変数αiを算出するステップS133。その後、特徴量算出部127は、固体相モデルの熱収支式例えば、上記式142を主に利用して、固体相温度式142における変数TSを算出するステップS135。このような処理を実施することで、特徴量算出部127は、固体相の高炉操業特徴量を算出することができる。○液体相の演算処理について次に、図12を参照しながら、図10におけるステップS121で実施される液体相の演算処理の流れについて、簡単に説明する。特徴量算出部127は、演算制御部121から液体相の演算処理の実施を要請されると、まず、液体相モデルの物質収支式例えば、上記式145と、移動速度決定式例えば、上記式147及び式148とを主に利用して、液体相速度式145における変数vLを算出するステップS141。続いて、特徴量算出部127は、液体相モデルの物質収支式例えば、上記式145と、移動速度決定式例えば、上記式147及び式148とを主に利用して、各液体成分の密度式145等における変数γiを算出するステップS143。その後、特徴量算出部127は、液体相モデルの熱収支式例えば、上記式146を主に利用して、液体相温度式146における変数TLを算出するステップS145。このような処理を実施することで、特徴量算出部127は、液体相の高炉操業特徴量を算出することができる。○気体相の演算処理について次に、図13を参照しながら、図10におけるステップS123で実施される気体相の演算処理の流れについて、簡単に説明する。特徴量算出部127は、演算制御部121から気体相の演算処理の実施を要請されると、まず、収束計算の計算回数カウント用のパラメータjを初期化して、j=1に設定するステップS151。その後、特徴量算出部127は、後述する収束計算の収束判定パラメータである許容誤差EPSを、例えば10-6に設定するステップS153。次に、特徴量算出部127は、流れ関数φの初期化を行い、φ0=0と設定するステップS155。続いて、特徴量算出部127は、気体相モデルの物質収支式例えば、上記式149や移動速度決定式例えば、上記式151~式153等を主に利用して、気体相流れの速度ベクトルを流れ関数で表すように変形したうえで、流れ関数φjを算出するステップS157。その後、特徴量算出部127は、流れ関数φj-1の値と、算出した流れ関数φjの値とを利用して、以下の式161で表される流れ関数の差分Δφjを算出するステップS159。式161から明らかなように、流れ関数の差分Δφjは、流れ関数φjとφj-1とのユークリッドノルムとして定義される値である。続いて、特徴量算出部127は、算出した流れ関数の差分Δφjが許容誤差EPS未満となったか否かΔφj<EPSを判定するステップS161。算出した流れ関数の差分ΔφjがEPS以上であった場合には、特徴量算出部127は、流れ関数φjの算出結果は収束していないと判断し、パラメータjの値を1増加させたうえでステップS163、ステップS157の処理を継続する。他方、算出した流れ関数の差分ΔφjがEPS未満であった場合には、特徴量算出部127は、算出した流れ関数の算出結果は収束していると判断し、得られた流れ関数φを用いて、気体相速度式149における変数vGを算出するステップS165。続いて、特徴量算出部127は、気体相モデルの移動速度決定式であるErgun式上記式151~式153を利用して、気体相圧力式151における変数Pを算出するステップS167。次に、特徴量算出部127は、気体相モデルの物質収支式例えば、上記式149を主に利用して、各気体成分の密度式149等における変数βiを算出するステップS169。その後、特徴量算出部127は、気体相モデルの熱収支式例えば、上記式150を主に利用して、気体相温度式150における変数TGを算出するステップS171。このような処理を実施することで、特徴量算出部127は、気体相の高炉操業特徴量を算出することができる。以上、図10~図13を参照しながら、本実施形態に係る高炉操業特徴量演算部103により実施される高炉操業特徴量の演算方法の流れの一例について、詳細に説明した。<高炉操業の安定性評価について>次に、図14~図18を参照しながら、本実施形態に係る操業制御部105の構成と、これにより実施される高炉3の操業状態の安定性評価処理と高炉3の制御方法の決定について説明する。[高炉操業状態制御システムの構成概念]操業制御部105に関する説明をするに先立ち、本実施形態に係る高炉操業状態制御システムの構成概念を、図14及び図15に基づいて説明する。高炉3は、上述したように、伝熱や流体、化学反応や相変化を伴う大規模複雑な動的システムである。このような動的システムの制御は、従来、モデルの計算負荷やシステム実装上の課題から線形定常モデルに基づき行われていた。この場合、動的システムに対する過去の入力やシステムの内部状態等は考慮されずにシステム出力が決定されるため、操業目標や操業動作点の変更や外乱に対する高炉3の動的挙動は計算することができなかった。すなわち、線形定常モデルに基づき制御した場合には、システム全体の制御性を動的安定性の観点から捉えることができない。一方で、大規模複雑な動的システムの制御のアプローチとして、ローカルなアクション計測、予測、制御で制御対象のグローバルな機能を実現するという提案がなされている上記非特許文献1。この提案によれば、図14に示すように、動的システムを制御する制御系は、制御対象を記述するための機能物理ネットワークの階層構造と、制御対象に所望の機能を実現するための機能動的情報ネットワークとして捉えられる。そして、動的システムの局所的な計測を、個々の各種センサを用いて物理ネットワークにて行い、局所的な計測による計測値に基づいて大局的なシステム全体に対して安定性の評価を動的情報ネットワークにて行うという、ネットワークの相互連携によって動的システムを制御する。本実施形態では、物理ネットワークで制御対象である高炉3を記述し、動的情報ネットワークで高炉3の操業状態の安定性を評価して安定するように制御する。物理ネットワーク、すなわち制御対象である動的システムは、当該制御対象の操業データや制御結果を動的情報ネットワークに出力し、動的情報ネットワークより入力された操作情報に基づいて制御対象を操業する。ここで、動的システムである高炉3は、上述したように、複数のサブシステムが相互に作用し合う相互作用系として捉えることができるとともに、異なるスケール例えば、高分解能な秒単位~低分解能な週単位での進行、高分解能な化学反応上の分子単位~低分解能のメートル単位での進行などで複数の非定常現象が起きているマルチスケールシステムでもある。換言すれば、図14に示すように、動的システムは、当該システムを構成する異なる分解能のサブシステムが相互に作用し合っているものといえる。動的情報ネットワークでは、制御対象から入力された操業データや制御結果である高炉操業情報に基づいて高炉3の操業状態を表す高炉操業特徴量を演算し、当該特徴量及び高炉操業情報に基づき高炉3という動的システム全体の安定性を定量的に評価する。そして、動的情報ネットワークは、動的システム全体の安定性評価の結果より、システム全体が安定して操業するための操作量を算出し、制御対象に出力する。すなわち、動的情報ネットワークは、制御対象を制御する動的制御システムとして考えることができる。本実施形態に係る動的情報ネットワークは、高炉3の操業状態をシミュレートする非定常物理モデルシミュレータと、高炉3の操業安定性を大局的に評価する操業安定性評価モデルとからなる。非定常物理モデルシミュレータでは、物理ネットワークから入力された高炉操業情報に基づいて、高炉3を非線形非定常分布定数モデルとして定式化し、高炉3の操業状態をシミュレートする。その結果、高炉3の操業状態を特徴づける高炉操業特徴量を取得し、操業安定性評価モデルへ出力する。操業安定性評価モデルは、高炉3の操業安定性を評価するためのモデルを同定し、高炉3を安定して操業させるために必要な各種制御手段7の操作量を決定する。操業安定性評価モデルは、高炉3内で進行している反応間での熱バランスやマスバランスといったダイナミクス動的視点に基づき同定することができる。なお、操業安定性評価モデルの具体的な説明は後述する。ここで、操業安定性評価モデルにおける動的システムの安定性判別処理の概要を図15に示す。図15では、高炉3内で進行している発熱反応と吸熱反応との熱バランスに起因する操業不調に関し、システムの安定特性を安定性特性曲線により表している。なお、安定性特性曲線は、正確には動的システムの多次元平面特性を示す。図15の安定性特性曲線は、左側から右側へ、上に凸、下に凸、上に凸の形状となっており、右側に向かうほど発熱反応が大きくなり、左側に向かうほど吸熱反応が大きくなるものとする。安定性特性曲線において、上に凸となっている部分の2つの頂点の間の領域下に凸となっている部分を含む。が発熱反応と吸熱反応との熱バランスのとれた操業安定領域である。また、上に凸となっている部分の頂点より下に凸となっている部分と反対側の領域が、発熱反応と吸熱反応との熱バランスが崩れる操業不安定領域である。このような安定性特性曲線を用いることにより、その曲線上に置かれたボールすなわち、動作点の動きによってシステムの安定不安定を考えることができる。例えば、図15に示すように、システムが動作点P1で操業されているとき、動作点P1は左右どちら側に僅かに動かされても安定性特性曲線の窪み部分に位置するように動く。すなわち、動作点P1は操業安定領域から外れないことから、高炉3内で進行している発熱反応と吸熱反応との熱バランスがとれており、システム高炉3は安定して操業されていると考えることができる。一方、システムが動作点P2で操業されているとき、動作点P2は右側に動かされると安定性特性曲線を下り落ち、再び同じ位置に戻ることはできない。すなわち、動作点P2の位置から熱バランスが右側に傾くと発熱反応が大きくなり、また、不可逆性が強いために、熱バランスが崩れてしまった状態となる。この場合、高炉3が熱を持ち過ぎてしまい、高炉3の操業が不安定操業不調となってしまう。システムが動作点P3で操業されているときも同等であり、動作点P3は左側に動かされると安定性特性曲線を下り落ち、再び同じ位置には戻らない。これは、動作点P3の位置から熱バランスが左側に傾くと吸熱反応が大きくなり、また、不可逆性が強いために、熱バランスが崩れてしまった状態となる。この場合、高炉3が冷え過ぎてしまい、高炉3の操業が不安定操業不調となってしまう。図15では、高炉3内で進行している発熱反応と吸熱反応との熱バランスに起因する操業不調に関して安定性特性曲線を示したが、熱バランスだけでなくマスバランス等のダイナミクスについても同等に考えることができる。本実施形態の操業制御部105は、このような動的視点に基づき、高炉3という動的システム全体の安定性を判断している。実際には、後述するように、操業制御部105の安定性評価部151の演算によって取得される操業安定性評価情報を用いて、図16に示した概念で高炉3の操業状態の安定性が判別される。なお、安定性特性曲線の操作は、操業条件の操作に相当する。したがって、操業条件を操作すると、安定性特性曲線が変化し、高炉3が操業安定となる操業安定領域も変化することになる。図14に戻り、動的情報ネットワークは、操業安定性評価モデルより現在のシステムの操業状態の安定性を判別すると、判別結果より不安定あるいは不安定となる方向へ動いていると認識された場合には、システムが安定状態となるように制御する制御方法を決定する。制御方法は、非定常物理モデルシミュレータシミュレータ、操業安定性評価モデル、および制御理論に基づき逐次操作量を算出し、高炉3が操業安定領域で操業する操作量を特定することにより決定される。このように、動的情報ネットワークでは、制御対象の操業状態をシミュレーションし、シミュレーション結果から、システム全体のマクロなレベルでの安定性を評価する。動的情報ネットワークでのシステムの安定性評価結果より、システムの動作点における操業安定性を判定することができる。そして、動的情報ネットワークは、システムが安定して操業できる操業安定領域内でシステムが動作するための操作量を算出し、操作情報として制御対象に出力する。以上、高炉操業状態制御システムの構成概念について説明した。高炉操業状態制御システムは、実際に制御対象を操業する階層の物理ネットワークと、システム全体が安定して操業するための計測、予測、制御を行う階層の動的情報ネットワークとの相互連携によって、大きく展開していく枠組みを構築している。そして、高炉3というシステム全体の安定性を評価することにより操作量を操作し、所望の状態を実現する。[操業制御部の構成について]次に、図16に基づいて、本実施形態に係る操業制御部105の構成について説明する。本実施形態に係る操業制御部105は、図16に示すように、安定性評価部151と、制御方法決定部153とを備える。安定性評価部151は、例えば、CPU、ROM、RAM等により実現される。安定性評価部151は、高炉操業特徴量演算部103の演算制御部121から入力された演算結果すなわち、高炉操業特徴量に関する情報に基づいて、高炉3の操業安定性を評価する。安定性評価部151による高炉3の操業安定性の評価は、マスバランスやエネルギバランスの動特性により記述された操業安定性評価モデル非線形システムに対する制御理論的視点でのシステム的な安定性評価条件に基づいて行われる。これにより、高炉3の操業安定性を定量的な指標操業安定性評価情報を用いて評価することができる。ここで、高炉3のダイナミクスを表現する状態空間表現モデルは、下記式200及び式201で表される。ここで、xtは状態量ベクトル、utは入力ベクトル、ytは出力ベクトルである。状態量ベクトルxtは、各種センサ5から入力される情報、例えば、炉内の熱エネルギや炉内の質量の時間変化量、炉内の固体相の充填レベル高さストックレベル等の変数からなる。また、入力ベクトルutは、高炉操業状態制御装置10から各種制御手段7に出力する情報、高炉3に対する操作量である、例えば、炉頂からの固体相流入質量流量や羽口からの気体相流入質量流量、羽口からの気体相流入温度、炉頂からの固体相流入質量流量における鉄鉱石とコークスの比率、炉体冷却水量等の変数からなる。そして、出力ベクトルytは、各種センサ5から入力される情報であって高炉3の操業状態を監視している量、例えば、出銑口からの液体相流出質量流量、出銑口からの液体相流出温度、炉頂からの気体相流出質量流量、炉頂からの気体相流出温度等の変数からなる。なお、状態量ベクトルxt、入力ベクトルut、及び出力ベクトルytの各変数はかかる例に限定されず、上記以外の値とすることもできる。式200及び式201のAt、Bt、Ct、Dtは行列関数である。行列関数は、高炉3の操業動作点において、高炉操業特徴量演算部103と実操業データとから同定される。同定されたAt、Bt、Ct、Dtは、高炉3の大局的な非定常システム特性を表している。そこで、安定性評価部151は、行列関数At、Bt、Ct、Dtを操業安定性評価モデルとして用いることにより、ある操業動作点における高炉3の操業安定性を行列関数At、Bt、Ct、Dtより制御理論に基づき逐次評価する。制御理論に基づく安定性評価条件としては、例えば、a可安定性可制御性に基づく条件可制御性行列のランク、固有値評価bナイキストの安定条件周波数特性に基づく安定特性評価cリアプノフの安定条件大局漸近安定特性評価内部安定性等がある。以下、これらの安定性評価条件について説明する。これらは、例えば、非特許文献9等でまとめられる如く広く知られており、本発明の実施形態において利用することが可能である。a可安定性可制御性に基づく条件上記式200及び式201の状態空間表現モデルで表されるシステムにおいて、入力ベクトルutを特定の時間関数を選択したとき、状態量ベクトルxtを任意の値にもっていくことができるならば、そのシステムは可制御である。システムが可制御であるためには、可制御性行列Ucのランクがnとなることが必要十分条件となる。なお、可制御性行列は、下記式202で表される。すなわち、式202の可制御性行列Ucのランクがnであれば、システムは可制御であることがわかる。したがって、安定性評価部151は、可制御性行列Ucのランクがnであるか否かを判別することにより、高炉3を制御可能か否かを判断することができる。なお、システムにおいて、ある有限な時刻sがあり、その間0≦t≦sの出力ベクトルytと入力ベクトルutとを測定することで、初期状態量x0を唯一に決定できるとき、そのシステムは可観測である。A、Bが可制御であるm入力の下記式203で表される制御対象において、全状態変数x1t、、xntが直接測定可能とし、制御入力を下記式204としたとする。式204を状態フィードバック制御といい、Fをフィードバック係数行列という。式203に式204を代入することより式205を導出することができ、式205の解は式206となる。ここで、x0は初期外乱を表し、t=0以前の区間で存在した外乱によって原点からずれたt=0での初期ベクトルを表す。したがって、Fを適切に選択し行列A-BFを安定行列にすることができれば、すべてのx0≠0に対して、xt→0t→∞とすることができ、状態変数を漸近的に原点に帰すことができる。換言すると、行列A-BFの固有値の値や配置を評価することにより、システムの安定性を評価することができる。bナイキストの安定条件ナイキストの安定判別法は、フィードバック系の内部安定性を判別する手法であり、開ループ伝達関数の周波数応答に基づき、ナイキスト線図を用いて図的に安定性を判別するものである。ナイキスト線図は、複素平面上にωを変化させたときの開ループ伝達関数の絶対値と位相角との関係を表したベクトル軌跡ナイキスト軌跡である。ナイキスト軌跡が点-1,0の周りを反時計方向にまわる回数が、開ループ伝達関数の不安定極の個数と等しいとき、制御系は安定であると判定することができる。ナイキストの安定判別法は、複雑な計算の必要がなく、高次の系やむだ時間系にも容易に適用できる。また、実データに基づいて判定することができ、直感的にもわかりやすい。cリアプノフの安定条件非線形システムである下記式207について、fx=0を満たす平衡点がx=0にあるものとする。リアプノフの安定性の定義は、x0が原点0すなわち、平衡点から少しずれたとき、xtが依然として原点近傍に留まり得るか否かを述べており、図16に示した安定性特性曲線のイメージで考えることができる。ここで、式207で表されるシステムのリアプノフ関数Vxは、下記式208で表される。リアプノフ関数Vxは、ベクトルxtについてのスカラー関数であり、かつ正定関数である。また、式208は準負定関数であり、∂V/∂xは連続であるとする。このとき、原点近傍のある範囲内でリアプノフ関数が存在すれば、原点は安定である。さらに式208が負定関数ならば、原点は漸近安定であることが知られている。特に、大域的な漸近安定性については、1xの全域でのリアプノフ関数Vxが存在する2∥x∥→∞のとき、Vx→∞となる3Vdotx≡0の解xtが式207の原点以外の開と恒等的に一致しないの3つの条件が成立するならば、原点は大域的に漸近安定であることが知られている。なお、Vdotxは、式208に示すリアプノフ関数Vxの時間微分を示す。このように、原点近傍のある範囲内でリアプノフ関数が存在するか否かを調べることにより、システムの安定性を判別することができる。上記の手法以外にも、例えば、スモールゲイン定理受動定理の安定性条件を用いてシステムの安定性を評価、判別することも可能である。また、安定性評価部151は、高炉3の非定常物理モデルシミュレータが算出する高炉操業特徴量に対して、部分空間同定法SubspacebasedStateSpaceIdentification:以後、4SID法と称する。を用いて、高炉3のダイナミクスを表現する上記式200及び式201の状態空間表現モデルを導出することも可能である。4SID法については、非特許文献10及び非特許文献11に、また、高炉の非定常物理シミュレータ算出結果に対する4SID法の適用方法の一例が非特許文献12に例示されており、これらの結果やその応用を本発明の実施形態に適用することが可能である。基本的な考え方は、まず、作成した高炉3の非定常物理モデルシミュレータの境界条件に外乱を加えたものを、高炉操業特徴量を算出するための高炉非定常物理モデルへの入力情報とする。次いで、入力情報に基づいて高炉3の非定常物理モデルシミュレータによる非定常計算を実行し、取得された高炉操業特徴量から外乱の影響を求める。これを高炉システムの出力情報とする。そして、安定性評価部151は、高炉操業情報、高炉操業特徴量、入力情報および出力情報に基づいて、4SID法を用いて、入力情報および出力情報のダイナミクスを近似する状態空間表現モデルを導出する。この状態空間表現モデルを操業安定性評価モデルとして利用することができる。このように同定された状態空間表現モデルの評価は、例えば、次式209の適合率FIT[]の値で判断することができる。式209では、一例として、入力として羽口における気体相温度を設定したときの、出力として設定した炉内の任意のメッシュにおける固体相温度TSに対する適合率を表している。また、Nは入出力情報のデータ長、kは離散化時間である。上述のとおり、高炉3は時間的にも空間的にも多分解能の複数のサブシステムが相互連携した動的なシステムと考えられる。したがって、4SID法の適用において、各メッシュにおける分解能を変化させ個々の分解能に対応する同定モデルを導出することができる。具体的には、1時間分解能のみを変化させた場合、2空間分解能のみを変化させた場合、3時間分解能と空間分解能の両方を同時に変化させた場合の3パターンが手法として提案できる。1~3の各パターンにおける各同定モデルの同定結果の一例を、固体相温度の適合率FIT[]の炉内2次元分布コンタ図として、図17A~図17Cに示す。図17A~図17Cの各図は、高炉3を12×18のメッシュに区分したときの各メッシュにおける式209により算出される同定モデルの適合率を表している。適合率が高いほど、高い数値を示す。図17Aは、時間分解能のみを変化させたパターン1の同定モデルの同定結果を示し、図17Bは、空間分解能のみを変化させたパターン2の同定モデルの同定結果を示す。そして、図17Cは、時間分解能を粗く設定した状態図17Aの状態bにおいて、空間分解能を変化させたパターン3の同定モデルの同定結果を示す。図17A~図17Cでは、分解能を変化させた3つの状態a、b、cを示しており、状態a、b、cの順に分解能が粗くなっている。なお、図17Aおよび図17Bの状態aは同一の状態を表している。図17Aおよび図17Bの状態aに示すように、羽口から離れるほど同定モデルの適合率は低くなっている。これは、羽口から離れたメッシュの出力は遅い時定数を持ち、また空間的にも高い周波数を持たないことによる。このため、図17Aおよび図17Bに示すように、羽口から離れたメッシュ近傍では時空間分解能を粗くした方が同定モデルの適合率が上がることがわかる。したがって、時間的にも空間的にも多分解能の複数のサブシステムが相互連携した動的なシステムである高炉のダイナミクスを、比較的少ないシステムパラメータで精度よく再現するには、時空間分解能の異なる同定モデル例えば、4SID法による状態空間表現モデルを組み合わせることが有効である。これは、図17Cに示すように、時間分解能を粗くした図17Aの状態bから空間分解能を粗くしていくと同定モデルの適合率が著しく向上していることからも示されている。例えば、組み合わせ法による伝達関数は、式211のように表現することができる。ここで、tは時間分解能、sは空間分解能、iは炉高方向メッシュ番号、jは炉半径方向メッシュ番号、wは重み係数、fは同定モデル、zは離散系伝達関数における演算子である。なお、ここでは、一例として、式211の空間メッシュ毎の伝達関数fi,jを例示したが、本発明はかかる例に限定されない。実高炉での計測データや高炉非定常物理モデルシミュレータによるシミュレーション結果等を通じて得られる、高炉という動的システムが有する時間分解能、空間分解能に関する知見に基づき、近傍に連続する複数のメッシュを一つの領域として取り扱い、当該領域毎の伝達関数fを定義することによって、更にシステムパラメータを少なくした形で精度のよい伝達関数を取得できることは言うまでもない。以上のような制御理論に基づき、安定性評価部151は、ある操業動作点における高炉3の操業安定性を逐次評価する。なお、本発明はかかる例に限定されず、上記以外の制御理論に基づき安定性を評価することもできる。安定性評価部151は、可安定性可制御性の状態を示す値や、ナイキスト線図のようなグラフから取得できる定量的な安定性判別情報を取得すると、制御方法決定部153へ出力する。制御方法決定部153は、例えば、CPU、ROM、RAM等により実現される。制御方法決定部153は、安定性評価部151から入力された安定性判別情報に基づいて、動的システムである高炉3の安定状態を判断し、高炉3の制御方法を決定する。すなわち、安定性判別情報より高炉3の操業状態が安定であると判断したとき、制御方法決定部153は、現状と同等、高炉3の安定操業を維持するように各種制御手段7の操作量を決定する。一方、安定性判別情報より高炉3の操業状態が不安定であると判断したとき、制御方法決定部153は、高炉3の操業状態が不安定から安定となるように、各種制御手段7の操作量を決定する。これにより、高炉3をシステム全体として安定して操業させることができる。制御方法決定部153により高炉3の制御方法が決定されると、各種制御手段7へ操作量が出力される。[高炉操業の安定性評価]図4に示す操業制御部105による高炉操業の安定性評価の処理ステップS15の詳細を図18に示す。図18は、本実施形態に係る操業制御部105による高炉操業の安定性評価の処理を示すフローチャートである。まず、操業制御部105は、高炉操業情報及び高炉操業特徴量を取得するステップS201。高炉操業情報は、例えば記憶部11のデータベースから取得することができ、高炉操業特徴量は、高炉操業特徴量演算部103から取得することができる。次いで、操業制御部105の安定性評価部151は、高炉3の状態空間表現モデルの同定を行うステップS203。高炉3の状態空間表現モデルは、ステップS201で取得した高炉操業情報及び高炉操業特徴量より、上記式200及び式201で表される高炉3の大局的な非定常システム特性を表す行列関数At、Bt、Ct、Dtを同定することにより取得される。ここで、行列関数At、Bt、Ct、Dtは、高炉3という動的システムの操業安定性を評価する操業安定性評価モデルに相当する。そして、安定性評価部151は、ステップS203にて同定された行列関数At、Bt、Ct、Dtを用いて、ある時刻における操業動作点における高炉操業安定性を評価するステップS205。高炉操業安定性の評価は、上述した制御理論に基づき、定量的に行われる。安定性評価部151は、評価結果を制御方法決定部153へ出力する。このように、操業制御部105は、高炉操業特徴量演算部103により演算された高炉操業特徴量や、高炉操業情報等を利用して、高炉操業の安定性評価を実施することができる。そして、操業制御部105の制御方法決定部153は、高炉操業特徴量、高炉操業情報、安定性評価の結果等を利用して、高炉3の制御方法を決定するステップS17。ステップS15及びS17の処理は、周期的に行われ、時々刻々の操業動作点において高炉3の操業安定性が評価されるとともに、高炉3の操業を安定させるための制御出力が決定される。以上、操業制御部105における高炉3の操業状態の安定性判別と制御出力の決定方法について説明した。高炉操業状態制御装置のハードウェア構成について次に、図19を参照しながら、本発明の実施形態に係る高炉操業状態制御装置10のハードウェア構成について、詳細に説明する。図19は、本発明の実施形態に係る高炉操業状態制御装置10のハードウェア構成を説明するためのブロック図である。高炉操業状態制御装置10は、主に、CPU901と、ROM903と、RAM905と、を備える。また、高炉操業状態制御装置10は、更に、バス907と、入力装置909と、出力装置911と、ストレージ装置913と、ドライブ915と、接続ポート917と、通信装置919とを備える。CPU901は、演算処理装置および制御装置として機能し、ROM903、RAM905、ストレージ装置913、またはリムーバブル記録媒体921に記録された各種プログラムに従って、高炉操業状態制御装置10内の動作全般またはその一部を制御する。ROM903は、CPU901が使用するプログラムや演算パラメータ等を記憶する。RAM905は、CPU901が使用するプログラムや、プログラムの実行において適宜変化するパラメータ等を一次記憶する。これらはCPUバス等の内部バスにより構成されるバス907により相互に接続されている。バス907は、ブリッジを介して、PCIPeripheralComponentInterconnect/Interfaceバスなどの外部バスに接続されている。入力装置909は、例えば、マウス、キーボード、タッチパネル、ボタン、スイッチおよびレバーなどユーザが操作する操作手段である。また、入力装置909は、例えば、赤外線やその他の電波を利用したリモートコントロール手段いわゆる、リモコンであってもよいし、高炉操業状態制御装置10の操作に対応したPDA等の外部接続機器923であってもよい。さらに、入力装置909は、例えば、上記の操作手段を用いてユーザにより入力された情報に基づいて入力信号を生成し、CPU901に出力する入力制御回路などから構成されている。高炉操業状態制御装置10のユーザは、この入力装置909を操作することにより、高炉操業状態制御装置10に対して各種のデータを入力したり処理動作を指示したりすることができる。出力装置911は、取得した情報をユーザに対して視覚的または聴覚的に通知することが可能な装置で構成される。このような装置として、CRTディスプレイ装置、液晶ディスプレイ装置、プラズマディスプレイ装置、ELディスプレイ装置およびランプなどの表示装置や、スピーカおよびヘッドホンなどの音声出力装置や、プリンタ装置、携帯電話、ファクシミリなどがある。出力装置911は、例えば、高炉操業状態制御装置10が行った各種処理により得られた結果を出力する。具体的には、表示装置は、高炉操業状態制御装置10が行った各種処理により得られた結果を、テキストまたはイメージで表示する。他方、音声出力装置は、再生された音声データや音響データ等からなるオーディオ信号をアナログ信号に変換して出力する。ストレージ装置913は、高炉操業状態制御装置10の記憶部の一例として構成されたデータ格納用の装置である。ストレージ装置913は、例えば、HDDHardDiskDrive等の磁気記憶部デバイス、半導体記憶デバイス、光記憶デバイス、または光磁気記憶デバイス等により構成される。このストレージ装置913は、CPU901が実行するプログラムや各種データ、および外部から取得した各種のデータなどを格納する。ドライブ915は、記録媒体用リーダライタであり、高炉操業状態制御装置10に内蔵、あるいは外付けされる。ドライブ915は、装着されている磁気ディスク、光ディスク、光磁気ディスク、または半導体メモリ等のリムーバブル記録媒体921に記録されている情報を読み出して、RAM905に出力する。また、ドライブ915は、装着されている磁気ディスク、光ディスク、光磁気ディスク、または半導体メモリ等のリムーバブル記録媒体921に記録を書き込むことも可能である。リムーバブル記録媒体921は、例えば、CDメディア、DVDメディア、Blu-rayメディア等である。また、リムーバブル記録媒体921は、コンパクトフラッシュ登録商標CompactFlash:CF、フラッシュメモリ、または、SDメモリカードSecureDigitalmemorycard等であってもよい。また、リムーバブル記録媒体921は、例えば、非接触型ICチップを搭載したICカードIntegratedCircuitcardまたは電子機器等であってもよい。接続ポート917は、機器を高炉操業状態制御装置10に直接接続するためのポートである。接続ポート917の一例として、USBUniversalSerialBusポート、IEEE1394ポート、SCSISmallComputerSystemInterfaceポート、RS-232Cポート等がある。この接続ポート917に外部接続機器923を接続することで、高炉操業状態制御装置10は、外部接続機器923から直接各種のデータを取得したり、外部接続機器923に各種のデータを提供したりする。通信装置919は、例えば、通信網925に接続するための通信デバイス等で構成された通信インターフェースである。通信装置919は、例えば、有線または無線LANLocalAreaNetwork、Bluetooth登録商標、またはWUSBWirelessUSB用の通信カード等である。また、通信装置919は、光通信用のルータ、ADSLAsymmetricDigitalSubscriberLine用のルータ、または、各種通信用のモデム等であってもよい。この通信装置919は、例えば、インターネットや他の通信機器との間で、例えばTCP/IP等の所定のプロトコルに則して信号等を送受信することができる。また、通信装置919に接続される通信網925は、有線または無線によって接続されたネットワーク等により構成され、例えば、インターネット、家庭内LAN、赤外線通信、ラジオ波通信または衛星通信等であってもよい。以上、本発明の実施形態に係る高炉操業状態制御装置10の機能を実現可能なハードウェア構成の一例を示した。上記の各構成要素は、汎用的な部材を用いて構成されていてもよいし、各構成要素の機能に特化したハードウェアにより構成されていてもよい。従って、本実施形態を実施する時々の技術レベルに応じて、適宜、利用するハードウェア構成を変更することが可能である。第1実施例<高炉操業特徴量の演算に関する実施例>以下では、まず、図20~図25を参照しながら、本発明の第1の実施形態に係る高炉操業特徴量の演算に関する実施例について説明する。本実施例では、本発明の第1の実施形態に係る高炉操業特徴量の演算方法を用いることで得られる演算の高速化効果について、具体的に説明する。まず、本実施例で用いた演算モデル等について説明する。本実施例では、第1の実施形態において説明した高炉2次元非定常物理モデルを計算モデルとした。この計算モデルについて、図8に示した12個の化学反応式について、固体相液体相気体相の3相を考慮した高炉操業特徴量の演算を実施した。なお、かかる演算処理において、着目した成分は、Fe2O3S、Fe3O4S、FeOS、FeOL、FeL、CS、COG、CO2G、H2G、H2OGの10成分である。また、高炉操業特徴量を算出する際に用いた非定常方程式として、第1の実施形態において式141~式153で示した方程式すなわち、物質収支式、熱収支式、移動速度決定式を利用した。なお、高炉操業特徴量の演算にあたって、非定常演算の条件を以下のように設定した。1初期状態t=0以前において:3相の温度分布373[K]、3相の流速分布0[m/s]2t=0[sec]において:a:羽口熱風条件をステップ入力:熱風温度1500[K]、熱風流量6300[Nm3/min]、送風圧力4[atm]b:炉頂装入物鉄鉱石コークス条件をステップ入力:装入物温度373[K]、装入流量375[ton/hour]以上のような条件設定のもと、本発明の第1の実施形態に係る高炉操業特徴量演算方法の高速化効果を検証するために、演算リソース等の条件が同一のコンピュータを利用して、以下のようにして演算を行った。○従来法による演算固体相液体相気体相の各非定常方程式を差分法及び一次オイラー法を用いて、t=0~200[hour]の範囲で求解した。ここで、演算の経過時間更新のための時間刻みΔtは、3相共通に設定した。具体的には、Δt=0.0001[sec]である。○第1の実施形態に係る演算方法時間分解能の差に着目し、固体相液体相については非定常方程式を利用し、気体相については定常方程式を利用したすなわち、時間tによる偏微分の項を0とみなした。。ここで、上記方程式の解は、t=0~200[hour]の範囲で算出した。また、具体的な演算は、図11~図14に示したような流れで実施するものとし、Δt1すなわち、液体相及び気体相の時間刻み=100[sec]、Δt2すなわち、固体相の時間刻み=1[sec]とした。○演算に要した時間に関して上述のような条件で高炉操業特徴量の演算を実施し、演算に要した時間を測定した。その上で、第1の実施形態に係る演算方法の演算時間/従来法による演算時間を算出することで、本演算方法による高速化の度合いを算出した。このような演算の結果、第1の実施形態に係る演算方法は、従来法による演算結果に比べて、73.8倍の高速化を実現することができた。○高炉操業特徴量の演算結果に関して上述のような条件のもとで、第1の実施形態に係る演算方法を実施することで算出された高炉操業特徴量の一例を、図20~図25に示した。ここで、図20は、炉内に存在するFe2O3の存在比すなわち、体積分率の演算結果を示したものである。また、図21は、炉内の気体相温度[K]の演算結果を示したものであり、図22は、炉内の固体相温度[K]の演算結果を示したものである。更に、図23は、炉内のCO濃度すなわち、体積分率の演算結果を示したものであり、図24は、炉内のH2濃度すなわち、体積分率の演算結果を示したものである。また、図25は、炉内圧力気体相圧力[atm]の演算結果を示したものである。なお、図20~図25において、縦軸は、高炉の高さを表しており、横軸は、高炉の半径を表している。図20~図25より明らかなように、本発明の第1の実施形態に係る演算方法を用いることで、200時間後までの長期的な高炉操業特徴量を、高速に算出することが可能であることがわかる。以上、添付図面を参照しながら本発明の好適な実施形態について詳細に説明したが、本発明はかかる例に限定されない。本発明の属する技術の分野における通常の知識を有する者であれば、特許請求の範囲に記載された技術的思想の範疇内において、各種の変更例または修正例に想到し得ることは明らかであり、これらについても、当然に本発明の技術的範囲に属するものと了解される。"