特許 6483571 性能推定装置及び性能推定方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6483571

(24)【登録日】2019年2月22日

(45)【発行日】2019年3月13日

(54)【発明の名称】性能推定装置及び性能推定方法

(51)【国際特許分類】

   G05B 13/04 20060101AFI20190304BHJP

   F02C 7/00 20060101ALI20190304BHJP

   F02C 9/00 20060101ALI20190304BHJP

   F01D 25/00 20060101ALI20190304BHJP

【ＦＩ】

   G05B13/04

   F02C7/00 A

   F02C9/00 A

   F01D25/00 V

【請求項の数】3

【全頁数】15

(21)【出願番号】特願2015-161242(P2015-161242)

(22)【出願日】2015年8月18日

(65)【公開番号】特開2017-16616(P2017-16616A)

(43)【公開日】2017年1月19日

【審査請求日】2018年3月9日

(31)【優先権主張番号】特願2015-133540(P2015-133540)

(32)【優先日】2015年7月2日

(33)【優先権主張国】JP

(73)【特許権者】

【識別番号】000000099

【氏名又は名称】株式会社ＩＨＩ

(73)【特許権者】

【識別番号】899000079

【氏名又は名称】学校法人慶應義塾

(74)【代理人】

【識別番号】100088155

【弁理士】

【氏名又は名称】長谷川 芳樹

(74)【代理人】

【識別番号】100113435

【弁理士】

【氏名又は名称】黒木 義樹

(74)【代理人】

【識別番号】100176245

【弁理士】

【氏名又は名称】安田 亮輔

(74)【代理人】

【識別番号】100171583

【弁理士】

【氏名又は名称】梅景 篤

(72)【発明者】

【氏名】中村 恵子

(72)【発明者】

【氏名】木村 麻衣

(72)【発明者】

【氏名】垣内 大紀

(72)【発明者】

【氏名】若林 優一

(72)【発明者】

【氏名】足立 修一

【審査官】 牧 初

(56)【参考文献】

【文献】 特開２０１４−２１５８２２（ＪＰ，Ａ）

【文献】 特開２０１１−１８００８５（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０５Ｂ １３／００−１３／０４

(57)【特許請求の範囲】

【請求項1】

推定対象の装置である対象装置を制御するための制御入力値及び前記対象装置に設けられたセンサによって計測されたセンサ計測値に基づいて、センサによる計測が不可能な前記対象装置の状態量の推定値である状態推定値を出力する性能推定装置であって、
前記対象装置の非線形シミュレーションモデルを用いて、前記制御入力値に基づいて、前記状態推定値の推定誤差を調整する前の推定値である事前推定値、及び前記センサ計測値の推定値であるセンサ推定値を算出する推定部と、
前記推定部によって算出された前記センサ推定値と前記センサ計測値との誤差を算出する誤差算出部と、
前記誤差算出部によって算出された前記誤差とロバストフィルタゲインとに基づいて、前記推定部に前記事前推定値を調整させるための調整値を生成する調整部と、
を備え、
前記推定部は、前記調整部によって生成された前記調整値に基づいて、前記事前推定値を調整して前記状態推定値を算出し、算出した前記状態推定値を出力する、性能推定装置。

【請求項2】

前記ロバストフィルタゲインは、推定誤差を予め定められた範囲に収めるように算出された固定値である、請求項１に記載の性能推定装置。

【請求項3】

推定対象の装置である対象装置に設けられたセンサによる計測が不可能な状態量の推定値である状態推定値を出力する性能推定方法であって、
前記対象装置を制御するための制御入力値を受信するステップと、
前記対象装置に設けられたセンサによって計測されたセンサ計測値を受信するステップと、
前記対象装置の非線形シミュレーションモデルを用いて、前記制御入力値に基づいて、前記状態推定値の推定誤差を調整する前の推定値である事前推定値、及び前記センサ計測値の推定値であるセンサ推定値を算出するステップと、
前記センサ推定値と前記センサ計測値との誤差を算出するステップと、
前記誤差とロバストフィルタゲインとに基づいて、前記事前推定値を調整するための調整値を生成するステップと、
前記調整値に基づいて、前記事前推定値を調整して前記状態推定値を算出するステップと、
前記状態推定値を出力するステップと、
を含む性能推定方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、性能推定装置及び性能推定方法に関する。

【背景技術】

【0002】

ガスタービンエンジン等では、センサによる計測が不可能な内部状態量及び性能パラメータが存在する。このようなセンサによる計測が不可能な内部状態量及び性能パラメータの値を推定する手法として、カルマンフィルタを用いた手法がある。例えば、特許文献１には、エンジンの物理モデルを含むカルマンフィルタを用いて、センサ計測値から観測不能変数推定値を演算するガスタービンエンジンの性能推定システムが記載されている。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特許第５０４６１０４号公報

【発明の概要】

【発明が解決しようとする課題】

【0004】

特許文献１に記載の性能推定システムでは、カルマンフィルタゲインが一定とされている。このカルマンフィルタゲインは、ガスタービンエンジンのある作動点において求められるので、カルマンフィルタゲインを求めた作動点から外れた作動点では、センサによる計測が不可能な状態量の推定精度が低下し、推定ができない場合もある。

【0005】

本発明は、センサによる計測が不可能な状態量の推定精度を向上可能な性能推定装置及び性能推定方法を提供する。

【課題を解決するための手段】

【0006】

本発明の一側面に係る性能推定装置は、推定対象の装置である対象装置を制御するための制御入力値及び対象装置に設けられたセンサによって計測されたセンサ計測値に基づいて、センサによる計測が不可能な対象装置の状態量の推定値である状態推定値を出力する装置である。この性能推定装置は、対象装置の非線形シミュレーションモデルを用いて、制御入力値に基づいて、状態推定値の推定誤差を調整する前の推定値である事前推定値、及びセンサ計測値の推定値であるセンサ推定値を算出する推定部と、推定部によって算出されたセンサ推定値とセンサ計測値との誤差を算出する誤差算出部と、誤差算出部によって算出された誤差とロバストフィルタゲインとに基づいて、推定部に事前推定値を調整させるための調整値を生成する調整部と、を備える。推定部は、調整部によって生成された調整値に基づいて、事前推定値を調整して状態推定値を算出し、算出した状態推定値を出力する。

【0007】

この性能推定装置では、対象装置の非線形シミュレーションモデルを用いて、事前推定値、及びセンサ推定値が算出され、センサ推定値及びセンサ計測値の誤差とロバストフィルタゲインとに基づいて生成された調整値に基づいて、事前推定値が調整され状態推定値が算出される。このように、対象装置の非線形シミュレーションモデルを用いることにより、線形シミュレーションモデルを用いるよりも推定精度が向上する。また、ロバストフィルタでは、カルマンフィルタと比較して、シミュレーションモデルと推定対象の対象装置との誤差、及び、対象装置の作動点が変化すること等の外乱要素による対象装置の状態の変化によらず、状態推定値を算出することができ、性能推定の安定性を向上することが可能となる。その結果、センサによる計測が不可能な対象装置の性能の推定精度を向上することが可能となる。

【0008】

ロバストフィルタゲインは、推定誤差を予め定められた範囲に収めるように算出された固定値であってもよい。この場合、ロバストフィルタゲインが固定値であるので、ロバストフィルタゲインの更新が不要となる。これにより、演算負荷を低減することが可能となる。

【0009】

本発明の別の側面に係る性能推定方法は、推定対象の装置である対象装置に設けられたセンサによる計測が不可能な状態量の推定値である状態推定値を出力する方法である。この性能推定方法は、対象装置を制御するための制御入力値を受信するステップと、対象装置に設けられたセンサによって計測されたセンサ計測値を受信するステップと、対象装置の非線形シミュレーションモデルを用いて、制御入力値に基づいて、状態推定値の推定誤差を調整する前の推定値である事前推定値、及びセンサ計測値の推定値であるセンサ推定値を算出するステップと、センサ推定値とセンサ計測値との誤差を算出するステップと、誤差とロバストフィルタゲインとに基づいて、事前推定値を調整するための調整値を生成するステップと、調整値に基づいて、事前推定値を調整して状態推定値を算出するステップと、状態推定値を出力するステップと、を含む。

【0010】

この性能推定方法では、対象装置の非線形シミュレーションモデルを用いて、事前推定値、及びセンサ推定値が算出され、センサ推定値及びセンサ計測値の誤差とロバストフィルタゲインとに基づいて生成された調整値に基づいて、事前推定値が調整され状態推定値が算出される。このように、対象装置の非線形シミュレーションモデルを用いることにより、線形シミュレーションモデルを用いるよりも推定精度が向上する。また、ロバストフィルタでは、カルマンフィルタと比較して、シミュレーションモデルと推定対象の対象装置との誤差、及び、対象装置の作動点が変化すること等の外乱要素による対象装置の状態の変化によらず、状態推定値を算出することができ、性能推定の安定性を向上することが可能となる。その結果、センサによる計測が不可能な対象装置の性能の推定精度を向上することが可能となる。

【発明の効果】

【0011】

本発明によれば、センサによる計測が不可能な状態量の推定精度を向上できる。

【図面の簡単な説明】

【0012】

【図1】一実施形態に係る性能推定装置を含む性能推定システムの概略構成を示す図である。

【図2】図１のガスタービンエンジンの各要素を模式的に示す図である。

【図3】図１の性能推定装置のハードウェア構成図である。

【図4】図１の性能推定装置が行う性能推定方法の一連の処理を示すフローチャートである。

【図5】比較例の性能推定装置を含む性能推定システムの概略構成を示す図である。

【発明を実施するための形態】

【0013】

以下、本発明の実施形態について、図面を参照しながら説明する。なお、図面の説明において同一要素には同一符号を付し、重複する説明は省略する。

【0014】

図１は、一実施形態に係る性能推定装置を含む性能推定システムの概略構成を示す図である。図１に示される性能推定システム１は、推定対象の装置である対象装置に設けられる複数のセンサによって計測されたセンサ計測値に基づいて、対象装置の状態を推定するシステムである。以下、対象装置の一例として、ガスタービンエンジンを用いて説明する。なお、時間変化する変数のある時刻ｋにおける値を表す場合に、当該変数を表す符号に（ｋ）を付して説明を行う場合があるが、特定の時刻に限られず、任意の時刻において成り立つ。また、性能推定システム１の単位時間を「１」としており、括弧内の値が大きいほど時間が経過していることを示している。また、図及び数式において、各変数を表す符号にΔが付されている場合があるが、これは後述のロバストフィルタゲインＫを算出した作動点における変数の値からの変化量を意味する。以下において、値及び変化量のいずれでも説明が成り立つ場合には、値を用いて説明を行うこととする。

【0015】

性能推定システム１は、ガスタービンエンジン１０と、性能推定装置２０と、を備えている。ガスタービンエンジン１０は、例えば、車両、船舶、及び航空機等のエンジンとして用いられる。

【0016】

図２は、ガスタービンエンジン１０の各要素を模式的に示す図である。図２に示されるように、ガスタービンエンジン１０は、ファン（ＦＡＮ）１１と、低圧圧縮機（ＬＰＣ）１２と、高圧圧縮機（ＨＰＣ）１３と、燃焼器（ＣＯＭＢ）１４と、高圧タービン（ＨＰＴ）１５と、低圧タービン（ＬＰＴ）１６と、ロータ１７と、ロータ１８と、を備えている。

【0017】

ファン１１は、ガスタービンエンジン１０の外部から空気を吸い込んで、吸い込んだ空気の一部を低圧圧縮機１２に供給する。低圧圧縮機１２は、ファン１１から供給された空気を圧縮して圧縮空気を生成し、圧縮空気を高圧圧縮機１３に供給する。高圧圧縮機１３は、低圧圧縮機１２から供給された圧縮空気をさらに圧縮して圧力をさらに高くして高圧空気を生成し、高圧空気を燃焼器１４に供給する。燃焼器１４は、高圧圧縮機１３から供給された高圧空気に燃料を混ぜ合わせて燃焼する。燃焼器１４は、燃焼して得られた高温高圧の燃焼ガスを高圧タービン１５に供給する。

【0018】

高圧タービン１５は、燃焼器１４から供給された高温高圧の燃焼ガスによって回転し、ロータ１８を回転させる。低圧タービン１６は、高圧タービン１５を通過した燃焼ガスによって回転し、ロータ１７を回転させる。ロータ１７は、ファン１１、低圧圧縮機１２及び低圧タービン１６を一体的に回転可能に接続している。ロータ１８は、高圧圧縮機１３及び高圧タービン１５を一体的に回転可能に接続している。

【0019】

図１に戻って、ガスタービンエンジン１０には、複数のセンサ（不図示）が設けられている。ガスタービンエンジン１０は、エンジン制御器（不図示）によって制御される。ガスタービンエンジン１０は、ある時刻ｋにおいて、エンジン制御器から出力される制御入力値ｕ（ｋ）を入力し、ガスタービンエンジン１０に設けられているセンサによって検知されるセンサ計測値ｙ（ｋ）を出力する。制御入力値ｕ（ｋ）は、ガスタービンエンジン１０を制御するための値であり、エンジン動作条件に基づいて出力される。ガスタービンエンジン１０が航空機に設けられる場合、エンジン動作条件は、高度、機速及び温度日等である。制御入力値ｕ（ｋ）は、燃料流量及び抽気バルブの開度等を含む。センサ計測値ｙ（ｋ）は、ｒ次元であり、ガスタービンエンジン１０のロータ１７，１８の回転数、各部における温度及び圧力を含む。

【0020】

性能推定装置２０は、ガスタービンエンジン１０の非線形モデル及び固定のロバストフィルタゲインＫを用いたロバストフィルタによって、ガスタービンエンジン１０の状態量ｘ（ｋ）を推定する。ロバストフィルタは、例えばＨ∞フィルタである。非線形モデルは、ガスタービンエンジン１０ごとに予め定められた非線形のシミュレーションモデルであり、非線形関数ｆ及び非線形関数ｇを用いて、後述の式（１）〜式（３）で表される。非線形関数ｆは、対象装置（ガスタービンエンジン１０）ごとに予め定められており、内部状態量ｘ_ｅ（ｋ−１）、性能パラメータ値ｑ（ｋ−１）及び制御入力値ｕ（ｋ）の関数である。非線形関数ｇは、対象装置ごとに予め定められており、内部状態量ｘ_ｅ（ｋ）及び性能パラメータ値ｑ（ｋ）の関数である。内部状態量ｘ_ｅは、センサ計測が不可能な状態量を含むガスタービンエンジン１０内部の状態の量である。内部状態量ｘ_ｅは、ｎ_ｅ次元のベクトルであり、各要素の回転数、各要素の内部エネルギー、各要素の空気流量、各要素のエンタルピー等を含む。性能パラメータ値ｑは、センサ計測が不可能なエンジン性能パラメータの値である。性能パラメータ値ｑは、ｐ次元のベクトルであり、ファンの流量、ファンの効率変化、各圧縮機の流量、各圧縮機の効率変化等を含む。

【0021】

具体的には、性能推定装置２０は、制御入力値ｕ（ｋ）及びセンサ計測値ｙ（ｋ）に基づいて、ガスタービンエンジン１０の状態推定値ｘ_ｈ（ｋ）を出力する。状態推定値ｘ_ｈ（ｋ）は、ガスタービンエンジン１０の状態量ｘ（ｋ）の推定値である。状態量ｘ（ｋ）は、ガスタービンエンジン１０の状態量のうち、センサによる計測が不可能な状態量を含む状態量であり、後述するように、内部状態量ｘ_ｅ（ｋ）と性能パラメータ値ｑ（ｋ）とを合わせたｎ（ｎ＝ｎ_ｅ＋ｐ）次元のベクトルである。状態推定値ｘ_ｈ（ｋ）は、内部状態量ｘ_ｅ（ｋ）の推定値と性能パラメータ値ｑ（ｋ）の推定値とを合わせたｎ（ｎ＝ｎ_ｅ＋ｐ）次元のベクトルである。

【0022】

図３は、性能推定装置２０のハードウェア構成図である。図３に示されるように、性能推定装置２０は、物理的には、１又は複数のプロセッサ２０１、主記憶装置であるＲＡＭ（Random Access Memory）及びＲＯＭ（Read Only Memory)等の記憶装置２０２、ハードディスク装置等の補助記憶装置２０３、キーボード等の入力装置２０４、ディスプレイ等の出力装置２０５、並びに、データ送受信デバイスである通信装置２０６等のハードウェアを備えるコンピュータとして構成され得る。性能推定装置２０の図１に示される各機能は、プロセッサ２０１、記憶装置２０２等のハードウェアに１又は複数の所定のコンピュータプログラムを読み込ませることにより、１又は複数のプロセッサ２０１の制御のもとで各ハードウェアを動作させるとともに、記憶装置２０２及び補助記憶装置２０３におけるデータの読み出し及び書き込みを行うことで実現される。

【0023】

図１に戻って、性能推定装置２０の詳細について説明する。性能推定装置２０は、誤差算出部２１と、調整部２２と、推定部３０と、を備えている。

【0024】

誤差算出部２１は、ガスタービンエンジン１０から出力されたセンサ計測値ｙ（ｋ）と推定部３０によって算出されたセンサ推定値ｙ_ｈ（ｋ）との誤差ｅ（ｋ）を算出する。この例では、誤差算出部２１は、センサ計測値ｙ（ｋ）からセンサ推定値ｙ_ｈ（ｋ）を減算し、その減算結果を誤差ｅ（ｋ）として調整部２２に出力する。なお、センサ推定値ｙ_ｈ（ｋ）は、センサ計測値ｙ（ｋ）の推定値である。

【0025】

調整部２２は、誤差算出部２１によって算出された誤差ｅ（ｋ）とロバストフィルタゲインＫとに基づいて、誤差ｅ（ｋ）が小さくなるように、推定部３０に事前推定値ｘ_ｈ⁻（ｋ）を調整させるための調整値を生成する。事前推定値ｘ_ｈ⁻（ｋ）は、状態推定値ｘ_ｈ（ｋ）の推定誤差を調整する前の推定値である。具体的には、調整部２２は、誤差算出部２１から出力された誤差ｅ（ｋ）とロバストフィルタゲインＫとを乗算し、その乗算結果Ｋｅ（ｋ）を調整値として推定部３０に出力する。ロバストフィルタゲインＫは、推定誤差を予め定められた範囲に収めるように算出された固定値のゲインである。ロバストフィルタゲインＫは、オフラインで算出され、調整部２２に予め記憶されている。

【0026】

ここで、ロバストフィルタゲインＫの算出方法について説明する。まず、ガスタービンエンジン１０の非線形モデルが準備される。この非線形モデルは、対象装置ごとに予め定められており、ガスタービンエンジン１０の非線形モデルは、以下の式（１）〜式（３）によって表される。具体的に説明すると、内部状態量ｘ_ｅ（ｋ）は、式（１）に示されるように、非線形関数ｆ及びシステム雑音ｖ_ｘｅ（ｋ）の和として表される。システム雑音ｖ_ｘｅ（ｋ）は、ｎ_ｅ次元のベクトルである。性能パラメータ値ｑ（ｋ）は、式（２）に示されるように、性能パラメータ値ｑ（ｋ−１）及び性能パラメータ変化量ｖ_ｑ（ｋ）の和として表される。性能パラメータ変化量ｖ_ｑ（ｋ）は、ｐ次元のベクトルである。センサ計測値ｙ（ｋ）は、式（３）に示されるように、非線形関数ｇ及び観測雑音ｗ（ｋ）の和として表される。観測雑音ｗは、ｒ次元のベクトルである。

【数1】

【0027】

続いて、ガスタービンエンジン１０の非線形モデルが線形化される。この線形化は、ある作動点を基準として行われる。この作動点としては、例えば、地上静止状態における最大出力時が用いられる。この作動点において線形化された非線形モデルは、作動点における状態量ｘからの変化量Δｘ（ｋ）、作動点におけるセンサ計測値ｙからの変化量Δｙ（ｋ）、作動点における制御入力値ｕからの変化量Δｕ（ｋ）、雑音ｖ（ｋ）、観測雑音ｗ（ｋ）を用いて、例えば式（４）及び式（５）によって表される。行列Ａ，Ｂ，Ｃは、対象装置及び作動点の組み合わせによって、予め定められている。

【数2】

【0028】

なお、式（６）に示されるように、状態量ｘは、内部状態量ｘ_ｅと性能パラメータ値ｑとを合わせたｎ（ｎ＝ｎ_ｅ＋ｐ）次元のベクトルである。雑音ｖは、システム雑音ｖ_ｘｅと性能パラメータ変化量ｖ_ｑとを合わせたｎ（ｎ＝ｎ_ｅ＋ｐ）次元のベクトルである。また、雑音ｖは、式（８）が成立する有限の雑音である。観測雑音ｗは、式（９）が成立する有限の雑音である。

【数3】

【数4】

【0029】

式（６）及び式（７）を用いて式（１）及び式（２）を変形し、式（１）’が得られる。また、式（６）を用いて式（３）を変形し、式（３）’が得られる。

【数5】

【0030】

ロバストフィルタでは、ガスタービンエンジン１０の状態量ｘ（ｋ）と推定部３０において推定された状態推定値ｘ_ｈ（ｋ）との推定誤差は、一定範囲内に収まる必要がある。このため、推定誤差の指標γを用いて、式（１０）が成り立つ必要がある。指標γは、所望の推定精度が得られ、かつ、式（１０）の解が得られる程度の値に設定される。

【数6】

【0031】

なお、ｘ_ｂ（０）は、状態量ｘの初期分布の平均値である。行列Ｑは雑音ｖに対する重み行列であり、式（１１）を満たす。行列Ｒは観測雑音ｗに対する重み行列であり、式（１２）を満たす。

【数7】

【0032】

ここで、式（１３）に示される評価関数Ｊを用いて、式（１０）の解法を検討する。

【数8】

【0033】

行列Ｑ及び行列Ｒは正の値を取るから、Ｑ＞０、Ｒ＞０のもとで、式（１０）を満たす状態推定値の変化量Δｘ_ｈを求める問題は、状態推定値の変化量Δｘ_ｈを固定し、任意の状態量ｘの初期値の変化量Δｘ（０）、雑音ｖ及び観測雑音ｗに対して、評価関数Ｊ＜０が成立することと等価であり、式（１４）として表される。

【数9】

【0034】

さらに、式（１０）は、状態量ｘの初期値の変化量Δｘ（０）、雑音ｖ及び観測雑音ｗを固定し、評価関数Ｊを最小化する状態推定値の変化量Δｘ_ｈを求めることと等価であり、式（１５）として表される。

【数10】

【0035】

式（１５）を満たすために評価関数Ｊが極値をもつためには、ｄＪ＝０となることが必要であることから、推定誤差の共分散Ｐ及び単位行列Ｉを用いて、以下のリカッチ方程式（１６）が導出される。

【数11】

【0036】

式（１６）に基づいて、推定誤差の共分散Ｐを求めるために、具体的には以下の方程式（１７）〜（１９）が用いられる。

【数12】

【0037】

方程式（１７）〜（１９）を解くことによって、式（２０）に示されるように、ロバストフィルタゲインＫが求められる。

【数13】

【0038】

推定部３０は、ガスタービンエンジン１０の非線形モデルを用いて、制御入力値ｕ（ｋ）に基づいて、事前推定値ｘ_ｈ⁻（ｋ）及びセンサ推定値ｙ_ｈ（ｋ）を算出する。推定部３０は、調整部２２によって生成された調整値に基づいて、事前推定値ｘ_ｈ⁻（ｋ）を調整して状態推定値ｘ_ｈ（ｋ）を算出し、算出した状態推定値ｘ_ｈ（ｋ）を出力する。推定部３０は、算出したセンサ推定値ｙ_ｈ（ｋ）を誤差算出部２１に出力し、算出した状態推定値ｘ_ｈ（ｋ）を外部に出力する。推定部３０は、時間更新部３１と、推定値更新部３２と、推定値更新部３３と、遅延部３４と、を備えている。

【0039】

時間更新部３１は、制御入力値ｕ（ｋ）と、遅延部３４から出力された状態推定値ｘ_ｈ（ｋ−１）と、に基づいて、事前推定値ｘ_ｈ⁻（ｋ）を算出する。具体的には、時間更新部３１は、式（１）で用いられた非線形関数ｆを用い、式（２１）によって事前推定値ｘ_ｈ⁻（ｋ）を算出する。時間更新部３１は、算出した事前推定値ｘ_ｈ⁻（ｋ）を推定値更新部３２及び推定値更新部３３に出力する。

【数14】

【0040】

推定値更新部３２は、時間更新部３１から出力された事前推定値ｘ_ｈ⁻（ｋ）に基づいて、センサ推定値ｙ_ｈ（ｋ）を算出する。具体的には、推定値更新部３２は、式（３）で用いられた非線形関数ｇを用い、式（２２）によってセンサ推定値ｙ_ｈ（ｋ）を算出する。推定値更新部３２は、算出したセンサ推定値ｙ_ｈ（ｋ）を誤差算出部２１に出力する。

【数15】

【0041】

推定値更新部３３は、時間更新部３１から出力された事前推定値ｘ_ｈ⁻（ｋ）と、調整部２２から出力された調整値とに基づいて、状態推定値ｘ_ｈ（ｋ）を算出する。具体的には、推定値更新部３３は加算器であり、式（２３）に示されるように、事前推定値ｘ_ｈ⁻（ｋ）と調整値とを加算し、その加算結果を状態推定値ｘ_ｈ（ｋ）として算出する。推定値更新部３３は、算出した状態推定値ｘ_ｈ（ｋ）を遅延部３４及び外部に出力する。

【数16】

【0042】

遅延部３４は、推定値更新部３３から出力された状態推定値ｘ_ｈ（ｋ）を単位時間だけ遅延させ、次の時刻における推定に対し、状態推定値ｘ_ｈ（ｋ−１）として時間更新部３１に出力する。

【0043】

次に、図４を参照して、性能推定装置２０が行う性能推定方法の一連の処理を説明する。図４は、性能推定装置２０が行う性能推定方法の一連の処理を示すフローチャートである。図４に示される処理は、単位時間ごとに繰り返し行われる。ここでは、ある時刻ｋにおける処理を説明する。

【0044】

まず、エンジン制御器が、エンジン動作条件に基づいて、制御入力値ｕ（ｋ）をガスタービンエンジン１０及び性能推定装置２０に出力する。そして、性能推定装置２０の時間更新部３１は、制御入力値ｕ（ｋ）を受信する（ステップＳ０１）。また、ガスタービンエンジン１０は、エンジン制御器から出力された制御入力値ｕ（ｋ）を入力し、ガスタービンエンジン１０に設けられているセンサによって計測されたセンサ計測値ｙ（ｋ）を出力する。そして、性能推定装置２０の誤差算出部２１は、センサ計測値ｙ（ｋ）を受信する（ステップＳ０２）。

【0045】

続いて、時間更新部３１は、ガスタービンエンジン１０の非線形モデルを用いて、制御入力値ｕ（ｋ）に基づいて、事前推定値ｘ_ｈ⁻（ｋ）を算出する（ステップＳ０３）。具体的には、時間更新部３１は、制御入力値ｕ（ｋ）と、遅延部３４から出力された状態推定値ｘ_ｈ（ｋ−１）と、を用いて、式（２１）によって事前推定値ｘ_ｈ⁻（ｋ）を算出する。そして、時間更新部３１は、算出した事前推定値ｘ_ｈ⁻（ｋ）を推定値更新部３２及び推定値更新部３３に出力する。

【0046】

続いて、推定値更新部３２は、ガスタービンエンジン１０の非線形モデルを用いて、センサ推定値を算出する（ステップＳ０４）。具体的には、推定値更新部３２は、時間更新部３１から出力された事前推定値ｘ_ｈ⁻（ｋ）を用いて、式（２２）によってセンサ推定値ｙ_ｈ（ｋ）を算出する。そして、推定値更新部３２は、算出したセンサ推定値ｙ_ｈ（ｋ）を誤差算出部２１に出力する。

【0047】

続いて、誤差算出部２１は、センサ計測値ｙ（ｋ）とセンサ推定値ｙ_ｈ（ｋ）との誤差ｅ（ｋ）を算出する（ステップＳ０５）。具体的には、誤差算出部２１は、センサ計測値ｙ（ｋ）からセンサ推定値ｙ_ｈ（ｋ）を減算し、その減算結果を誤差ｅ（ｋ）として調整部２２に出力する。

【0048】

続いて、調整部２２は、誤差算出部２１によって算出された誤差ｅ（ｋ）とロバストフィルタゲインＫとに基づいて、事前推定値ｘ_ｈ⁻（ｋ）を調整するための調整値を生成する（ステップＳ０６）。具体的には、調整部２２は、誤差ｅ（ｋ）とロバストフィルタゲインＫとを乗算し、その乗算結果Ｋｅ（ｋ）を調整値として推定部３０に出力する。

【0049】

続いて、推定値更新部３３は、調整値に基づいて、事前推定値ｘ_ｈ⁻（ｋ）を調整して状態推定値ｘ_ｈ（ｋ）を算出する（ステップＳ０７）。具体的には、推定値更新部３３は、式（２３）に示されるように、事前推定値ｘ_ｈ⁻（ｋ）と調整値とを加算し、その加算結果を状態推定値ｘ_ｈ（ｋ）として算出する。そして、推定値更新部３３は、算出した状態推定値ｘ_ｈ（ｋ）を外部及び遅延部３４に出力する（ステップＳ０８）。以上のようにして、時刻ｋにおける性能推定方法の一連の処理が終了する。

【0050】

このように、性能推定装置２０は、ガスタービンエンジン１０の非線形モデル及び固定のロバストフィルタゲインＫを用いたロバストフィルタによって、実機（ガスタービンエンジン１０）の状態量ｘ（ｋ）（内部状態量ｘ_ｅ（ｋ）及び性能パラメータ値ｑ（ｋ））を推定する。そして、性能推定装置２０は、センサ計測値ｙ（ｋ）と、推定部３０から出力されるセンサ推定値ｙ_ｈ（ｋ）との誤差ｅ（ｋ）が小さくなるように、誤差ｅ（ｋ）にロバストフィルタゲインＫを乗算し、事前推定値ｘ_ｈ⁻（ｋ）の調整を繰り返し行うことで、ガスタービンエンジン１０の状態量ｘ（ｋ）を推定、つまり、状態推定値ｘ_ｈ（ｋ）を算出する。

【0051】

次に、比較例の性能推定装置１２０を含む性能推定システム１００と比較しながら性能推定装置２０の作用効果を説明する。図５は、比較例の性能推定装置１２０を含む性能推定システム１００の概略構成を示す図である。図５に示されるように、性能推定システム１００は、性能推定システム１と比較して、性能推定装置２０に代えて、性能推定装置１２０を備える点で相違している。

【0052】

性能推定装置１２０は、性能推定装置２０と比較して、ロバストフィルタゲインＫ（ｋ）が固定でない点、ガスタービンエンジン１０の非線形モデルではなく線形モデルを用いる点で主に相違している。つまり、性能推定装置１２０は、ガスタービンエンジン１０の線形モデル及びロバストフィルタゲインＫ（ｋ）を用いて、ガスタービンエンジン１０の状態推定値ｘ_ｈ（ｋ）を推定する装置であり、誤差算出部１２１と、調整部１２２と、推定部１３０と、を備えている。

【0053】

線形モデルは、行列Ａ，Ｂ，Ｃを用いて、後述の式（２４）及び式（２５）で表される。行列Ａ，Ｂ，Ｃは、対象装置（ガスタービンエンジン１０）ごとに、各作動点に対して予め定められている。この行列Ａ，Ｂ，Ｃは、式（４）及び式（５）の行列Ａ，Ｂ，Ｃと同じである。

【0054】

性能推定装置１２０が行う性能推定方法の一連の処理を説明する。まず、エンジン制御器が、エンジン動作条件に基づいて、制御入力値ｕ（ｋ）をガスタービンエンジン１０及び性能推定装置１２０に出力する。そして、性能推定装置１２０の時間更新部１３１は、制御入力値ｕ（ｋ）を受信する。また、ガスタービンエンジン１０は、エンジン制御器から出力された制御入力値ｕ（ｋ）を入力し、ガスタービンエンジン１０に設けられているセンサによって計測されたセンサ計測値ｙ（ｋ）を出力する。そして、性能推定装置１２０の誤差算出部１２１は、センサ計測値ｙ（ｋ）を受信する。

【0055】

続いて、時間更新部１３１は、ガスタービンエンジン１０の線形モデルを用いて、制御入力値ｕ（ｋ）に基づいて、事前推定値ｘ_ｈ⁻（ｋ）を算出する。具体的には、時間更新部１３１は、演算部１３１ａと、演算部１３１ｂと、加算器１３１ｃと、を備えている。まず、演算部１３１ａは、行列Ａと遅延部１３４から出力された状態推定値ｘ_ｈ（ｋ−１）との乗算を行い、乗算結果を加算器１３１ｃに出力する。また、演算部１３１ｂは、行列Ｂと制御入力値ｕ（ｋ）との乗算を行い、乗算結果を加算器１３１ｃに出力する。そして、加算器１３１ｃは、演算部１３１ａの乗算結果と演算部１３１ｂの乗算結果とを加算し、加算結果を事前推定値ｘ_ｈ⁻（ｋ）として算出する。つまり、時間更新部１３１は、式（２４）に示される演算を行って、事前推定値ｘ_ｈ⁻（ｋ）を算出している。そして、時間更新部１３１は、算出した事前推定値ｘ_ｈ⁻（ｋ）を推定値更新部１３２及び推定値更新部１３３に出力する。

【数17】

【0056】

続いて、推定値更新部１３２は、ガスタービンエンジン１０の線形モデルを用いて、センサ推定値を算出する。具体的には、推定値更新部１３２は、式（２５）に示されるように、行列Ｃと時間更新部１３１から出力された事前推定値ｘ_ｈ⁻（ｋ）との乗算を行い、乗算結果をセンサ推定値ｙ_ｈ（ｋ）として算出する。そして、推定値更新部１３２は、算出したセンサ推定値ｙ_ｈ（ｋ）を誤差算出部１２１に出力する。

【数18】

【0057】

続いて、誤差算出部１２１は、誤差算出部２１と同様に、センサ計測値ｙ（ｋ）とセンサ推定値ｙ_ｈ（ｋ）との誤差ｅ（ｋ）を算出する。具体的には、誤差算出部１２１は、センサ計測値ｙ（ｋ）からセンサ推定値ｙ_ｈ（ｋ）を減算し、その減算結果を誤差ｅ（ｋ）として調整部１２２に出力する。

【0058】

続いて、調整部１２２は、誤差算出部１２１によって算出された誤差ｅ（ｋ）とロバストフィルタゲインＫ（ｋ）とに基づいて、事前推定値ｘ_ｈ⁻（ｋ）を調整するための調整値を生成する。具体的には、調整部１２２は、まず、式（２６）に示されるように、単位行列Ｉ、行列Ｃ、観測雑音ｗに対する重み行列Ｒ、推定誤差の指標γ、及び推定誤差の共分散Ｐ（ｋ−１）を用いて行列Ψ（ｋ）を算出する。そして、調整部１２２は、式（２７）に示されるように、行列Ａ、推定誤差の共分散Ｐ（ｋ−１）、算出した行列Ψ（ｋ）、及び雑音ｖに対する重み行列Ｑを用いて推定誤差の共分散Ｐ（ｋ）を算出する。さらに、調整部１２２は、式（２８）に示されるように、算出した推定誤差の共分散Ｐ（ｋ）、行列Ｃ、及び重み行列Ｒを用いてロバストフィルタゲインＫ（ｋ）を算出する。そして、調整部１２２は、誤差算出部１２１から出力された誤差ｅ（ｋ）とロバストフィルタゲインＫ（ｋ）とを乗算し、その乗算結果Ｋ（ｋ）ｅ（ｋ）を調整値として推定部１３０に出力する。

【数19】

【数20】

【数21】

【0059】

続いて、推定値更新部１３３は、推定値更新部３３と同様に、時間更新部１３１から出力された事前推定値ｘ_ｈ⁻（ｋ）と、調整部１２２から出力された調整値とに基づいて、状態推定値ｘ_ｈ（ｋ）を算出する。具体的には、推定値更新部１３３は加算器であり、式（２９）に示されるように、事前推定値ｘ_ｈ⁻（ｋ）と調整値とを加算し、その加算結果を状態推定値ｘ_ｈ（ｋ）として算出する。そして、推定値更新部１３３は、算出した状態推定値ｘ_ｈ（ｋ）を遅延部１３４及び外部に出力する。以上のようにして、性能推定装置１２０が行う時刻ｋにおける性能推定方法の一連の処理が終了する。

【数22】

【0060】

このように、性能推定装置１２０は、ガスタービンエンジン１０の線形モデルを用いているので、非線形なガスタービンエンジン１０では推定精度が低下する。また、性能推定装置１２０では、作動点に応じた線形モデルを準備する必要があり、さらにロバストフィルタゲインＫ（ｋ）を単位時間ごとに毎回更新する必要があるので、演算負荷が高い。

【0061】

一方、性能推定装置２０では、ガスタービンエンジン１０の非線形モデルを用いて、事前推定値ｘ_ｈ⁻（ｋ）、及びセンサ推定値ｙ_ｈ（ｋ）が算出され、センサ推定値ｙ_ｈ（ｋ）及びセンサ計測値ｙ（ｋ）の誤差ｅ（ｋ）とロバストフィルタゲインＫとに基づいて生成された調整値に基づいて、事前推定値ｘ_ｈ⁻（ｋ）が調整され状態推定値ｘ_ｈ（ｋ）が算出される。このように、性能推定装置２０では、ガスタービンエンジン１０の非線形モデルが用いられることにより、線形モデルが用いられるよりも推定精度が向上する。また、ロバストフィルタゲインＫは、推定誤差を予め定められた範囲に収めるように算出された固定値である。このため、ロバストフィルタゲインＫの更新が不要となる。これにより、演算負荷を低減することが可能となる。

【0062】

また、ロバストフィルタでは、カルマンフィルタと比較して、シミュレーションモデルとガスタービンエンジン１０との誤差、及び、ガスタービンエンジン１０の作動点が変化すること等の外乱要素によるガスタービンエンジン１０の状態の変化によらず、状態推定値ｘ_ｈ（ｋ）を算出することができ、性能推定の安定性を向上することが可能となる。つまり、フィルタ設計点以外の作動点での性能推定の安定性が向上する。

【0063】

以上、本発明の実施形態について説明したが、本発明は上記実施形態に限定されない。例えば、上記実施形態では、対象装置の一例としてガスタービンエンジンを用いているがこれに限られない。対象装置は、制御入力値ｕ（ｋ）に対して非線形な動作を行う装置であって、制御入力値ｕ（ｋ）を入力とし、センサ計測値ｙ（ｋ）及び状態量ｘ（ｋ）を出力とする非線形シミュレーションモデルで表現し得る装置であればよい。

【0064】

また、上記実施形態では、調整部２２は、固定値のロバストフィルタゲインＫを用いているが、調整部２２は、単位時間ごとにロバストフィルタゲインＫ（ｋ）を算出してもよい。この場合、調整部２２は、各時刻において、式（４）及び式（５）に示されるように、その作動点で非線形モデルを線形化し、式（１６）に基づいて（具体的には式（１７）〜式（１９）を用いて）、推定誤差の共分散Ｐを算出する。そして、調整部２２は、式（２０）を用いてロバストフィルタゲインＫ（ｋ）を算出する。

【符号の説明】

【0065】

１０ ガスタービンエンジン
２０ 性能推定装置
２１ 誤差算出部
２２ 調整部
３０ 推定部
ｅ 誤差
Ｋ ロバストフィルタゲイン
ｕ 制御入力値
ｘ 状態量
ｘ_ｈ 状態推定値
ｘ_ｈ⁻ 事前推定値
ｙ センサ計測値
ｙ_ｈ センサ推定値

【図1】

【図2】

【図3】

【図4】

【図5】