(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2022183711
(43)【公開日】2022-12-13
(54)【発明の名称】外乱推定装置、外乱推定方法およびプログラム
(51)【国際特許分類】
G05B 13/02 20060101AFI20221206BHJP
F23G 5/50 20060101ALI20221206BHJP
【ＦＩ】
G05B13/02 C
F23G5/50 A
【審査請求】未請求
【請求項の数】9
【出願形態】ＯＬ
(21)【出願番号】P 2021091155
(22)【出願日】2021-05-31
(71)【出願人】
【識別番号】000006208
【氏名又は名称】三菱重工業株式会社
(74)【代理人】
【識別番号】100149548
【弁理士】
【氏名又は名称】松沼 泰史
(74)【代理人】
【識別番号】100162868
【弁理士】
【氏名又は名称】伊藤 英輔
(74)【代理人】
【識別番号】100161702
【弁理士】
【氏名又は名称】橋本 宏之
(74)【代理人】
【識別番号】100189348
【弁理士】
【氏名又は名称】古都 智
(74)【代理人】
【識別番号】100196689
【弁理士】
【氏名又は名称】鎌田 康一郎
(72)【発明者】
【氏名】広江 隆治
(72)【発明者】
【氏名】井手 和成
(72)【発明者】
【氏名】佐瀬 遼
【テーマコード（参考）】
3K062
5H004
【Ｆターム（参考）】
3K062AA24
3K062AC01
5H004GA07
5H004GB01
5H004GB02
5H004GB04
5H004GB14
5H004HA01
5H004HA02
5H004HA04
5H004HB01
5H004HB02
5H004HB04
5H004JB07
5H004JB18
(57)【要約】
【課題】制御対象に生じる外乱を推定する装置を提供する。
【解決手段】外乱推定装置は、制御対象が備えるセンサが計測した計測値を取得する取得部と、計測値を要素とする計測ベクトルの分散共分散行列を算出し、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、その特異ベクトルに基づいて制御対象に生じる外乱を推定する推定部と、を備える。
【選択図】図１
【解決手段】外乱推定装置は、制御対象が備えるセンサが計測した計測値を取得する取得部と、計測値を要素とする計測ベクトルの分散共分散行列を算出し、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、その特異ベクトルに基づいて制御対象に生じる外乱を推定する推定部と、を備える。
【選択図】図１
【特許請求の範囲】
【請求項1】
制御対象が備えるセンサが計測した計測値を取得する取得部と、
前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定する推定部と、
を備える外乱推定装置。
前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定する推定部と、
を備える外乱推定装置。
【請求項2】
前記取得部は、前記制御対象の制御の目的となる変数である制御量の計測値を取得し、 前記推定部は、前記計測ベクトルにおける前記最大特異値の特異ベクトルに基づいて、前記外乱を前記制御量の変動として推定する、
請求項１に記載の外乱推定装置。
請求項１に記載の外乱推定装置。
【請求項3】
前記推定部は、推定した前記外乱の分散の大きさに基づいて、推定した前記外乱の確からしさを判定し、推定した前記外乱とともに前記判定の結果を出力する、
請求項１または請求項２に記載の外乱推定装置。
請求項１または請求項２に記載の外乱推定装置。
【請求項4】
前記推定部は、前記外乱が前記計測値の変動として現れるまでの遅れ時間を補正する補正手段、を備え、
前記推定部は、前記取得部が取得した前記計測値の遅れ時間を前記補正手段によって補正し、補正後の前記計測値を要素とする前記計測ベクトルを用いて、前記外乱を推定する、請求項１から請求項３の何れか１項に記載の外乱推定装置。
前記推定部は、前記取得部が取得した前記計測値の遅れ時間を前記補正手段によって補正し、補正後の前記計測値を要素とする前記計測ベクトルを用いて、前記外乱を推定する、請求項１から請求項３の何れか１項に記載の外乱推定装置。
【請求項5】
前記推定部は、推定した前記外乱に基づいて、前記計測値のうちの前記遅れ時間が最も長い前記計測値の推定値を推定し、前記計測値の推定値と前記計測値の実測値の差の分散に基づいて、推定した前記外乱の確からしさを判定し、推定した前記計測値の推定値とともに前記判定の結果を出力する、
請求項４に記載の外乱推定装置。
請求項４に記載の外乱推定装置。
【請求項6】
前記推定部を複数備え、前記外乱の確からしさが最も高い前記推定部が推定した外乱を選択する、
請求項５に記載の外乱推定装置。
請求項５に記載の外乱推定装置。
【請求項7】
複数の前記推定部のそれぞれは、他の前記推定部と異なる前記計測値に基づいて前記外乱を推定するか、または、他の前記推定部と異なる前記遅れ時間の補正を行って前記外乱を推定するか、または、他の前記推定部と異なる周期で前記分散共分散行列および前記特異ベクトルを更新して前記外乱を推定する
請求項６に記載の外乱推定装置。
請求項６に記載の外乱推定装置。
【請求項8】
制御対象が備えるセンサが計測した計測値を取得するステップと、
前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定するステップと、
を有する外乱推定方法。
前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定するステップと、
を有する外乱推定方法。
【請求項9】
コンピュータに、
制御対象が備えるセンサが計測した計測値を取得するステップと、
前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定するステップと、
を実行させるプログラム。
制御対象が備えるセンサが計測した計測値を取得するステップと、
前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定するステップと、
を実行させるプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本開示は、外乱推定装置、外乱推定方法およびプログラムに関する。
【背景技術】
【0002】
ゴミ焼却炉にボイラを設置し、ゴミ焼却の際に発生する熱を回収し、発生した蒸気により発電を行なうゴミ発電は、ゴミを単に廃棄物としてではなく、ゴミに燃料としての付加価値を生じせしめる点で経済的に重要である。ゴミの燃料としての付加価値を向上するには発生する蒸気量を安定化させ、計画したとおりの発電ができるようにすることが効果的である。特許文献１には、ゴミ発電に関し、ゴミ焼却炉の発熱量の変動原因として廃棄物が有する水分に注目し、廃棄物の水分率の変動に応じて、単位時間あたりの焼却炉への廃棄物の供給量を調節する制御方法が開示されている。つまり、特許文献１には、蒸気流量に対する外乱を廃棄物の水分含有率から推定し、外乱が制御量（蒸気流量）に影響する前にゴミ供給量などを調整することによって、安定した発電を実現する技術が開示されている。
【0003】
特許文献２には、廃棄物の単位供給量当たりの発熱量を推定し、ゴミ焼却炉の燃焼制御を行う方法が開示されている。しかし、廃棄物の単位供給量当たりの発熱量の推算には、数時間のデータが必要であり、推算した値は数時間を平均化した値であるため、特に廃棄物の性質が時間的に変動する場合には、現時点での「廃棄物の単位時間当たりの発熱量」をタイムリーに推算することができない。このため、焼却炉に供給する廃棄物や燃焼空気の調節に使う推定ボイラ蒸発量（蒸気流量）が不正確となり、発電が変動となる可能性がある。特許文献２の技術では、（１）排ガス中の酸素および水分の成分濃度をセンサによって計測して、廃棄物の発熱量を算出し、（２）算出された廃棄物の発熱量に基づいてボイラ蒸発量を算出し、（３）算出したボイラ蒸発量を基に焼却炉に投入する廃棄物、燃焼空気などの供給量を制御する。つまり、特許文献２では、廃棄物の発熱量の変動によって外乱を推定し、外乱の表れである発熱量から算出したボイラ蒸発量に基づいてゴミ焼却炉の燃焼制御を行っている。
【先行技術文献】
【特許文献】
【0004】
【特許文献1】特開２０１９－１７８８５０号公報
【特許文献2】特許第５９９６７６２号公報
【発明の概要】
【発明が解決しようとする課題】
【0005】
特許文献１、２に開示の制御方法は、ゴミ焼却炉の特性を基として、蒸気流量などの制御量の変動に対して、制御量が安定するよう制御する技術であり、ゴミ焼却炉に固有の技術である。最近の機械学習が、同じ技術を他の対象物に適用することを可能とするような汎用性が無い。制御対象に生じる外乱を推定する汎用的な方法が求められている。
【0006】
本開示は、上記課題を解決することができる外乱推定装置、外乱推定方法およびプログラムを提供する。
【課題を解決するための手段】
【0007】
本開示の外乱推定装置は、制御対象が備えるセンサが計測した計測値を取得する取得部と、前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定する推定部と、を備える。
【0008】
また、本開示の外乱推定方法は、制御対象が備えるセンサが計測した計測値を取得するステップと、前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定するステップと、を有する。
【0009】
また、本開示のプログラムは、コンピュータに、制御対象が備えるセンサが計測した計測値を取得するステップと、前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定するステップと、を実行させる。
【発明の効果】
【0010】
上述の外乱推定装置、外乱推定方法およびプログラムによれば、外乱を推定することができる。
【図面の簡単な説明】
【0011】
【図1】各実施形態に係る制御システムの概略図である。
【図2】第一実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
【図3】第一実施形態に係る外乱推定処理の一例を示す図である。
【図4】第一実施形態に係る処理の更新タイミングの一例を示す図である。
【図5】第二実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
【図6】第二実施形態に係る外乱推定処理の一例を示す図である。
【図7】第三実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
【図8】第三実施形態に係る外乱推定処理の一例を示す図である。
【図9】第四実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
【図10】第四実施形態に係る外乱の発生時刻とその影響が計測値に現れるまでの時間差の一例を示す図である。
【図11】第四実施形態に係る外乱推定処理の一例を示す図である。
【図12】第五実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
【図13】第五実施形態に係る外乱推定処理の一例を示す図である。
【図14】第六実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
【図15】各実施形態に係る外乱推定装置のハードウェア構成の一例を示す図である。
【発明を実施するための形態】
【0012】
以下、実施形態の外乱推定装置を、図面を参照して説明する。以下の説明では、同一または類似の機能を有する構成に同一の符号を付す。そして、それら構成の重複する説明は省略する場合がある。「ＸＸまたはＹＹ」とは、ＸＸとＹＹのうちいずれか一方の場合に限定されず、ＸＸとＹＹの両方の場合も含み得る。これは選択的要素が３つ以上の場合も同様である。「ＸＸ」および「ＹＹ」は、任意の要素（例えば任意の情報）である。
【0013】
＜第一実施形態＞
（構成）
図１は、各実施形態に係る制御システムの概略図である。
制御システム１００は、外乱推定装置１０と、制御装置２０と、制御対象３０と、を含む。外乱推定装置１０は、制御対象３０の制御にあたり、その制御量に影響を与える外乱を推定する。外乱推定装置１０は、取得部１１と、推定部１２と、出力部１３と、を備える。取得部１１は、制御対象３０に設けられたセンサが計測する計測値を取得する。制御対象３０に外乱が生じると、その外乱は、取得部１１が取得する計測値の変動として表れる。推定部１２は、取得部１１が取得した計測値を用いて、制御量に生じた外乱ｑの大きさを推定する。出力部１３は、外乱ｑの推定値（以下、外乱ｑ又は推定値ｑと記載する。）を制御装置２０へ出力する。
（構成）
図１は、各実施形態に係る制御システムの概略図である。
制御システム１００は、外乱推定装置１０と、制御装置２０と、制御対象３０と、を含む。外乱推定装置１０は、制御対象３０の制御にあたり、その制御量に影響を与える外乱を推定する。外乱推定装置１０は、取得部１１と、推定部１２と、出力部１３と、を備える。取得部１１は、制御対象３０に設けられたセンサが計測する計測値を取得する。制御対象３０に外乱が生じると、その外乱は、取得部１１が取得する計測値の変動として表れる。推定部１２は、取得部１１が取得した計測値を用いて、制御量に生じた外乱ｑの大きさを推定する。出力部１３は、外乱ｑの推定値（以下、外乱ｑ又は推定値ｑと記載する。）を制御装置２０へ出力する。
【0014】
制御装置２０は、外乱推定装置１０から外乱ｑを取得し、制御対象３０のセンサが計測した計測値を取得する。制御装置２０は、外乱ｑと計測値に基づいて、制御対象３０を制御する。制御対象３０とは、例えば、ゴミ焼却炉、発電プラント、化学プラントなどの各種プラント、船舶、ガスタービン、蒸気タービン、コンプレッサーなどの各種機械などである。以下、ゴミ焼却炉を制御対象３０の一例として取り上げて本開示に係る外乱推定処理の説明を行うが、各実施形態の適用対象はゴミ焼却炉に限定されない。
【0015】
制御装置２０が取得する計測値には制御量が含まれる。制御装置２０は、例えば、制御量が一定となるように制御対象３０を制御する。例えば、ゴミ焼却炉では発生する蒸気の流量を一定にして運転することが理想である。発生する蒸気流量が一定であれば、焼却炉として最大能力で蒸気を持続的に発生させることができるので、ゴミの焼却量、すなわち処理量と、発電による売電収入を最大化することができる。しかし、市中から回収したゴミは多様であり、例えば、時間的に同じ割合で炉にゴミを供給しても蒸気流量を一定にできない。特許文献１の技術では、ゴミの水分を計測し、特許文献２の技術では、ゴミの単位質量当たりの発熱量を推定し、炉の発熱量の変動原因を推定して、廃棄物や燃焼空気の調節に利用している。蒸気流量などの制御の目的となる変数は一般に制御量と呼ばれている。そして、制御量を含む制御対象の計測値を意図に反し変動させるものは一般に外乱と呼ばれる。ゴミ焼却炉では、ゴミ中の水分の変動や、ゴミの単位質量あたりの発熱量の変動は代表的な外乱である。さらに、これらをマクロに捉えるならば、炉全体での発熱量の変動も外乱である。特許文献１、２に記載の技術は、どちらもゴミ焼却炉に関わる熱学的な知識に基づき燃焼速度の変動に関する特定の外乱を推定する技術である。これに対し、本実施形態の技術は制御対象についての例えば熱学のような事前知識が無くても、外乱ｑを推定することを可能とする。具体的には、制御対象に設けられたセンサから得られる計測値から支配的な外乱を推定する。以下では、ゴミ焼却炉を例に、推定部１２による外乱ｑの推定手順を説明する。
【0016】
図２は、第一実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
図２に外乱推定装置１０のうち、推定部１２の要部の構成を示す。推定部１２は、制御対象３０のセンサが計測した計測値を要素とするｍ行１列の計測ベクトルｙを構成する手段１２１と、計測ベクトルｙの分散共分散行列を算出する手段１２２と、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出する手段１２３と、最大特異値の特異ベクトルに基づき制御対象３０の外乱ｑを推定する手段１２４と、を備える。
図２に外乱推定装置１０のうち、推定部１２の要部の構成を示す。推定部１２は、制御対象３０のセンサが計測した計測値を要素とするｍ行１列の計測ベクトルｙを構成する手段１２１と、計測ベクトルｙの分散共分散行列を算出する手段１２２と、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出する手段１２３と、最大特異値の特異ベクトルに基づき制御対象３０の外乱ｑを推定する手段１２４と、を備える。
【0017】
（外乱の推定手順）
外乱をｑ∈Ｒ１で表す。ゴミ焼却炉では外乱ｑとして支配的なものは燃焼速度である。外乱ｑが変動すると計測値ｙ∈Ｒｍも変動する。両者の変動を、式（１）のように、１次式で近似する。
ｙ＝ｃ１×ｑ ・・・（１）
式（１）のｃ１はｍ行１列の係数ベクトルである。ｃ１は、外乱ｑが増加するときの計測値ベクトルｙの応答を表している。ゴミ焼却炉では、燃焼速度が増加すると、蒸気流量は増加し、燃焼室温度は増加し、排ガスの酸素濃度は減少する。ｃ１は増加や減少を定量化する係数ベクトルであり、以下のようにして値を決定する。まず、計測値を列要素とする計測値ベクトルｙの分散共分散行列Ｑ０∈Ｒｍ×ｍを、式（２）のように計算する。式（２）のＶａｒは分散の記号である。
□Ｑ０＝Ｖａｒ（ｙ）・・・（２）
次に、分散共分散行列Ｑ０を特異値分解（Ｓｉｎｇｌａｒ Ｖａｌｕｅ Ｄｅｃｏｍｐｏｓｉｔｉｏｎ，ＳＶＤ）し、式（３）の特異ベクトルｕｉ（ｉ＝１，２，・・・，ｍ）∈Ｒｍと特異値σ２ ｉ（ｉ＝１，２，…，ｍ）∈Ｒ＋を求める。ここに特異値分解の慣例に従って特異値は大きさの順にソートする。すなわちσ２ １が最大特異値、σ２ ｍが最小特異値である。右肩の記号Ｔは、行列の転置を表す。
外乱をｑ∈Ｒ１で表す。ゴミ焼却炉では外乱ｑとして支配的なものは燃焼速度である。外乱ｑが変動すると計測値ｙ∈Ｒｍも変動する。両者の変動を、式（１）のように、１次式で近似する。
ｙ＝ｃ１×ｑ ・・・（１）
式（１）のｃ１はｍ行１列の係数ベクトルである。ｃ１は、外乱ｑが増加するときの計測値ベクトルｙの応答を表している。ゴミ焼却炉では、燃焼速度が増加すると、蒸気流量は増加し、燃焼室温度は増加し、排ガスの酸素濃度は減少する。ｃ１は増加や減少を定量化する係数ベクトルであり、以下のようにして値を決定する。まず、計測値を列要素とする計測値ベクトルｙの分散共分散行列Ｑ０∈Ｒｍ×ｍを、式（２）のように計算する。式（２）のＶａｒは分散の記号である。
□Ｑ０＝Ｖａｒ（ｙ）・・・（２）
次に、分散共分散行列Ｑ０を特異値分解（Ｓｉｎｇｌａｒ Ｖａｌｕｅ Ｄｅｃｏｍｐｏｓｉｔｉｏｎ，ＳＶＤ）し、式（３）の特異ベクトルｕｉ（ｉ＝１，２，・・・，ｍ）∈Ｒｍと特異値σ２ ｉ（ｉ＝１，２，…，ｍ）∈Ｒ＋を求める。ここに特異値分解の慣例に従って特異値は大きさの順にソートする。すなわちσ２ １が最大特異値、σ２ ｍが最小特異値である。右肩の記号Ｔは、行列の転置を表す。
【0018】
【数1】
【0019】
次に、外乱ρｉ（ｉ＝１，２，…，ｍ）∈Ｒ１があるものとして、式（４）のように、計測値ベクトルｙの変動を特異ベクトルｕと未知の外乱ρとで表す。未知の外乱ρの要素は線形独立、すなわちｉ≠ｊならばＣｏｖ（ρｉ，ρｊ）＝０である。Ｃｏｖは共分散の記号である。ｕは計測値ベクトルｙの分散共分散行列Ｑ０を特異値分解すると値が分るので、未知の外乱ρは計測値ベクトルｙから値を計算することができる。
【0020】
【数2】
【0021】
分散共分散行列Ｑ０の対称性から、特異ベクトルｕには、式（５）の性質がある。
【0022】
【数3】
【0023】
従って、式（４）の両辺に左からｕＴを掛けて外乱ρを陽に定義することができる。
【0024】
【数4】
【0025】
外乱ρの分散共分散行列を計算すると式（７）を得る。
【0026】
【数5】
【0027】
未知の外乱の第１要素であるρ１の分散は、式（７）が示すように、最大特異値σ１
２である。したがって、計測値ベクトルｙの分散はρ１に起因する成分が最も大きいといえる。なぜなら、特異値の性質から、
Ｖａｒ（ｙ１）＋Ｖａｒ（ｙ２）＋・・・＋Ｖａｒ（ｙｍ）＝
σ１ ２＋σ２ ２＋・・・＋σｍ ２ ・・・（８）
が成り立ち、特にσ１ ２＞＞σ２ ２＋σ３ ２＋・・・＋σｍ ２であるならば、次式（８Ａ）のように近似され、計測値ベクトルｙの変動はρ１に支配される。
Ｖａｒ（ｙ１）＋Ｖａｒ（ｙ２）＋・・・＋Ｖａｒ（ｙｍ）＝
σ１ ２＋σ２ ２＋・・・＋σｍ ２ ・・・（８）
が成り立ち、特にσ１ ２＞＞σ２ ２＋σ３ ２＋・・・＋σｍ ２であるならば、次式（８Ａ）のように近似され、計測値ベクトルｙの変動はρ１に支配される。
【0028】
【数6】
【0029】
式（６）からρ１に関係する部分を取り出して外乱qの推定式として式（９）を得る。
【0030】
【数7】
【0031】
ゴミ焼却炉の例では、外乱ｑとして燃焼速度の変動が支配的であることが分かっている。取得部１１が計測値ｙを取得し、上述の手順によって、推定部１２が式（９）により支配的な外乱ｑを推定すると、それは燃焼速度の変動である。出力部１３は、外乱ｑを制御装置２０へ出力する。制御装置２０は、外乱ｑを燃焼速度の変動と見なして、それを相殺するように燃焼空気やごみ供給を調整する。これにより、蒸気流量を一定にしてゴミ焼却炉を運転することができる。ゴミ焼却炉に限らず、どのような制御対象であっても、支配的な外乱について見当をつけることはできる。例えば、船舶の自動操舵では潮流が支配的な外乱であり、自動車の速度制御では路面の勾配などが支配的な外乱である。これらは、運動方程式などの解析的な手段ではなく、経験から得られる知見である。
【0032】
（動作）
上記の手順を図３に示す。図３は、第一実施形態に係る外乱推定処理の一例を示す図である。まず、取得部１１が、ゴミ焼却炉に設けられたセンサが計測した蒸気流量、燃焼室温度、排ガスの酸素濃度などの計測値を取得する（ステップＳ１）。推定部１２の手段１２１は、取得部１１が取得した計測値を用いて計測ベクトルｙを構成する（ステップＳ２）。例えば、手段１２１は、蒸気流量、燃焼室温度、排ガスの酸素濃度の各計測値を要素とする計測ベクトルｙを構成する。次に推定部１２の手段１２２は、式（２）により分散共分散行列Ｑ０を算出する（ステップＳ３）。次に推定部１２の手段１２３は、式（３）により、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出する（ステップＳ４）。次に推定部１２の手段１２４は、式（９）により、外乱ｑを推定する（ステップＳ５）。推定部１２は、外乱ｑを制御装置２０へ出力する。
上記の手順を図３に示す。図３は、第一実施形態に係る外乱推定処理の一例を示す図である。まず、取得部１１が、ゴミ焼却炉に設けられたセンサが計測した蒸気流量、燃焼室温度、排ガスの酸素濃度などの計測値を取得する（ステップＳ１）。推定部１２の手段１２１は、取得部１１が取得した計測値を用いて計測ベクトルｙを構成する（ステップＳ２）。例えば、手段１２１は、蒸気流量、燃焼室温度、排ガスの酸素濃度の各計測値を要素とする計測ベクトルｙを構成する。次に推定部１２の手段１２２は、式（２）により分散共分散行列Ｑ０を算出する（ステップＳ３）。次に推定部１２の手段１２３は、式（３）により、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出する（ステップＳ４）。次に推定部１２の手段１２４は、式（９）により、外乱ｑを推定する（ステップＳ５）。推定部１２は、外乱ｑを制御装置２０へ出力する。
【0033】
上記の外乱推定処理において計測ベクトルｙ、外乱ｑ、Ｑ０、ｕ１の更新タイミングの関係を図４に示す。例えば、計測値ベクトルｙは外乱推定装置１０に新たな計測値が到来する周期ＴＡで更新し、それに合わせて、外乱ｑも周期ＴＡで更新する。また、式（２）による分散共分散行列Ｑ０や式（３）の特異値分解によるｕ１の更新は、周期ＴＢで更新する。特異値ベクトルｕ１は、分散共分散行列から値が決まるので、分散共分散行列Ｑ０を更新する周期ＴＢで計算する。更新周期は対象の特性に基づいて決定するが、一般に、ＴＡ＜＜ＴＢである。特殊なケースとして、対象の特性が一定であるならば、特異値ベクトルｕ１は事前に決めた値に固定しても良い。
【0034】
本実施形態によれば、制御対象３０にて計測された計測値に基づいて、制御量に影響する外乱ｑを推定することができる。多くの場合、制御対象３０のセンサは、制御量や制御量に影響する物理量を計測するために設けられているので、新たなセンサの追加を必要とせずに既設のセンサの計測値を用いて、外乱ｑを推定することができる。また、図４に例示したように、計測ベクトルｙが得られる周期ＴＡで、最新の計測値に基づくタイムリーな外乱ｑを推定することができる。また、制御対象３０に設けられたセンサによる計測値と上記の手順を実行すればよいので、制御対象３０の特性に依存せず、さまざまな制御対象３０の外乱推定に汎用的に適用することができる。
【0035】
＜第二実施形態＞
図５、図６を用いて第二実施形態の外乱推定装置について説明する。
第一実施形態では、支配的な外乱（例えば、燃焼速度）が事前に明らかであることを前提とした。ゴミ焼却炉では、燃焼速度が支配的な外乱であることが分かっているので、式（９）で算出した外乱ｑが燃焼速度の変動であり、燃焼速度の変動に基づき、燃焼空気やゴミの供給量を調節した。第二実施形態では、外乱を制御量の変動に換算して推定する。これにより、支配的な外乱が事前に分からない場合にも第一実施形態と同様の手法を適用できるようにする。例えば、ゴミ焼却炉の場合、制御量とは、蒸気流量であり、外乱を制御量の変動に換算して推定するとは、燃焼速度の変動をそれが及ぼす蒸気流量の変動に換算して推定するという事である。制御量に換算して推定することにより、燃焼速度の変動が支配的な外乱であるという知識が無くても外乱ｑを推定することができるばかりではなく、ゴミ焼却炉を制御するうえで、外乱ｑの大きさ（燃焼速度）を制御量である蒸気流量に換算する必要が無く、蒸気流量のまま扱って制御に流用することができるという利点が得られる。
図５、図６を用いて第二実施形態の外乱推定装置について説明する。
第一実施形態では、支配的な外乱（例えば、燃焼速度）が事前に明らかであることを前提とした。ゴミ焼却炉では、燃焼速度が支配的な外乱であることが分かっているので、式（９）で算出した外乱ｑが燃焼速度の変動であり、燃焼速度の変動に基づき、燃焼空気やゴミの供給量を調節した。第二実施形態では、外乱を制御量の変動に換算して推定する。これにより、支配的な外乱が事前に分からない場合にも第一実施形態と同様の手法を適用できるようにする。例えば、ゴミ焼却炉の場合、制御量とは、蒸気流量であり、外乱を制御量の変動に換算して推定するとは、燃焼速度の変動をそれが及ぼす蒸気流量の変動に換算して推定するという事である。制御量に換算して推定することにより、燃焼速度の変動が支配的な外乱であるという知識が無くても外乱ｑを推定することができるばかりではなく、ゴミ焼却炉を制御するうえで、外乱ｑの大きさ（燃焼速度）を制御量である蒸気流量に換算する必要が無く、蒸気流量のまま扱って制御に流用することができるという利点が得られる。
【0036】
（構成）
図５は、第二実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
第二実施形態に係る外乱推定装置１０Ａは、推定部１２に代えて推定部１２Ａを備える。第二実施形態の取得部１１が取得する計測値には、制御量が含まれる。以下では、制御量を計測値ベクトルｙの第一要素に配置するとして説明する。
第二実施形態に係る推定部１２Ａは、手段１２４に代えて、最大特異値の特異ベクトルに基づき対象物に加わる外乱を制御量の変動として推定する手段１２５を備える。推定部１２Ａは、第一実施形態と同様に、計測値ベクトルｙの分散共分散行列Ｑ０を式（２）で算出し、分散共分散行列Ｑ０を式（３）のように特異値分解する。そして、推定部１２Ａは、特異値ベクトルｕｉ（ｉ＝１，２，…，ｍ）と特異値σ２ ｉ（ｉ＝１，２，…，ｍ）を算出する。式（６）の第一行目を取り出したものが式（６Ａ）である。
図５は、第二実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
第二実施形態に係る外乱推定装置１０Ａは、推定部１２に代えて推定部１２Ａを備える。第二実施形態の取得部１１が取得する計測値には、制御量が含まれる。以下では、制御量を計測値ベクトルｙの第一要素に配置するとして説明する。
第二実施形態に係る推定部１２Ａは、手段１２４に代えて、最大特異値の特異ベクトルに基づき対象物に加わる外乱を制御量の変動として推定する手段１２５を備える。推定部１２Ａは、第一実施形態と同様に、計測値ベクトルｙの分散共分散行列Ｑ０を式（２）で算出し、分散共分散行列Ｑ０を式（３）のように特異値分解する。そして、推定部１２Ａは、特異値ベクトルｕｉ（ｉ＝１，２，…，ｍ）と特異値σ２ ｉ（ｉ＝１，２，…，ｍ）を算出する。式（６）の第一行目を取り出したものが式（６Ａ）である。
【0037】
【数8】
【0038】
ここに、ｕ１，ｊ（ｊ＝１，２，…，ｍ）は最大特異値についての特異値ベクトルｕ１の第ｊ要素である。第一要素（制御量）に対応する最大特異値が支配的であるならば、即ちσ２
１＞＞σ２
２＋σ２
３＋・・・＋σ２
ｍであるならば、計測値ベクトルｙは、式（１０Ａ）のように近似される。
【0039】
【数9】
【0040】
式（６Ａ）のξは、支配的なρ１により式（１１）のように表される。
【0041】
【数10】
【0042】
ここにξ｜ρ１は、ρ１を入力条件としたときのξであることを示す。
同様に、制御量は、ρ１により次のように表される。
ｙ１｜ρ１＝ｕ１１ρ１ ・・・（１２）
ここにｙ１｜ρ１は、ρ１を入力条件としたときのｙ１であることを示す。
式（１２）により、ξの値が分かると、制御量の変動に換算した外乱をｑｙ１｜ξと記すと、それは次式（１３）で表される。
同様に、制御量は、ρ１により次のように表される。
ｙ１｜ρ１＝ｕ１１ρ１ ・・・（１２）
ここにｙ１｜ρ１は、ρ１を入力条件としたときのｙ１であることを示す。
式（１２）により、ξの値が分かると、制御量の変動に換算した外乱をｑｙ１｜ξと記すと、それは次式（１３）で表される。
【0043】
【数11】
【0044】
（動作）
上記の手順を図６に示す。図６は、第二実施形態に係る外乱推定処理の一例を示す図である。まず、取得部１１が、ゴミ焼却炉に設けられたセンサが計測した計測値を取得する（ステップＳ１）。計測値には制御量が含まれる。次に推定部１２Ａの手段１２１は、計測ベクトルｙを構成する（ステップＳ２）。手段１２１は、制御量を第一要素として、計測ベクトルｙを構成する。次に推定部１２Ａの手段１２２は、式（２）により分散共分散行列Ｑ０を算出する（ステップＳ３）。次に推定部１２Ａの手段１２３は、式（３）により、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出する（ステップＳ４）。次に推定部１２Ａの手段１２５は、式（１３）により、制御量の変動に換算した外乱ｑｙ１｜ξを推定する（ステップＳ６）。推定部１２Ａは、外乱ｑｙ１｜ξを制御装置２０へ出力する。
上記の手順を図６に示す。図６は、第二実施形態に係る外乱推定処理の一例を示す図である。まず、取得部１１が、ゴミ焼却炉に設けられたセンサが計測した計測値を取得する（ステップＳ１）。計測値には制御量が含まれる。次に推定部１２Ａの手段１２１は、計測ベクトルｙを構成する（ステップＳ２）。手段１２１は、制御量を第一要素として、計測ベクトルｙを構成する。次に推定部１２Ａの手段１２２は、式（２）により分散共分散行列Ｑ０を算出する（ステップＳ３）。次に推定部１２Ａの手段１２３は、式（３）により、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出する（ステップＳ４）。次に推定部１２Ａの手段１２５は、式（１３）により、制御量の変動に換算した外乱ｑｙ１｜ξを推定する（ステップＳ６）。推定部１２Ａは、外乱ｑｙ１｜ξを制御装置２０へ出力する。
【0045】
本実施形態によれば、第一実施形態の効果に加え、外乱についての知識が無くても、制御対象３０にて計測された計測値に基づいて、外乱ｑを制御量の変動に換算した推定値ｑｙ１｜ξを算出することができる。例えば、第一実施形態と同様、計測値ベクトルｙと外乱ｑｙ１｜ξは周期ＴＡで更新する。一方、特異値ベクトルｕ１は、分散共分散行列Ｑ０から値が決まるので、分散共分散行列Ｑ０を更新する周期ＴＢで計算する。また、対象の特性が一定であるならば、特異値ベクトルｕ１は事前に決めた値に固定しても良い。
【0046】
＜第三実施形態＞
図７、図８を用いて第三実施形態の外乱推定装置について説明する。
第三実施形態では、外乱を制御量の変動として推定した推定値ｑｙ１｜ξを制御量の実測値ｙ１と比較して外乱推定の精度を判定する。制御装置２０は、精度が悪ければ外乱を相殺する調節を取りやめる。ゴミ焼却炉の例では、推定した蒸気流量の変動と実際の蒸気流量との差が小さければ、推定した蒸気流量の変動に基づき燃焼空気やごみ供給を調節して変動を相殺する。一方、推定した蒸気流量の変動と実際の蒸気流量との差が大きければ、調節を取りやめる。
図７、図８を用いて第三実施形態の外乱推定装置について説明する。
第三実施形態では、外乱を制御量の変動として推定した推定値ｑｙ１｜ξを制御量の実測値ｙ１と比較して外乱推定の精度を判定する。制御装置２０は、精度が悪ければ外乱を相殺する調節を取りやめる。ゴミ焼却炉の例では、推定した蒸気流量の変動と実際の蒸気流量との差が小さければ、推定した蒸気流量の変動に基づき燃焼空気やごみ供給を調節して変動を相殺する。一方、推定した蒸気流量の変動と実際の蒸気流量との差が大きければ、調節を取りやめる。
【0047】
（構成）
図７は、第三実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
第三実施形態に係る外乱推定装置１０Ｂは、推定部１２に代えて推定部１２Ｂを備える。推定部１２Ｂは、第二実施形態で説明した処理により推定値ｑｙ１｜ξを算出し、推定値ｑｙ１｜ξの精度を判定する。推定部１２Ｂは、第二実施形態に係る推定部１２Ａの構成に加え、推定した制御量の変動と実際の制御量との誤差を算出する手段１２６と、計算された誤差の分散を算出する手段１２７と、誤差の分散から外乱推定の確からしさを判定する手段１２８と、を備える。第二実施形態の取得部１１が取得する計測値には、制御量が含まれる。以下では、制御量を計測値ベクトルｙの第一要素に配置するとして説明する。第三実施形態の出力部１３は、推定値ｑｙ１｜ξに加え、外乱推定精度の判定結果（調整制限指令）を出力する。
図７は、第三実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
第三実施形態に係る外乱推定装置１０Ｂは、推定部１２に代えて推定部１２Ｂを備える。推定部１２Ｂは、第二実施形態で説明した処理により推定値ｑｙ１｜ξを算出し、推定値ｑｙ１｜ξの精度を判定する。推定部１２Ｂは、第二実施形態に係る推定部１２Ａの構成に加え、推定した制御量の変動と実際の制御量との誤差を算出する手段１２６と、計算された誤差の分散を算出する手段１２７と、誤差の分散から外乱推定の確からしさを判定する手段１２８と、を備える。第二実施形態の取得部１１が取得する計測値には、制御量が含まれる。以下では、制御量を計測値ベクトルｙの第一要素に配置するとして説明する。第三実施形態の出力部１３は、推定値ｑｙ１｜ξに加え、外乱推定精度の判定結果（調整制限指令）を出力する。
【0048】
推定部１２Ｂは、推定した制御量（例えば、蒸気流量）の変動の推定値ｑｙ１｜ξと実際の制御量ｙ１とを取得し、変動ｑｙ１｜ξと制御量ｙ１の差の分散を次式（１４）で計算する。
Ｊ＝Ｖａｒ（ｙ１－ｑｙ１｜ξ） ・・・（１４）
推定部１２Ｂは、分散Ｊが予め定めた閾値よりも小さければ調整制限指令にＯＦＦを設定し、分散Ｊが予め定めた閾値よりも大きければ調整制限指令にＯＮを設定する。出力部１３は、推定部１２Ｂが算出した推定値ｑｙ１｜ξと調整制限指令を制御装置２０へ出力する。制御装置２０は、調整制限指令がＯＦＦ（分散Ｊが閾値よりも小さい）であれば外乱を相殺する調節を実施する。例えば、ゴミ焼却炉の場合、蒸気流量の変動（推定値ｑｙ１｜ξ）を抑えるようにゴミ供給や燃焼空気の供給量を調整する。調整制限指令がＯＮ（分散Ｊが閾値よりも大きい）であれば外乱を相殺する調整を実施しない。
Ｊ＝Ｖａｒ（ｙ１－ｑｙ１｜ξ） ・・・（１４）
推定部１２Ｂは、分散Ｊが予め定めた閾値よりも小さければ調整制限指令にＯＦＦを設定し、分散Ｊが予め定めた閾値よりも大きければ調整制限指令にＯＮを設定する。出力部１３は、推定部１２Ｂが算出した推定値ｑｙ１｜ξと調整制限指令を制御装置２０へ出力する。制御装置２０は、調整制限指令がＯＦＦ（分散Ｊが閾値よりも小さい）であれば外乱を相殺する調節を実施する。例えば、ゴミ焼却炉の場合、蒸気流量の変動（推定値ｑｙ１｜ξ）を抑えるようにゴミ供給や燃焼空気の供給量を調整する。調整制限指令がＯＮ（分散Ｊが閾値よりも大きい）であれば外乱を相殺する調整を実施しない。
【0049】
（動作）
上記の手順を図８に示す。図８は、第三実施形態に係る外乱推定処理の一例を示す図である。まず、第二実施形態で説明した処理により、推定部１２Ｂは、制御量の変動に換算した外乱ｑｙ１｜ξを推定する（ステップＳ１０）。次に手段１２６が、推定した制御量の変動、つまり外乱ｑｙ１｜ξと制御量ｙ１の誤差を算出する（ステップＳ１１）。次に手段１２７が、ステップＳ１１で算出された誤差の分散Ｊを算出する（ステップＳ１２）。次に手段１２８が、ステップＳ１２で算出された誤差の分散Ｊに基づいて、外乱ｑｙ１｜ξの推定の確からしさを判定する（ステップＳ１３）。手段１２８は、分散Ｊが所定の閾値よりも大きい場合、推定は不確であると判定し、分散Ｊが所定の閾値よりも小さい場合、推定は確かであると判定する。図７に示すように、この判定にはヒステリシス幅が設けられていてもよい。ヒステリシス幅を設けることで、制御量ｙ１の計測誤差や変動を吸収し、安定した制御を行うことができる。手段１２８は、推定が不確であると判定すると、調整制御指令にＯＮを設定し、推定は確かであると判定すると、調整制御指令にＯＦＦを設定する。推定部１２Ｂは、外乱の推定値ｑｙ１｜ξと調整制限指令（ＯＮ又はＯＦＦ）を制御装置２０へ出力する（ステップＳ１４）。
上記の手順を図８に示す。図８は、第三実施形態に係る外乱推定処理の一例を示す図である。まず、第二実施形態で説明した処理により、推定部１２Ｂは、制御量の変動に換算した外乱ｑｙ１｜ξを推定する（ステップＳ１０）。次に手段１２６が、推定した制御量の変動、つまり外乱ｑｙ１｜ξと制御量ｙ１の誤差を算出する（ステップＳ１１）。次に手段１２７が、ステップＳ１１で算出された誤差の分散Ｊを算出する（ステップＳ１２）。次に手段１２８が、ステップＳ１２で算出された誤差の分散Ｊに基づいて、外乱ｑｙ１｜ξの推定の確からしさを判定する（ステップＳ１３）。手段１２８は、分散Ｊが所定の閾値よりも大きい場合、推定は不確であると判定し、分散Ｊが所定の閾値よりも小さい場合、推定は確かであると判定する。図７に示すように、この判定にはヒステリシス幅が設けられていてもよい。ヒステリシス幅を設けることで、制御量ｙ１の計測誤差や変動を吸収し、安定した制御を行うことができる。手段１２８は、推定が不確であると判定すると、調整制御指令にＯＮを設定し、推定は確かであると判定すると、調整制御指令にＯＦＦを設定する。推定部１２Ｂは、外乱の推定値ｑｙ１｜ξと調整制限指令（ＯＮ又はＯＦＦ）を制御装置２０へ出力する（ステップＳ１４）。
【0050】
ステップＳ１３で、調整制限指令がＯＮとなる状況が続くような場合、推定部１２Ｂは、図４を用いて説明したＱ０やＵ１の更新頻度を上げる等して、精度の向上を試みてもよい。それでも精度が向上しない場合には、推定値ｑｙ１｜ξの算出に用いた計測値ベクトルｙの第二要素以下として用いる計測値を選定し直してもよい。
【0051】
本実施形態によれば、第二実施形態の効果に加え、外乱推定値ｑｙ１｜ξの精度を確認しながら、制御対象３０の制御を行うことができる。また、調整制限指令の値に基づいて、外乱の推定値ｑｙ１｜ξに基づく調整の実行と停止を自動的に切り替える機能を制御装置２０に組み込むことにより、制御装置２０の制御精度を確保することができる。
【0052】
＜第四実施形態＞
図９～図１１を用いて第四実施形態の外乱推定装置１０Ｃについて説明する。
外乱が起きてから制御量や計測値の変動となって表れるまでには、一般的に、時間遅れがある。例えば、ゴミ焼却炉では、外乱として燃焼速度が変わったあと、その影響が炉内の温度に表れるには、例えば１０秒の時間がかかり、蒸気流量の変動として表れるには、例えば３００秒の時間がかかるとする。すなわち、例えば時刻ｔに燃焼速度が変わっても、その影響が炉内温度に表れる時刻はｔ＋１０であり、蒸気流量の変動に表れる時刻はｔ＋３００である。したがって、この場合、燃焼速度の変動について炉内温度と蒸気流量には２９０秒の応答の時間差がある。それが分かっているならば、両者には２９０秒の時間差をつけて計測値ベクトルｙを構成すべきである。例えば、計測値ベクトルｙが蒸気流量と炉内温度を含むとすると、時刻ｔの計測値ベクトルの要素には、時刻ｔの蒸気流量と、時刻ｔ－２９０の炉内温度とを用いる。遅れ時間の値が不明な場合は、遅れ時間の値を変えて複数設定してもよい。例えば、前述の例では、計測値ベクトルｙの要素は、時刻ｔの蒸気流量と、時刻ｔ－２９０の炉内温度とした。これに、時刻ｔ－３５０の炉内温度、時刻ｔ－３２０の炉内温度、時刻ｔ－２６０の炉内温度、などを計測値ベクトルｙの要素として加えてもよい。
図９～図１１を用いて第四実施形態の外乱推定装置１０Ｃについて説明する。
外乱が起きてから制御量や計測値の変動となって表れるまでには、一般的に、時間遅れがある。例えば、ゴミ焼却炉では、外乱として燃焼速度が変わったあと、その影響が炉内の温度に表れるには、例えば１０秒の時間がかかり、蒸気流量の変動として表れるには、例えば３００秒の時間がかかるとする。すなわち、例えば時刻ｔに燃焼速度が変わっても、その影響が炉内温度に表れる時刻はｔ＋１０であり、蒸気流量の変動に表れる時刻はｔ＋３００である。したがって、この場合、燃焼速度の変動について炉内温度と蒸気流量には２９０秒の応答の時間差がある。それが分かっているならば、両者には２９０秒の時間差をつけて計測値ベクトルｙを構成すべきである。例えば、計測値ベクトルｙが蒸気流量と炉内温度を含むとすると、時刻ｔの計測値ベクトルの要素には、時刻ｔの蒸気流量と、時刻ｔ－２９０の炉内温度とを用いる。遅れ時間の値が不明な場合は、遅れ時間の値を変えて複数設定してもよい。例えば、前述の例では、計測値ベクトルｙの要素は、時刻ｔの蒸気流量と、時刻ｔ－２９０の炉内温度とした。これに、時刻ｔ－３５０の炉内温度、時刻ｔ－３２０の炉内温度、時刻ｔ－２６０の炉内温度、などを計測値ベクトルｙの要素として加えてもよい。
【0053】
図９は、第四実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
図９に第二実施形態と組み合わせた場合の第四実施形態の構成を示す。第四実施形態の推定部１２Ｃは、第二実施形態の構成に加え、遅れ時間補正手段１２９を備える。遅れ時間補正手段１２９は、計測値ベクトルｙの各要素について、過去に計測された計測値と、遅れ時間とを対応付けて記憶している。例えば、計測値ベクトルｙが蒸気流量と炉内温度を含む場合、遅れ時間補正手段１２９は、時刻ｔにおいて、取得部１１によって取得された２９０秒前の炉内温度の計測値と、炉内温度が蒸気流量よりも２９０秒の遅れがあることと、を対応付けて記憶している。遅れ時間補正手段１２９は、計測値ベクトルｙを取得し、内部に記憶した遅れ時間の値により計測ベクトルの一つ一つの要素の遅れ時間を補正して、遅れ時間補正後の計測値ベクトルｙ～を出力する。第四実施形態では、推定部１２Ｃは、計測ベクトルｙの代わりに遅れ時間補正後の計測ベクトルｙ～に基づいて外乱を推定する。
図９に第二実施形態と組み合わせた場合の第四実施形態の構成を示す。第四実施形態の推定部１２Ｃは、第二実施形態の構成に加え、遅れ時間補正手段１２９を備える。遅れ時間補正手段１２９は、計測値ベクトルｙの各要素について、過去に計測された計測値と、遅れ時間とを対応付けて記憶している。例えば、計測値ベクトルｙが蒸気流量と炉内温度を含む場合、遅れ時間補正手段１２９は、時刻ｔにおいて、取得部１１によって取得された２９０秒前の炉内温度の計測値と、炉内温度が蒸気流量よりも２９０秒の遅れがあることと、を対応付けて記憶している。遅れ時間補正手段１２９は、計測値ベクトルｙを取得し、内部に記憶した遅れ時間の値により計測ベクトルの一つ一つの要素の遅れ時間を補正して、遅れ時間補正後の計測値ベクトルｙ～を出力する。第四実施形態では、推定部１２Ｃは、計測ベクトルｙの代わりに遅れ時間補正後の計測ベクトルｙ～に基づいて外乱を推定する。
【0054】
時間補正について、図１０を用いてさらに具体的に説明する。図１０は、外乱の発生時刻とその影響が計測値に現れるまでの時間差の一例を示す図である。計測ベクトルｙにはｍ個の要素があり、要素ごとの遅れ時間をτｉ（ｉ＝１，２，…，ｍ）とすると、外乱が発生してから計測ベクトルｙの要素に応答が表れるまでの時間は図１０のように表される。ここで、説明を簡単にするために、計測ベクトルｙ～の要素は遅れ時間の大きい順に並んでいるものとする。制御量は、ゴミ焼却炉などの設備の最終的な出力であるので、遅れ時間が計測ベクトルｙ～の要素のなかで最大であることは一般的である。制御量の変動を補償するための外乱推定だから、外乱の推定に制御量より応答が遅いものを使うことは意味がない。したがって、計測ベクトルｙの要素のなかで、制御量が最も遅れ時間が大きくなるのは当然である。時刻ｔにおいて式（６Ａ）のξを算出することを考える。ξの算出に用いるのは、図１０に示すように時刻ｔ以前の情報である。一方、制御量ｙ１に外乱の影響が表れるのは、なお、τΔ＝τ２－τ１後になる。遅れ時間補正後の計測ベクトルｙ～は、遅れ時間を用いて次式（１５）のように表すことができる。
【0055】
【数12】
【0056】
この、遅れ時間補正後の計測ベクトルｙ～を用いて分散共分散行列Ｑ０を求め、さらにその特異ベクトルｕを計算する。推定部１２Ｃは、時刻ｔの時点で｛ｙ～
２，ｙ～
３，…，ｙ～
ｍ｝は現在値または過去値であるので、これらを用いて式（１１）により、ξ｜ρ１を算出し、式（１３）を用いて外乱ｑｙ＾１｜ξを算出する。外乱ｑｙ＾１｜ξは、時刻ｔの時点で予測した時刻ｔ＋τΔにおける制御量（例えば、蒸気流量）の変動の予測である。将来の値が分かるので、ゴミ焼却炉の運転操作盤などに表示すれば運転操作の助けになる。
【0057】
（動作）
上記の手順を図１１に示す。図１１は、第四実施形態に係る外乱推定処理の一例を示す図である。まず、取得部１１が、ゴミ焼却炉に設けられたセンサが計測した蒸気流量、燃焼室温度、排ガスの酸素濃度などの計測値を取得する（ステップＳ１）。次に推定部１２Ｃの手段１２１が、計測ベクトルｙを構成する（ステップＳ２）。次に、推定部１２Ｃの遅れ時間補正手段１２９が、計測ベクトルｙを取得し、各要素の遅れ時間を補正した、遅れ時間補正後の計測ベクトルｙ～を出力する（ステップＳ７）。以降は、第二実施形態と同様である。つまり、推定部１２Ｃの手段１２２は、式（２）により分散共分散行列Ｑ０を算出する（ステップＳ３）。次に推定部１２Ｃの手段１２３は、式（３）により、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出する（ステップＳ４）。次に推定部１２Ｃの手段１２５は、式（１３）により、制御量の変動に換算した外乱ｑｙ＾１｜ξを推定する（ステップＳ６）。推定部１２Ａは、外乱の推定値（予測値）ｑｙ＾１｜ξを制御装置２０へ出力する。
上記の手順を図１１に示す。図１１は、第四実施形態に係る外乱推定処理の一例を示す図である。まず、取得部１１が、ゴミ焼却炉に設けられたセンサが計測した蒸気流量、燃焼室温度、排ガスの酸素濃度などの計測値を取得する（ステップＳ１）。次に推定部１２Ｃの手段１２１が、計測ベクトルｙを構成する（ステップＳ２）。次に、推定部１２Ｃの遅れ時間補正手段１２９が、計測ベクトルｙを取得し、各要素の遅れ時間を補正した、遅れ時間補正後の計測ベクトルｙ～を出力する（ステップＳ７）。以降は、第二実施形態と同様である。つまり、推定部１２Ｃの手段１２２は、式（２）により分散共分散行列Ｑ０を算出する（ステップＳ３）。次に推定部１２Ｃの手段１２３は、式（３）により、分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出する（ステップＳ４）。次に推定部１２Ｃの手段１２５は、式（１３）により、制御量の変動に換算した外乱ｑｙ＾１｜ξを推定する（ステップＳ６）。推定部１２Ａは、外乱の推定値（予測値）ｑｙ＾１｜ξを制御装置２０へ出力する。
【0058】
第四実施形態によれば、第二実施形態の効果に加え、一つ一つの計測値において外乱の影響が表れるまでの遅れ時間を補正する為、外乱推定精度を向上することができる。第四実施形態は、第二実施形態だけではなく、第一実施形態や第三実施形態と組み合わせることができる。また、第二実施形態、第三実施形態と組み合わせた場合、第四実施形態によれば、将来の制御量を予測することができる。
【0059】
＜第五実施形態＞
図１２～図１３を用いて第五実施形態の外乱推定装置１０Ｄについて説明する。
第四実施形態では、外乱による制御量の変動の予測値を推定した。第五実施形態では、第四実施形態と第三実施形態を組み合わせて予測の精度を判定する。予測精度が悪いことがわかれば、誤差が悪影響を及ぼさないように、制御装置２０は、外乱を相殺する調節を取りやめる。例えば、ゴミ焼却炉の例では、予測した蒸気流量の変動と実際の蒸気流量との差が小さければ、制御装置２０は、予測した蒸気流量の変動に基づき燃焼空気やごみ供給を調節して変動を相殺する。一方、予測した蒸気流量の変動と実際の蒸気流量との差が大きければ、制御装置２０は、調節を取りやめる。
図１２～図１３を用いて第五実施形態の外乱推定装置１０Ｄについて説明する。
第四実施形態では、外乱による制御量の変動の予測値を推定した。第五実施形態では、第四実施形態と第三実施形態を組み合わせて予測の精度を判定する。予測精度が悪いことがわかれば、誤差が悪影響を及ぼさないように、制御装置２０は、外乱を相殺する調節を取りやめる。例えば、ゴミ焼却炉の例では、予測した蒸気流量の変動と実際の蒸気流量との差が小さければ、制御装置２０は、予測した蒸気流量の変動に基づき燃焼空気やごみ供給を調節して変動を相殺する。一方、予測した蒸気流量の変動と実際の蒸気流量との差が大きければ、制御装置２０は、調節を取りやめる。
【0060】
（構成）
図１２は、第五実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
第五実施形態に係る外乱推定装置１０Ｄは、推定部１２に代えて推定部１２Ｄを備える。推定部１２Ｄは、以下で述べる第四実施形態の説明に類似の処理により推定値ｑｙ＾１｜ξを算出し、実際の制御量（例えば、蒸気流量）ｙ１との差の分散に基づいて、推定値ｑｙ＾１｜ξの精度を判定する。推定部１２Ｄは、第四実施形態に係る推定部１２Ｃの構成と、第三実施形態に係る手段１２６～１２８と、に加え、制御量に外乱の影響が表れるまでの遅れ時間と外乱の予測値とに基づく制御量の変動から制御量を推定する手段１３０と、を備える。第二実施形態の取得部１１が取得する計測値には、制御量が含まれる。第五実施形態の出力部１３は、制御量の予測値ｙ＾ １と制御量の予測値の予測精度の判定結果（調整制限指令）を出力する。
図１２は、第五実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
第五実施形態に係る外乱推定装置１０Ｄは、推定部１２に代えて推定部１２Ｄを備える。推定部１２Ｄは、以下で述べる第四実施形態の説明に類似の処理により推定値ｑｙ＾１｜ξを算出し、実際の制御量（例えば、蒸気流量）ｙ１との差の分散に基づいて、推定値ｑｙ＾１｜ξの精度を判定する。推定部１２Ｄは、第四実施形態に係る推定部１２Ｃの構成と、第三実施形態に係る手段１２６～１２８と、に加え、制御量に外乱の影響が表れるまでの遅れ時間と外乱の予測値とに基づく制御量の変動から制御量を推定する手段１３０と、を備える。第二実施形態の取得部１１が取得する計測値には、制御量が含まれる。第五実施形態の出力部１３は、制御量の予測値ｙ＾ １と制御量の予測値の予測精度の判定結果（調整制限指令）を出力する。
【0061】
第５の実施形態では、第４の実施形態で説明した遅れ時間補正後の計測値ベクトルを式（１６）のように構成する。
【0062】
【数13】
【0063】
計測ベクトルｚ～は、計測ベクトルｙ～の第2要素に制御量ｙ１（ｔ）を追加している。計測ベクトルｚ～の分散共分散行列Ｑ０を求め、さらに、その特異ベクトルｕを計算する。時刻ｔにおいて、｛ｚ～２，ｚ～３，…，ｚ～ｍ＋１｝は現在値または過去値であるので、これらを用いて式（１７）を用いてξを計算する。
【0064】
【数14】
【0065】
ξから、式（１８）のように、ｚ～１すなわち時刻ｔ＋ｙ１の予測値が得られる。
【0066】
【数15】
【0067】
外乱の予測値ｑｚ～１|ξ（ｔ）はｑｙ＾１|ξ（ｔ）であり、時刻ｔにおける制御量ｙ１（ｔ＋τΔ）の予測値である。ｙ１（ｔ＋τΔ）は予測値だから、実際の計測値と区別するためにｙ＾１（ｔ＋τΔ）と記し、時刻ｔと時刻ｔ＋τΔの間での予測値の増分をτΔ×予測値の時間微分で近似すると（１９）式が導かれる。
【0068】
【数16】
【0069】
式（１９）から、式（２０）のように、予測値の時間的な変化を表す微分方程式が得られる。
【0070】
【数17】
【0071】
式（２０）を時間について数値積分することにより、外乱の予測値ｑｙ＾１｜ξ（ｔ）を基に現時刻ｔの制御量の推定値ｙ＾１（ｔ）が分る。実際の計算では、式（２０）は次式（２１）のように時定数がτΔ，ゲインが１の１次遅れフィルタにより簡便に演算できる。
【0072】
【数18】
【0073】
手段１３０は、式（２１）により、時刻ｔにおいて時刻ｔ＋τΔにおける外乱を予測した値ｑｙ＾１｜ξ（ｔ）に基づき時刻ｔの制御量を推定する。手段１２７は、式（２１）によって算出されたｙ＾１（ｔ）と実際の制御量の計測値ｙ１（ｔ）とに基づいて、第三実施形態と同様に次式（２２）を用いて推定値と計測値（実測値）の差の分散を計算する。
Ｊ ＝ Ｖａｒ（ｙ１（ｔ）－ｙ＾１（ｔ）） ・・・（２２）
Ｊ ＝ Ｖａｒ（ｙ１（ｔ）－ｙ＾１（ｔ）） ・・・（２２）
【0074】
（動作）
図１３は、第五実施形態に係る外乱推定処理の一例を示す図である。
第四実施形態で説明した処理により、推定部１２Ｄは、時刻ｔ＋τΔにおける制御量に換算した外乱ｑｙ＾１｜ξを予測する（ステップＳ２０）。次に手段１３０が、時刻ｔにおける制御量を推定する（ステップＳ２１）。上述のように、外乱の予測値ｑｙ＾１|ξ（ｔ）は、時刻ｔ＋τΔの制御量ｙ＾１（ｔ＋τΔ）を表している。手段１３０は、式（２１）によって、時間を巻き戻す計算を行い、制御量ｙ＾１（ｔ＋τΔ）から時刻ｔにおける制御量ｙ＾１（ｔ）を推定する。次に手段１２６が、時刻ｔにおける制御量の推定値と実際の制御量との誤差を算出する（ステップＳ２２）。次に手段１２７が、ステップＳ２２で算出された誤差の分散Ｊを算出する（ステップＳ２３）。手段１２７は、式（２２）によって分散Ｊを算出する。次に手段１２８が、ステップＳ２３で算出された誤差の分散Ｊに基づいて、制御量ｙ＾１（ｔ）の推定の確からしさを判定する（ステップＳ２４）。手段１２８は、分散Ｊが所定の閾値よりも大きい場合、推定は不確であると判定し、分散Ｊが所定の閾値よりも小さい場合、推定は確かであると判定する。図１２に示すように、この判定にはヒステリシス幅が設けられていてもよい。ヒステリシス幅を設けることで、制御量ｙ１の計測誤差や変動を吸収し、安定した制御を行うことができる。手段１２８は、予測が不確であると判定すると、調整制御指令にＯＮを設定し、推定は確かであると判定すると、調整制御指令にＯＦＦを設定する。推定部１２Ｄは、制御量の予測値ｙ＾１と調整制限指令（ＯＮ又はＯＦＦ）を制御装置２０へ出力する（ステップＳ２５）。
図１３は、第五実施形態に係る外乱推定処理の一例を示す図である。
第四実施形態で説明した処理により、推定部１２Ｄは、時刻ｔ＋τΔにおける制御量に換算した外乱ｑｙ＾１｜ξを予測する（ステップＳ２０）。次に手段１３０が、時刻ｔにおける制御量を推定する（ステップＳ２１）。上述のように、外乱の予測値ｑｙ＾１|ξ（ｔ）は、時刻ｔ＋τΔの制御量ｙ＾１（ｔ＋τΔ）を表している。手段１３０は、式（２１）によって、時間を巻き戻す計算を行い、制御量ｙ＾１（ｔ＋τΔ）から時刻ｔにおける制御量ｙ＾１（ｔ）を推定する。次に手段１２６が、時刻ｔにおける制御量の推定値と実際の制御量との誤差を算出する（ステップＳ２２）。次に手段１２７が、ステップＳ２２で算出された誤差の分散Ｊを算出する（ステップＳ２３）。手段１２７は、式（２２）によって分散Ｊを算出する。次に手段１２８が、ステップＳ２３で算出された誤差の分散Ｊに基づいて、制御量ｙ＾１（ｔ）の推定の確からしさを判定する（ステップＳ２４）。手段１２８は、分散Ｊが所定の閾値よりも大きい場合、推定は不確であると判定し、分散Ｊが所定の閾値よりも小さい場合、推定は確かであると判定する。図１２に示すように、この判定にはヒステリシス幅が設けられていてもよい。ヒステリシス幅を設けることで、制御量ｙ１の計測誤差や変動を吸収し、安定した制御を行うことができる。手段１２８は、予測が不確であると判定すると、調整制御指令にＯＮを設定し、推定は確かであると判定すると、調整制御指令にＯＦＦを設定する。推定部１２Ｄは、制御量の予測値ｙ＾１と調整制限指令（ＯＮ又はＯＦＦ）を制御装置２０へ出力する（ステップＳ２５）。
【0075】
ステップＳ２４で、調整制限指令がＯＮとなる状況が続くような場合、推定部１２Ｄは、図４を用いて説明したＱ０やｕ１の更新頻度を上げる等して、精度の向上を試みてもよい。それでも精度が向上しない場合には、推定値ｑｙ＾１｜ξの算出に用いた計測値ベクトルｚ～の第二要素以下として用いる計測値を選定し直したり、その遅れ時間を設定し直したりしてもよい。
【0076】
本実施形態によれば、第四実施形態の効果に加え、外乱の予測値ｑｙ＾１｜ξの予測精度を確認しながら制御対象３０の制御を行うことができる。また、調整制限指令の値に基づいて、外乱の予測値ｑｙ＾１｜ξに基づく調整の実行と停止を自動的に切り替える機能を制御装置２０に組み込むことにより、制御装置２０の制御精度を確保することができる。
【0077】
＜第六実施形態＞
図１４を用いて第六実施形態の外乱推定装置１０Ｅについて説明する。
本開示の外乱の推定では、特異ベクトルｕ１∈Ｒｍが重要である。特異ベクトルｕ１は、予め指定した遅れ時間｛τ１，τ２，・・・，τｍ｝に従って周期ＴＢごとに更新する。特異ベクトルの値は更新の度に変化する。変化した値をそのまま利用しても良いが、例えば、特異ベクトルを複数通り計算して、そのなかで最も望ましいものを利用する多数決方式にすれば、外乱推定の信頼度は一つだけのときよりも高まると期待される。遅れ時間委ついても同様である。また、例えば、制御対象３０の運転形態（起動時、定格運転時、低出力運転時）に応じて、外乱の影響が反映される計測値の種類が変化すること等が考えられる。従って、第六実施形態では、特異ベクトルの更新タイミング、遅れ時間、計測ベクトルｙを構成する計測値などを様々に異ならせて、それぞれについて第五実施形態の方法で制御量の予測精度を判定し、最も精度の高い制御量を用いて、制御対象３０の制御を行う。
図１４を用いて第六実施形態の外乱推定装置１０Ｅについて説明する。
本開示の外乱の推定では、特異ベクトルｕ１∈Ｒｍが重要である。特異ベクトルｕ１は、予め指定した遅れ時間｛τ１，τ２，・・・，τｍ｝に従って周期ＴＢごとに更新する。特異ベクトルの値は更新の度に変化する。変化した値をそのまま利用しても良いが、例えば、特異ベクトルを複数通り計算して、そのなかで最も望ましいものを利用する多数決方式にすれば、外乱推定の信頼度は一つだけのときよりも高まると期待される。遅れ時間委ついても同様である。また、例えば、制御対象３０の運転形態（起動時、定格運転時、低出力運転時）に応じて、外乱の影響が反映される計測値の種類が変化すること等が考えられる。従って、第六実施形態では、特異ベクトルの更新タイミング、遅れ時間、計測ベクトルｙを構成する計測値などを様々に異ならせて、それぞれについて第五実施形態の方法で制御量の予測精度を判定し、最も精度の高い制御量を用いて、制御対象３０の制御を行う。
【0078】
（構成）
図１４は、第六実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
第六実施形態に係る外乱推定装置１０Ｅは、第五実施形態に係る推定部１２Ｄを複数と、複数の推定部１２Ｄがそれぞれに算出した分散Ｊの中から最も小さなものを選ぶ選択手段１３１と、選択手段１３１が選択した分散Ｊに対応する制御量の予測値ｙ＾１を選択する選択手段１３２と、選択手段１３１が選択した分散Ｊに対応する調整制御指令を選択する選択手段１３３と、を備える。第六実施形態の出力部１３は、選択手段１３２が選択した制御量の予測値ｙ＾１と選択手段１３３が選択した調整制限指令を出力する。
図１４は、第六実施形態に係る外乱推定装置の要部の機能構成の一例を示す図である。
第六実施形態に係る外乱推定装置１０Ｅは、第五実施形態に係る推定部１２Ｄを複数と、複数の推定部１２Ｄがそれぞれに算出した分散Ｊの中から最も小さなものを選ぶ選択手段１３１と、選択手段１３１が選択した分散Ｊに対応する制御量の予測値ｙ＾１を選択する選択手段１３２と、選択手段１３１が選択した分散Ｊに対応する調整制御指令を選択する選択手段１３３と、を備える。第六実施形態の出力部１３は、選択手段１３２が選択した制御量の予測値ｙ＾１と選択手段１３３が選択した調整制限指令を出力する。
【0079】
例えば、図１４のように、推定部１２Ｄ－１～１２Ｄ－２を備え、それぞれが推定した外乱による制御量の変動を考慮した制御量の予測値と実際の制御量の差の分散［Ｊ］１、［Ｊ］２のうち最小のものを選択手段１３１によって選出し、分散が最小となる番号を選出しｉ＊とする（式（２３））。
【0080】
【数19】
【0081】
そして、番号ｉ＊の推定部１２Ｄからの出力を、選択手段１３２、１３３がそれぞれ制御量の変動の推定値ｙ＾ｉ＊、調整制限指令＊として選択し、出力部１３がこれらを制御装置２０へ出力する。
【0082】
本実施形態によれば、外乱推定の信頼度を向上することができる。推定部１２Ｅは、計測ベクトルｙの分散共分散行列Ｑ０の特異ベクトルを更新する周期ＴＢについて異なる値を設定した複数の推定部１２Ｄを備えていてもよいし、外乱が制御量の値に変動として現れるまでの制御量の遅れ時間について異なる値を設定した複数の推定部１２Ｄを備えていてもよいし、異なる種類の計測ベクトルｙに基づいて外乱推定を行う複数の推定部１２Ｄを備えていてもよいし、これら３つのパラメータのうちの２つまたは３つを様々に異ならせた複数の推定部１２Ｄを備えていてもよい。
【0083】
図１５は、各実施形態に係る外乱推定装置のハードウェア構成の一例を示す図である。
コンピュータ９００は、ＣＰＵ９０１、主記憶装置９０２、補助記憶装置９０３、入出力インタフェース９０４、通信インタフェース９０５を備える。
上述の外乱推定装置１０～１０Ｅは、コンピュータ９００に実装される。そして、上述した各機能は、プログラムの形式で補助記憶装置９０３に記憶されている。ＣＰＵ９０１は、プログラムを補助記憶装置９０３から読み出して主記憶装置９０２に展開し、当該プログラムに従って上記処理を実行する。また、ＣＰＵ９０１は、プログラムに従って、記憶領域を主記憶装置９０２に確保する。また、ＣＰＵ９０１は、プログラムに従って、処理中のデータを記憶する記憶領域を補助記憶装置９０３に確保する。
コンピュータ９００は、ＣＰＵ９０１、主記憶装置９０２、補助記憶装置９０３、入出力インタフェース９０４、通信インタフェース９０５を備える。
上述の外乱推定装置１０～１０Ｅは、コンピュータ９００に実装される。そして、上述した各機能は、プログラムの形式で補助記憶装置９０３に記憶されている。ＣＰＵ９０１は、プログラムを補助記憶装置９０３から読み出して主記憶装置９０２に展開し、当該プログラムに従って上記処理を実行する。また、ＣＰＵ９０１は、プログラムに従って、記憶領域を主記憶装置９０２に確保する。また、ＣＰＵ９０１は、プログラムに従って、処理中のデータを記憶する記憶領域を補助記憶装置９０３に確保する。
【0084】
なお、外乱推定装置１０～１０Ｅの全部または一部の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより各機能部による処理を行ってもよい。ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、ＣＤ、ＤＶＤ、ＵＳＢ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。また、このプログラムが通信回線によってコンピュータ９００に配信される場合、配信を受けたコンピュータ９００が当該プログラムを主記憶装置９０２に展開し、上記処理を実行しても良い。また、上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであってもよい。
【0085】
以上のとおり、本開示に係るいくつかの実施形態を説明したが、これら全ての実施形態は、例として提示したものであり、発明の範囲を限定することを意図していない。これらの実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で種々の省略、置き換え、変更を行うことができる。これらの実施形態及びその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれる。
【0086】
＜付記＞
各実施形態に記載の外乱推定装置１０～１０Ｅ、外乱推定方法およびプログラムは、例えば以下のように把握される。
各実施形態に記載の外乱推定装置１０～１０Ｅ、外乱推定方法およびプログラムは、例えば以下のように把握される。
【0087】
（１）第１の態様に係る外乱推定装置１０～１０Ｅは、制御対象が備えるセンサが計測した計測値を取得する取得部１１と、前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定する推定部１２と、を備える。
これにより、制御対象が備えるセンサが計測した計測値に基づいて、制御対象に生じる外乱を推定することができる（第一実施形態）。
これにより、制御対象が備えるセンサが計測した計測値に基づいて、制御対象に生じる外乱を推定することができる（第一実施形態）。
【0088】
（２）第２の態様に係る外乱推定装置１０Ａ～１０Ｅは、（１）の外乱推定装置１０Ａ～１０Ｅであって、前記取得部は、前記制御対象の制御の目的となる変数である制御量の計測値を取得し、前記推定部は、前記計測ベクトルにおける前記最大特異値の特異ベクトルに基づいて（式（１３））、前記外乱を前記制御量の変動として推定する。
これにより、外乱を制御量の変動に換算して推定することができ、支配的な外乱が事前に分からない場合でも外乱の大きさを推定することができる（第二実施形態）。
これにより、外乱を制御量の変動に換算して推定することができ、支配的な外乱が事前に分からない場合でも外乱の大きさを推定することができる（第二実施形態）。
【0089】
（３）第３の態様に係る外乱推定装置１０Ｂ～１０Ｃは、（１）～（２）の外乱推定装置１０Ｂ～１０Ｃであって、前記推定部は、推定した前記外乱の分散の大きさに基づいて、推定した前記外乱の確からしさを判定し、推定した前記外乱とともに前記判定の結果を出力する。
これにより、外乱推定の確からしさを判定することができる。例えば、確からしさが低い場合には、その推定結果を用いないで制御対象を制御することにより、制御の精度を確保できる（第三実施形態）。
これにより、外乱推定の確からしさを判定することができる。例えば、確からしさが低い場合には、その推定結果を用いないで制御対象を制御することにより、制御の精度を確保できる（第三実施形態）。
【0090】
（４）第４の態様に係る外乱推定装置１０Ｃ～１０Ｄは、（１）～（３）の外乱推定装置１０Ｃ～１０Ｄであって、前記推定部は、前記外乱が前記計測値に変動として現れるまでの遅れ時間を補正する補正手段、を備え、前記推定部は、前記取得部が取得した前記計測値の遅れ時間を前記補正手段によって補正し、補正後の前記計測値を要素とする前記計測ベクトルを用いて、前記外乱を推定する。
制御対象に外乱が起きてから制御量や計測値の変動となって表れるまでには一般的に時間遅れがある。第４の態様によれば、時間遅れを考慮して外乱を推定することによって、外乱の推定精度を向上することができる（第四実施形態）。
制御対象に外乱が起きてから制御量や計測値の変動となって表れるまでには一般的に時間遅れがある。第４の態様によれば、時間遅れを考慮して外乱を推定することによって、外乱の推定精度を向上することができる（第四実施形態）。
【0091】
（５）第５の態様に係る外乱推定装置１０Ｄは、（４）の外乱推定装置１０Ｄであって、前記推定部は、推定した前記外乱に基づいて、前記計測値のうちの前記遅れ時間が最も長い前記計測値の推定値を推定し、前記計測値の推定値と前記計測値の実測値の差の分散に基づいて、推定した前記外乱の確からしさを判定し、推定した前記計測値の推定値とともに前記判定の結果を出力する。
これにより、外乱を制御量の変動に換算し、制御量の実測値と比較して外乱推定の精度を判定する。これにより、外乱を用いた制御の精度を確保できる（第五実施形態）。
これにより、外乱を制御量の変動に換算し、制御量の実測値と比較して外乱推定の精度を判定する。これにより、外乱を用いた制御の精度を確保できる（第五実施形態）。
【0092】
（６）第６の態様に係る外乱推定装置１０Ｅは、（５）の外乱推定装置１０Ｄの前記推定部を複数備え、前記外乱の確からしさが最も高い前記推定部が推定した外乱を選択する。
これにより、最も精度よく推定された外乱を制御に利用することができる（第六実施形態）。
これにより、最も精度よく推定された外乱を制御に利用することができる（第六実施形態）。
【0093】
（７）第７の態様に係る外乱推定装置１０Ｅは、（６）の外乱推定装置１０Ｅであって、複数の前記推定部のそれぞれは、他の前記推定部と異なる前記計測値に基づいて前記外乱を推定するか、または、他の前記推定部と異なる前記遅れ時間の補正を行って前記外乱を推定するか、または、他の前記推定部と異なる周期で前記分散共分散行列および前記特異値ベクトルを更新して前記外乱を推定する。
様々な条件を与えて外乱を推定することで、精度よく外乱推定ができる可能性を高めることができる。
様々な条件を与えて外乱を推定することで、精度よく外乱推定ができる可能性を高めることができる。
【0094】
（８）第８の態様に係る外乱推定方法は、制御対象が備えるセンサが計測した計測値を取得するステップと、前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定するステップと、を有する。
【0095】
（９）第９の態様に係るプログラムは、コンピュータに、制御対象が備えるセンサが計測した計測値を取得するステップと、前記計測値を要素とする計測ベクトルの分散共分散行列を算出し、前記分散共分散行列を特異値分解して最大特異値の特異ベクトルを算出し、前記特異ベクトルと前記計測ベクトルとに基づいて、前記制御対象に生じる外乱を推定するステップと、を実行させる。
【符号の説明】
【0096】
１００・・・制御システム、１０～１０Ｅ・・・外乱推定装置、１１・・・取得部、１２～１２Ｅ・・・推定部、１３・・・出力部、２０・・・制御装置、３０・・・制御対象、９００・・・コンピュータ、９０１・・・ＣＰＵ、９０２・・・主記憶装置、９０３・・・補助記憶装置、９０４・・・入出力インタフェース、９０５・・・通信インタフェース
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】