WO2016152618A1 - 制御モデルのパラメータと外乱との同時推定方法、及びこの同時推定方法を用いた制御対象の制御方法 - Google Patents

Abstract

　制御モデル内のパラメータと制御モデルが表す制御対象へ入力される外乱とを同時に推定する方法であって、外乱を時系列変数値として、制御モデルに与えることとし、制御対象における実績データの時系列値と、制御モデルをシミュレーションすることで得られたシミュレーション結果のデータの時系列値を含む評価関数が、最大又は最小となるように、前記外乱と前記パラメータとを同時に推定することを特徴とする。

Description

制御モデルのパラメータと外乱との同時推定方法、及びこの同時推定方法を用いた制御対象の制御方法

　本発明は、外乱が加わる動的システムに関し、制御モデルのパラメータと外乱とを同時に推定することのできる推定方法と、この推定方法を用いた制御対象の制御方法に関する。

　自然法則に従う動的な対象をモデル化し、モデルに基づいた制御（モデルベースド制御）を行う際には、その制御モデルの精度が重要である。制御モデルの精度は、例えば、フィードフォワード量の決定や、フィードバックゲインの決定に必要である。
　また、アクチュエータの飽和を避けるために、制御モデルに加わる外乱も同時に精度よく知る必要がある。その理由としては、例えば、外乱が大きいと、フィードバックゲインを大きくすることができず、アクチュエータの飽和を避けることができないため、外乱も精度よく推定する必要がある。

　さらに、制御系に異常状態（例えば、ハンチング、定常偏差が生じる等）が発生しているときには、その異常の原因が制御系に起因するものか、又は外乱に起因するものかによっては、講じる対策が異なるため、制御モデルの制御パラメータと、制御モデルに加わる外乱を同時に推定する必要がある。例えば、制御系による起因であれば、制御ゲインを小さくする必要があり、外乱による起因であれば、制御ゲインを大きくすることが考えられるからである。

　また、プロセスの最適化や異常診断等の他の用途にも、制御モデル（制御パラメータ）の精度や、制御モデルに加わる外乱の精度が重要となる。
　このような、制御モデルに基づいて制御系を制御する技術としては、特許文献１～特許文献１２、非特許文献１～非特許文献５に開示された技術がある。
　特許文献１には、外乱推定オブザーバによる衝突検出方法において、速度ループに対して構成された外乱推定オブザーバを用いて外乱を推定して、この推定した外乱に基づいて産業用ロボット等の衝突を検出する方法が開示されている。すなわち、特許文献１は、外乱のみを推定する技術である。

　特許文献２には、自律的制御方法および制御システムにおいて、制御モデルのパラメータをチューニングする方法が開示されている。すなわち、特許文献２は、制御モデルのパラメータのみを推定して（外乱は推定せず）、制御モデルを調整する技術である。
　特許文献３には、連続鋳造の鋳型内湯面レベル制御方法において、外乱オブザーバで外乱を推定し、推定した外乱を打ち消すようにオブザーバ出力を計算し、オブザーバ出力を加算した開度指令値で、タンディッシュのスライディングノズルの開度を制御する方法が開示されている。すなわち、特許文献３は、外乱のみを推定して（パラメータは推定せず）、スライディングノズルの開度指令値を算出する技術である。

　特許文献４には、制御対象同定方法において、予め定められた所定の同定用入力信号及び、制御対象の出力信号から制御対象の特性を推定し、推定された制御対象の特性などから、未来の制御対象の状態を予測し、予測された制御対象の状態が制約に対して満足するまで修正して制御対象モデルを作成する方法が開示されている。すなわち、特許文献４は、所定の入力信号を入力し、制御モデルのパラメータのみを推定して（外乱は推定せず）、制御モデルを作成する技術である。

　特許文献５には、プロセスモデルのパラメータ調整装置、調整支援装置及び方法において、プロセス外部入力データが入力されたプロセスモデルにより計算されるプロセス出力シミュレーション値と、収集されているプロセス出力値との誤差を最小にするように、パラメータの値を同定する方法が開示されている。すなわち、特許文献５は、制御モデルのパラメータのみを推定して（外乱は推定せず）、パラメータの値を同定する技術である。

　なお、一般に、同文献の段落［０００２］～［００１０］に記載されているように、現実のプロセスを詳細に模擬することのできるプロセスシミュレータがあれば、設計支援、運転管理支援、シナリオシミュレーション、ソフトウエアセンサー・オブザーバ、モデルベース制御、短期のプロセス最適化、長期のプロセス最適化、教育／運転訓練等が可能である。

　特許文献６には、ΔΣ変調アルゴリズムを用いてプラントを制御する制御装置において、同文献の式（２）に表されている、制御対象の制御モデルに印加されている外乱推定値ｃ１を、部分モデルパラメータ同定器により、モデル化誤差がなくなるように、逐次的に同定する方法が開示されている。
　ところで、特許文献６は、逐次的に同定するため、同文献の段落［００４４］に記載されているように、収束時間が短縮されるとはいえ、実際にはモデルパラメータの最適値への収束時間がかかるものとなっている。すなわち、同文献で得られる外乱やパラメータは、時間的に遅れた値となる。また、同文献の式（２）は線形の式であり、非線形の場合は開示されていない。また、同文献の式（２）のような１つの式で表せない場合、例えば、むだ時間を有する場合や、モデルパラメータに対して、線形式で表せない場合には、対応するものとはなっていない。

　特許文献７には、シミュレーションモデルの同定方法およびそのプログラムにおいて、モデル出力と実機の出力との残差に基づいて、複数パラメータの外乱の大きさを同定する方法が開示されている。すなわち、特許文献７は、同文献の段落［００２１］や図３に記載されているように、複数の外乱モデルを用意して、その外乱モデルの出力に乗算する大きさの最適値を求めるものである（外乱は推定せず）。

　特許文献８には、位置制御装置において、駆動手段への位置指令等から外乱補償量を算出し、算出された外乱補償量を駆動手段への入力に足し合わせて外乱を補償する方法が開示されている。この特許文献８は、例として、同文献の段落［００２８］に記載されている式（５）用いるため、任意の時系列パターンをもつ外乱を推定することはできないものとなっている。

　特許文献９には、ΔΣ変調アルゴリズムを用いてプラントを制御する制御装置おいて、同文献の式（２）に表されている、制御対象の制御モデルに印加されている外乱推定値ｃ１を、部分モデルパラメータ同定器により、モデル化誤差がなくなるように、逐次的に同定する方法が開示されている。ところで、特許文献９は、上記した特許文献６を発展させた技術であり、同文献で得られる外乱やパラメータは、時間的に遅れた値となっている。

　特許文献１０には、外乱推定装置、制御対象モデル推定装置、フィードフォワード量推定装置および制御装置において、複数の運転条件でそれぞれ計測される複数の外乱波形に基づいて、その複数の運転条件とは異なる他の運転条件の外乱波形を推定する方法が開示されている。すなわち、特許文献１０は、予め記憶された外乱波形に基づいて、異なる他の運転条件の外乱波形を推定するものであり、任意の時系列パターンをもつ外乱を推定することはできないものとなっている。

　特許文献１１には、連続鋳造機の湯面レベル制御系における周波数推定装置、方法及びプログラムにおいて、逐次最小２乗法（ＲＬＳ法）により、外乱の周波数を推定し、非線形最適化手法により、外乱の周波数成分を推定する方法が開示されている。すなわち、特許文献１１は、外乱のみを推定する方法である。
　特許文献１２には、ＰＩＤ調整装置およびＰＩＤ調整プログラムにおいて、ＰＳＯ（Particle Swarm Optimization）により、ＰＩＤ調節器の比例ゲイン、積分時間、微分時間を調整する方法が開示されている。ところで、特許文献１２は、外乱を用いているものの、その外乱は推定されたものではない。

　非特許文献１には、周期外乱の影響を受ける線形時不変システムの同定と外乱の推定において、周期外乱の影響を受ける線形時不変システムの同定と外乱の推定に関する技術が開示されている。ところで、非特許文献１は、外乱が周期的であり、対象となるシステムが線形時不変なものに限定されるなど、実用上の制限がある。
　非特許文献２には、独立成分分析によるパラメータと外乱の同時推定において、独立成分分析によるパラメータと外乱の同時推定に関する技術が開示されている。ところで、非特許文献２は、信号の統計的性質に差があること等、独立成分分析を適用するための条件を満たす必要がある。

　非特許文献３には、オブザーバを用いたイナーシャ誤差及び摩擦外乱の同時推定に関する研究において、同時推定の技術が開示されている。ところが、非特許文献３は、外乱のモデルを仮定して行っているが、実際の外乱とは異なることが考えられる。
　非特許文献４には、外乱によって生成された入出力データを用いた外乱抑制ＦＲＩＴ法において、大きさが未知の外乱のパラメータと制御器のパラメータを同時推定している。ところで、非特許文献４の外乱は、インパルス状外乱であり、制限があるので、実際の外乱と異なることが考えられる。

　非特許文献５には、ＰＩＤ制御において、多変数制御系のＰＩＤ調節計調整法が開示されている。ところが、非特許文献５では、外乱を用いているものの、その外乱は推定されたものではない。

日本国特開平３－１９６３１３号公報 日本国特開平１０－２５４５０４号公報 日本国特開平１０－３２８８０１号公報 日本国特開平１１－３１２００３号公報 日本国特開２００２－３２８７０２号公報 日本国特開２００５－１１５４８８号公報 日本国特開２００８－９６８２号公報 日本国特開２００８－２１７２７９号公報 日本国特開２００８－２５７７４１号公報 日本国特開２０１０－２１８００７号公報 日本国特開２０１４－１１１２６９号公報 日本国特開２０１０－２６６９６７号公報

金子修、木村和暉、林祐樹、山本茂、「周期外乱の影響をうける線形時不変システムの同定と外乱の推定」、電気学会論文誌Ｃ、一般社団法人電気学会 、2014年7月1日、第134巻、第7号、P.917-923 杉本謙二、宮浦恭弘、「独立成分分析によるパラメータと外乱の同時推定」、第４４回自動制御連合講演会前刷、2001年11月22日、P.610-611 岡部弘佑、小黒龍一、蘭林峰、「オブザーバを用いたイナーシャ誤差及び摩擦外乱の同時推定に関する研究」、電気学会産業応用部門大会(CD-ROM)、2011年9月6日、論文番号2-32 増田士朗、武田郷平、「外乱によって生成された入出力データを用いた外乱抑制ＦＲＩＴ法」、電気学会論文誌Ｃ、一般社団法人電気学会、2011年11月、第131巻、第4号、P.788-793 システム制御情報学会編、須田信英著者代表「ＰＩＤ制御」朝倉書店、1992年7月20日、P.144

　ところで、制御系の評価や設計等においては、以下に述べる問題が生じていた。
　例えば、ハンチングが生じたり、定常偏差が生じたり、収束が遅かったりするなどの制御異常が発生している場合、異常の原因が制御系に起因するものか、又は外乱に起因するものかを知る必要がある。なぜならば、異常の原因によって、講じる対策が異なるためである。例えば、制御系起因ならゲインを下げる必要があり、外乱起因ならゲインを上げる必要がある。

　アナログ機器や、古い設備などでは、制御系のパラメータの値が不確かな場合がある。すなわち、制御系のパラメータの値が高精度ではないことがある。不確かなパラメータの場合、制御モデル自体が得られなかったり、不精確となったりして、外乱の推定が不可能、または、不精確となってしまう虞がある。それ故、その制御モデルに基づいた制御系の設計も不可能、または、不精確となってしまう。

　また、制御系を最初に適用させた時などは、得られるデータの測定時間が短いため、測定されるデータの量が少ないという場合がある。
　さらに、制御系では、コントローラの目標値（ＳＶ）、制御量（ＰＶ）、操作量（ＭＶ）は、操業データとして保存されていることが多いが、その他の値は保存されていないことが多い。

　例えば、制御対象が、ゲイン、一次遅れ、積分器、ローパスフィルタ、むだ時間の要素で順に構成されている場合、例えば一次遅れと積分器の間の値は保存されていない。このため、各要素については、入出力データが得られていないことになり、各要素を個別に推定（あるいは同定）することはできない。そこで、限定された入出力データから推定するようにしておくことが望ましい。

　ここで、同定とは、［相良節夫、秋月影雄、中溝高好、片山徹、「システム同定」社団法人計測自動制御学会、1987年11月10日発行、第3版、P.4］に記載されているように、「実システムの入出力データに基づいて、与えられたモデルのパラメータ、あるいはシステムの周波数応答、インパルス応答を決定すること」をいうものとする。例えば、制御モデルの構造が決定しており、その制御モデルのパラメータを同定する場合は、パラメータ推定と同じとなる。

　また、推定とは、「未確定の値を決定すること」をいうものとする。
　そして、パラメータ推定に関し、ステップ信号、Ｍ系列信号、インパルス信号などの信号を、対象に加えることは、生産設備等においては製品の品質維持等の観点から、難しい場合がある。また、不安定系など、開ループでのパラメータ推定（あるいは同定）は困難である。

　それ故、操業上、閉ループでのパラメータ推定が望ましい。さらに、通常の操業条件でのパラメータ推定がより望ましい。
　操業の初期状態は、定常状態になるまでの前であったり、外乱が印加した状態であったりするので、得られた操業データ（制御量（ＰＶ）、操作量（ＭＶ）等）の初期値は、定常値とは異なる値となっている。これらの初期値も一致させなければ、制御モデルを用いたシミュレーションで、実績値を再現することはできない。

　また、制御モデルを、定常値からの偏差系で構築した場合、操作量（ＭＶ）の定常値が不明となっている。しかし、操作量（ＭＶ）の定常値を知っておかなければ、制御モデルを用いたシミュレーションで、実績値を再現することはできない。
　ここで、偏差系とは、［長縄明大、矢田晴義、「最適サーボ設計における重み行列決定の一手法」電気学会論文誌Ｃ、１２０巻、１２号、平成１２年、P.2101］に記載されているように、「定常時との偏差を新たに変数とした系」をいうものとする。

　例えば、線形時不変系の場合、操作量の定常値との偏差を新たな操作量に、状態量の定常値との差を新たな状態量に、制御量の目標値との差を新たな目標値とすることとなる。なお、状態量の定義については、後述することとする。
　ところで、制御対象は、線形系とは限らない。例えば、ゲインが変化する場合や、切り換えを含む場合等、非線形な場合がある。また、シミュレーションが可能であっても、線形系に関する推定（あるいは同定）理論が使えない場合もある。

　また、制御系を設計する場合、制御モデルのパラメータが不精確であれば、得られるコントローラが不適切になる場合がある。また、例えば、外乱の大きさを評価せずに、ハイゲインコントローラを設計してしまうと、操作量の飽和が生じることとなり、制御の不安定化につながる虞がある。
　そこで、本発明は、上記のような課題に鑑みてなされたものであり、外乱とパラメータとを同時且つより精確に推定することができる制御モデルのパラメータと外乱との同時推定方法、及び、この同時推定方法で推定された結果を用いて制御系の制御を行うことのできる制御方法を提供することを目的とする。

　上記課題を解決するため、本発明においては以下の技術的手段を講じている。
　即ち、本発明の制御モデルのパラメータと外乱との同時推定方法は、制御モデル内のパラメータと当該制御モデルが表す制御対象へ入力される外乱とを同時に推定する方法であって、前記外乱を時系列変数値として、前記制御モデルに与えることとし、前記制御対象における実績データの時系列値と、前記制御モデルをシミュレーションすることで得られたシミュレーション結果のデータの時系列値を含む評価関数が、最大又は最小となるように、前記外乱と前記パラメータとを同時に推定することを特徴とする。

　本発明によれば、シミュレーション可能な制御モデルを用いることで、外乱と制御モデル内のパラメータを同時に推定することが可能である。そしてシミュレーション結果のデータと実績データとが一致することとなり、精確な外乱を得ることが可能である。また、本発明は、制御対象が線形系でも、非線形でも、分布定数系でも、切り換え系（例えば、炉頂圧制御）でも、外乱とパラメータを推定することが可能である。本発明は、同定用の特別な信号を用いずに、外乱とパラメータを推定することが可能である。

　本発明によれば、制御異常の原因が、制御起因か、若しくは、外乱起因かを判断することができる。また、本発明は、外乱を遅れなく推定することができる。なお、本発明はオフライン計算であるので、外乱推定オブザーバのようなフィルタは不要である。また、制御モデルに加わる外乱は、特定の関数やパターンに限られないのではあるが、本発明を用いることで、任意の時系列パターンを持つ外乱を推定することができる。

　本発明においては、パラメータ推定のために、外乱を与える必要はないので、容易に外乱とパラメータを推定することが可能である。本発明によれば、工場内の設備立上げ後のような、実績データが数少ない場合や、実績データが非定常データであった場合においても、パラメータと外乱を推定することが可能である。
　本発明においては、得られる実績データの項目が、目標値（ＳＶ）、制御量（ＰＶ）、操作量（ＭＶ）のように限られていても、パラメータと外乱を推定することが可能である。本発明は、設計支援、運転管理支援、シナリオシミュレーション、ソフトウエアセンサー・オブザーバ、モデルベース制御、短期のプロセス最適化、長期のプロセス最適化、教育／運転訓練等に用いることが可能である。本発明によれば、独立成分分析を適用するための条件を満足することなく、パラメータと外乱を推定することが可能である。

　本発明の制御モデルのパラメータと外乱との同時推定方法は、制御モデル内のパラメータと当該制御モデルが表す制御対象へ入力される外乱とを同時に推定する方法であって、前記外乱を時系列変数値として、前記制御モデルに与えることとし、前記制御対象における実績データの時系列値と、前記制御モデルをシミュレーションすることで得られたシミュレーション結果のデータの時系列値と、前記外乱の時系列変数値とを含む評価関数が、最大又は最小となるように、前記外乱と前記パラメータとを同時に推定することを特徴とする。

　本発明によれば、外乱の時系列変数値が含まれていないときに、外乱の誤差とパラメータの誤差とが大きくなる場合があるが、この構成を有する本発明によれば、誤差が大きくなることを防止することができる。
　好ましくは、前記評価関数には、実績データおよびシミュレーション結果のデータとして、制御量（ＰＶ）、操作量（ＭＶ）の少なくとも一つが与えられているとよい。

　この構成を有する本発明によれば、保存されている実績データが目標値（ＳＶ）、制御量（ＰＶ）、操作量（ＭＶ）といった数少ない場合でも、パラメータと外乱を推定することが可能である。
　好ましくは、前記外乱に周波数成分が存在する場合には、当該周波数成分をカットするフィルタを適用した後の外乱が前記評価関数へ含まれるとよい。

　本発明のフィルタに、ローパスフィルタやバンドストップフィルタを適用することで、外乱を適切に推定することができる。なお、パラメータと外乱の組み合わせは無数にある場合があるが、本発明を用いることで、例えば、制御の変動がすべて外乱により生じているという解を算出することを避けることができる。また、外乱が高周波に存在する場合、本発明を用いることで、外乱の定常偏差を小さくすることができる。なお、外乱の定常偏差とは、推定された外乱の時系列値の平均値のことをいうものとする。

　好ましくは、前記実績データの時系列値が複数群、存在する場合においては、前記パラメータを前記実績データの全ての群にて共通する値をもつパラメータと、前記実績データの各群に連動する値をもつパラメータとに分離した上で、シミュレーション結果のデータの時系列値を評価関数へ導入すると共に、前記外乱を実績データの各群に連動する時系列変数値として評価関数へ導入するとよい。

　この構成を有する本発明によれば、状態量の変動に含まれる周波数成分が変化する等により、最適解が得られやすく、精確に推定することが可能である。なお、本発明における線形系の同定では、ＰＥ性（Persistently Exciting）が必要ということと類似している。
　好ましくは、前記パラメータとして、未知の初期値が採用されるとよい。

　この構成を有する本発明によれば、定常状態からシミュレーションを始める場合は勿論のこと、非定常な状態が初期状態の場合においても、初期値を推定することができる。
　好ましくは、前記パラメータとして、未知の定常値が採用されるとよい。
　この構成を有する本発明によれば、制御モデルが偏差系の場合であっても、推定することができる。

　好ましくは、前記外乱がランプ状である場合には、前記外乱の微分値を評価関数へ導入すると共に、前記外乱の微分値を積分した上で、当該積分値を前記制御モデルに与えるとよい。
　この構成を有する本発明によれば、制御モデルに加わる外乱がランプ状の場合であっても、より精確に推定することができる。

　本発明の制御対象の制御方法は、上記の制御モデルのパラメータと外乱との同時推定方法で推定されたパラメータと外乱とのうち少なくともパラメータを用いて、前記制御対象を制御することを特徴とする。
　本発明によれば、制御対象を制御する際に、精確なパラメータや外乱を用いているので、制御性能を向上させることが可能である。すなわち、本発明においては、フィードバック補償、フィードフォワード補償、ゲインテーブル等を制御モデルや外乱に基づいて適正化しているので、制御性能を向上させることが可能である。

　本発明によれば、制御対象に加わる外乱を定量的に見積もることができるため、操作量（ＭＶ）が飽和しないような制御系を設計することができる。

　本発明によれば、外乱とパラメータとを同時、且つより精確に推定することができる。

連続系の制御モデルの一例を示した図である。 離散系の制御モデルの一例を示した図である。 本発明にかかる制御モデルのパラメータと外乱との同時推定方法のフローチャート図である。 外乱の一例を示した図である。 高炉の炉頂圧発電設備（高炉ガスのフロー）の概略を示した図である。 ＴＲＴの内部構造を示した断面模式図である。 バイパス弁の内部構造を示した断面模式図である。 ＲＳＥの内部構造を示した断面模式図である。 ＴＲＴのＣＶ値の一例を示した図である。 バイパス弁のＣＶ値の一例を示した図である。 ＲＳＥのＣＶ値の一例を示した図である。 ＴＲＴ緊急停止時の圧力変化と開度変化を示した図である。 本発明のバイパス弁の制御のブロック図である。 ＴＲＴ翼弁開度からフィードフォワード補償量を求めるフローチャート図である。 ＴＲＴ緊急停止時の圧力変化と開度変化を示した図である。 炉頂圧および前圧の制御モデルの概略を示した図である。 外乱に適用するフィルタの特性の一例を示した図である。 推定された外乱とシミュレーション結果であり、時間（ｓ）に対する炉頂圧の制御量（ＰＶ）の関係を示した図である。 推定された外乱とシミュレーション結果であり、時間（ｓ）に対する第１バイパス弁の開度の関係を示した図である。 推定された外乱とシミュレーション結果であり、時間（ｓ）に対する前圧の制御量（ＰＶ）の関係を示した図である。 推定された外乱とシミュレーション結果であり、時間（ｓ）に対する第２バイパス弁の操作量（ＭＶ）の関係を示した図である。 推定された外乱とシミュレーション結果であり、時間（ｓ）に対する炉頂圧の操作量（ＭＶ）の関係を示した図である。 推定された外乱とシミュレーション結果であり、時間（ｓ）に対する高炉ガス流量の外乱の関係を示した図である。 推定された外乱とシミュレーション結果であり、時間（ｓ）に対するＴＲＴ前圧の操作量（ＭＶ）の関係を示した図である。 推定結果に基づきフィードフォワード補償を追加した場合のシミュレーション結果であり、時間（ｓ）に対する炉頂圧の制御量（ＰＶ）の関係を示した図である。 推定結果に基づきフィードフォワード補償を追加した場合のシミュレーション結果であり、時間（ｓ）に対する第１バイパス弁の開度の関係を示した図である。 推定結果に基づきフィードフォワード補償を追加した場合のシミュレーション結果であり、時間（ｓ）に対する前圧の制御量（ＰＶ）の関係を示した図である。 推定結果に基づきフィードフォワード補償を追加した場合のシミュレーション結果であり、時間（ｓ）に対する第２バイパス弁の操作量（ＭＶ）の関係を示した図である。 推定結果に基づきフィードフォワード補償を追加した場合のシミュレーション結果であり、時間（ｓ）に対する炉頂圧の操作量（ＭＶ）の関係を示した図である。 推定結果に基づきフィードフォワード補償を追加した場合のシミュレーション結果であり、時間（ｓ）に対する高炉ガス流量の外乱の関係を示した図である。 推定結果に基づきフィードフォワード補償を追加した場合のシミュレーション結果であり、時間（ｓ）に対するＴＲＴ前圧の操作量（ＭＶ）の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、時間（ｓ）に対する炉頂圧の制御量（ＰＶ）の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、時間（ｓ）に対する第１バイパス弁の開度の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、時間（ｓ）に対する前圧の制御量（ＰＶ）の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、時間（ｓ）に対する第２バイパス弁の操作量（ＭＶ）の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、時間（ｓ）に対する炉頂圧の操作量（ＭＶ）の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、時間（ｓ）に対する高炉ガス流量の外乱の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、時間（ｓ）に対するＴＲＴ前圧の操作量（ＭＶ）の関係を示した図である。 連続鋳造機と湯面レベルの制御系の概略を示した図である。 湯面レベルの制御系のブロック図である。 制御テスト時における湯面レベルのオフセットを持ったハンチング現象を示した図である。 連続鋳造機の湯面レベル制御系の制御モデルと外乱の概略を示した図である。 外乱に適用するフィルタの特性の一例を示した図である。 ランプ状の外乱の場合における、その外乱の微分値に適用するフィルタの特性の一例を示した図である。 外乱が主にランプ状である場合の処理の仕方の一例を示した図である。 本実施形態の２例目における最適化のフローチャート図である。 推定された外乱とシミュレーション結果を示した図であり、コントローラＡで制御時の、時間（ｓ）に対するバルブ開度の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、コントローラＢで制御時の、時間（ｓ）に対するバルブ開度の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、コントローラＡで制御時の、時間（ｓ）に対する湯面レベルの関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、コントローラＢで制御時の、時間（ｓ）に対する湯面レベルの関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、コントローラＡで制御時の、時間（ｓ）に対する湯面レベルの外乱の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、コントローラＢで制御時の、時間（ｓ）に対する湯面レベルの外乱の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、コントローラＡで制御時の、時間（ｓ）に対するコントローラの外乱の関係を示した図である。 推定された外乱とシミュレーション結果を示した図であり、コントローラＢで制御時の、時間（ｓ）に対するコントローラの外乱の関係を示した図である。 制御系改善後のレベル変動を示した図である。 外乱の評価なし、フィルタなし、外乱の微分値不使用の場合における、ゲイン、一次遅れ、むだ時間、ＰＩ制御からなる系の制御モデルのパラメータと外乱の概略を示した図である。 真値のパラメータと、真値の外乱の場合の実績データを示した図であり、時間（ｓ）に対する目標値（ＳＶ）実績の関係を示した図である。 真値のパラメータと、真値の外乱の場合の実績データを示した図であり、時間（ｓ）に対する制御量（ＰＶ）実績の関係を示した図である。 真値のパラメータと、真値の外乱の場合の実績データを示した図であり、時間（ｓ）に対する操作量（ＭＶ）実績の関係を示した図である。 真値のパラメータと、真値の外乱の場合の実績データを示した図であり、時間（ｓ）に対する外乱実績の関係を示した図である。 外乱の評価なし、フィルタなし、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する目標値（ＳＶ）の関係を示した図である。 外乱の評価なし、フィルタなし、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する制御量（ＰＶ）の関係を示した図である。 外乱の評価なし、フィルタなし、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する操作量（ＭＶ）の関係を示した図である。 外乱の評価なし、フィルタなし、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する外乱の関係を示した図である。 外乱の評価あり、フィルタなし、外乱の微分値不使用の場合における、ゲイン、一次遅れ、むだ時間、ＰＩ制御からなる系の制御モデルのパラメータと外乱の概略を示した図である。 外乱の評価あり、フィルタなし、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する目標値（ＳＶ）の関係を示した図である。 外乱の評価あり、フィルタなし、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する制御量（ＰＶ）の関係を示した図である。 外乱の評価あり、フィルタなし、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する操作量（ＭＶ）の関係を示した図である。 外乱の評価あり、フィルタなし、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する外乱の関係を示した図である。 外乱の評価あり、フィルタあり、外乱の微分値不使用の場合における、ゲイン、一次遅れ、むだ時間、ＰＩ制御からなる系の制御モデルのパラメータと外乱の概略を示した図である。 外乱の評価あり、フィルタあり、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する目標値（ＳＶ）の関係を示した図である。 外乱の評価あり、フィルタあり、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する制御量（ＰＶ）の関係を示した図である。 外乱の評価あり、フィルタあり、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する操作量（ＭＶ）の関係を示した図である。 外乱の評価あり、フィルタあり、外乱の微分値不使用の場合における、実績データとシミュレーション結果を示した図であり、時間（ｓ）に対する外乱の関係を示した図である。

　以下、本発明に係る制御モデルのパラメータと外乱との同時推定方法の実施形態を、図を基に説明する。なお、本発明においては、制御対象が線形系とは限らないので、線形系に限らないこととする。
　まず、制御モデルについて、説明する。
　図１に、対象が連続系の場合の制御モデルを示す。

　ここで、制御モデルとは、自然法則に従う動的な対象をシミュレーション可能なように表したものをいうものとする。
　例えば、後述する具体例のように、高炉１７の炉頂圧の制御系を、状態方程式（ガスの質量を求める時間積分計算を含む。）・弁の流量特性・制御の切り換え機能等で表したり、連続鋳造機３０の湯面レベルの制御系を、積分器・一次遅れ等で表したりしたものである。この制御モデルは、以下に示す連続系で表現されたり、後述する離散系で表現されたりする。また、制御モデルに、微分方程式だけでなく、微分代数方程式を含んでもよい。さらに、制御モデルは、集中定数系だけでなく、分布定数系でもよい。

　このような制御モデルは、例えば、公用ソフトウェアであるThe MathWorks,Inc.製のSimulink（Simulinkは、The MathWorks,Inc.の登録商標）を用いて作成することができる。そして、連続系の制御モデルの一例を、式（１）～式（４）に示す。

　ここで、各変数や各関数は次の通りである。
　ｔは、時刻を示す。
　ｘ（ｔ）は、時系列関数のｍ_ｘ行１列のベクトルである。成分は、シミュレーション可能な制御モデルの中のすべての状態量の関数を示す。状態量とは、任意の時点でシステム全体の状態を表せるシステム変数群の最小の部分集合である。例えば線形系の場合は、積分器の出力信号となる。また、状態変数と呼ばれることもある。

　これら状態量は、自動制御の分野では、制御対象の状態量を表すことが多いが、本発明では、制御モデルの構成要素すべての状態量を表すこととする。すなわち、本発明においては、制御対象の状態量に加えて、必要に応じて、コントローラ７～９，４０の状態量や、最適化時に付加するフィルタや積分器等の状態量を含むこととする。
　ｘ０は、ｍ_ｘ行１列のベクトルであり、ｘ（ｔ）の初期値である。

　ｒ（ｔ）は、時系列関数のｍ_ｒ行１列のベクトルである。なお、本発明においては、制御モデルに加えられる入力のうち、既知のものを示す。各成分は、例えば、目標値や、既知の外乱、制御モデルに影響を与える既知の値（例えば、高炉１７の炉頂圧制御での定常時の高炉ガス流量、すなわち、公称高炉ガス流量や、連続鋳造機３０の鋳造速度など）等の関数である。

　ｄ（ｔ）は、時系列関数のｍ_ｄ行１列のベクトルである。各成分は、制御モデルに加えられる入力のうち、未知のもの、すなわち、外乱を示す変数の関数となる。
　ｐは、ｍ_ｐ行１列のベクトルであり、調整すべきパラメータである。本発明においては、調整の必要に応じて、制御モデルに含まれるゲイン、時定数、むだ時間、定常値、状態量の初期値等を含むこととする。なお、コントローラ７～９，４０のゲインが未知のときは、このゲインを含めることもできる。また、ｐを、時系列の値とすれば、時間とともに変化するパラメータとすることもできる。さらに、偏差に応じた折れ線ゲインをパラメータとすること等も可能である。

　ｙ（ｔ）は、時系列関数のｍ_ｙ行１列のベクトルである。各成分は、時系列の実績データが得られ、かつ、時系列のシミュレーション結果のデータと一致させる変数の関数である。変数として、例えばプロセス制御の場合、操作量、制御量が挙げられる。また、他にも、時系列の実績データが得られ、かつ、時系列のシミュレーション結果のデータと一致させたい変数があれば、加えればよい。

　ｄｄ（ｔ）は、時系列関数のｍ_ｄ行１列のベクトルである。各成分は、外乱ｄ（ｔ）に、フィルタや積分等の操作を行った後の変数の関数である。なお、フィルタ等の操作が不要な成分については、必要に応じて補間等を行い、ｄ（ｔ）をそのままｄｄ（ｔ）として用いればよい。
　ｆｃ（）は、連続系において、状態量を計算する制御モデルである。例えば、線形系の場合には、３つの要素から構成され、その３つの要素はそれぞれ、微分方程式（微分代数方程式を含む）、出力方程式と、外乱に対してフィルタや積分等の操作する式となる。

　ｃｃ_ｉｎｅｑ（）は、ｍ_{ｃｃｉｎｅｑ}行１列の関数のベクトルであり、不等式制約式である。
　ｃｃ_ｅｑ（）は、ｍ_ｃｃｅｑ行１列の関数のベクトルであり、等式制約式である。
　なお、上記の制御モデルは、ｒ（ｔ）、ｄ（ｔ）を入力とし、ｙ（ｔ）、ｄｄ（ｔ）を出力とするシミュレーション可能なモデルであれば、特に限定はしない。

　図２に、離散系の場合の制御モデルを示す。そして、離散系の制御モデルの一例を、式（５）～式（８）に示す。

　ここで、各変数や各関数は次の通りである。
　ｉ（ｉ＝０,…,Ｌ－１）は、実績データのサンプリング時間毎のカウントを示す。本実施形態では、シミュレーション結果を実績データと同じ時間間隔で保存するものとする。
　ｊ（ｊ＝０,…,Ｍ－１）は、外乱のサンプリング時間毎のカウントを示す。なお、（ｉ＝０）と（ｊ＝０）は、初期時刻に相当し、（ｉ＝Ｌ－１）と（ｊ＝Ｍ－１）は、最終時刻に相当するものとする。

　ｘ_ｉは、ｍ_ｘ行１列のベクトル（ｉ＝０,…,Ｌ－１）であり、成分は、カウントｉでの、シミュレーション可能な制御モデルの中のすべての状態量である。状態量は、任意の時点でシステム全体の状態を表せるシステム変数群の最小の部分集合である。例えば、線形系の場合は、サンプリング周期の遅れを示す遅延要素ｚ^―１の出力信号となる。
　これら状態量は、自動制御の分野では、制御対象の状態量を表すことが多いが、本実施形態では、制御モデルの構成要素すべての状態量を表すこととする。すなわち、本実施形態においては、制御対象の状態量に加えて、必要に応じて、コントローラ７～９，４０の状態量や、最適化時に付加するフィルタや積分器等の状態量を含むこととする。

　なお、連続系の場合のｍ_ｘと、離散系の場合のｍ_ｘとは、離散化方法の違い等により同じ値でない場合がある。
　ｘ０は、ｍ_ｘ行１列のベクトルであり、ｘ_ｉの初期値である。
　ｒ_ｉは、ｍ_ｒ行１列のベクトル（ｉ＝０,…,Ｌ－１）であり、制御モデルに加えられる入力のうち、既知のものを示す。各成分は、カウントｉでの、例えば、目標値や、既知の外乱、制御モデルに影響を与える既知の値（例えば、高炉１７の炉頂圧制御での定常時の高炉ガス流量、すなわち、公称高炉ガス流量や、連続鋳造機３０の鋳造速度。）等である。

　ｄ_ｊは、ｍ_ｄ行１列のベクトル（ｊ＝０,…,Ｍ－１）である。各成分は、カウントｉでの、制御モデルに加えられる入力のうち、未知のもの、すなわち、外乱を示す。
　ｐは、ｍ_ｐ行１列のベクトルであり、調整すべきパラメータである。本発明においては、調整の必要に応じて、制御モデルに含まれるゲイン、時定数、むだ時間、定常値、状態量の初期値等を含むこととする。なお、コントローラ７～９，４０のゲインが未知のときは、このゲインを含めることもできる。また、ｐを、時系列の値とすれば、時間とともに変化するパラメータとすることもできる。さらに、偏差に応じた折れ線ゲインをパラメータとすること等も可能である。

　ｙ_ｉは、ｍ_ｙ行１列のベクトル（ｉ＝０,…,Ｌ－１）である。各成分は、カウントｉでの、時系列の実績データと、時系列のシミュレーション結果のデータとを一致させる変数である。例えば、プロセス制御の場合、操作量、制御量が挙げられる。また、他にも、時系列の実績データが得られ、かつ、時系列のシミュレーション結果のデータと一致させたい変数があれば、加えればよい。

　ｄｄ_ｉは、ｍ_ｄ行１列のベクトル（ｉ＝０,…,Ｌ－１）である。各成分は、カウントｉでの、外乱ｄ_ｊに、必要に応じて、フィルタや積分等の操作を行った後の変数である。なお、フィルタ等の操作が不要な成分については、必要に応じて、時刻を合わせる補間等を行い、ｄ_ｊの対応する成分を、ｄｄ_ｉの成分として用いればよい。
　ｆｄ（）は、離散系において、状態量を計算する制御モデルであり、上記した連続系のｆｃ（）に対応している。

　ｃｄ_ｉｎｅｑ（）は、ｍ_{ｃｄｉｎｅｑ}行１列の関数のベクトルであり、不等式制約式である。
　ｃｄ_ｅｑ（）は、ｍ_ｃｄｅｑ行１列の関数のベクトルであり、等式制約式である。
　なお、上記した制御モデルは、ｒ_ｉ、ｄ_ｊを入力とし、ｙ_ｉ、ｄｄ_ｉを出力とするシミュレーション可能なモデルであれば、特に限定しない。

　図３に、本発明にかかる制御モデルのパラメータと外乱との同時推定方法のフローチャート図を示す。
　後の説明のため、離散系の制御モデルの場合を例に挙げてパラメータと外乱の最適化を説明するが、連続系の制御モデルの場合も同様の処理過程を有する。
　なお、連続系の制御モデルから離散系の制御モデルへの変換は、公知の手法、例えば、双一次変換や０次ホールドを用いて変換することができる。

　まず、外乱ｄ_ｊの初期値と、パラメータｐの初期値とを設定する。
　ここで、図４に示すように、外乱ｄ_ｊ（ｊ＝０,…,Ｍ－１）は、時系列のデータの各点を変数とした時系列変数値である。すなわち、外乱ｄ_ｊ（ｊ＝０,…,Ｍ－１）の値が、推定すべき変数である。なお、（ｊ＝０）が初期時刻、（ｊ＝Ｍ－１）が、最終時刻に対応する。

　この場合、外乱の基底関数として、ｊのとき１であり、ｊでないときは０であるＭ個の関数を選んでいることとなる。他に、外乱ｄ_ｊを表す方法として、離散フーリエ変換や離散ウェーブレット変換のような直交する関数を基底として選んでもよい。また、直交していなくても、一次独立な関数を、基底として選んでもよい。
　次に、ｉ＝０（ｊ＝０と同時刻、初期時刻に相当）から、ｉ＝Ｌ－１（ｊ＝Ｍ－１と同時刻、最終時刻に相当）までのシミュレーションを、式（５）～式（８）を用いて行う。シミュレーションには、例えば、公用ソフトウェアであるThe MathWorks,Inc.製のSimulink（Simulinkは、The MathWorks,Inc.の登録商標）を用いて行うことができる。

　続いて、評価関数を計算する。

　ここで、各変数は次の通りである。
　ｙ_{ａｃｔｕａｌ ｉ}は、ｍ_ｙ行１列のベクトル（ｉ＝０,…,Ｌ－１）である。各成分は、時系列の実績データと、時系列のシミュレーション結果のデータとを一致させる場合の、カウントｉでの、実績データの変数である。なお、ｙ_ｉが、シミュレーション結果のデータの変数となる。

　Ｗ_ｙは、ｍ_ｙ行ｍ_ｙ列の行列であり、時系列の実績データと、時系列のシミュレーション結果のデータとを一致させる場合の重みを示す。本実施形態においては、対称な準正定行列とする。
　Ｗ_ｄｄは、ｍ_ｄ行ｍ_ｄ列の行列であり、外乱にフィルタや積分等の操作を行った後の変数の重みを示す。本実施形態においては、対称な準正定行列とする。

　次に、評価関数Ｊが収束したか否かを判断する。収束の判断の一例は、評価関数Ｊが、ｄ_ｊ、ｐを変更する１ループの間の変化量が所定の値より小さくなった場合である。
　ところで、最適化計算は、公知の方法、例えば逐次二次計画法を用いて行うこととしている。例えば、公用ソフトウェア、The MathWorks,Inc.製のMatlab（Matlabは、The MathWorks,Inc.の登録商標）のOptimization Toolbox（Optimization Toolboxは、The MathWorks,Inc.の登録商標）を用いて、最適化計算を行うとよい。

　そして、最適化計算の結果、評価関数Ｊが収束していなければ、ｄ_ｊ、ｐを変更して、再び、シミュレーションと評価関数Ｊの計算を繰り返す。
　最適化計算の結果、評価関数Ｊが収束していれば、得られたｄ_ｊを外乱の最適値ｄ_{ｊ ｏｐｔ}とし、ｐをパラメータの最適値ｐ_ｏｐｔとして、計算を終了する。
　なお、ローカルミニマムに陥る場合には、外乱ｄ_ｊの初期値と、パラメータｐの初期値との組み合わせを複数個設定し、図３のフローを実行し、評価関数の値が最も小さい場合の外乱の最適値ｄ_{ｊ ｏｐｔ}とパラメータの最適値ｐ_ｏｐｔとを採用すればよい。

　なお、図３は、公知の逐次二次計画法等を想定したフローチャート図としているが、他の最適化手法、例えば、公知のＧＡ（Genetic algorithm）、ＰＳＯ（Particle Swarm Optimization）でもよい。
　さらに、得られた外乱の最適値ｄ_{ｊ ｏｐｔ}とパラメータの最適値ｐ_ｏｐｔを利用して、制御系の設計が可能である。

　例えば、上記の特許文献１２や非特許文献５の手法に、本実施形態の制御モデルのパラメータと外乱との同時推定方法で得られた結果を用いることで、ＰＩＤ制御の制御パラメータを調整することができる。なお、ＰＩＤ制御以外の制御器であっても、同様に、本実施形態の制御モデルのパラメータと外乱との同時推定方法で得られた結果を用いることで、制御パラメータを調整することができる。

　また、The MathWorks,Inc.製のMatlab（Matlabは、The MathWorks, Inc.の登録商標）のSimulink Control Design（Simulink Control Designは、The MathWorks,Inc.の登録商標）等の公用ソフトウェアに、本実施形態の制御モデルのパラメータと外乱との同時推定方法を適用させることで、制御パラメータの調整が可能である。
　また、本実施形態においては、精確なモデルが得られているため、例えば、公知の方法を用いてフィードフォワード補償が可能となる。

　まとめると、本発明は、制御対象へ入力される外乱を時系列変数値として、制御モデルに与えることとし、制御対象における実績データの時系列値と、制御モデルをシミュレーションすることで得られたシミュレーション結果のデータの時系列値を含む評価関数が、最大又は最小となるように、外乱と制御モデルのパラメータとを同時に推定する方法である。

　また、本発明は、制御対象へ入力される外乱を時系列変数値として、制御モデルに与えることとし、制御対象における実績データの時系列値と、制御モデルをシミュレーションすることで得られたシミュレーション結果のデータの時系列値と、外乱の時系列変数値とを含む評価関数が、最大又は最小となるように、外乱と制御モデルのパラメータとを同時に推定する方法である。

　ここで、上記の評価関数には、実績データおよびシミュレーション結果のデータとして、制御量（ＰＶ）、操作量（ＭＶ）の少なくとも一つが与えられている。また、上記の外乱に周波数成分が存在する場合には、当該周波数成分をカットするフィルタを適用した後の外乱が評価関数へ含まれることとしている。
　さらに、上記の実績データの時系列値が複数群、存在する場合においては、パラメータを実績データの全ての群にて共通する値をもつパラメータと、実績データの各群に連動する値をもつパラメータとに分離した上で、評価関数へ導入すると共に、外乱を実績データの各群に連動する時系列変数値として評価関数へ導入することとしている。

　また、上記のパラメータとして、未知の初期値が採用されている。さらに、上記のパラメータとして、未知の定常値が採用されている。
　また、上記の外乱がランプ状である場合には、外乱の微分値を評価関数へ導入すると共に、外乱の微分値を積分した上で、当該積分値を制御モデルに与えることとしている。
　以上述べた方法に従って、制御モデルのパラメータ及び、制御対象へ入力される外乱を推定することで、同時、且つより精確に推定することができる。

　本発明にかかる制御モデルのパラメータと外乱との同時推定方法で推定されたパラメータと外乱とのうち少なくともパラメータを用いることで、制御対象を制御することができる。
　続いて、本実施形態の制御モデルのパラメータと外乱との同時推定方法を用いて制御対象を制御する方法を、３つの具体例を挙げて説明する。
［１例目］
　１例目は、鉄鋼プロセスの１つである高炉１７の炉頂圧発電設備１におけるバイパス弁６ａ，６ｂの操作方法について、外乱とパラメータとを同時推定し、フィードフォワード補償量を決定する方法について説明する。

　まず、高炉１７の炉頂圧発電設備１の概略について、図を基に説明する。
　図５に示すように、高炉１７の炉頂圧発電設備１は、高炉１７の炉頂で排出されたガス（高炉ガス）を回収し、回収した高炉ガスを用いて発電を行うもので、高炉ガスを回収する高炉ガス流路４を備えている。
　高炉ガス流路４には、高炉ガスに含まれるダスト等を除去する集塵設備が接続されている。この集塵設備は、ダストキャッチャー（図示省略）やリングスリットウォッシャー（ＲＳＷ）であり、これらの内部に高炉ガスを通すことにより、ダスト等を除去することができる。このうち、リングスリットウォッシャー（ＲＳＷ）は、複数（例えば、３つ）のリングスリットエレメント２（ＲＳＥ）を並列に接続することにより構成されている。

　また、ＲＳＥ２の下流側の高炉ガス流路４には、ミストセパレータ（図示省略）が設けられ、少なくとも２つの流路に枝分かれしている。１つ目の流路は、高炉ガスにより発電を行う発電部３（ＴＲＴ）を通過するメイン流路であり、２つ目の流路は、ＴＲＴ３のガス入側とガス出側との間をメイン流路とは別に結ぶバイパス流路である。このバイバス流路５には、並列して、２つのバイパス弁６ａ，６ｂが設けられている。バイパス流路であって、バイパス弁６ａ，６ｂを設置した下流側には、高炉ガスの消音のためのサイレンサ１８が設けられている。

　次に、上記した各装置（ＲＳＥ２、ＴＲＴ３、バイバス弁）について、詳しく説明する。
　図８に示すように、ＲＳＥ２は、高炉ガスが流れる円錐管１０と、この円錐管１０に挿入された移動自在なコーン状（円錐状）の可動部１１とを備えている。可動部１１は、図５に示す炉頂圧コントローラ７により制御される。具体的には、炉頂圧コントローラ７は、炉頂圧の実績値と、炉頂圧の目標値とに基づき、可動部１１の操作量を求め、この操作量に基づいて可動部１１を移動させる。即ち、炉頂圧コントローラ７は、円錐管１０と可動部１１との間隔を操作し、炉頂圧が目標値になるように、通過する高炉ガス流量を調整する。

　つまり、ＲＳＥ２では、可動部１１と円錐管１０（ケース）とで構成された弁機構１６ａによって、両者間の間隙が変化することで、高炉ガスの流量や圧力を調整する。なお、流入ガスには細かな水滴が含まれており、間隙を通過する際に、流速が増えることにより、細かなダストと水滴が付着して、除塵することができる。
　図６に示すように、ＴＲＴ３（発電部）は、発電機１３に直結するタービンの軸１４（ロータ）に設けられた動翼１５と、ケーシング１２側に固定された静翼１６とを備えている。

　また、静翼１６の一部は、角度を変更することができる翼弁１６ａ（弁機構）で構成されており、図５に示すＴＲＴ前圧コントローラ８によって翼の角度（弁の角度）が変更できるようになっている。具体的には、ＴＲＴ前圧コントローラ８は、前圧の実績値と前圧の目標値とに基づき、翼弁１６ａ等の操作量を求め、この操作量に基づいて翼弁１６ａの角度を変更し、これにより、前圧が目標値になるように、ＴＲＴ３を通過するガス流量を調整する。

　つまり、流入ガスである高炉ガスがタービンに吸気され、静翼１６と動翼１５を通過する際に、動翼１５がタービンの軸１４を回転させる。タービンの軸１４は発電機１３に接続されており、電力を発生させる。最前段の静翼１６は、前述のように、角度が変えられる翼弁１６ａとなっており、高炉ガスの流量や圧力を調整することができる。
　図７に示すように、バイパス弁６ａ，６ｂは、バタフライ弁と呼ばれる弁であり、円状のディスクを軸芯周りで回動させて角度を変更することにより、当該バイバス弁の開度（弁開度という）を変更することができる。バイパス弁６ａ，６ｂでは、図５に示すバイパス弁前圧コントローラ９が、前圧の実績値と前圧の目標値とに基づき、バイパス弁６ａ，６ｂの操作量を求め、この操作量に基づいて弁開度を変更することによって、バイバス流路５を通過するガス流量を調整する、即ち、前圧を目標値に制御する。

　なお、第１バイパス弁６ａは、フィードフォワード補償（ＦＦ補償）により、弁開度が制御され、第２バイパス弁６ｂは、フィードフォワード補償（ＦＦ補償）とフィードバック補償の（ＦＢ補償）の組み合わせにより、弁開度が制御される。また、図５では、バイパス弁６ａ，６ｂをバイバス流路５に２つ設置した例を示しているが、この例に限定されることなく、バイパス流路に、バイバス弁を１つ、または、３つ以上設置してもよい。また、バイパス弁６ａ，６ｂは、バタフライ弁以外の弁であってもよい。また、ＲＳＥ２等の炉頂圧の制御系が存在しない場合には、ＴＲＴ３やバイパス弁６ａ，６ｂが炉頂圧を制御することとなる。

　以上、炉頂圧発電設備１によれば、高炉１７の炉頂から排出される高炉ガスは、通常はＴＲＴ３を流れ、その後、高炉ガスホルダへと回収されることとなる。ＴＲＴ３では高炉ガスによりタービンが回転し発電が行われる。
　ここで、ＴＲＴ３の保守作業を行うために、ＴＲＴ３に流す高炉ガスを停止する状況が発生する。その場合には、徐々に時間をかけ、高炉ガスをバイパス流路へ流すようにする。このとき、バイパス弁６ａ，６ｂは通常のＦＦ補償及びＦＢ補償によりコントロールされる。

　このように、ＴＲＴ３に高炉ガスを流すことにより、発電を続けることができるが、何らかの事情で、ＴＲＴ３を緊急停止する必要がある。緊急停止する場合は、保守作業の場合とは異なり、即座にＴＲＴ３を停止することから、前圧が急激に上昇して、炉頂圧も大きくなってしまう。炉頂圧の上昇は、高炉１７の安定操業の妨げとなる。
　そのため、本発明では、ＴＲＴ３を緊急停止した場合でも、前圧、即ち、炉頂圧の上昇を抑えて、緊急停止による炉頂圧の変動が発生しない対策を講じている。詳しくは、本発明では、緊急停止する前の高炉ガスの流量を計算する。そのうえで、緊急停止時には、バイパス弁６ａ，６ｂのフィードフォワード補償による流量が、計算して得られた発電部３内を流れる高炉ガスの流量と等しくなるように、バイパス弁６ａ，６ｂのＣＶ値を用いて、バイパス弁６ａ，６ｂの開度のフィードフォワード補償量を算出し、算出したフィードフォワード補償量をバイパス弁６ａ，６ｂに適用することで、緊急停止時の高炉ガスをバイパス弁６ａ，６ｂ側に逃がすこととしている。

　本発明では、良好な制御精度を得るために、制御モデル内のパラメータと、当該制御モデルが表す制御対象（バイパス弁６ａ，６ｂ）へ入力される外乱とを同時に推定して、制御モデルの精度を向上させて、精確なフィードフォワード補償量を得ることとしている。
　続いて、緊急停止時における処理について、ＴＲＴ緊急停止前に、バイパス弁６ａ，６ｂに流れる流量が０（ゼロ）の場合を詳しく説明する。

　まず、Ｃ_Ｖ値について説明する。
　図９～１１は、それぞれ、ＴＲＴ３、バイパス弁６ａ，６ｂ、ＲＳＥ２のＣ_Ｖ値またはＣ_Ｖ値相当の値の例を示している。ここで、開度は、翼弁１６ａの角度、バタフライ弁の角度、ＲＳＥ２における弁機構１６ａ（可動部１１及び円錐管１０）における間隙（変位量）である。なお、開度は、ＴＲＴ３、バイパス弁６ａ，６ｂ、ＲＳＥ２のそれぞれの角度や変位を正規化し、パーセント表示で表してもよい。

　簡便に説明するために、バイパス弁６ａ，６ｂのそれぞれのＣ_Ｖ値の特性は、同じとする。なお、バイパス弁６ａ，６ｂのそれぞれのＣ_Ｖ値の特性は、異なっていてもよい。
　Ｃ_Ｖ値は、公知文献の「改訂第２版，工業プロセス用調節弁の実技ハンドブック,（株）山武，調節弁ハンドブック編纂委員会，日本工業出版，２０１２年」のＰ．５２に記されているように、容量係数と呼ばれる値であり、比重１の水を弁差圧１［ｐｓｉ］として通過する体積流量を［ＵＳｇａｌ／ｍｉｎ］であらわした数値である。なお、本発明でのＣＶ 値は、定義にかかわらず、少なくとも弁機構（翼弁１６ａ、バイパス弁６ａ，６ｂ、ＲＳＥ２等）の開度と、入側の圧力と、弁機構１６ａを流れる流体の流量との関係を示し、開度から流量を求めることができ、あるいは、流量から開度を求めることができる値であればよい。

　すなわち、本発明でのＣ_Ｖ値は、Ｃ_Ｖ値相当の値を含むものとする。弁機構１６ａの出側の圧力や、流体の温度等は、弁機構１６ａを流れる流体の流量との関係を示す式の中で、用いられても、用いられなくてもよい。以後の説明で、用いられない値の入力は不要となる。
　また、弁機構１６ａの入側の圧力または出側の圧力は、流路を構成する配管やサイレンサ１８を間に含めた場合の圧力を含むこととする。理由は、入側の圧力または出側の圧力の箇所にセンサがなく、配管の上流側や下流側にセンサがあったり、さらに途中の配管に弁やサイレンサ１８が設置されている場合があったりするからである。

　続いて、従来の場合のＴＲＴ緊急停止時の制御の例を説明する。フィードフォワード補償を１つのバイパス弁６にのみ適用する場合である。
　ＴＲＴ３を緊急停止する場合は、当該ＴＲＴ３の上流側に設けられた遮断弁（非表示）が閉じ、同時に、翼弁１６ａの開度が０（ゼロ）まで閉じ、ＴＲＴ３を流れる高炉ガス流量は０となる。翼弁１６ａが閉じると同時に、バイパス弁６ａは、フィードフォワード補償により開く。バイパス弁６ｂは、前圧目標値と前圧実績値が一致するようなフィードバック補償がなされ、ＴＲＴ３に流れていたガスと同流量のガスがバイパス弁６ａ、６ｂに流れるようになるまで、バイパス弁６ｂを開けていくこととなる。

　このとき、バイパス弁６の開度が定常値に達するまでは、バイパス弁６の流量は、元のＴＲＴ３へ流れる流量より少ないこととなる。このため、図１２に示すように、ＲＳＥ２とバイパス弁の間の高炉ガスの質量が増加することとなり、この部分の圧力、すなわち、前圧が一時的に増加する。さらに、前圧（言い換えれば、ＲＳＥ２の背圧）が増加することにより、ＲＳＥ２を流れる流量が一時的に減少し、高炉１７とＲＳＥ２との間の高炉ガスの質量が増加し、炉頂圧も増加することとなる。

　炉頂圧の増加は、高炉ガスの放散等を生じ、安定操業上好ましくない。また、前圧の一時的な増加により除塵能力が低下するため、前圧設定値を高くすること、すなわち、発電量を増加することができない。このように従来では、フィードフォワード補償量が不足しており、フィードフォワード補償量を適切に決定し、バイパス弁に適切にフィードフォワード補償できれば、炉頂圧や前圧の変動を抑制できる。

　次に、本発明おける緊急停止について説明する。
　図１３は、本発明のバイパス弁６ａ，６ｂの制御のブロック図を示している。
　本発明では、従来の手法とは異なり、バイパス弁６ａ，６ｂへのフィードフォワード補償量を、ＴＲＴ３のＣ_Ｖ値と、バイパス弁６ａ，６ｂのＣ_Ｖ値とを用いて定量的に定めている。また、従来ではフィードフォワード補償量を１つのバイパス弁のみに適用していたが、本発明では、複数のバイパス弁６ａ，６ｂに適用している。

　図１４は、フィードフォワード補償量を求めるフローチャートである。
　ＴＲＴ３の緊急停止前は、バイパス弁６ａ，６ｂは閉鎖していて、高炉ガスの全量がＴＲＴ３側に流れる状況である。本発明では、翼弁１６ａの開度（翼弁開度という）を取得し、図９により、翼弁１６ａのＣ_Ｖ値（翼弁Ｃ_Ｖ値という）を計算する。ここが図１４のＳＴＥＰ１である。

　次に、翼弁Ｃ_Ｖ値（Ｃ_ＶＴ）、入側圧力（前圧）（Ｐ_１［－］）、出側圧力（背圧）（Ｐ_２［－］）、温度（Ｔ［℃］）、比重（Ｇ［－］）から、ＴＲＴ３を流れる高炉ガスの流量（Ｑ_Ｔ［－］）を計算する。ここで、［－］は、無次元化したことを示す。ここが図１４のＳＴＥＰ２である。
　具体的には、式（１０）を用いて、高炉ガスの流量（Ｑ_Ｔ［－］）を求める。

　なお、各パラメータＰ_１，Ｐ_２，Ｔ，Ｇ，ａ_Ｔは、センサ等から得られた実績値、シミュレーションから得られた値、バイパス弁６ａ，６ｂのメーカーから提供された値等そのもの、あるいは、無次元化に合わせて単位変換した値である。
　次に、翼弁Ｃ_Ｖ値（Ｃ_ＶＴ）及び高炉ガスの流量（Ｑ_Ｔ）を算出した後、バイパス弁６ａ，６ｂの開度を計算する。ここが図１４のＳＴＥＰ３である。

　さて、バイパス弁６ａ，６ｂの開度の計算方法について説明する。この方法では、一方のバイパス弁（第２バイパス弁６ｂ）の開度を仮決めしたうえで、他方のバイパス弁（第１バイパス弁６ａ）の開度を求めることとしている。例えば、第２バイパス弁６ｂの開度を、４２％とした場合について説明する。この説明では、第１バイパス弁６ａの特性と第２バイパス弁６ｂの特性とは同じであって、図１０に示す特性を有しているとする。

　バイパス弁６ａ，６ｂのＣ_Ｖ値（Ｃ_ＶＢ）、バイパス弁６ａ，６ｂの流量（Ｑ_Ｂ［－］）、バイパス弁６ａ，６ｂの入側圧力（前圧）（Ｐ_１［－］）、バイパス弁６ａ，６ｂに流れる高炉ガスの温度（Ｔ［℃］）について、式（１１）が成り立つ。

　なお、各パラメータＰ_１，Ｔ，Ｇ，ａ_Ｂは、センサ等から得られた実績値、シミュレーションから得られた値、バイパス弁６ａ，６ｂのメーカーから提供された値等そのもの、あるいは、無次元化に合わせて単位変換した値である。
　第２バイパス弁６ｂの開度が４２％であるため、図１０に示す開度及びＣＶ値（Ｃ_ＶＢ）の関係から、第２バイパス弁６ｂのＣ_Ｖ値であるＣ_ＶＢ２を得ることができる。そして、図１０により得られたＣＶＢ２と、Ｔ，Ｐ_１，Ｇ，ａ_Ｂを、式(１１)に代入して、第２バイパス弁６ｂの流量であるＱ_Ｂ２を求めることができる。

　ここで、高炉ガスの流量（Ｑ_Ｔ）については、式（１２）が成立することから、第２バイパス弁６ｂの流量であるＱ_Ｂ２と、予め得られた高炉ガスの流量Ｑ_Ｔとを用いて、第１バイパス弁６ａの流量であるＱ_Ｂ１を求める。

　第１バイパス弁６ａの特性は、第２バイパス弁６ｂと同じ特性としており、Ｑ_Ｂ１，Ｔ，Ｐ_１，Ｇ，ａ_Ｂを式（１１）に適用することにより、第１バイパス弁６ａに対応するＣ_Ｖ値であるＣ_ＶＢ１を得ることができる。式（１１）により得られたＣ_ＶＢ１と、図１０とを用いて、第１バイパス弁６ａの開度を得ることができる。
　なお、この実施形態では、各流量は、標準状態での流量［Ｎｍ^３／ｈ］を無次元化したもの、あるいは、質量の流量［ｋｇ／ｓ］を無次元化したものとしている。このように、緊急停止前の高炉ガス流量と等しくすることによって、同じ状態での流量が等しいことを意味する。

　以上、開度の計算法の例を説明したが、この例に限られるものではない。バイパス弁６ａ，６ｂが２つでなく、３つ以上の場合でも同様に計算することができる。
　なお、２つ以上のバイパス弁６ａ，６ｂを１つのバイパス弁６として扱うスプリットレンジ制御（例えば、改訂第２版 工業プロセス用調節弁の実技ハンドブック,（株）山武，調節弁ハンドブック編纂委員会Ｐ．２３７）を行っている場合には、１つのバイパス弁６とみなした場合のＣＶ値を求めることで、同様に計算することができる。

　さて、バイパス弁６ａ，６ｂ（第１バイパス弁６ａ、第２バイパス弁６ｂ）の開度を求めた後は、図１４に示すように、バイパス弁６ａ，６ｂのフィードフォワード補償量を計算する。ここが図１４のＳＴＥＰ４である。
　第１バイパス弁６ａ、第２バイパス弁６ｂそれぞれについて、ＴＲＴ緊急停止前の開度が０（ゼロ）の場合、「フィードフォワード補償量＝計算された開度」となる。また、ＴＲＴ緊急停止前に、オフセット等により、バイパス弁６ａ，６ｂの開度に０でない開度初期値が存在する場合には、バイパス弁６ａ，６ｂそれぞれについて、「バイパス弁のフィードフォワード補償量＝計算された開度－ＴＲＴ緊急停止前の開度初期値」とすればよい。

　フィードフォワード補償量の設定が終了すると、バイパス弁前圧コントローラ９がフィードフォワード補償量を出力する。ここが図１４のＳＴＥＰ５である。
　以上、上記した実施形態では、発電部３（ＴＲＴ）が、当該発電部３内を流れる高炉ガスの流量（Ｑ_Ｔ）を調整可能とする弁機構１６ａ（翼部）を備えており、発電部３の弁機構１６ａの開度とＣＶ値とを基に、発電部３を緊急停止する前の高炉ガスの流量を計算している。そして、バイパス弁６ａ，６ｂのフィードフォワード補償による流量が、計算して得られた発電部３内を流れる高炉ガスの流量と等しくなるように、バイパス弁６ａ，６ｂのＣ_Ｖ値を用いて、バイパス弁６ａ，６ｂの開度のフィードフォワード補償量を算出し、算出したフィードフォワード補償量をバイパス弁６ａ，６ｂに適用している。

　図１５は、本発明を実施した場合の結果を示しており、理想的な制御状況である。
　本発明では、第２バイパス弁６ｂのフィードフォワード補償量（図１５中の最下段）と、第１バイパス弁６ａのフィードフォワード補償量（図１５中の下から２つ目）が適正化され、ＴＲＴ３側に流れていた高炉ガスの全量が、フィードフォワード補償により速やかに、バイパス流路に流れることになる。このため、図１５に示すように、前圧の変動（図１５中の上から２つ目）が抑制され、さらに、炉頂圧の変動（図１５中の最上段）を抑制することができる。

　ところで、発電部３のＣ_Ｖ値およびバイパス弁６ａ，６ｂのＣ_Ｖ値の精度は、フィードフォワード補償量の精度に影響を与える。この発電部３のＣ_Ｖ値は、配管４，５の影響や、上流側・下流側の設備（サイレンサ１８等）の影響や、炉頂圧発電設備１の経時変化の影響（ガス成分、配管変更など）等を受ける。このため、フィードフォワード補償量を算出する際に用いる発電部３のＣ_Ｖ値、およびバイパス弁６ａ，６ｂのＣ_Ｖ値については、公称の値（メーカー提供の値、設備設置時の計算値等）から補正する必要がある。

　以下、発電部３のＣ_Ｖ値、およびバイパス弁６ａ，６ｂのＣ_Ｖ値の補正方法について説明する。
　図１６に、炉頂圧および前圧の制御モデルを示す。なお、炉頂圧および前圧の制御モデルは、連続系の制御モデルである式（１）、式（２）に相当するが、実際には、公用ソフトウェアであるThe MathWorks,Inc.製のSimulink（Simulinkは、The MathWorks, Inc.の登録商標）を用いてシミュレーションを行うと、離散化されるので、式（５）、式（６）に対応していることとなる。

　高炉－ＲＳＷ間の状態方程式は、次の通りである。なお、高炉ガスの質量は、「炉頂での高炉ガス流量－ＲＳＥでの高炉ガス流量」を時間積分して求める。

　ＲＳＷ－入口間の状態方程式は、次の通りである。なお、高炉ガスの質量は、「ＲＳＥでの高炉ガス流量－ＴＲＴでの高炉ガス流量－第１バイパス弁での高炉ガス流量－第２バイパス弁での高炉ガス流量」を時間積分して求める。

　炉頂圧コントローラ７、ＴＲＴ前圧コントローラ８、バイパス弁前圧コントローラ９に関しては、ＰＩＤコントローラを採用しており、このコントローラ部分の制御モデルは次の通りである。なお、ゲインは既知とする。

　式（１５）においては、炉頂圧コントローラ７からの出力が、炉頂圧ＭＶ、あるいは、ＲＳＥ２の開度に相当する。

　式（１６）においては、ＴＲＴ前圧コントローラ８からの出力が、ＴＲＴ前圧ＭＶ、あるいは、翼弁開度に相当する。

　式（１７）においては、バイパス弁前圧コントローラ９からの出力が、バイパス弁前圧ＭＶ、あるいは、バイパス弁開度に相当する。
　なお、炉頂圧ＳＶは炉頂圧目標値であり、前圧ＳＶは前圧目標値である。
　一方、ＲＳＥ公称流量特性として、図１１に示す公称Ｃ_Ｖ値を用いる。なお、ＲＳＥ２では、図１１の縦軸のＣ_Ｖ値は流量［１／ｓ］そのものを指すものとする。なお、流量は、標準状態での流量、または、質量の流量としているが、前述のように質量は無次元化している。炉頂圧ＭＶを入力すれば、ＲＳＥ２の流量が出力される。

　ＴＲＴ公称流量特性として、図９に示す公称Ｃ_Ｖ値を用いる。ＴＲＴ前圧ＭＶを入力し、Ｃ_Ｖ値を求め、このＣ_Ｖ値と式（１０）とから、ＴＲＴ３の流量が出力される。
　第２バイパス弁流量特性として、図１０に示す公称Ｃ_Ｖ値を用いる。バイパス弁前圧ＭＶを入力し、Ｃ_Ｖ値を求め、このＣ_Ｖ値と式（１１）とから、第２バイパス弁６ｂの流量が出力される。

　第１バイパス弁流量特性として、図１０に示す公称Ｃ_Ｖ値を用いる。ＦＦ補償量を入力し、Ｃ_Ｖ値を求め、このＣ_Ｖ値と式（１１）とから、第１バイパス弁６ａの流量が出力される。
　なお、公称流量は、既知の公称高炉ガス流量である。この公称流量は、高炉ガス流量の定常値とし、一定の値とする。この定常値に対して、原料装入時等に外乱が印加することとなる。

　ところで、本具体例のフィルタは、図１７に示すような低周波域でゲインが大きく、高周波域でゲインが小さくなるようなフィルタを採用している。このフィルタは、外乱が存在する高周波域で、ゲインを小さくするようにしている。本例では、時定数１００[s]の一次遅れのフィルタとしている。折点周波数は、１．５９×１０^－３[Hz]（周期６２８[s]）である。

　このように、推定すべき外乱に上記のフィルタを適用することで、フィルタ適用後の信号で高周波外乱が小さくなり、その小さくなった高周波外乱は後述の評価関数で評価されないこととなる。
　すなわち、上記のフィルタを用いることで、最適化の際に、フィルタ適用前の高周波外乱の値は大きくなることができ、推定可能となる。

　なお、フィルタ適用後においても、低周波外乱の信号は減衰しないので、後述の評価関数で評価されることとなり、フィルタ適用前の低周波外乱の値は、最適化により小さな値となる。
　以上をまとめると、外乱として、主に高周波域の外乱が、大きな値として推定されることとなる。

　一方、Ｋ_ＣＲ，Ｋ_ＣＢ２，Ｋ_Ｒ，Ｋ_Ｔ，Ｋ_Ｂは、推定されるべき、制御モデルのパラメータである。
　なお、Ｋ_ＣＲ，Ｋ_ＣＢ２は、ＰＩＤコントローラ７，９に未知のゲインがある場合を想定している。また、Ｋ_Ｒ，Ｋ_Ｔ，Ｋ_Ｂは、ＲＳＥ２、ＴＲＴ３、バイパス弁６ａ，６ｂの特性に補正係数が必要な場合を想定している。

　ここで、高炉ガス流量外乱は、主に、鉄鉱石やコークスを高炉１７に投入（装入）する際に、高炉ガス流量が減少するという変動を示した外乱である。なお、高炉ガス流量外乱には、高炉１７内のガスの流れの乱れ等により、高炉ガス流量が変化することも含まれる。
　図４に示すように、高炉ガス流量外乱を、時系列の各時刻の値に対応したカウントでの値を変数とする。

　制御モデル内のパラメータの評価関数としては、次の式（１８）を用いた。

　ここで、式（１８）において、添え字の「実績」は既に得られている実績値を意味し、「シミュレーション」は、シミュレーション値を意味する。
　式（１８）の第１項～第５項は、ｙ_ｉについて、シミュレーション値を実績値に一致させるための項であり、具体的には、制御量（ＰＶ）と、操作量（ＭＶ）を、実績とシミュレーションで一致させるための項である。なお、式（１８）の第６項は、フィルタ適用後の外乱を、０（ゼロ）に近づけるための項である。

　制約条件については、次の通りであり、式（７）に対応する。

　なお、各制約条件については、パラメータの取り得る値を考慮して、広めに設定している。
　以上をまとめると、本具体例における炉頂圧制御は、一般形を具体化した形となっている。
　次に、本具体例における炉頂圧制御について、具体的に説明する。

　ｘ_ｉ（状態量）の成分は、高炉－ＲＳＷ間の高炉ガスの質量と、ＲＳＷ－入口間の高炉ガスの質量とを求めるための２つの積分器、各ＰＩＤコントローラ７～９の積分器を離散化した場合の状態量となる。例えば、［１／ｓ＝サンプリング時間／（１－ｚ^－１）］と離散化した場合には、ｍ_ｘ＝５であるが、離散化方法により異なる。なお、Ｌ＝１０１とした。

　ｘ０（初期値）については、この実施例では、変数とせず、固定値とする。炉頂とＲＳＷとの間の高炉ガスの質量の初期値、及びＲＳＷとＴＲＴ３とバイパス弁６ａ，６ｂとの間の高炉ガスの質量の初期値は、時刻が０のときの、圧力、体積、温度が既知であるため、状態方程式から得られる計算値とする。
　また、各ＰＩＤコントローラ７～９のＩゲインの積分器の初期値は、時刻が０のときの操作量（ＭＶ）が得られているため、Ｐゲイン、Ｄゲインを近似的に無視すると、操作量（ＭＶ）から逆算できるため、この値を採用する。例えば、「Ｉゲイン×積分器＝操作量（ＭＶ）」であり、Ｉゲイン、積分器の順に掛け算がなされて操作量（ＭＶ）が得られる場合には、その得られた操作量（ＭＶ）が初期値となる。

　ｒ_ｉ（既知の入力）の成分は、公称高炉ガス流量、炉頂圧の目標値（ＳＶ）、前圧の目標値（ＳＶ）、温度を用いる箇所での高炉ガス温度、ＴＲＴ３・バイパス弁６ａ，６ｂの背圧、（従来制御での）ＦＦ補償量とする。この場合、ｍ_ｒ＝６である。
　ｄ_ｊ（未知の入力）は、推定すべき外乱であり、公称高炉ガス流量に印加する、高炉ガス流量外乱である。この場合、ｍ_ｄ＝１である。なお、Ｍ＝４１とした。

　ｐ（推定すべきパラメータ）は、図１６中のＫ_ＣＲ，Ｋ_ＣＢ２，Ｋ_Ｒ，Ｋ_Ｔ，Ｋ_Ｂである。この場合、ｍ_ｐ＝５である。
　ｙ_ｉ（一致すべき変数）の成分は、炉頂圧の制御量（ＰＶ）と、前圧の制御量（ＰＶ）、炉頂圧の操作量（ＭＶ）、ＴＲＴ前圧の操作量（ＭＶ）、バイパス弁前圧の操作量（ＭＶ）としている。ｍ_ｙ＝５である。

　ｄｄ_ｉ（外乱にフィルタや積分等の操作を行った後の値）は、上記の高炉ガス流量外乱にフィルタを適用した後のフィルタ後高炉ガス流量外乱である。上述のように、ｍ_ｄ＝１である。
　ｆｄ（）は、離散系において、状態量を計算する制御モデルであり、図１６に示すモデルが離散化されたモデルである。ここでは、公用ソフトウェアであるThe MathWorks,Inc.製のSimulink（Simulinkは、The MathWorks,Inc.の登録商標）を用いて、離散化と計算を行っている。

　なお、離散系の方程式は、固定サンプリング周期で解くことを前提として説明しているが、上記ソフトウェアで可能なように、サンプリング周期は可変でもよい。また、モデル内の積分器は、上記ソフトウェアでは、離散化されて計算されることとなる。
　ｃｄ_ineq（）は、不等式制約式であり、上記の式（１９）である（ｍ_cdineq＝９２）。
　ｃｄ_eq（）は、等式制約式であるが、本具体例では用いていない（ｍ_cdeq＝０）。

　評価関数の重みは、それぞれ、以下の式（２０）、式（２１）のようになる。

　図３のフローチャートに基づき、最適化の結果、次のパラメータを得た。ここで、第１バイパス弁６ａのフィードフォワード量は実績値を用いている。また、翼弁開度は、トリップ後に実績と同じ値の０（ゼロ）としている。さらに、トリップ後には、ＴＲＴ３を流れる流量が０となるようにＴＲＴ３の前の遮断弁が閉じられるので、ＴＲＴ３を流れる流量を０としている。

　また、このときの高炉ガス流量外乱は、図１８Ｆに示すグラフ（太破線）のように、推定される。外乱の定常偏差はほぼ０（ゼロ）であり、周期が約２０秒の高周波域の外乱が推定されていることがわかる。なお、図１８Ａ～Ｅ、図１８Ｇの各グラフに、他の実績データ（実線）とシミュレーション結果のデータ（破線）を示す。
　このように、式（２２）に基づいて、ＣＶ値を補正し、精確なフィードフォワード補償量を得ることができる。

　図１９に、精確なフィードフォワード補償を適用した場合のシミュレーション結果のデータを示す。ここで、図１９Ｆのグラフ（太破線）に示すように、外乱として、既に推定されたものをシミュレーションに使用する。また、既に推定されたパラメータを使用する。
　なお、図１９Ａ～Ｅ、図１９Ｇの各グラフに、他の実績データ（実線）とシミュレーション結果のデータ（破線）を示す。図１９Ａの炉頂圧の制御量（ＰＶ）の実績データ、及びシミュレーション結果のデータを示すグラフを参照すると、目標の炉頂圧が１．３２５［－］であるので、炉頂圧変動が、図１８Ａ中の０．０３６から、図１９Ａ中の０．０１１へと低減されていることがわかる。

　なお、他に、上記した同時推定方法を用いて、ＰＩＤコントローラ７～９のゲイン調整等に使用することができる。例えば、外乱の影響を受ける操作量（ＭＶ）の大きさや速度について、飽和しないようなゲインを選ぶことができる。また、推定された外乱を、設備診断等に用いることもできる。さらに、Ｈ∞制御やオブザーバ等のコントローラ設計に用いることができる。

　以上述べた具体例では、制御系を切り換えるという非定常時のデータを用いて、パラメータと外乱を同時に推定したが、制御系の切り換えを伴わない非定常時の場合でも推定ができることを、変形例として以下に示す。
[変形例]
　図２０に、第１バイパス弁６ａを閉じたまま、第２バイパス弁６ｂを徐々に閉めていく場合の例の結果を示す。図２０Ｆのグラフ（太破線）には、推定された外乱が示されている。外乱の定常偏差はほぼ０（ゼロ）であり、また、周期が約４０秒の高周波域の外乱が推定されていることがわかる。なお、図２０Ａ～Ｅ、図２０Ｇの各グラフには、他の実績値（実線）とシミュレーション結果（破線）が示されている。

　ここで、評価関数と不等式制約条件は、式（１８）、（１９）と同じである。
　なお、第２バイパス弁６ｂの操作量（ＭＶ）には、実績値を与えているため、Ｋ_ＣＢ２は、推定すべきパラメータに含まれない。このことから、Ｋ_ＣＢ２をｃｄ_ｅｑ（）（等式制約式）で次のように設定しておけばよい。

　そして、図３のフローチャート図に基づき、最適化を行う。その結果、以下に示す式（２４）のようなパラメータを得た。ここで、なお、第１バイパス弁６ａ、第２バイパス弁６ｂの開度は、実績値をそのまま用いている。

　このように、操業条件に合わせて、シミュレーション条件を変え、また、推定すべきパラメータを変更することで、パラメータや外乱を適切に推定することができる。
［２例目］
　続いて、具体例の２例目について、説明する。２例目は、鉄鋼プロセスの１つである連続鋳造機３０の鋳型３１内にある溶鋼Ｘの湯面レベルを制御する方法について、外乱とパラメータとを同時推定し、制御偏差が生じる原因を特定し、さらに、制御系を調整する方法について説明する。

　図２１に示す如く、タンディッシュ３３の底部には、タンディッシュ３３内の溶鋼Ｘを鋳型３１へ排出する溶鋼孔３４が形成されている。この溶鋼孔３４に連通するように浸漬ノズル３５が設けられ、浸漬ノズル３５の先端部は鋳型３１内に挿入されるものとなっている。タンディッシュ３３の底部であって、浸漬ノズル３５の根本には、溶鋼孔３４に隣接するように、浸漬ノズル３５内を流下する溶鋼Ｘの流量を調整するスライドバルブ３６が設けられている。

　スライドバルブ３６は、定形耐火物で構成された固定の板状の上プレート３６ａ、下プレート３６ｂと、定形耐火物で構成された移動する板状の中間プレート３６ｃとを有している。 上プレート３６ａの下面と中間プレート３６ｃの上面とは面接触していて、中間プレート３６ｃの上面は、上プレート３６ａの下面を水平方向に摺動するようになっている。なお、下プレート３６ｂと中間プレート３６ｃとの関係は、上下の位置関係が異なる点を除いて、上プレート３６ａと中間プレート３６ｃとの関係と同じである。上プレート３６ａ、下プレート３６ｂと中間プレート３６ｃとには、上下に貫通していてタンディッシュ３３から排出された溶鋼Ｘを通す流下孔３７（バルブ）が設けられている。

　中間プレート３６ｃは、ロッドが水平方向に伸縮する電気油圧シリンダ３９（アクチュエータ）に連結されていて、この電気油圧シリンダ３９のロッドを伸縮させることで、中間プレート３６ｃは、上プレート３６ａと下プレート３６ｂとの間を水平方向に摺動し、流下孔３７の断面積（バルブ断面積）が変化し、流下する溶鋼Ｘの流量が可変とされる。さらに、図２１に示す如く、鋳型３１の上方には、鋳型３１内の湯面レベルを計測する湯面レベル計３８が設けられている。なお、湯面レベルとは、鋳型３１内の溶鋼表面（湯面）の高さのことをいう。加えて、この湯面レベル計３８の計測結果などを基に、鋳型３１内の湯面レベルを目標値に一致させるべく、スライドバルブ３６の位置、言い換えれば電気油圧シリンダ３９の伸縮量を算出し、算出された量をスライドバルブ３６へ適用する制御部４０（コントローラ）が、連続鋳造機３０には備えられている。

　図２２は、制御部４０で行われる鋳型３１内の湯面レベル制御系（単に、レベル制御系と呼ぶこともある）のブロック図を示したものである。
　図２２において、Ｃ_ＦＢ（ｚ）がフィードバックコントローラであり、遅延要素ｚ－１の関数である。なお、遅延要素がない場合（動特性がない場合）も含むものとする。アクチュエータ３９や湯面レベル計３８は、それぞれ「一次遅れ＋むだ時間」で表され、シミュレーション時には、離散化される。

　シミュレーションには、既述のように、公用ソフトウェアであるThe MathWorks,Inc.製のSimulink（Simulinkは、The MathWorks,Inc.の登録商標）を用いることができる。スライドバルブ３６は、流量係数Ｋ_ｆのゲインで表される。鋳型３１は、鋳型断面積Ａ_moldの逆数をゲインとする積分器で表される。
　図２３に、湯面レベルが定常偏差を持ってハンチングしている異常状態の一例を示す。

　スライドバルブ開度のグラフ（実線）は、実績値を示す。湯面レベルのグラフは、実線が実績値（湯面レベル制御のＰＶ）を示し、一点鎖線が、目標値（湯面レベル制御のＳＶ）を示す。
　時刻０～約５０秒の間は、コントローラＡ（４０）が制御している。時刻約５０～１６０秒の間は、コントローラＢ（４０）が制御している。湯面レベルのグラフを見てわかるように、周期が約４０秒のハンチングを生じている。さらに、湯面レベル目標値（湯面レベルＳＶ）が、－１．０２［－］（一点鎖線）であるのに対して（［－］は無次元化した値であることを示す。）、湯面レベル実績値は、約－０．９５［－］を中心に振動しており、レベル偏差にオフセットを生じている。これらのハンチングと定常偏差が制御上の問題である。

　ハンチングと定常偏差の原因の候補として、外乱（例えば、バルジング等による外乱）、制御によるハンチング（例えば、制御ゲインが高すぎる、時定数が長すぎる、むだ時間が長すぎる、把握していない制御ゲインがループ内に存在する等）、ハードウェアの異常動作（想定外の動作）等が挙げられる。しかし、外乱と制御系パラメータの両方が不明であるため、異常原因の特定が困難であった。

　そこで、制御モデルのパラメータと、その制御モデルに加わる外乱の推定により、このレベル制御異常の原因を特定することとした。
　図２４に、連続鋳造機３０の湯面レベル制御系の制御モデルと外乱を示す。
　図２４中のフィルタｓは、図２５に示すようなバンドストップフィルタとしている。具体的には、（ｓ^２＋ω_ｂ ^２）／（ｓ^２＋２ζ_ｂω_ｂｓ＋ω_ｂ ^２），ω_ｂ＝２π／１３．５，ζ_ｂ＝０．８としている。すなわち、湯面レベル外乱の周期は、１３．５秒近傍としている。湯面レベルの外乱の主な成分は、バルジングと呼ばれる凝固殻が、鋳型３１下方のロール３２の間で膨らんだり凹んだりすることで、鋳型３１の湯面高さが変わることによる外乱（バルジング外乱）であり、この外乱の周波数は、鋳造速度／ロール３２間隔である。

　ロール３２間隔が一定でない場合は、バルジングを生じる部分のロール３２間隔の最大・最小から、周波数の最小・最大を求め、この最小・最大の間の周波数の信号をストップするようなフィルタとすればよい。
　図２４中のフィルタｃは、図２６に示すように、コントローラ４０に加わる外乱の微分値（コントローラ外乱の微分値）に適用するフィルタである。フィルタの形は、図２５と同じとしている。外乱の微分値に周期が１３．５秒近傍の成分が含まれていると予想して推定することとなる。なお、後述のように、結果として、この周期のコントローラ外乱は存在しないと推定された。すなわち、周期１３．５秒近傍の外乱は、湯面レベル外乱であることが確認された。コントローラ４０に加わる外乱は、図２７に示すように、積分され、ランプ状の外乱を表現することができる。

　続いて、本具体例の制御モデルについて、説明する。
　コントローラＡ，Ｂ（４０）は、ｋ＝０，１に対応する。ここで、コントローラＡ（４０）で制御している場合と、コントローラＢ（４０）で制御している場合を、それぞれ別のデータとして切出すこととする。すなわち、実績データの時系列値が複数群存在する場合である。この例では、２つの群が存在することとなる。コントローラＡ，Ｂ（４０）は、次の式で表される。

　なお、ｚ^－１は、既述のように、サンプリング周期の時間の遅延要素を意味している。また、ｚ^－ｎｏｄは、ｎｏｄサンプリング周期分の時間の遅延要素を意味している。
　ここで、各係数は既知であり、次の値としている。

　各ｋの場合に、分母分子のそれぞれの遅延要素の初期値として、式（２５）の分子のｚ^－１，ｚ^－２，ｚ^－３，ｚ^－４にそれぞれ、ｕ_{1init_k}，ｕ_{2init_k}，ｕ_{3init_k}，ｕ_{４init_k}（ｋ＝０，１）を代入し、式（２５）の分母のｚ^－１，ｚ^－２，ｚ^－３，ｚ^－４にそれぞれ、ｙ_{1init_k}，ｙ_{2init_k}，ｙ_{3init_k}，ｙ_{４init_k}（ｋ＝０，１）を代入する。これにより、コントローラ４０の出力の初期値を計算することができる。

　コントローラ４０の未知のゲインを、ｋ_conとする。例えば、コントローラ４０の仕様が明確でない場合等を想定している。
　制御対象のゲインである流量係数／鋳型断面積は、次のように決定した。

　次の積分器は、１／ｓで表され、出力の初期値をｌｅｖｅｌ_{init_k}（ｋ＝０，１）とする。
　次の一次遅れは、制御対象の一次遅れ成分を意味し、一次遅れ「１／（Ｔ_lag・ｓ＋１）」で表され、時定数Ｔ_lagがパラメータであり、状態量を１つ含んでいるため、出力の初期値をｌｅｖｅｌ_{laginit_k}（ｋ＝０，１）とする。

　この制御モデルでは、アクチュエータ３９の時定数は小さいとして無視し、湯面レベル真値と湯面レベル観測値との間に一次遅れがあるものとしている。この一次遅れには、センサの特性と、センサに追加されたフィルタの特性とを含むものとする。なお、シミュレーション時には、適宜、離散化され解かれることとなる。
　次のむだ時間要素は、制御対象のむだ時間の成分を意味し、むだ時間Ｔ_delay［ｓ］前の値を現在の値とする要素である。シミュレーションの初期には、過去の値を初期値として与える必要があり、過去のＴ_delay［ｓ］前までの初期値をｌｅｖｅｌ_{delayinit_k}（ｋ＝０，１）とするものとする。

　この制御モデルでは、アクチュエータ３９のむだ時間は小さいとして無視し、湯面レベル真値と湯面レベル観測値との間にむだ時間があるものとしている。このむだ時間には、センサの特性、フィルタの特性、制御のサンプリング周期の特性等を含むものとする。なお、シミュレーション時には、適宜、離散化され解かれることとなる。
　ゲインｋ_plantは、制御対象の未知のゲインを示し、本具体例では、流量係数Ｋ_ｆの変動分や他の動特性による変動分を合わせた制御対象全体の未知ゲインを表す。ｋ_plantは、制御モデルと一致していれば１である。

　開度定常値ＭＶ_{steady_k}（ｋ＝０，１）は、バルブ開度の定常値である。式（２５）から得られる操作量（ＭＶ）は、偏差系に対する操作量で、０（ゼロ）が基準になるため、初期の状態での定常値を加えた値が、シミュレーションで得られるスライドバルブ開度となる。
　本具体例では、図２４などに示されている「湯面レベル目標値」、「湯面レベル」、「バルブ開度」が、それぞれ、レベル制御系の「湯面レベルＳＶ」、「湯面レベルＰＶ」、「湯面レベルＭＶ」となる。

　以上より、推定すべきパラメータは、ｕ_{1init_k}，ｕ_{2init_k}，ｕ_{3init_k}，ｕ_{4init_k}，ｙ_{1init_k}，ｙ_{2init_k} ，ｙ_{3init_k}，ｙ_{4init_k}(ｋ＝０，１），ｋ_con，ｌｅｖｅｌ_{init_k}(ｋ＝０，１），Ｔ_lag，ｌｅｖｅｌ_{laginit_k}(ｋ＝０，１），Ｔ_delay，ｌｅｖｅｌ_{delayinit_k}(ｋ＝０，１），ｋ_plant，ＭＶ_{steady_k}(ｋ＝０，１）となる。
　この中で、ｋ_con，Ｔ_lag，Ｔ_delay，ｋ_plantが、共通パラメータであり、残りが個別のパラメータとなる。ここで、共通パラメータとは、複数の実績データを用いるときに、モデル内で共通の値として用いられるパラメータである。すなわち、実績データのすべての群にて共通する値をもつパラメータである。個別のパラメータとは、複数の実績データを用いるときに、モデル内で個別の値として用いられるパラメータである。すなわち、実績データの各群に連動する値をもつパラメータである。

　また、ｕ_{1init_k}，ｕ_{2init_k}，ｕ_{3init_k}，ｕ_{4init_k}，ｙ_{1init_k}，ｙ_{2init_k}，ｙ_{3init_k}，ｙ_{4init_k}(ｋ＝０，１），ｌｅｖｅｌ_{init_k}(ｋ＝０，１），ｌｅｖｅｌ_{laginit_k}(ｋ＝０，１），ｌｅｖｅｌ_{delayinit_k}(ｋ＝０，１），ＭＶ_{steady_k}(ｋ＝０，１）が、未知の初期値となる。ここで、未知の初期値とは、シミュレーションの初期に与えるべき値のうち、未知のものをいう。

　さらに、ＭＶ_{steady_k}(ｋ＝０，１）は、未知の定常値である。未知の定常値とは、制御系として偏差系を構成し、シミュレーションをする際に与えるべき定常値のうち、未知のものをいう。
　続いて、本具体例において、推定すべき外乱は、湯面レベルに加わる外乱と、スライドバルブ開度に加わるランプ状の外乱との２つである。

　１つ目の湯面レベルに加わる外乱は、よく知られたバルジングに起因する外乱ｄ_level k j (ｋ＝０，１、ｊ＝０,…,Ｍ_ｋ－１）であり、図４に示すように、時系列の各時刻の値に対応したカウントでの値を変数とする。この１つ目の外乱に、既述のように、フィルタｓを適用した後のフィルタ後レベル外乱ｄ_{levelfilterd k j} (ｋ＝０，１、ｊ＝０,…,Ｍ_ｋ－１）を、後述のように評価する。

　２つ目の外乱は、制御用のコンピュータが操作量（ＭＶ）を出力し、その後、実際にスライドバルブ開度の変化となるまでに印加する外乱である。
　図２３に示すように、２つ目の外乱は、湯面レベルが定常偏差を持っていることから、ランプ状の外乱であると推定されたものである。なお、この２つ目の外乱は、調査により、操作信号の伝達経路で信号が欠落することにより発生したことが後程確認された。

　本具体例では、スライドバルブ開度に加わる外乱の微分値（コントローラ外乱の微分値）、ｄ（ドット）_MV k j（ｋ＝０，１、ｊ＝０，…，Ｍ_ｋ－１）を推定することとしている。図４に示すように、外乱の微分値を、時系列の各時刻の値に対応したカウントでの値を変数とした。この外乱の微分値に、既述のように、フィルタｃを適用した後のフィルタ後コントローラ外乱微分値ｄ（ドット）_{MVfilterd k j}（ｋ＝０，１、ｊ＝０，…，Ｍ_ｋ－１）を、後述のように評価する。

　なお、本具体例では、２つの実績データ、すなわち、コントローラＡ（４０）で制御している場合の実績データと、コントローラＢ（４０）で制御している場合の実績データとを用いているため、外乱は各実績データに対して別々であり且つ個別の外乱を用いている。個別の外乱とは、複数の実績データを用いてシミュレーションを行うときに、実績データ毎に個別の時系列の値として用いられる外乱である。すなわち、実績データの各群に連動する時系列変数値である。

　本具体例の評価関数は、次のように定めた。すなわち、２つの実績データに対応するモデルで共通なパラメータは共通の変数とし、別々のパラメータは別の変数として推定することとしている。外乱については、別々なので、それぞれ別に推定することとなる。
　そして、制御モデル内のパラメータの評価関数としては、次の式（２８）を用いた。

　第１項、第２項は、湯面レベルとバルブ開度とを、実績とシミュレーションで一致させるための項である。第３項、第４項は、フィルタ後レベル外乱と、フィルタ後コントローラ外乱微分値とを０（ゼロ）に近づけるための項である。なお、Ｎ＝２，Ｌ_０＝５２，Ｌ_１＝１０１である。実績データと、シミュレーション結果のデータのサンプリング周期は、１秒である。

　制御モデルの制約条件は、次の通りである。

　各制約条件については、パラメータの取り得る値を考慮して、広めに設定している。
　以上をまとめると、本具体例における湯面レベル制御は、データ数が１つの場合を、データ数が複数の場合に拡張した一般形を、具体化した形となっている。
　次に、本具体例における湯面レベル制御について、具体的に説明する。
　ｘ_ｋｉ（状態量）は、各ｋ（ｋ＝０，１）について、コントローラ４０の遅延要素の８個（分子４＋分母４）、鋳型３１の積分器１個に対応する遅延要素（離散化の方法に依存する）、一次遅れの積分器１個に対応する遅延要素（離散化の方法に依存する）、むだ時間要素の遅延要素（数は、離散化のサンプリング時間に依存する）、コントローラ外乱に対する積分器１個に対応する遅延要素である（状態量の数はシミュレーション方法に依存する。）。

　ｘ０_ｋ（初期値）は、コントローラ４０の遅延要素の８個（分子４＋分母４）の初期値、鋳型３１の積分器に対応する遅延要素の初期値、一次遅れの積分器に対応する遅延要素の初期値、むだ時間要素の遅延要素の初期値である。なお、コントローラ外乱の微分値に適用する積分器に対応する遅延要素の初期値は０とした。
　ｒ_ｋｉ（既知の入力）は、各ｋに対応した湯面レベルＳＶ（レベル目標値）である。

　ｄ_ｋｊ（未知の入力）は、推定すべき外乱であり、コントローラ外乱の微分値と、レベル外乱との２つである（ｋ毎にｍ_ｄｋ＝２）。
　パラメータとしては、図４のように、時系列の各時刻の値を変数とし、Ｍ_０＝２０，Ｍ_１＝４０とした。外乱の変数の時間間隔が、実績データやシミュレーション結果のデータのサンプリング周期１秒より長いのは、外乱の変数を減らすことを可能とするためである。例えば、外乱の時間的変化が滑らかと予想される場合には、外乱の変数の時間間隔を長くとることができ、最適化のための計算量を減らすことができる。

　ｐ（推定すべきパラメータ）は、既述のように、ｕ_{1init_k}，ｕ_{2init_k}，ｕ_{3init_k}，ｕ_{4init_k}，ｙ_{1init_k}，ｙ_{2init_k}，ｙ_{3init_k}，ｙ_{4init_k}(ｋ＝０，１），ｋ_con，ｌｅｖｅｌ_{init_k}(ｋ＝０，１），Ｔ_lag，ｌｅｖｅｌ_{laginit_k}(ｋ＝０，１），Ｔ_delay，ｌｅｖｅｌ_{delayinit_k}(ｋ＝０，１），ｋ_plant，ＭＶ_{steady_k}(ｋ＝０，１）となる（ｍ_ｐ＝２８）。

　ｙ_ｋｉは、一致すべき変数であり、湯面レベルＰＶ（湯面レベル）、湯面レベルＭＶ（バルブ開度）としている（ｋ毎にｍ_ｙｋ＝２）。
　ｄｄ_ｋｉは、外乱にフィルタや積分等の操作を行った後の値であり、コントローラ外乱の微分値にローパスフィルタを適用した後の値（フィルタ後コントローラ外乱微分値）と、レベル外乱にバンドストップフィルタを適用した後の値（フィルタ後レベル外乱）である。

　ｆｄ（）は、離散系において、状態量を計算する制御モデルであり、図２４に示すモデルから離散化して得られる制御モデルである。本具体例では、公用ソフトウェアであるThe MathWorks,Inc.製のSimulink（Simulinkは、The MathWorks,Inc.の登録商標）を用いて行っている。
　なお、本具体例の離散系の方程式は、固定サンプリング周期で解くことを前提として説明しているが、上記ソフトウェアで可能なように、サンプリング周期は可変でよい。また、制御モデル内の積分器は、上記ソフトウェアにおいて、離散化されて計算されることとなる。

　ｃｄ_ineq（）は、不等式制約式であり、上記の式（２９）である（ｍ_cdineq＝２９６）。
　ｃｄ_eq（）は、等式制約式であるが、本具体例では用いていない（ｍ_cdeq＝０）。
　制御モデル内のパラメータの評価関数は、式（９）を、複数の実績データに対応して、拡張した次の式（３０）となる。

　ここで、各変数は次の通りである。
　Ｎは、実績データの数である。
　ｙ_actual k iは、ｋ番目のデータの、ｍ_ｙｋ行１列のベクトル（ｉ＝０,・・・,Ｌ_ｋ－１）。各成分は、カウントｉでの、時系列の実績データと、シミュレーション結果のデータとを一致させる場合の、実績データである。

　Ｗ_ｙｋは、ｍ_ｙｋ行ｍ_ｙｋ列の行列。シミュレーション結果のデータの重みを示す。本具体例では、対称な準正定行列とする。
　Ｗ_ｄｄｋは、ｍ_ｄｋ行ｍ_ｄｋ列の行列。外乱にフィルタや積分等の操作を行った後の値の重みを示す。本具体例では、対称な準正定行列とする。
　具体的な数値は、次のようになる。

　図２８のフローチャート図に基づいた最適化計算を行い、その結果、以下に示すパラメータを得た。

　また、このときの２つの実績データに対応する外乱は、図２９Ｅ～図２９Ｈのグラフ（太破線）のように推定される。なお、図２９Ｇ、図２９Ｈに示すコントローラ外乱については、積分した後の値を表示している。図２９Ｅ、図２９Ｆに示すように、推定されたレベル外乱が小さいことから、連続鋳造機３０において、小さなバルジングが発生していることが推測される。バルジング外乱は、周期が約１３．５秒と推定されている。また、湯面レベル外乱の定常偏差はほぼ０（ゼロ）と推定されている。また、コントローラ外乱として、ランプ状の外乱が印加されていることが推測される。細かくは、コントローラ外乱として、ランプ状の外乱と、重畳した周期約４０秒の小さな外乱が印加していると推定されている。これらの周波数帯域の外乱は、外乱微分値が存在すると予想した周波数帯域外にあり、外乱微分値が小さな値として推定されていることとなる。なお、存在すると予想した周期約１３．５秒近傍の外乱は、推定されておらず、周期１３．５秒近傍の外乱は、コントローラ外乱ではなく、レベル外乱であると推測される。

　ちなみに、外乱の時系列値を評価しない場合、すなわち、Ｗ_ｄｄ０，Ｗ_ｄｄ１が零行列の場合には、レベル外乱がランプ状であるような想定外の推定結果が得られた。したがって、外乱の時系列値を評価しない場合、外乱が大きくなる解がある場合には、外乱の大きさに制限を加えることが望ましい。すなわち、外乱の時系列値を評価する方が望ましい。
　なお、図２９Ａ～図２９Ｄは、実績データ（実線）とシミュレーション結果のデータ（破線）である。一点鎖線は、目標値を示している。

　以上の推定結果から、制御コンピュータからスライドバルブ開度までの信号経路に異常があることが推測される。なお、既述のように、実際に異常が存在したことが確認されている。
　また、周期約４０秒の振動は、主に、制御対象の位相遅れによる制御起因であると推定される。

　そこで、信号経路上の異常に対策を実施し、制御対象の位相遅れを除く処理を行った。
　その結果、図３０に示すように、この周期約４０秒の大きな振動はなくなり、レベル精
度を向上させることができた。
　なお、制御モデル内のパラメータと、その制御モデルに加わる外乱が得られているので、従来のＰＩＤ制御を適用する場合に、ゲイン調整を行い、制御することも、既述のように可能である。
［３例目］
　さて、具体例の３例目について、３つ例を挙げて説明する。この３例目は、評価関数に外乱の時系列変数値を含まない例である。
（３例目の１）
　まず、３例目の１は、図３１に示すように、ゲインＫ、時定数Ｔの一次遅れ、むだ時間Ｌで構成された制御対象である。なお、コントローラとしてＰＩ制御が使用され、比例ゲインをＫ_Ｐ、積分ゲインをＫ_Ｉとする。また、一次遅れの初期値をｌａｇｉｎｉｔ＿_ｔｖ、むだ時間の初期値をｄｅｌａｙｉｎｉｔ＿_ｔｖ、コントローラの積分要素の初期値をＩｉｎｉｔ＿_ｔｖとする。また、連続系の要素は、シミュレーション時に離散化されることとなる。また、外乱をｄ（ｔ）とおく。

　次の式（３６）のように、制御モデルのパラメータと、外乱の真値を決定する。なお、＿ｔｖは、真値を示す添え字である。

　図３２Ａ～Ｄに、真値のパラメータと、真値の外乱の場合の実績データを示す。なお、ここでは、真値、または、真値を用いたシミュレーション結果を、実績データとしている。本具体例では、この実績データから、制御モデルのパラメータと外乱とを逆算することを考えている。なお、実績データの外乱は、ランプ状外乱やオフセットを含んでいることとしている。

　評価関数としては、次の式（３７）を用いた。

　ここでは、式（３７）において、添え字の「実績」は、既に得られている実績値を意味し、「シミュレーション」は、シミュレーション値を意味する。
　また、式（３７）の第１項～第２項は、ｙ_ｉについて、シミュレーション値を実績値に一致させるための項であり、具体的には、制御量（ＰＶ）と、操作量（ＭＶ）を、実績とシミュレーションで一致させるための項である。

　なお、式（３７）では、外乱を評価していない。式（９）で、Ｗ_ｄｄ＝０とした場合に相当し、評価関数が外乱の時系列変数値を含まないこととなる。
　制約条件については、次の式（３８）の通りであり、式（７）に対応する。

　なお、各制約条件については、パラメータの取り得る値を考慮して、広めに設定している。
　以上をまとめると、本具体例は、一般形を具体化した形となっている。以下に、本具体例の詳細を説明する。
　ｘ_ｉ（状態量）の成分は、まず、時定数の積分器、ＰＩＤコントローラの積分器を離散化した場合の状態量である。さらに、むだ時間を離散化した場合の状態量、例えば、むだ時間１秒をサンプリング時間０．１秒で離散化した場合には過去の１０点の１０個の状態量がある。ｍｘは、離散化方法やむだ時間により異なることとなる。なお、Ｌ＝５１とした。

　ｘ０（初期値）については、時定数の積分器とＰＩＤコントローラの積分器を離散化した場合の状態量の初期値、むだ時間を離散化した場合の状態量の初期値である。
　ｒ_ｉ（既知の入力）の成分は、目標値（ＳＶ）である。この場合、ｍ_ｒ＝１である。
　ｄ_ｊ（未知の入力）は、推定すべき外乱である。この場合、ｍ_ｄ＝１である。なお、Ｍ＝２６とした。

　ｐ（推定すべきパラメータ）は、Ｋ，Ｔ，Ｌ，Ｋ_Ｐ，Ｋ_Ｉ，ｌａｇｉｎｉｔ，ｄｅｌａｙｉｎｉｔ，Ｉｉｎｉｔとである。この場合、ｍ_ｐ＝８である。
　ｙ_ｉ（一致すべき変数）の成分は、ＰＶとＭＶとしている。ｍ_ｙ＝２である。
　ｆｄ（）は、離散系において、状態量を計算する制御モデルであり、図３１に示すモデルが離散化されたモデルである。ここでは、公用ソフトウェアである、The MathWorks,Inc.製のSimulink（Simulinkは、The MathWorks,Inc.の登録商標）を用いて、離散化と計算を行っている。

　なお、離散系の方程式は、固定サンプリング周期で解くことを前提として説明しているが、上記ソフトウェアで可能なように、サンプリング周期は可変でもよい。また、モデル内の積分器は、上記ソフトウェアでは、離散化されて計算されることとなる。
　ｃｄ_ineq（）は、不等式制約式であり、上記の式（１２）である（ｍ_cdineq＝６８）。
　ｃｄ_eq（）は、等式制約式であるが、本具体例では用いていない（ｍ_cdeq＝０）。

　評価関数の重みは、それぞれ、以下の式（３９）、式（４０）のようになる。

　図３に示しているフローチャートに基づき、最適化を行う。なお、最適化の初期値は、制約条件を満足する範囲で与えて、複数回最適化を行い、最小の評価関数の値０．００を持つ場合を採用した。
　なお、このときの初期値は、次の式（４１）の通りである。

　図３３Ａ～Ｄに、最適化後のシミュレーション結果を示す。なお、ＳＶは、シミュレーション時に実績と同じＳＶを与えている。
　図３３Ａ～Ｄに示すように、評価関数値が０．００となっていることからわかるように、ＰＶ，ＭＶについては、実績とシミュレーションの時系列値が一致している。得られた制御モデルのパラメータは、次の式（４２）の通りである。

　なお、外乱については、図３３Ｄのグラフにあるように、推定された外乱は、実績の外乱と一致している。
（３例目の２）
　続いて、３例目の１の例において、さらに評価関数が外乱の時系列変数値を含む場合の推定の例、すなわち３例目の２について説明する。なお、検討する理由としては、２例目の場合のように、外乱の時系列変数値を評価しないと、精確なパラメータや外乱を推定できない場合があるからである。

　図３４に従って、外乱を評価する場合を説明する。
　なお、制御モデルのパラメータ、外乱の真値は、式（３６）と同じとする。図３２Ａ～Ｄが、真値のパラメータと、真値の外乱の場合の実績データを示す。なお、実績データの外乱は、ランプ状外乱やオフセットを含んでいる。
　評価関数としては、次の式（４３）を用いた。

　なお、式（４３）の第３項は、外乱を０（ゼロ）に近づけるための項である。また、制約条件については、式（３８）と同じである。
　以上をまとめると、本具体例は、一般形を具体化した形となっている。既述の３例目の１との違いのみを、以下に説明する。
　ｄｄ_ｉは、外乱ｄ_ｊと、同時刻で同じ値を持っている。なお、ｄ_ｊからｄｄ_ｉを求めるには、補間等を行うとよい。

　評価関数の重みは、それぞれ、式（３９）と、以下の式（４４）のようになる。

　図３のフローチャートに基づき、最適化を行う。なお、最適化の初期値は、制約条件を満足する範囲で与えて、複数回最適化を行い、最小の評価関数の値、１．１９×１０^４を持つ場合を採用した。
　なお、このときの初期値は、次の式（４５）の通りである。

　図３５Ａ～Ｄに、最適化後のシミュレーション結果を示す。なお、ＳＶは、シミュレーション時に実績と同じＳＶを与えている。図３５Ａ～Ｄに示すように、ＰＶ，ＭＶについては、実績とシミュレーションの時系列値が一致している。得られたパラメータは次の式（４６）の通りである。

　以上より、ゲインで３３％、時定数で２９％、むだ時間で６％の誤差がある。
　また、外乱についても、図３５Ｄのグラフに示すように、推定された外乱は、実績の外乱と誤差があることがわかる。
　なお、外乱に対する重みＷ_ｄｄを小さくすると、外乱の推定値の精度が良くなるが、最適化の収束に時間がかかる傾向にある。

　したがって、本具体例では、外乱に単純に重みを掛けて評価した場合、若干の誤差を持つが、制御モデルのパラメータと外乱の推定は可能である。
（３例目の３）
　続いて、３例目の２の例において、外乱に周波数が存在する場合に、当該周波数成分をカットするフィルタを適用し、評価関数が、フィルタ適用後の外乱の時系列変数値を含む場合の推定の例、すなわち３例目の３について説明する。

　本具体例は、フィルタを適用することにより、外乱をより精確に推定できる。
　図３６に示すように、フィルタ適用後の外乱を評価する場合を説明する。フィルタは、時定数３０[s]の一次遅れ要素である。
　制御モデルのパラメータ、外乱の真値は、式（３６）と同じとする。図３２Ａ～Ｄが、真値のパラメータと、真値の外乱の場合の実績データを示す。なお、実績データの外乱は、ランプ状外乱やオフセットを含んでいる。

　評価関数としては、次の式（４７）を用いた。

　なお、式（４７）の第３項は、フィルタ後の外乱を０（ゼロ）に近づけるための項である。外乱の成分はフィルタにより低減され、評価されないため、外乱自体は大きな値を取ることができる。このため、より精確な外乱を推定可能となる。
　また、制約条件については、式（３８）と同じである。
　以上をまとめると、本具体例は、一般形を具体化した形となっている。既述の３例目の２との違いのみを、以下に説明する。

　ｄｄ_ｉは、外乱ｄ_ｊにフィルタを適用した後の外乱である。
　評価関数の重みは、それぞれ、式（３９）と、以下の式（４８）のようになる。

　図３のフローチャートに基づき、最適化を行う。なお、最適化の初期値は、制約条件を満足する範囲で与えて、複数回最適化を行い、最小の評価関数の値３．３３×１０^２を持つ場合を採用した。
　なお、このときの初期値は、次の式（４９）の通りである。

　図３７Ａ～Ｄに、最適化後のシミュレーション結果を示す。なお、ＳＶは、シミュレーション時に実績と同じＳＶを与えている。図３７Ａ～Ｄに示すように、ＰＶ，ＭＶについては、実績とシミュレーションの時系列値が一致している。得られたパラメータは次の式（５０）の通りである。

　以上より、ゲインで４．４％、時定数で４．４％、むだ時間で０．７％だけの誤差しかない。
　また、外乱についても、図３７Ｄのグラフにあるように、推定された外乱は、ランプ状であっても、ほぼ一致している。
　したがって、本具体例では、外乱にフィルタを適用して評価した場合、外乱を精度よく推定可能であることがわかる。

　以上、今回開示された実施形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は上記した説明ではなくて特許請求の範囲によって示され、特許請求の範囲と均等の意味及び範囲内でのすべての変更が含まれることが意図される。

　本出願は、２０１５年３月２４日出願の日本特許出願（特願２０１５－０６１０５９）、及び２０１６年１月１２日出願の日本特許出願（特願２０１６－００３７７６）に基づくものであり、その内容はここに参照として取り込まれる。

　１　炉頂圧発電設備
　２　リングスリットエレメント（ＲＳＥ）
　３　発電部（ＴＲＴ、炉頂圧発電機）
　４　高炉ガス流路
　５　バイバス流路
　６　バイパス弁
　６ａ　第１バイパス弁
　６ｂ　第２バイパス弁
　７　炉頂圧コントローラ
　８　ＴＲＴ前圧コントローラ
　９　バイパス弁前圧コントローラ
　１０　円錐管
　１１　可動部
　１２　ケーシング
　１３　発電機
　１４　タービンの軸（ロータ）
　１５　動翼
　１６　静翼
　１６ａ　翼弁（弁機構）
　１７　高炉
　１８　サイレンサ
　３０　連続鋳造機
　３１　鋳型
　３２　ロール
　３３　タンディッシュ
　３４　溶鋼孔
　３５　浸漬ノズル
　３６　スライドバルブ
　３６ａ　上プレート
　３６ｂ　下プレート
　３６ｃ　中間プレート
　３７　流下孔（バルブ）
　３８　湯面レベル計
　３９　電気油圧シリンダ（アクチュエータ）
　４０　制御部（コントローラ）
　Ｘ　溶鋼

Claims (9)

1. 　制御モデル内のパラメータと当該制御モデルが表す制御対象へ入力される外乱とを同時に推定する方法であって、前記外乱を時系列変数値として、前記制御モデルに与えることとし、前記制御対象における実績データの時系列値と、前記制御モデルをシミュレーションすることで得られたシミュレーション結果のデータの時系列値を含む評価関数が、最大又は最小となるように、前記外乱と前記パラメータとを同時に推定することを特徴とする制御モデルのパラメータと外乱との同時推定方法。
2. 　制御モデル内のパラメータと当該制御モデルが表す制御対象へ入力される外乱とを同時に推定する方法であって、前記外乱を時系列変数値として、前記制御モデルに与えることとし、前記制御対象における実績データの時系列値と、前記制御モデルをシミュレーションすることで得られたシミュレーション結果のデータの時系列値と、前記外乱の時系列変数値とを含む評価関数が、最大又は最小となるように、前記外乱と前記パラメータとを同時に推定することを特徴とする制御モデルのパラメータと外乱との同時推定方法。
3. 　前記評価関数には、実績データおよびシミュレーション結果のデータとして、制御量（ＰＶ）、操作量（ＭＶ）の少なくとも一つが与えられていることを特徴とする請求項１又は２に記載の制御モデルのパラメータと外乱との同時推定方法。
4. 　前記外乱に周波数成分が存在する場合には、当該周波数成分をカットするフィルタを適用した後の外乱が前記評価関数へ含まれることを特徴とする請求項１～３のいずれか一項に記載の制御モデルのパラメータと外乱との同時推定方法。
5. 　前記実績データの時系列値が複数群、存在する場合においては、前記パラメータを前記実績データの全ての群にて共通する値をもつパラメータと、前記実績データの各群に連動する値をもつパラメータとに分離した上で、シミュレーション結果のデータの時系列値を評価関数へ導入すると共に、前記外乱を実績データの各群に連動する時系列変数値として評価関数へ導入することを特徴とする請求項１～４のいずれか一項に記載の制御モデルのパラメータと外乱との同時推定方法。
6. 　前記パラメータとして、未知の初期値が採用されることを特徴とする請求項１～５のいずれか一項に記載の制御モデルのパラメータと外乱との同時推定方法。
7. 　前記パラメータとして、未知の定常値が採用されることを特徴とする請求項１～６のいずれか一項に記載の制御モデルのパラメータと外乱との同時推定方法。
8. 　前記外乱がランプ状である場合には、前記外乱の微分値を評価関数へ導入すると共に、前記外乱の微分値を積分した上で、当該積分値を前記制御モデルに与えることを特徴とする請求項２～７のいずれか一項に記載の制御モデルのパラメータと外乱との同時推定方法。
9. 　請求項１～８のいずれか一項に記載された制御モデルのパラメータと外乱との同時推定方法で推定されたパラメータと外乱とのうち少なくともパラメータを用いて、前記制御対象を制御することを特徴とする制御対象の制御方法。

Images