JP2023025646A - 制御装置、方法及びプログラム - Google Patents

Abstract

【課題】プロセスの制御において、外乱発生時に生じ得る操業変動を速やかに、精度良く抑制する。【解決手段】高炉プロセス１の自動制御系は、制御器２と、外乱オブザーバ３とを備え、高炉プロセス１を制御対象として、操作量ｕに対する状態変数ｘの時間変化を表す状態空間モデルを利用したモデル予測制御を実行する。高炉プロセス１は、操作量ｕを入力として、制御量ｚ及び制御量ｚを含む観測量ｙを出力し、操作量ｕに入力外乱ｄ1が重畳し、観測量ｙに出力外乱ｄ2が重畳するものとする。外乱オブザーバ３は、操作量ｕ及び出力外乱ｄ2が重畳した観測量ｙを入力として、状態変数の推定値ｘ＾、入力外乱の推定値ｄ1＾、出力外乱の推定値ｄ2＾を出力する。制御器２は、目標値ｒ及び外乱オブザーバ３の出力を入力として、モデル予測制御を実行することにより、操作量ｕを出力する。【選択図】図１

Description

本発明は、プロセスを制御対象とする制御装置、方法及びプログラムに関する。

高炉操業においては、安価な鉄鉱石の利用や低還元材比条件下での生産が近年指向されている。しかしながら、安価で低品質の原料等を用いることによって高炉への装入原料の品質に変動が生じると、高炉操業に対する外乱となり得る。その結果として、生産量の変動や炉熱の低下といった、高炉操業の変動（以下では単に「操業変動」と称する場合がある。）を引き起こし、炉況の不安定につながるおそれがある。安定的な高炉操業の実現のためには、外乱発生時に生じ得る操業変動を速やかに抑制することが求められる。

高炉操業に関する技術として、特許文献１及び特許文献２には、非定常状態における高炉内の状態を計算可能な物理モデルを用いて高炉における溶銑温度を予測する手法が開示されている。

特開２０１８－２４９３５号公報 特開２０１９－１９３８５号公報

足立：ＭＡＴＬＡＢによる制御のための上級システム同定、東京電機大学出版局、２００４ 萩原：ディジタル制御入門、コロナ社、ｐ．４６、１９９９ B.Moore：Principal component analysis in linear systems: Controllability, observability, and model reduction , IEEE Transactions on Automatic Control ( Volume: 26, Issue: 1, February 1981) 安藤他：数値解析手法による制御系設計、コロナ社、１９８６ J. Maciejowski：モデル予測制御、東京電機大学出版局、２００５ K. Murota: Matrices and Matroids for Systems Analysis, Springer-Verlag, Berlin,2000.

特許文献１では、外乱によるプロセス変化を補償するため、過去所定の時間に渡り遡及し、モデル内温度分布の初期値とモデルパラメータの調整により外乱に対処するとしている。しかしながら、外乱発生に伴うプロセス変化をモデル初期値とモデルパラメータの調整だけで補償するため、初期値とモデルパラメータの変化で調整できない外乱（調整対象としないモデルパラメータの影響やモデルで考慮されない反応の影響、測定誤差等）が発生した場合に、却って制御性能が悪化するおそれがある。また、非定常計算が可能な高炉数学モデルの利用を前提としているが、通常、数値計算モデルは極めて計算負荷が高く、短い制御周期での計算が難しいため、この制御周期がネックとなり、外乱発生時に生じる操業変動を速やかに抑制できなくなるおそれがある。

特許文献２では、非定常モデルで計算された溶銑温度予測値を補正するために、溶銑温度以外のソルロスカーボン量やＲＡＲ（還元材比）等の操業指標の時間変化率を説明変数とする回帰式を利用している。しかしながら、この手法では、高炉プロセスの操業変動を引き起こす外乱が、高炉プロセスの操作端側に重畳する外乱入力か、高炉プロセスの出力端側に重畳する外乱出力であるかの区別がされておらず、これらの外乱発生時に生じる操業変動を精度良く抑制できなくなるおそれがある。

本発明は上記のような点に鑑みてなされたものであり、プロセスの制御において、外乱発生時に生じ得る操業変動を速やかに、精度良く抑制することを目的とする。

本発明の制御装置は、プロセスを制御対象として、状態変数の時間変化を表す状態空間モデルを利用したモデル予測制御を実行する制御装置であって、
前記プロセスに操作量を入力したときの前記プロセスから出力される制御量を含む観測量を取得し、
前記操作量に入力外乱が重畳可能な設計とし、前記観測量に出力外乱が重畳可能な設計とするものとして、
前記入力外乱及び出力外乱のうち少なくともいずれかの外乱が重畳した前記操作量及び観測量を入力として、前記状態変数の推定値と、前記入力外乱の推定値及び前記出力外乱の推定値のうち少なくともいずれかとを出力する推定手段を有することを特徴とする。
本発明の制御方法は、プロセスを制御対象として、状態変数の時間変化を表す状態空間モデルを利用したモデル予測制御を実行する制御方法であって、
前記プロセスに操作量を入力したときの前記プロセスから出力される制御量を含む観測量を取得し、
前記操作量に入力外乱が重畳可能な設計とし、前記観測量に出力外乱が重畳可能な設計とするものとして、
前記入力外乱及び出力外乱のうち少なくともいずれかの外乱が重畳した前記操作量及び観測量を入力として、前記状態変数の推定値と、前記入力外乱の推定値及び前記出力外乱の推定値のうち少なくともいずれかとを出力するステップを有することを特徴とする。
本発明のプログラムは、プロセスを制御対象として、状態変数の時間変化を表す状態空間モデルを利用したモデル予測制御を実行するためのプログラムであって、
前記プロセスに操作量を入力したときの前記プロセスから出力される制御量を含む観測量を取得し、
前記操作量に入力外乱が重畳可能な設計とし、前記観測量に出力外乱が重畳可能な設計とするものとして、
前記入力外乱及び出力外乱のうち少なくともいずれかの外乱が重畳した前記操作量及び観測量を入力として、前記状態変数の推定値と、前記入力外乱の推定値及び前記出力外乱の推定値のうち少なくともいずれかとを出力する推定手段としてコンピュータを機能させる。

本発明によれば、入力外乱と出力外乱を区別して重畳可能な設計としているので、プロセスの制御において、外乱発生時に生じ得る操業変動を速やかに、精度良く抑制できる。

高炉プロセスの概要を説明するための図である。 高炉プロセスと、操作量、制御量及び観測量との関係を示す図である。 第１の実施形態に係る高炉プロセスの自動制御系の構成を示す図である。 時刻と制御量及び操作量との関係の例を示す図である。 情報処理装置の構成例を示す図である。 実施例１の結果を示す図である。 第２の実施形態に係る高炉プロセスの自動制御系の構成を示す図である。 実施例２の結果を示す図である。 第３の実施形態に係る高炉プロセスの自動制御系における更新装置を説明するための図である。 基準値からの観測量の摂動量と、基準値からのパラメータの摂動量との関係の例を示す特性図である。 計算ステップと、間接還元反応のパラメータ及びその観測量であるガス利用率との関係の例を示す特性図である。 送風量の変更操作に対する出銑量及び溶銑温度の時間応答の計算結果を示す特性図である。

以下、添付図面を参照して、本発明の好適な実施形態について説明する。
＜第１の実施形態＞
図１は、高炉プロセスの概要を説明するための図である。高炉プロセスでは、図１に示すように、高炉１０１の炉上部から焼結鉱とコークスを交互に層を造るように装入し、炉下部にある送風羽口から熱風、還元材として微粉炭、等が吹込まれる。この熱風で微粉炭及びコークスがガス化し、一酸化炭素、水素等の高温の還元ガスが炉内を吹き昇り、焼結鉱を溶かしながら酸素を奪い、溶けた鉄分はコークスの炭素と接触して還元され、炭素５％弱を含む溶銑となって炉底の湯溜まり部に溜まる。この溶銑は、炉底横に設けられた出銑口から取り出され、次の製鋼プロセスへと運ばれる。なお、送風羽口（或いはその他の羽口）からは、コークス炉ガス（ＣＯＧ）や天然ガス（ＮＧ）等のガスが吹き込まれることがある。

また、高炉操業においては、操業管理指標が設定され、当該操業管理指標の時間推移の監視、及びそれに対応した操作量での高炉プロセスに対する操業操作の実行によって、炉況を制御している。そのため、適切な操業管理指標の設定は、高炉の安定的な操業のために極めて重要である。操業管理指標の詳細は後述するが、設定された目標値に保つべき指標（本実施形態では、制御量）や、主として炉況変化を把握するために利用する指標（本実施形態では、観測量）等がある。

図２は、高炉プロセス１と、操作量ｕ、制御量ｚ及び観測量ｙとの関係を示す図である。図２に示すように、高炉プロセス１は、操作量ｕを入力し、制御量ｚ及び観測量ｙを出力とする。また、高炉プロセス１の状態をｘとする。
操作量ｕは、送風量、微粉炭量に関する値（微粉炭吹込み量ＰＣＩ又は微粉炭比ＰＣＲ）、酸素富化に関する値（酸素富化量又は酸素富化率）、送風温度、湿分量（送風湿分）、及び液化天然ガス量に関する値（ＬＮＧ吹込み量）、コークス比、Ｏ／Ｃ（鉱石とコークスの重量比）のいずれか一つを含む（複数選択可）。
制御量ｚは、出銑量及び炉熱指標を含む。炉熱指標は、溶銑温度、溶銑中Ｓｉ、カーボンソリューションロス量、直接還元反応量、羽口上熱バランス、及び融着帯レベルのいずれか一つを含む（複数選択可）。また、制御量ｚは、後述する観測量ｙの線形結合により表現されたものであってよい。
観測量ｙは、制御量ｚ（出銑量及び炉熱指標）を含み、かつ、その他の炉熱指標（制御量ｚとして扱われなかった炉熱指標）、ガス組成（例えばＣＯ濃度、ＣＯ₂濃度、Ｎ₂濃度、Ｈ₂濃度）、ガス利用率（例えばＣＯ利用率、Ｈ₂利用率）、Ｔｆ（羽口先理論燃焼温度）、炉頂温度、熱流比、及び送風圧力のいずれかを含む（複数選択可）。

なお、上述した炉熱指標のうち、羽口上熱バランスとは、送風羽口（以下、羽口とも称する。）よりも上部の熱バランスを指す指標である。具体的には、羽口から吹込むガスが有する送風顕熱（単に「送風顕熱」ともいう。）、レースウエイで燃焼するカーボンの燃焼熱（単に「カーボン燃焼熱」ともいう。）、及び一酸化炭素による間接還元反応による反応熱（単に「間接還元熱」ともいう。）を含む羽口上の（すなわち羽口レベル及びこれより上部で生じる）インプット熱量と、羽口から吹込むガスが有する送風湿分の分解熱（単に「送風湿分分解熱」ともいう。）、微粉炭の分解熱（単に「微粉炭分解熱」ともいう。）、吸熱反応であるソリューションロス反応による反応熱（単に「ソリューションロス反応熱」ともいう。）及び直接還元反応による反応熱（単に「直接還元熱」ともいう。）、並びに水素還元反応（すなわち水素による鉄鉱石の還元反応）による反応熱（単に「水素還元熱」ともいう。）を含む羽口上のアウトプット熱量と、の差分として定義される指標である。

［状態空間モデル］
本実施形態では、高炉プロセス１の制御に、高炉操業を表現する（１)式の離散時間線形状態空間モデルを利用する。（１）式の離散時間線形状態空間モデルには、現時刻をｋで表し、現時刻において操作量ｕ（ｋ）を加えたときに高炉プロセス１の状態（状態変数）がｘ（ｋ）からｘ（ｋ＋１）に変化する離散時間の変化を表す式が含まれる。また、（１）式の離散時間線形状態空間モデルには、ここで推定した状態変数ｘと制御量ｚ及び観測量ｙとの関係式を含む。なお、説明の簡単のため、（１）式では直達項を割愛するが、必要に応じ考慮してもよい。

ここで、（１）式の離散時間線形状態空間モデルの作成方法について述べる。
各操作量ｕと各観測量ｙとの間で成立する（２）式の伝達関数行列Ｇ（ｓ）が得られているものとする。ｓはラプラス演算子であり、Ｇ（ｓ）はｐ行ｍ列の伝達関数行列である。ｐは観測量ｙの次元、ｍは操作量ｕの次元である。

伝達関数行列Ｇ（ｓ）は、例えば、高炉操業中に操業水準（動作点）の変更中期間を含むデータを利用し、非特許文献１等に記載されている開ループ同定手法又は閉ループ同定手法により作成することができる。操業中のデータを利用しない方法として、非定常計算が可能な高炉数学モデル（高炉非定常モデル）が利用できる場合には、高炉非定常モデルの動作点をいくつか変更し得られた時系列データを用いて、開ループ同定により伝達関数行列を求めてもよい。いずれの場合も、定常操業時における動作点の摂動に関し、通常、操作量の加算性が概ね成立するので、線形モデルを仮定してよい。したがって、例えば、操作量毎に伝達関数行列を求めればよい。

また、このようにして得られた伝達関数行列Ｇ（ｓ）に対して、いわゆる状態実現によって可制御正準形等の適当な連続時間の状態空間表現に変換した後、ゼロ次ホールド（制御周期間で入力値を一定とする）等の手法で時間離散化することにより、（１）式の離散時間線形状態空間モデルを作成すること、すなわち、行列Ａ，Ｂ，Ｃ，Ｃ_zを求めることができる。また、後述する外乱オブザーバと呼ばれるタイプの状態推定オブザーバとモデル予測制御を適用するために、この離散時間線形状態空間モデルが可到達性及び可観測性を満たすことを前提とする。なお、離散時間線形状態空間モデルにおいては、非特許文献２等に記載されているように、連続時間状態空間モデルの可制御性に対応する概念としては、離散時間系においては可到達性がより正確な概念となるが、便宜上、この性質を可制御性とも呼ぶこととする。なお、このとき得られた離散時間線形状態空間モデルが可制御性及び可観測性を有さない場合であっても、例えば、非特許文献３に記載されている平衡実現を適用すれば、可制御可観測なモデルを構成することができる。平衡実現では、非特許文献４の５章により計算の手続きが記載されている。所与の状態空間モデルに対して定義される可制御グラミアン及び可観測グラミアンを利用し、可制御可観測性の程度を示すＨａｎｋｅｌ特異値を算出する。Ｈａｎｋｅｌ特異値が所定の閾値を下回るものを打ち切ることで、可制御可観測性を満たす状態空間モデルを得ることができる。
観測量の数が多い場合、伝達関数行列Ｇ（ｓ）を単純に時間離散化して状態空間モデルを作成すれば、そのままでは可制御可観測な状態空間モデルが得られない場合がある。上述した平衡実現により再構成された状態空間モデルは、このような観測量の次元が多い場合に、特に顕著な効果を有する。

なお、離散時間線形状態空間モデルの作成方法は、上述したような伝達関数行列の作成を経由する方法に限定されるものではない。非特許文献１に記載されている部分空間同定法（ＭＯＥＳＰ法やＮ４ＳＩＤ法）等のシステム同定手法を用いて、伝達関数行列の作成を経由せず、直接、入力と出力の時系列データから離散時間線形状態空間モデルを求めてもよい。
また、離散時間線形状態空間モデルを利用しているが、連続時間状態空間モデルが得られている場合には、それを時間離散化して利用してもよい。また、非線形の状態空間モデルが得られている場合であっても、例えば、動作点周りで線形近似して線形状態空間モデルを作成してもよいし、Ｋｏｏｐｍａｎモード分解又はＣａｒｌｅｍａｎ線形化のような手法で、線形化されたモデルを利用してもよい。

［外乱オブザーバを含む高炉プロセスの自動制御系］
図３は、第１の実施形態に係る高炉プロセス１の自動制御系の構成を示す図である。高炉プロセス１の自動制御系は、制御器２と、上述の方法で求めた離散時間線形状態空間モデル（Ａ，Ｂ，Ｃ，Ｃ_z行列）を利用し構成した外乱オブザーバ３とを備え、高炉プロセス１を制御対象として、制御器２ではモデル予測制御を実行する。高炉プロセス１において、操作量ｕに入力外乱ｄ₁が重畳可能な設計とし、観測量ｙに出力外乱ｄ₂が重畳可能な設計とする。
制御器２は、目標値ｒ及び外乱オブザーバ３の出力である推定値ｘ＾，ｄ₁＾，ｄ₂＾（例えばｘ＾の表記は、ｘの上に＾が付されているものとする）を入力として、モデル予測制御を実行することにより、操作量ｕを出力する。
外乱オブザーバ３は、本発明でいう推定手段として機能し、高炉プロセス１に操作量ｕを入力したときの高炉プロセス１から出力される制御量ｚを含む観測量ｙを取得し、入力外乱ｄ₁及び出力外乱ｄ₂のうち少なくともいずれかの外乱が重畳した操作量ｕ及び観測量ｙを入力として、状態変数の推定値ｘ＾と、入力外乱の推定値ｄ₁＾及び出力外乱の推定値ｄ₂＾のうち少なくともいずれかとを出力する。図３に示す例では、外乱オブザーバ３は、操作量ｕ及び出力外乱ｄ₂が重畳した観測量ｙを入力として、状態変数の推定値ｘ＾、入力外乱の推定値ｄ₁＾、出力外乱の推定値ｄ₂＾を出力する。

図３では、高炉プロセス１に重畳する外乱の位置を模式的に示す。入力外乱ｄ₁は、高炉プロセス１の入側（操作端）に重畳する外乱又は計測誤差であり、出力外乱ｄ₂は、高炉プロセス１の出側（出力端）に重畳する外乱又は計測誤差である。実際の高炉プロセスでは、正確には図３に示す位置で外乱が重畳するとは限らないが、実用上、入力外乱ｄ₁及び出力外乱ｄ₂のいずれか近いものとみなしてよい。例えば、装入原料に起因する還元性の低下や性状ばらつきによる微粉炭分解熱の変化等は、入力外乱ｄ₁とみなすことができる。また、溶銑温度の測温ばらつき等は、出力外乱ｄ₂とみなすことができる。また、高炉現象を離散時間線形状態空間モデルで近似した際に生じるモデル誤差も、入力外乱ｄ₁又は出力外乱ｄ₂とみなすことができる。

次に、Ｌｕｅｎｂｅｒｇｅｒ型の状態推定オブザーバに、高炉プロセス１に重畳する入力外乱ｄ₁及び出力外乱ｄ₂を状態変数として考慮した外乱オブザーバ３の構築法について述べる。外乱オブザーバ３の構築法については、例えば非特許文献５等に記載されている。ここでは、より一般的に入力外乱と出力外乱とをそれぞれ備えた（３）式のように構成した場合について説明する。以下では、推定する外乱が定数値外乱である場合を例として説明するが、後述するように定数値外乱に限らない。

なお、オブザーバゲイン（Ｌ_x，Ｌ_d）については、例えばＬＱＲ設計法により離散型Ｒｉｃｃａｔｉ方程式を計算して求めた値を使用すればよい。

非負整数ｎ_d1，ｎ_d2は、それぞれ入力外乱の推定値ｄ₁＾の次元、出力外乱の推定値ｄ₂＾の次元であり、その合計をｎ_dとする（ｎ_d1＋ｎ_d2＝ｎ_d）。制御量ｚの次元ｒ≦ｎ_d≦観測量ｙの次元ｐであってよいが、ｎ_d＝ｐとして説明する。

ここで、Ｂ_d1（（４）式を参照）は、行列Ｂから入力外乱ｄ₁が重畳するとみなした操作量ｕに応じた列を抜き取って得られる行列である。また、Ｃ_d2（（５）式を参照）は、出力外乱ｄ₂が重畳するとみなした観測量ｙに応じて定まる行列である。例えば、外乱をすべて入力外乱ｄ₁として構成し、出力外乱ｄ₂を含めない場合、Ｂ_d1＝Ｂ,Ｃ_d2＝Ｏとなる。また、外乱をすべて出力外乱ｄ₂として構成し、入力外乱ｄ₁を含めない場合、Ｂ_d1＝Ｏ,Ｃ_d2＝Ｉとなる。
なお、Ｂ_d1とＣ_d2は、外乱が状態変数に与える影響を定める行列であり、上記で説明したような、入力外乱又は出力外乱に対応付けた形に限るものではないが、本実施形態では、この形で説明するものとする。

行列Ｂ_d1，Ｃ_d2の構成にあたり、具体的には、ｎ_d1＋ｎ_d2＝ｎ_d（＝ｐ）（０≦ｎ_d1≦操作量ｕの次元ｍ、０≦ｎ_d2≦観測量ｙの次元ｐ）を満たす入力外乱ｄ₁及び出力外乱ｄ₂の任意の組み合わせを選択することができる。このような組み合わせは（ｍ＋ｐ ｎ_d）（ｍ＋ｐが上段、ｎ_dが下段に記載されているものとする）通りある。例えばｍ＝２、ｐ＝４、ｎ_d＝４とすると、１５通りの構成法があることになる。

このように、行列Ｂ_d1，Ｃ_d2は、想定される外乱の重畳位置に応じて定まるものであり、すなわち、外乱の重畳位置の選択は、外乱オブザーバ３を構成する行列Ｂ_d1，Ｃ_d2を設計することに対応する。
行列Ｂ_d1，Ｃ_d2の設計指針について説明する。高炉プロセス１に重畳する外乱の位置、影響の程度に関し、事前知識が利用できる場合には、その知識に基づき、外乱オブザーバ３で考慮すべき外乱の重畳位置を選択し、行列Ｂ_d1，Ｃ_d2を設計すればよい。また、特段に事前知識がない場合は、（ｍ＋ｐ ｎ_d）通りの全ての組み合わせから、制御性能を表す評価関数の最悪性能が最も低い、すなわち、最もロバスト性の高い組み合わせのものを選択するようにしてもよい。

以下では、入力外乱の推定値ｄ₁＾及び出力外乱の推定値ｄ₂＾を（６）式のように一つのベクトルにまとめて、（３）式を（７）式のように表記する。

本実施形態では、制御量ｚ以外の観測量ｙを利用可能である場合も含めて、外乱オブザーバ３を構築する。例えば、非特許文献５に記載されるように、外乱オブザーバの構成においては、制御量ｚと観測量ｙを同一、すなわちＣ_z＝Ｃとして設定した上で、制御量ｚと同じ数の出力外乱のみを考慮したオブザーバを構築することが多いが、本発明では、入力外乱ｄ₁及び出力外乱ｄ₂の最適な組み合わせを考慮することによって、外乱発生時の操業変動の効果的な抑制効果を得ることができる。

また、外乱オブザーバ３で構成する外乱の定式化は定数値外乱に限るものではなく、ランプ状外乱であってもよい。また、一次遅れ系や正弦波等の周期性を持つダイナミクスを持たせてもよい。そのような外乱のダイナミクスは一般に（９）式のようにしてオブザーバに考慮すればよいが、例えば、（９）式において、特に、Ａ_ξ＝Ｏ,Ｂ_ξ＝Ｉとすれば、定数値外乱の場合に相当し、Ａ_ξ＝Ｉ,Ｂ_ξ＝Ｉとすれば、ランプ状外乱の場合に相当する。

［モデル予測制御］
外乱オブザーバ３とモデル予測制御とを組み合わせた制御について説明する。モデル予測制御では、オンラインで最適化問題を計算しながらフィードバック制御を実行する。モデル予測制御については、非特許文献５に詳しい。モデル予測制御では、状態推定オブザーバ（本実施形態では外乱オブザーバ３）により、現時刻ｋの状態変数を推定し、その推定値を予測軌道Ｒｅｃｅｄｉｎｇ Ｈｏｒｉｚｏｎの初期値とし、その初期値たる状態推定オブザーバの推定値に基づいて、将来の制御操作に対する予測軌道を計算する。
制御器２は、（１０）式のように、将来の制御操作に対する予測軌道と参照軌道（目標値ｒ）との誤差、及び将来時刻における操作量ｕの変化量の系列Δｕ＾（ｋ＋１｜ｋ），Δｕ＾（ｋ＋２｜ｋ），・・・，Δｕ＾（ｋ＋Ｈ_u｜ｋ）に関する２次形式の評価関数Ｖ（ｋ）を最小にする最適操作量の変化量の系列を計算するようにしてモデル予測制御を実行する。

ｘ＾（ｋ＋ｉ｜ｋ）は、現時刻ｋにおいて計算された、最適操作量の変化量の系列に対する将来時刻ｋ＋ｉにおける状態変数の推定値ｘ＾の値を表しており、ｚ＾（ｋ＋ｉ｜ｋ）も同様の意味である。また、ｒ（ｋ＋ｉ｜ｋ）は、現時刻ｋにおける計算で用いる将来時刻ｋ＋ｉにおける目標値ｒの値を意味する。また、Δｕ＾（ｋ＋ｉ｜ｋ）は、現時刻ｋにおける計算で用いる将来時刻ｋ＋ｉにおける最適操作量の変化量Δｕ＾の値を意味する。
図４は、時刻と制御量及び操作量との関係の例を示す図である。図４（ａ）に、時刻とｒ（ｋ＋ｉ｜ｋ），ｚ＾（ｋ＋ｉ｜ｋ）との関係を示す。また、図（ｂ）に、時刻とｕ（ｋ）＋Δｕ＾（ｋ＋ｉ｜ｋ）との関係を示す。

以下では、表記の簡単のため、（１０）式をベクトルと行列を用いて、（１１）式のように書き改める。

ここで、（１２）式のように、Ｚ（ｋ）は行列とベクトルの形でまとめられる。

外乱オブザーバ３で推定された入力外乱及び出力外乱を利用し、予測軌道と参照軌道との誤差である追従誤差ε（ｋ）を入力外乱と出力外乱の影響を考慮し、（１３）式のように定義する。このようにして、追従誤差を定義すれば、定常偏差を生じることなく目標値制御を行うことが可能となる。

さらに、適切に定義された行列とベクトルを用いて、（１０)式の評価関数Ｖ（ｋ）は、（１４）式のように２次形式で表現できる。

なお、モデル予測制御では、操作量ｕ及び制御量ｚのうち少なくともいずれか一方の制約条件が設定されるようにしてもよい。制約条件として、例えば操作量ｕの上下限制約、制御量ｚの上下限制約、操作量ｕの変化量の上下限制約、制御量ｚの変化量の上下限制約等がある。より具体的には、制御量ｚの制約条件としては、出銑量の上限値制約、炉熱指標の下限値制約等がある。また、操作量ｕの制約条件としては、送風量の上下限制約、送風量の変化量の上下限制約等がある。なお、制約条件はこれらに限るものではなく、線形不等式又は等式で表現できる制約であればよい。

以上より、モデル予測制御では、２次形式の評価関数と線形不等式制約からなる２次計画問題を求解する問題に帰着される。２次計画問題には各種ソルバーが整備され、高速求解が可能であるため、短い制御周期でのリアルタイム計算処理が実現できる。

なお、制約条件のうち一部に関しては制約違反を許容するようにしてもよい。制約条件を厳密に満たすべき制約（ハード制約）とそうでないもの（ソフト制約）とに分け、ソフト制約に関しては制約違反を許容する。例えば、操作量に関する制約はハード制約となるが、制御量に関する制約はソフト制約とする。また、ソフト制約の中で優先順位を設定することもできる。例えば、出銑量上限と溶銑温度下限の制約について、溶銑温度の下限制約をより優先する等を定式化すればよい。

また、観測量ｙとして、制御量ｚ以外の指標も含む場合に、（１５）式に示すように、評価関数Ｖ（ｋ）は、観測量の推定値の将来予測値の時間微分又は時間差分に関する二乗和の項をさらに追加してもよい。このようにすれば、観測量の時間変化を評価関数Ｖ（ｋ）に考慮することができ、応答の速応性向上が期待でき、外乱発生時の操業変動をより速やかに抑制することが可能になる。

図５は、高炉プロセス１の制御装置として機能する情報処理装置２００の構成例を示す図である。情報処理装置２００が、図３に示すように制御器２と外乱オブザーバ３とを備える高炉プロセス１の自動制御系を構成して、モデル予測制御を実行する。
情報処理装置２００は、ＣＰＵ２１１、ＲＯＭ２１２、ＲＡＭ２１３、補助記憶装置２１４、表示部２１５、操作部２１６、通信Ｉ／Ｆ（インターフェース）２１７、及びバス２１８を有する。
ＣＰＵ２１１は、ＲＯＭ２１２やＲＡＭ２１３に格納されているコンピュータプログラムやデータを用いて情報処理装置２００の全体を制御することで、制御器２及び外乱オブザーバ３の機能を実現して、モデル予測制御を実行する。なお、情報処理装置２００がＣＰＵ２１１とは異なる一又は複数の専用のハードウェアを有し、ＣＰＵ２１１による処理の少なくとも一部を専用のハードウェアが実行してもよい。専用のハードウェアの例としては、ＡＳＩＣ（Application Specific Integrated Circuit：特定用途向け集積回路）、ＦＰＧＡ（フィールドプログラマブルゲートアレイ）、及びＤＳＰ（デジタルシグナルプロセッサ）等がある。ＲＯＭ２１２は、変更を必要としないプログラム等を格納する。ＲＡＭ２１３は、補助記憶装置２１４から供給されるプログラムやデータ、及び通信Ｉ／Ｆ２１７を介して外部から供給されるデータ等を一時記憶する。補助記憶装置２１４は、例えばハードディスクドライブ等で構成され、画像データや音声データ等の種々のデータを記憶する。
表示部２１５は、例えば液晶ディスプレイやＬＥＤ等で構成され、ユーザが情報処理装置２００を操作するためのＧＵＩ（Graphical User Interface）等を表示する。操作部２１６は、例えばキーボードやマウス、ジョイスティック、タッチパネル等で構成され、ユーザによる操作を受けて各種の指示をＣＰＵ２１１に入力する。通信Ｉ／Ｆ２１７は、情報処理装置２００の外部の装置との通信に用いられる。例えば、情報処理装置２００が外部の装置と有線で接続される場合には、通信用のケーブルが通信Ｉ／Ｆ２１７に接続される。情報処理装置２００が外部の装置と無線通信する機能を有する場合には、通信Ｉ／Ｆ２１７はアンテナを備える。バス２１８は、情報処理装置２００の各部をつないで情報を伝達する。
なお、表示部２１５と操作部２１６が情報処理装置２００の内部に存在するものとして図示したが、表示部２１５と操作部２１６との少なくとも一方が情報処理装置２００の外部に別の装置として存在していてもよい。

以上述べたように、高炉プロセスの制御において、入力外乱ｄ₁と出力外乱ｄ₂を区別して重畳可能な設計としているので、外乱発生時に生じ得る操業変動を速やかに、精度良く抑制できる。

高炉操業においては、溶銑温度や溶銑中成分等の計測値又は分析値は、炉底の湯溜まり部を経て取得されるものであり、高炉内炉心部の反応タイミングに対して、ある程度の時間遅れが生じるとされる。一方、炉頂のガス組成や炉頂温度、送風圧力等は、溶銑中成分の分析値等と比較的すれば相対的には速やかに反応変化するとされる。そこで、観測量ｙとして、溶銑温度以外の炉頂で計測されるガス組成等の分析値も利用することで、操業中に発生した外乱（特に入力外乱）の影響を速やかに検知し、外乱発生時の操業変動を速やかに抑制することができる。観測量ｙとして利用するのが好ましい、比較的に速やかに反応変化する指標としては、主にガス組成等の炉頂で計測できる指標が考えられ、ガス組成以外にも、例えばガス利用率、カーボンソリューションロス量、直接還元反応量、羽口上熱バランス、炉頂温度、熱流比、及び送風圧力がある。

実施例では、制御量ｚとして出銑量及び溶銑温度を対象とし、操作量ｕとして送風量、酸素富化量及び微粉炭吹込み量を制御した。なお、酸素富化量は酸素富化率を一定として、送風量と連動し操作したため、操作端の実質の自由度は２となる。そのため、以下では、操作端として送風量及び微粉炭吹込み量のみ記載し、酸素富化量の記載を省略する。
また、制御量ｚ以外の観測量ｙとして、ガス組成（ＣＯ濃度、ＣＯ₂濃度）の観測量を利用するものとして説明する。この場合、観測量ｙの次元ｐは４、制御量ｚの次元ｒは２、操作量ｕの次元ｍは２である。
［制御量］出銑量、炉熱指標（溶銑温度）
［操作量］送風量、微粉炭吹込み量、酸素富化量 ※酸素富化量は酸素富化率一定として送風量と連動
［観測量］出銑量、溶銑温度、ガス組成（ＣＯ濃度、ＣＯ₂濃度）

（可制御可観測な状態空間モデルの構築）
２つの操作量ｕ（送風量、微粉炭吹込み量）と４つの観測量ｙ（出銑量、溶銑温度、ガス組成（ＣＯ濃度、ＣＯ₂濃度）との間の伝達関数行列を求め、Δｔ＝２０［ｍｉｎ］とし、０次ホールドで時間離散化した。このとき得られた状態空間モデルのＨａｎｋｅｌ特異値は０に近い値のものがあったので、８次で打ち切り、可制御可観測な状態空間モデルを構築した。

（外乱オブザーバの構築）
次に、外乱オブザーバ３を構築する。上述したように、外乱オブザーバ３で考慮可能な外乱の重畳位置には自由度があり、この例では、観測量ｙの次元ｐ＝４であるので、入力外乱ｄ₁及び出力外乱ｄ₂を合計４つ選択することができる。ここでは、外乱として、微粉炭分解熱の変化を想定した場合を説明する。この外乱は、操作量ｕに含まれる微粉炭の効果に影響を与えるため、入力外乱ｄ₁に相当する。したがって、ここでは、入力外乱ｄ₁を２つ選び、残る２つの外乱を出力外乱ｄ₂として、４つの観測量ｙから２つ選べばよい。事前検討から、ロバスト性の高いものとして、入力外乱ｄ₁を送風量、微粉炭吹込み量、出力外乱ｄ₂を出銑量、溶銑温度として、外乱オブザーバ３を構築した。

このようにして構築した外乱オブザーバ３とモデル予測制御とを組み合わせた制御において、外乱発生時の操業変動の抑制効果をシミュレーションにより検討した。比較例として、出銑量及び溶銑温度のみを観測量ｙとして外乱オブザーバ３を構築した場合の制御系において、外乱発生時の操業変動の抑制効果をシミュレーションにより検討した。
図６は、シミュレーション結果を示す図である。図６（ａ）は出銑量の時系列変化を表し、（ｂ）は溶銑温度の時系列変化を表す。点線が実施例の特性であり、実線が比較例の特性である。比較例では、制御量と観測量を同一とした。すなわち、出銑量及び溶銑温度のみを観測量ｙとし、出力外乱のみを考慮した外乱オブザーバを構築した。この場合も外乱抑制可能であり定常偏差を０にできるものの、実施例のようにより多くの観測量ｙを利用して外乱オブザーバ３を構築した場合には、顕著な外乱発生時の操業変動の抑制効果が得られる。この理由は、溶銑温度より比較的速やかに反応変化する指標（実施例ではガス組成）を利用し、入力外乱も含めて考慮することで、より速やかな外乱の検出と抑制が可能になるためである。

このように、高炉プロセスのような時定数の異なる物理現象が混在するプロセスでは、複数の観測量からそれらの時間変化の時定数の違いを利用することができる。すなわち、制御対象とする指標の時間変化に先立ち、別の指標を観測量として利用することで、速やかな外乱の検出と抑制が可能となる。本実施形態では、時定数の異なる物理現象を含むプロセスの例として高炉プロセスを対象として説明したが、これに限るものではない。例えば、一般に応答の遅い熱プロセスの温度制御において、より応答の速い流量や圧力が測定可能な場合等に適用できる。さらに、より一般には、マスターループ（メジャーループ）とスレーブループ（マイナーループ）を持つカスケード制御構造を有するプロセスにも適用可能である。

＜第２の実施形態＞
第２の実施形態では、外乱が測定可能である場合に、測定可能外乱に基づいて、制御操作のフィードフォワード補正機能を持たせるようにした例を説明する。なお、以下では、第１の実施形態との共通点についてはその詳細な説明を省略し、第１の実施形態との相違点を中心に説明する。
図７は、第２の実施形態に係る高炉プロセス１の自動制御系の構成を示す図である。図３に示すものと同様の構成要素には同一の符号を付し、その詳細な説明を省略する。本実施形態では、制御器２に、測定可能外乱ｄ_mを入力する。制御器２は、測定可能外乱ｄ_mに基づいて、制御操作のフィードフォワード補正を実行する補正手段として機能する。測定可能外乱ｄ_mとしては、装入原料（鉱石、コークス）の粒径や組成等がある。

以下、測定可能外乱ｄ_mに基づく最適操作量の計算について説明する。
図７では、測定可能外乱ｄ_mが高炉プロセス１に重畳するように図示するが、測定可能外乱ｄ_mを、高炉プロセス１の入側に重畳する入力外乱に相当するもの、高炉プロセス１の出側に重畳する出力外乱に相当するものに区別する。
ここでは、測定可能外乱ｄ_mが入力外乱であるものとして、その定式化を説明する。
この場合、離散時間線形状態空間モデルは、測定可能外乱ｄ_mを考慮し、（１６）式に示すようになる。なお、制御量ｚについては（１）式と同様であり、ここでは省略する。
また、高炉のインプット条件となる各種操業操作のうち、離散時間線形状態空間モデルに対する操作量uとして選択しなかった操作は、測定可能外乱ｄ_mとみなしてもよい。

この場合、状態変数ｘ及び制御量ｚの将来予測値は、（１７）式に示すように表される。

ここで、測定可能外乱の予測ベクトルＤ_m（ｋ）を、（１８）式のように定義する。

追従誤差ε（ｋ）については、（１８）式で求められた測定可能外乱を差し引き、計算すればよい。

このように、測定可能外乱ｄ_mに基づいて追従誤差ε（ｋ）を補正することで、最適操作量には、結果として、外乱に対するフィードフォワード補正の効果が付与される。
なお、外乱オブザーバ３で推定する入力外乱及び出力外乱も考慮する場合は、（１３）式で計算される追従誤差から、さらにΞＤ_m（ｋ）を差し引けばよい。

なお、測定可能外乱ｄ_mが出力外乱であった場合、離散時間線形状態空間モデルは、（２０）式に示すようになる。なお、制御量ｚについては（１）式と同様であり、ここでは省略する。

この場合、状態変数ｘ及び制御量ｚの将来予測値は、（２１）式に示すように表される。

追従誤差ε（ｋ）は、測定可能外乱を差し引き、（２２）式のように計算すればよい。

この場合にも、外乱オブザーバ３で推定する入力外乱及び出力外乱も考慮する場合は、同様にして（１３）式で計算される追従誤差から、さらに（Ｃ_m×Ｉ）ｄ_m（ｋ）（ここに
記載の「×」はクロネッカー演算子を表す。）を差し引けばよい。

上記では、入力外乱に相当する測定可能外乱と、出力外乱に相当する測定可能外乱それぞれを考慮する場合について、その定式化を説明したが、入力外乱に相当する測定可能外乱ｄ_m,1と出力外乱に相当する測定可能外乱ｄ_m,2を同時に考慮することもできる。

なお、入側位置であっても出側位置であっても、その位置に重畳する主要な外乱が測定できていれば、第１の実施形態で述べた外乱オブザーバ３において、外乱をその位置とは別の位置に設定したものを構築することができる。例えば入力外乱に相当する測定可能外乱があれば、入力外乱ｄ₁を減らし、出力外乱ｄ₂を増やして外乱オブザーバ３を設計することができる。このように第２の実施形態では、ある位置の外乱が測定可能であることで、外乱を測定するセンサを利用し、外乱オブザーバ３をより効果的に設計することができる。

なお、外乱の定式化は上記の定式化に限定されるものではなく、ダイナミクスを考慮してもよい。例えば、第１の実施形態と同様、測定可能外乱に対するプロセス変化を状態空間モデルにより定式化してもよい。操業中に測定された原料鉱石粒径の変化点前後の期間を含む実操業データ、又は非定常計算が可能な高炉数学モデルを利用し得られた時間応答データを利用し、これら時系列データから伝達関数行列を作成し、状態実現及び時間離散化により状態空間モデルを作成すればよい。
例えば、操作量と外乱との交互作用を無視し、制御操作の応答と測定可能外乱の応答との間に加算性が成り立つとすれば、（２６）式のような定式化が可能である。ここでは、外乱は出側位置に重畳するものとしているが、これに限らない。

このとき、外乱オブザーバ３は、（２７）式のように構築することができる。

上記の定式化で可観測性が成立しなくなるおそれがある場合には、オブザーバゲインＬ_m＝０として、より簡単な（２８）式の形で代用してもよい。

操作量と外乱との交互作用が無視できず、加算性が成り立たない場合には、例えば、交互作用効果を考慮した（２９）式のようなモデル、又はＫｏｏｐｍａｎモード分解やＣａｒｌｅｍａｎ線形化のような手法で線形化されたモデルを作成すればよい。

図８は、測定可能外乱として粒径（鉱石平均粒径）を測定した場合の計算結果を示す図である。図８（ａ）は粒径の時系列変化を表し、（ｂ）は微粉炭吹込み量の時系列変化を表し、（ｃ）は溶銑温度の時系列変化を表す。図８（ａ）に示すように、センサ情報から測定可能外乱である粒径変化を検知する。図８（ｂ）に示すように、操作量ｕである微粉炭吹込み量に対して、粒径変化に応じたフィードフォワード補正量を自動計算する（ＦＢ＋ＦＦ制御）。なお、比較例として、制御量と観測量を同一に設定し、第1の実施形態で述べたフィードバック制御のみを実行する例も載せる（ＦＢ制御）。この結果、図８（ｃ）に示すように、溶銑温度が変化する。なお、比較例として、ＦＢ制御、無制御（操作量の一定保持）の場合も載せる。制御操作のフィードフォワード補正機能を持たせることにより、顕著な外乱発生時の操業変動の抑制効果が得られている。

＜第３の実施形態＞
第３の実施形態では、高炉プロセス１の制御に利用される状態空間モデル（（１)式の離散時間線形状態空間モデル）を更新する機能を持たせるようにした例を説明する。第１、２の実施形態では、短期間での外乱変動を補償できるが、第３の実施形態では、それに加えて、比較的長期間での操業変動を補償できるようにすることを目的とする。なお、以下では、第１の実施形態との共通点についてはその詳細な説明を省略し、第１の実施形態との相違点を中心に説明する。

図９を参照して、第３の実施形態に係る高炉プロセスの自動制御系における更新装置４を説明する。図９（ａ）は、第３の実施形態に係る高炉プロセス１の自動制御系における更新装置４の機能構成例を示す図である。また、図９（ｂ）は、更新装置４が実行する処理例を示すフローチャートである。
更新装置４は、（１）式の離散時間線形状態空間モデル（以下、単に状態空間モデルと呼ぶ。）を更新する。例えば高炉プロセス１の制御装置として機能する情報処理装置２００（図５を参照）が更新装置４としても機能するように構成すればよい。更新装置４は、本発明でいう更新手段として機能する。
更新装置４は、判定部４１と、パラメータフィッティング計算部４２と、線形近似計算部４３とを有する。

ステップＳ１０１で、判定部４１は、状態空間モデルを更新するための所定の条件が成立したか否かを判定する。
本実施形態では、高炉プロセス１の所定の操業条件が変化したことが所定の条件に該当する（以下、条件１と呼ぶ。）。判定部４１は、各種センサで計測された操作量ｕ及び観測量ｙの測定値に基づいて、所定の操業条件が変化したことを検知する。所定の操業条件が変化するとは、操作量ｕ及び観測量ｙの瞬時値、或いは、所定の期間における時間平均値又は標準偏差等が所定の量を超えて変化したことをいう。なお、所定の期間は、実績データ等に基づいて予め設定されている。このような変化を引き起こす要因として、高炉プロセスの動作点（操業点）の変化や、鉱石の被還元性の変化（被還元性指数ＲＩの変化等）や、通気性変化等の高炉の内部状態の変化に相当するプロセスの性状の変化等がある。
また、所定の時間タイミングになったことが所定の条件に該当する（以下、条件２と呼ぶ。）。判定部４１は、例えば上位のプロセスコンピュータと通信を行うことで、所定の時間タイミングになったことを検知する。
また、高炉プロセス１で所定のイベント、例えば休風等のイベントが発生したことが所定の条件に該当する（以下、条件３と呼ぶ。）。判定部４１は、例えば上位のプロセスコンピュータと通信を行うことで、所定のイベントが発生したことを検知する。

判定部４１で条件１～３のいずれかが成立したと判定したとき、ステップＳ１０２で、パラメータフィッティング計算部４２は、状態空間モデルを作成するのに利用される数学モデルについてパラメータフィッティング計算を行う。
本実施形態では、判定部４１で条件１～３のいずれかが成立したと判定したとき、パラメータフィッティング計算部４２は、現時点での高炉数学モデルについてパラメータフィッティング計算を行う。条件１～３は、高炉数学モデルのパラメータの変化が発生していると想定されることから、パラメータフィッティング計算を行う。高炉数学モデルのパラメータについては後述する。
なお、条件１には、操業諸元の目標値変更等による高炉プロセスの動作点（操業点）の変化等も含まれるが、この場合、必ずしも高炉数学モデルのパラメータの変化は発生していないこともあり得る。このような場合には、後述するパラメータフィッティング計算処理における収束判定条件が満たされているため、結果として、パラメータフィッティング計算が省略されることがあり得る。

ステップＳ１０３で、線形近似計算部４３は、高炉数学モデルから状態空間モデルを作成する。その際、動作点周りでの線形近似計算を行う。
上述したように、判定部４１で条件１～３のいずれかが成立したと判定したとき、パラメータフィッティング計算部４２でパラメータフィッティング計算が行われるので、パラメータフィッティング計算後の高炉数学モデルから状態空間モデルを作成する。その際、動作点周りでの線形近似計算を行う。
なお、判定部４１で条件１～３のいずれも成立しない場合は、パラメータフィッティング計算及び線形近似計算はいずれも実行しない。

以上のようにした更新装置４は、線形近似計算部４３で作成した状態空間モデルを新しい状態空間モデルとして、更新を実行する（新しい状態空間モデルを現時点での状態空間モデルと置換する）。
なお、更新装置４がステップＳ１０１～Ｓ１０３を実行するようにしたが、その一部又は全ての処理を人が行うようにしてもよい。

以下、パラメータフィッティング計算部４２による処理を説明する。
まず、高炉数学モデルについて述べる。
高炉数学モデルは、高炉内に小領域を規定し、小領域内での挙動を塊状帯や融着帯での物質移動、反応、及び伝熱等の物理化学現象の計算式に基づいて、これらの時間発展をシミュレーションするものであり、高炉操業条件や原料性状から炉況を把握又は予測するために用いることができる。例えば三次元高炉数学モデルであれば、高炉の内部領域を高さ方向、径方向、周方向に分割することで複数のメッシュ（小領域）を規定し、各メッシュ内での挙動をシミュレーションする。高炉数学モデルについては種々の論文等が出されているが、例えばＫｏｊｉ ＴＡＫＡＴＡＮＩらの「Ｔｈｒｅｅ－ｄｉｍｅｎｓｉｏｎａｌ Ｄｙｎａｍｉｃ Ｓｉｍｕｌａｔｏｒ ｆｏｒ Ｂｌａｓｔ Ｆｕｒｎａｃｅ」，ＩＳＩＪ Ｉｎｔｅｒｎａｔｉｏｎａｌ，Ｖｏｌ．３９（１９９９），Ｎｏ．１，ｐ．１５－２２、又は、西岡らの「高炉数学モデルの開発」（新日鉄住金技報、第４１０号Ｐ．７３－７９、２０１８）等に記載の高炉数学モデルを好適に用いることができる。
例えば、上記のような非定常計算が可能な高炉数学モデルを利用すれば、高炉数学モデルで用いる各種モデルパラメータ、高炉原料の装入条件、及び送風条件等の各種境界条件の設定に対して、出銑量（トン／日）、コークス比ＣＲ（ｋｇ／ｔ）、微粉炭比ＰＣＲ（ｋｇ／ｔ）、送風圧力（ｈＰａ）、炉頂圧（ｈＰａ）、炉頂温度（℃）、ガス利用率ηＣＯ（％）、及び溶銑温度（℃）等の時間変化を計算可能である。

このような高炉数学モデルでは、高炉内で生じる様々な化学反応が考慮されているが、高炉数学モデルを好適に使用するためには、各種物理量の実測値とそのモデル計算値とが一致するように、これらの高炉内で生じる様々な化学反応の反応速度定数やその他モデル内パラメータを適切に調整しておくことが望ましい。このような物理量の実測値の例としては、送風圧力、Ｈ₂利用率、スラグ中Ｆｅｏ濃度、ＣＯ利用率、カーボンソリューション量（ＣＳＬ）、溶銑中Ｓｉ濃度、溶銑中Ｃ濃度、溶銑温度、出銑量、コークス比等がある。
また、調整対象の候補となるモデル内のパラメータの例としては、例えば西岡らの高炉数学モデルの開発」（新日鉄住金技報、第４１０号Ｐ．７３－７９、２０１８）に記載されているような、鉱石の還元反応、コークスのガス化反応、ＳｉＯ₂の還元反応、浸炭反応、水素ガスシフト反応等の化学反応の反応速度定数の補正係数や、固液熱交換係数の補正係数や炉内空隙率補正パラメータ等がある。

次に、パラメータフィッティング計算の前提条件について述べる。
パラメータフィッティング計算の前提条件として、例えば高炉の安定操業時であることがある。高炉の安定操業時においては、各操業諸元の変動は比較的小さく、所定の期間における各操業諸元の平均値や分散等の統計的性質がほぼ一定であるとみなせる。このような安定条件を具体的に判定する方法として、例えば出銑量や溶銑温度等の注目する諸元に関し、その瞬時値が所定の上下限内に収まるか否か、又はこれら諸元の時系列の時間平均値又は標準偏差又は変動係数等が所定の閾値内に収まるか否か等を判定するものが考えられる。

次に、パラメータフィッティング計算の計算対象となる、化学反応の反応速度定数の補正係数やそのモデル内パラメータについて述べる。
本実施形態では、間接還元反応の反応速度定数の補正係数（以降、単に間接還元反応と記載する。他の反応も同様）、コークスガス化反応、ＳｉＯ₂の還元反応、直接還元反応、浸炭反応、水素ガスシフト反応、固液熱交換係数、及び炉内の空隙率補正パラメータの８種類をパラメータフィッティング計算の計算対象とする。計算時間の短縮等の目的に応じて、８種類のすべてを計算対象とするのではなく、８種類の中から一部を計算対象として選択してもよい。
上記の８種類のパラメータのフィッティング対象となる観測量は、例えば、ＣＯ利用率（ηＣＯ）又はＣＳＬ、溶銑温度、溶銑中のＳｉ濃度、スラグＦｅＯ、溶銑中のＣ濃度、Ｈ₂利用率（ηＨ₂）、送風圧力等があるが、これに限らない。
パラメータフィッティング計算のために選定した観測量は、既に上に挙げた観測量と同じものでもよいし、別のものを選定してもよい。本実施形態では、別のものを選定している。

次に、パラメータフィッティング計算のアルゴリズムについて述べる。
パラメータフィッティング計算では、先に述べたような高炉プロセスがほぼ定常状態にある状況を想定する。このような場合に、現在時刻から所定の時間遡って時間平均された実績データを入力として高炉数学モデルの非定常計算を行い、定常状態まで計算させたモデル計算値（定常計算値と称する）と実績値が一致するように、パラメータをフィッティングする。なお、平均時間長については、例えば、操業諸元の実績値におけるマスバランスの一致を考慮するために、高炉に装入された原料が炉内を通過する時間を１単位として、その整数倍等の時間長を設定すればよい。例えば、８時間程度を１単位として、その整数倍等の時間長を設定すればよい。
すなわち、このパラメータフィッティング計算のステップでは、具体的には、（３０）式の連立非線形方程式を解く。

ここで、非線形関数Ｆのパラメータｐに関する勾配∂Ｆ／∂ｐが算出可能な場合には、Ｎｅｗｔｏｎ法等の反復法のアルゴリズムによって、少ない繰り返し回数で速やかにｐ*を求めることができるが、実際には、非線形関数Ｆは高炉数学モデルの計算手順に対応したものであり、勾配∂Ｆ／∂ｐを正確に算出することは困難である。そこで、以下で述べる手続きに従って、正確に算出された勾配∂Ｆ／∂ｐを利用することなく、（３０）式の非線形連立方程式を求解する。

パラメータフィッティング計算では、まず、計算対象とした８種類のパラメータの感度計算を実施する。感度は、各パラメータと各観測量の基準値周りでのスカラー量の摂動量の関係式で表す。特に、基準値からの観測量の摂動量Δｙと、基準値からのパラメータの摂動量Δｐとの関係は、概ね２次関数又は指数関数で良好に近似することができる。（３０）式は、２次関数で近似した場合の関係式を表している。図１０は、基準値からの観測量の摂動量Δｙと、基準値からのパラメータの摂動量Δｐとの関係を２次関数で近似した例を示す特性図である。これを全てのパラメータと観測量について計算しておく。すなわち、全てのパラメータと観測量について、（３１）式の関係式のa、ｂを求めておく。

パラメータフィッティングの具体的な計算手続きについて述べる。
本実施形態では、観測量の目標誤差δｙを解消するようパラメータ更新量を計算する。この段落でも、以下で定義するδｙとｐは、いずれもスカラー量であるものとして説明する。
ここで、観測量の目標誤差δｙは、（３２）式で定義される。ｐは、更新前のパラメータ値である。

例えば、（３３）式に示すように、感度の関係式を２次関数で近似した場合について、パラメータ更新量の計算式を説明する。感度計算実行時のパラメータ基準値ｐ₀からの摂動量δｘ（＝ｐ－ｐ₀）、観測量の目標誤差δｙを用いて、各ステップでのパラメータ更新量δｐを２次方程式の解の公式より算出する。なお、ａが小さい値をとる場合が多いので、桁落ちを防ぐため、根号を分母に移し、（３３）式のように、有理化しておくことが望ましい。根号の前の符号は、aの符号及び放物線の軸位置で場合分けする。

各パラメータと観測量には、おおむね以下の対応関係がある。間接還元反応に対応する観測量は、ガス利用率ηＣＯ又はＣＳＬである。コークスガス化反応に対応する観測量は、ガス利用率ηＣＯ又はＣＳＬのうち、間接還元反応に対応する観測量として選択しなかったものである。固液熱交換係数に対応する観測量は、溶銑温度である。ＳｉＯ₂還元反応に対応する観測量は、溶銑中Ｓｉである。直接還元反応に対応する観測量は、スラグＦｅＯである。浸炭反応に対応する観測量は、溶銑中のＣである。水素ガスシフト反応に対応する観測量は、水素利用率ηＨ₂である。空隙率補正パラメータに対応する観測量は、送風圧力である。
各パラメータ毎に、対応する観測量の目標誤差を解消するように、（３２）式のパラメータ更新計算を順番に行い、最終的に全ての観測量に関して、所定の収束判定条件を満たすまでモデル計算値が実績値に収束するまで、上記のパラメータ更新計算を繰り返す。

なお、上記の計算手続きのままでは収束が遅かったり、或いは、振動したり発散したりすることがあり、これを防ぐため、以下に述べる手順により、パラメータ更新計算の高速化及び安定化を実現することもできる。

まず、行方向に観測量、列方向にパラメータを並べた８×８サイズの行列を定義し、その各要素には、各観測量に与える各パラメータの影響が無い場合は０、影響が有る場合は１を格納する。ここで、パラメータの影響が無いとは、パラメータ摂動量に対する観測量の感度の絶対値が、所定の適当な閾値より小さい場合のことを指し、パラメータの影響が有るとは、パラメータ摂動量に対する観測量の感度の絶対値が、所定の適当な閾値以上である場合のことを指すものとする。さらに、この行列に対して、ＤＭ分解（Dulmage-Mendelsohn分解、非特許文献６に記載）することにより、ブロック下三角行列に並び替える。ＤＭ分解により並び替えることで、例えば、（３４）式のようなブロック下三角行列が得られる。ここで、行方向の観測量の並びは、スラグ中ＦｅＯ、溶銑中Ｃ濃度、Ｈ₂利用率、送風圧力、ＣＯ利用率、ＣＳＬ、溶銑温度、溶銑中Ｓｉであり、列方向のパラメータの並びは、直接還元反応、浸炭反応、水素ガスシフト反応、空隙率補正パラメータ、間接還元反応、コークスガス反応、ＳｉＯ２還元反応、固液熱交換係数となった。
この並び替えられた行列を、パラメータ影響行列と称するものとする。パラメータ影響行列の構造に基づいて、後述するパラメータ調整計算の計算順序を決定することができる。

（３４）式のパラメータ影響行列では、左上４×４行列の構造から、浸炭反応、直接還元反応、水素ガスシフト反応、及び空隙率補正パラメータの４つは互いに連関しておらず、その影響は独立したものと考えてよいことがわかる。このため、この４つは同時に更新してよいことがわかる。
また、右下４×４行列の構造から、間接還元反応、ＳｉＯ₂の還元反応、固液熱交換係数、コークスガス化反応の４つは互いに連関しており、それぞれの影響は独立したものと考えられない。このため、これら４つのパラメータはまとめずに順番に更新させる。後述する定式化に基づけば、互いの連関を考慮して、これら４つをまとめて更新することもできる。
以上、まとめると、例えば、以下の（１）～（５）を１サイクルとし、パラメータが収束するまでサイクルを回してパラメータ調整計算を進めていけばよい。
（１）浸炭反応、直接還元反応、水素ガスシフト反応、及び空隙率補正パラメータの４つ
（２）間接還元反応
（３）コークスガス化反応
（４）ＳｉＯ２の還元反応
（５）固液熱交換係数
なお、上記の（１）～（５）は順不同である。
図１１は、パラメータフィッティング計算の例を説明するための図であり、（ａ）が計算ステップと、間接還元反応の観測量であるガス利用率ηＣＯとの関係の例を示す特性図、（ｂ）が計算ステップと、間接還元反応のパラメータとの関係の例を示す特性図である。図１１（ａ）、（ｂ）に示すように、観測量の実績値（目標値）と計算値とが一致するようにパラメータが収束する。

なお、パラメータ影響行列の構造より、間接還元反応、ＳｉＯ₂の還元反応、固液熱交換係数、コークスガス化反応の４つは互いに連関しており、その影響は独立したものと考えられないため、これらの連関を考慮して、（３５）式の非線形計画問題を解くことにより、（２）、（３）、（４）、（５）はまとめて計算することができる。まとめて計算することで、計算時間の短縮が可能となる。
Given：重みＱ_i，Ｒ_i＞０、誤差δ_i、上下限Δ_pi,min，Δ_pi,max、２次関数近似パラメータａ_ij，ｂ_ij
決定変数（パラメータの更新量）：Δｐ_i

以上述べた手続きに従って、正確に算出された勾配∂Ｆ／∂ｐを利用することなく、２次関数や指数関数等で近似された感度を利用することにより、非線形方程式である連立方程式（（３０）式）を満たすモデルパラメータを正確に計算することができる。なお、この計算では、高炉数学モデルの定常計算値をフィッティングさせるものであり、その定常計算値は、高炉数学モデルの非定常計算開始時の初期状態によらないので、任意の初期状態で計算を開始してよい。特に、計算の高速化を図るために、前回計算終了時の計算結果を初期状態として設定してもよい。
また、パラメータフィッティング計算の最初に実行する感度計算については、パラメータフィッティング計算の都度毎回実施するようにしてもよいし、計算時間の短縮のため、都度の計算は省略し、予め算出しておいた感度計算結果で代用するようにしてもよい。

以下、線形近似計算部４３による処理を説明する。
前記パラメータフィッティング計算された高炉数学モデルを利用した非定常計算により、操作量ｕ∈Ｒ^mを高炉数学モデルの入力、制御量ｚ∈Ｒ^rを含む観測量ｙ∈Ｒ^pを出力とする適当な時間応答を計算する。例えば、各操作量を単独でステップ状に変更して時間応答（ステップ応答）を計算させればよい。この時間応答データを一次遅れ系又は二次遅れ系等の線形システムの伝達関数でフィッティングし、行方向に観測量を並べて、列方向に操作量を並べた、サイズp×ｍの伝達関数行列を作成する。なお、操作量uのほか、外乱ｄ∈Ｒ^mdを与えて同様の計算を行い、サイズp×(ｍ＋ｍｄ)の伝達関数行列に変換してもよい。また、上記フィッティングさせる伝達関数については、操作量の影響遅れに対応したむだ時間を考慮してもよい。
ここで、操作量は、例えば、送風量、微粉炭量に関する値、酸素富化に関する値、送風温度、湿分量、及び液化天然ガス量、コークス比又はＯ／Ｃ等である。外乱は、例えば、装入原料（鉱石、コークス）の粒径や組成等である。観測量は、出銑量及び炉熱指標（溶銑温度、溶銑中Ｓｉ、カーボンソリューションロス量、直接還元反応量、羽口上熱バランス、及び融着帯レベル等）を含む制御量のほか、ガス組成、ガス利用率、カーボンソリューションロス量、直接還元反応量、羽口上熱バランス、炉頂温度、熱流比、及び送風圧力等である。

図１２（ａ）は、送風量の変更操作に対する出銑量の時間応答の計算結果を示す特性図である。線形近似計算部４３では、図１２（ａ）に示すように、送風量の操作量変更操作に対する出銑量の時間応答を計算する。この時間応答は、図中の点線で示すように、一次遅れ系の伝達関数で近似することができる。

図１２（ｂ）は、送風量の変更操作に対する溶銑温度の時間応答の計算結果を示す特性図である。線形近似計算部４３では、図１２（ｂ）に示すように、送風量の操作量変更操作に対する溶銑温度の時間応答を計算する。この時間応答は、図中の点線で示すように、二次遅れ系の伝達関数で近似することができる。

このようにして得られた伝達関数行列を可制御正準形等の適当な正準形で状態実現したうえで、０次ホールド等の手法で適当なサンプリング時間で時間離散化すれば、離散時間の状態空間モデルが得られるが、通常、上記の伝達関数行列は冗長なダイナミクスを含むため、このままでは、可制御可観測性を満たす状態空間モデルを得ることができない。そこで、平衡実現法を適用して、可制御可観測な状態空間モデルを構成する。平衡実現法は、ＭＩＭＯ（多入力多出力系）線形モデルに対する低次元化法であり、元モデルを可制御可観測な線形状態空間モデルに縮約することができる（非特許文献３を参照）。

ここで、可制御可観測な線形状態空間モデルを作成する理由について述べる。一般に、可制御可観測な線形状態空間モデルであれば、例えば、Ｌｕｅｎｂｅｒｇｅｒ型オブザーバ等の、観測量を利用したモデル内変数の状態推定計算手法が適用できるが、高炉数学モデルは大規模非線形モデルであるとともに、線形システムにおける可観測性に相当する性質の成立が保証されていないため、そのままでは、Ｌｕｅｎｂｅｒｇｅｒ型オブザーバ等の状態推定計算手法が適用できない。また、同様に、線形システムにおける可制御性に相当する性質の成立も保証されていないため、モデル予測制御における最適操作量の計算もまた困難である。
高炉操業の場合、操業点の動作範囲を限定することによって、概ね線形性が成立するため、比較的安定状態における高炉数学モデルを線形近似モデルに置き換えることは有効といえる。また、上記手続きにより、可観測性や可制御性が担保された線形近似モデルを利用することにより、状態推定計算やモデル予測制御による最適操作量の演算が可能となり、好適に高炉自動制御を実施することが可能となる。

なお、本実施形態では、判定部４１で条件１～３が成立したと判定したとき、パラメータフィッティング計算部４２は、現時点での高炉数学モデルについてパラメータフィッティング計算を行うものとしたが、例えば高炉の安定状態が継続している場合、パラメータフィッティング計算を省略し、代表的なパラメータで固定したうえで、線形近似計算部４３で動作点周りでの線形近似計算を実行するようにしてもよい。
また、線形近似計算実行後に、バンプレスに状態空間モデルを更新するために、モデル更新直後に、スムージング等の計算処理を適宜行い、予め状態変数を適切に更新しておいてもよい。

上記実施形態は、本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその技術思想、又はその主要な特徴から逸脱することなく、様々な形で実施することができる。
また、本発明を適用した高炉プロセスの制御の機能は、ソフトウェア（プログラム）を、ネットワーク又は各種記憶媒体を介してシステム或いは装置に供給し、そのシステム或いは装置のコンピュータがプログラムを読み出して実行することによっても実現可能である。また、本発明を適用した高炉プロセスの制御の機能は、ＰＬＣ（Programmable Logic Controller）により実現されてもよいし、ＡＳＩＣ等の専用のハードウェアにより実現してもよい。

次のように記載することも可能である。
［請求項１］
プロセスを制御対象として、状態変数の時間変化を表す状態空間モデルを利用したモデル予測制御を実行する制御装置であって、
前記プロセスに操作量を入力したときの前記プロセスから出力される制御量を含む観測量を取得し、
前記操作量に入力外乱が重畳可能な設計とし、前記観測量に出力外乱が重畳可能な設計とするものとして、
前記入力外乱及び出力外乱のうち少なくともいずれかの外乱が重畳した前記操作量及び観測量を入力として、前記状態変数の推定値と、前記入力外乱の推定値及び前記出力外乱の推定値のうち少なくともいずれかとを出力する推定手段を有することを特徴とする制御装置。
［請求項２］
前記推定手段への入力は、前記制御量以外の観測量を含むことを特徴とする請求項１に記載の制御装置。
［請求項３］
前記プロセスは、時定数の異なる物理現象を含むプロセスであることを特徴とする請求項１又は２に記載の制御装置。
［請求項４］
目標値及び前記推定手段の出力を入力として、前記モデル予測制御を実行することにより、前記操作量を出力する制御手段を有することを特徴とする請求項１乃至３のいずれか１項に記載の制御装置。
［請求項５］
前記プロセスは、高炉プロセスであり、
前記制御量は、出銑量及び炉熱指標を含み、
前記炉熱指標は、溶銑温度、溶銑中Ｓｉ、カーボンソリューションロス量、直接還元反応量、羽口上熱バランス、及び融着帯レベルのいずれか一つを含むことを特徴とする請求項１乃至４のいずれか１項に記載の制御装置。
［請求項６］
前記プロセスは、高炉プロセスであり、
前記操作量は、送風量、微粉炭量に関する値、酸素富化に関する値、送風温度、湿分量、及び液化天然ガス量、コークス比、Ｏ／Ｃに関する値のいずれか一つを含むことを特徴とする請求項１乃至５のいずれか１項に記載の制御装置。
［請求項７］
前記プロセスは、高炉プロセスであり、
前記観測量は、前記制御量を含み、かつ、ガス組成、ガス利用率、カーボンソリューションロス量、直接還元反応量、羽口上熱バランス、炉頂温度、熱流比、及び送風圧力のいずれか一つであって前記制御量に含まれていないものを含むことを特徴とする請求項１乃至６のいずれか１項に記載の制御装置。
［請求項８］
前記操作量及び前記制御量のうち少なくともいずれか一方の制約条件が設定されていることを特徴とする請求項１乃至７のいずれか１項に記載の制御装置。
［請求項９］
前記制約条件のうち一部に関しては制約違反を許容することを特徴とする請求項８に記載の制御装置。
［請求項１０］
前記一部の前記制約条件の中で優先順位が設定されていることを特徴とする請求項９に記載の制御装置。
［請求項１１］
前記入力外乱の推定値の次元と前記出力外乱の推定値の次元との合計を、前記制御量の次元以上、前記観測量の次元以下とすることを特徴とする請求項１乃至１０のいずれか１項に記載の制御装置。
［請求項１２］
測定可能外乱に基づいて、制御操作のフィードフォワード補正を実行する補正手段を有することを特徴とする請求項１乃至１１のいずれか１項に記載の制御装置。
［請求項１３］
所定の条件が成立したとき、前記状態空間モデルを更新する更新手段を有することを特徴とする請求項１乃至１２のいずれか１項に記載の制御装置。
［請求項１４］
前記所定の条件は、前記プロセスの所定の操業条件が変化したこと、前記プロセスの所定の性状が変化したこと、所定の時間タイミングになったこと、及び前記プロセスで所定のイベントが発生したことのうち少なくともいずれか一つであることを特徴とする請求項１３に記載の制御装置。
［請求項１５］
前記更新手段は、前記状態空間モデルを作成するのに利用される数学モデルについてパラメータフィッティング計算を行うパラメータフィッティング計算手段を有することを特徴とする請求項１３又は１４に記載の制御装置。
［請求項１６］
前記パラメータフィッティング計算手段は、前記数学モデルを利用し、パラメータに関する感度計算を実行し、ブロック下三角構造に並び替えられたパラメータ影響行列の構造に基づいて、パラメータを更新することを特徴とする請求項１５に記載の制御装置。
［請求項１７］
前記プロセスは、高炉プロセスであり、
前記数学モデルは、高炉数学モデルであり、
前記パラメータフィッティング計算手段は、間接還元反応、コークスガス化反応、ＳｉО₂還元反応、直接還元反応、浸炭反応、水素ガスシフト反応のうちのそれぞれの反応速度定数の補正係数、固液熱交換係数、及び空隙率補正パラメータのうち少なくともいずれか一つを計算対象としてパラメータフィッティング計算を行うことを特徴とする請求項１５又は１６に記載の制御装置。
［請求項１８］
前記更新手段は、前記パラメータフィッティング計算手段でのパラメータフィッティング計算前の前記数学モデル、又は前記パラメータフィッティング計算手段でのパラメータフィッティング計算後の前記数学モデルから前記状態空間モデルを作成し、そのときに動作点周りでの線形近似計算を行う線形近似計算手段を有することを特徴とする請求項１５乃至１７のいずれか１項に記載の制御装置。
［請求項１９］
前記状態空間モデルは、離散時間線形状態空間モデルであり、
前記線形近似計算手段は、前記数学モデルを利用し、前記操作量を入力、前記制御量を出力とする時間応答を一次遅れ系又は二次遅れ系の伝達関数でフィッティングし、伝達関数行列を構成した後、平衡実現法を適用して、可制御可観測な前記離散時間線形状態空間モデルを構成することを特徴とする請求項１８に記載の制御装置。
［請求項２０］
プロセスを制御対象として、状態変数の時間変化を表す状態空間モデルを利用したモデル予測制御を実行する制御方法であって、
前記プロセスに操作量を入力したときの前記プロセスから出力される制御量を含む観測量を取得し、
前記操作量に入力外乱が重畳可能な設計とし、前記観測量に出力外乱が重畳可能な設計とするものとして、
前記入力外乱及び出力外乱のうち少なくともいずれかの外乱が重畳した前記操作量及び観測量を入力として、前記状態変数の推定値と、前記入力外乱の推定値及び前記出力外乱の推定値のうち少なくともいずれかとを出力するステップを有することを特徴とする制御方法。
［請求項２１］
プロセスを制御対象として、状態変数の時間変化を表す状態空間モデルを利用したモデル予測制御を実行するためのプログラムであって、
前記プロセスに操作量を入力したときの前記プロセスから出力される制御量を含む観測量を取得し、
前記操作量に入力外乱が重畳可能な設計とし、前記観測量に出力外乱が重畳可能な設計とするものとして、
前記入力外乱及び出力外乱のうち少なくともいずれかの外乱が重畳した前記操作量及び観測量を入力として、前記状態変数の推定値と、前記入力外乱の推定値及び前記出力外乱の推定値のうち少なくともいずれかとを出力する推定手段としてコンピュータを機能させるためのプログラム。

１：高炉プロセス、２：制御器、３：外乱オブザーバ、４：更新装置、４１：判定部、４２：パラメータフィッティング計算部、４３：線形近似計算部

Claims (21)

1. プロセスを制御対象として、状態変数の時間変化を表す状態空間モデルを利用したモデル予測制御を実行する制御装置であって、
   前記プロセスに操作量を入力したときの前記プロセスから出力される制御量を含む観測量を取得し、
   前記操作量に入力外乱が重畳可能な設計とし、前記観測量に出力外乱が重畳可能な設計とするものとして、
   前記入力外乱及び出力外乱のうち少なくともいずれかの外乱が重畳した前記操作量及び観測量を入力として、前記状態変数の推定値と、前記入力外乱の推定値及び前記出力外乱の推定値のうち少なくともいずれかとを出力する推定手段を有することを特徴とする制御装置。
2. 前記推定手段への入力は、前記制御量以外の観測量を含むことを特徴とする請求項１に記載の制御装置。
3. 前記プロセスは、時定数の異なる物理現象を含むプロセスであることを特徴とする請求項１又は２に記載の制御装置。
4. 目標値及び前記推定手段の出力を入力として、前記モデル予測制御を実行することにより、前記操作量を出力する制御手段を有することを特徴とする請求項１又は２に記載の制御装置。
5. 前記プロセスは、高炉プロセスであり、
   前記制御量は、出銑量及び炉熱指標を含み、
   前記炉熱指標は、溶銑温度、溶銑中Ｓｉ、カーボンソリューションロス量、直接還元反応量、羽口上熱バランス、及び融着帯レベルのいずれか一つを含むことを特徴とする請求項１又は２に記載の制御装置。
6. 前記プロセスは、高炉プロセスであり、
   前記操作量は、送風量、微粉炭量に関する値、酸素富化に関する値、送風温度、湿分量、及び液化天然ガス量、コークス比、Ｏ／Ｃに関する値のいずれか一つを含むことを特徴とする請求項１又は２に記載の制御装置。
7. 前記プロセスは、高炉プロセスであり、
   前記観測量は、前記制御量を含み、かつ、ガス組成、ガス利用率、カーボンソリューションロス量、直接還元反応量、羽口上熱バランス、炉頂温度、熱流比、及び送風圧力のいずれか一つであって前記制御量に含まれていないものを含むことを特徴とする請求項１又は２に記載の制御装置。
8. 前記操作量及び前記制御量のうち少なくともいずれか一方の制約条件が設定されていることを特徴とする請求項１又は２に記載の制御装置。
9. 前記制約条件のうち一部に関しては制約違反を許容することを特徴とする請求項８に記載の制御装置。
10. 前記一部の前記制約条件の中で優先順位が設定されていることを特徴とする請求項９に記載の制御装置。
11. 前記入力外乱の推定値の次元と前記出力外乱の推定値の次元との合計を、前記制御量の次元以上、前記観測量の次元以下とすることを特徴とする請求項１又は２に記載の制御装置。
12. 測定可能外乱に基づいて、制御操作のフィードフォワード補正を実行する補正手段を有することを特徴とする請求項１又は２に記載の制御装置。
13. 所定の条件が成立したとき、前記状態空間モデルを更新する更新手段を有することを特徴とする請求項１又は２に記載の制御装置。
14. 前記所定の条件は、前記プロセスの所定の操業条件が変化したこと、前記プロセスの所定の性状が変化したこと、所定の時間タイミングになったこと、及び前記プロセスで所定のイベントが発生したことのうち少なくともいずれか一つであることを特徴とする請求項１３に記載の制御装置。
15. 前記更新手段は、前記状態空間モデルを作成するのに利用される数学モデルについてパラメータフィッティング計算を行うパラメータフィッティング計算手段を有することを特徴とする請求項１３に記載の制御装置。
16. 前記パラメータフィッティング計算手段は、前記数学モデルを利用し、パラメータに関する感度計算を実行し、ブロック下三角構造に並び替えられたパラメータ影響行列の構造に基づいて、パラメータを更新することを特徴とする請求項１５に記載の制御装置。
17. 前記プロセスは、高炉プロセスであり、
    前記数学モデルは、高炉数学モデルであり、
    前記パラメータフィッティング計算手段は、間接還元反応、コークスガス化反応、ＳｉО₂還元反応、直接還元反応、浸炭反応、水素ガスシフト反応のうちのそれぞれの反応速度定数の補正係数、固液熱交換係数、及び空隙率補正パラメータのうち少なくともいずれか一つを計算対象としてパラメータフィッティング計算を行うことを特徴とする請求項１５に記載の制御装置。
18. 前記更新手段は、前記パラメータフィッティング計算手段でのパラメータフィッティング計算前の前記数学モデル、又は前記パラメータフィッティング計算手段でのパラメータフィッティング計算後の前記数学モデルから前記状態空間モデルを作成し、そのときに動作点周りでの線形近似計算を行う線形近似計算手段を有することを特徴とする請求項１５に記載の制御装置。
19. 前記状態空間モデルは、離散時間線形状態空間モデルであり、
    前記線形近似計算手段は、前記数学モデルを利用し、前記操作量を入力、前記制御量を出力とする時間応答を一次遅れ系又は二次遅れ系の伝達関数でフィッティングし、伝達関数行列を構成した後、平衡実現法を適用して、可制御可観測な前記離散時間線形状態空間モデルを構成することを特徴とする請求項１８に記載の制御装置。
20. プロセスを制御対象として、状態変数の時間変化を表す状態空間モデルを利用したモデル予測制御を実行する制御方法であって、
    前記プロセスに操作量を入力したときの前記プロセスから出力される制御量を含む観測量を取得し、
    前記操作量に入力外乱が重畳可能な設計とし、前記観測量に出力外乱が重畳可能な設計とするものとして、
    前記入力外乱及び出力外乱のうち少なくともいずれかの外乱が重畳した前記操作量及び観測量を入力として、前記状態変数の推定値と、前記入力外乱の推定値及び前記出力外乱の推定値のうち少なくともいずれかとを出力するステップを有することを特徴とする制御方法。
21. プロセスを制御対象として、状態変数の時間変化を表す状態空間モデルを利用したモデル予測制御を実行するためのプログラムであって、
    前記プロセスに操作量を入力したときの前記プロセスから出力される制御量を含む観測量を取得し、
    前記操作量に入力外乱が重畳可能な設計とし、前記観測量に出力外乱が重畳可能な設計とするものとして、
    前記入力外乱及び出力外乱のうち少なくともいずれかの外乱が重畳した前記操作量及び観測量を入力として、前記状態変数の推定値と、前記入力外乱の推定値及び前記出力外乱の推定値のうち少なくともいずれかとを出力する推定手段としてコンピュータを機能させるためのプログラム。

Images