特許 6260437 動的システムの推定装置及び方法

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

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

(11)【特許番号】6260437

(24)【登録日】2017年12月22日

(45)【発行日】2018年1月17日

(54)【発明の名称】動的システムの推定装置及び方法

(51)【国際特許分類】

   F02C 7/00 20060101AFI20180104BHJP

   F01D 25/00 20060101ALI20180104BHJP

   G05B 13/04 20060101ALI20180104BHJP

   G05B 23/02 20060101ALI20180104BHJP

【ＦＩ】

   F02C7/00 A

   F01D25/00 W

   F01D25/00 X

   G05B13/04

   G05B23/02 R

【請求項の数】4

【全頁数】19

(21)【出願番号】特願2014-92861(P2014-92861)

(22)【出願日】2014年4月28日

(65)【公開番号】特開2015-209828(P2015-209828A)

(43)【公開日】2015年11月24日

【審査請求日】2017年2月23日

(73)【特許権者】

【識別番号】000000099

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

(74)【代理人】

【識別番号】100083806

【弁理士】

【氏名又は名称】三好 秀和

(74)【代理人】

【識別番号】100100712

【弁理士】

【氏名又は名称】岩▲崎▼ 幸邦

(74)【代理人】

【識別番号】100101247

【弁理士】

【氏名又は名称】高橋 俊一

(74)【代理人】

【識別番号】100095500

【弁理士】

【氏名又は名称】伊藤 正和

(74)【代理人】

【識別番号】100098327

【弁理士】

【氏名又は名称】高松 俊雄

(72)【発明者】

【氏名】古川 洋之

(72)【発明者】

【氏名】垣内 大紀

(72)【発明者】

【氏名】木村 麻衣

【審査官】 瀬戸 康平

(56)【参考文献】

【文献】 特開２００９−１６２２３０（ＪＰ，Ａ）

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

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

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

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

【文献】 米国特許出願公開第２００９／０３２６８９０（ＵＳ，Ａ１）

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

Ｆ０１Ｄ ２５／００

Ｆ０２Ｃ ７／００

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

１７／００−１７／０２，２１／００−２１／０２，２３／０２

ＤＷＰＩ（Ｄｅｒｗｅｎｔ Ｉｎｎｏｖａｔｉｏｎ）

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

【請求項1】

動的システムの状態及び性能パラメータを推定する動的システムの推定装置であって、
動的システムから所定の観測値を検出する所定の数のセンサと、
前記センサにて検出した観測値に基づいて、前記動的システムの動的モデルを用いて前記動的システムの状態を推定するカルマンフィルタとを含み、
前記動的システムは、性能パラメータにより経時的な性能の変化が記述される要素を含み、前記カルマンフィルタは、前記要素に替えて一次遅れ要素を導入することにより、前記所定の数より多い性能パラメータの推定も可能とすることを特徴とする推定装置。

【請求項2】

前記性能パラメータは、前記一次遅れ要素の出力値の入力値に対する比として与えられることを特徴とする請求項１に記載の推定装置。

【請求項3】

前記動的システムはガスタービンエンジンであって、前記センサが検出する観測値は、ガスタービンエンジンにおける回転数、温度及び圧力を含むことを特徴とする請求項１又は２に記載の推定装置。

【請求項4】

動的システムの状態及び性能パラメータを推定する動的システムの推定方法であって、
所定の数のセンサにて動的システムから所定の観測値を検出するステップと、
前記センサにて検出した観測値に基づいて、前記動的システムの動的モデルを用いたカルマンフィルタにて前記動的システムの状態を推定するステップとを含み、
前記動的システムは、性能パラメータにより経時的な性能の変化が記述される要素を含み、前記カルマンフィルタは、前記要素に替えて一次遅れ要素を導入することにより、前記所定の数より多い性能パラメータの推定も可能とすることを特徴とする動的システムの推定方法。

【発明の詳細な説明】

【技術分野】

【0001】

この発明は、動的システムの状態及び性能パラメータを推定する動的システムの推定装置及び方法に関する。

【背景技術】

【0002】

従来、複数の要素から構成されるガスタービンエンジンの如き動的システムにおいて、この動的システムの状態を推定する推定装置が提供されている。この推定装置は、エンジンの回転数、温度、圧力などの観測値をセンサにて検出し、動的モデルを備えるカルマンフィルタを用いて動的システムの各種状態を推定することができる。このような推定装置としては、例えば引用文献１が知られている。

【0003】

動的システムを構成する要素においては、性能パラメータによってその性能が規定されるものがある。このような要素においては、経時的な変化を考慮しない限り性能パラメータは時間によらず一定である。カルマンフィルタによると、経時的な変化を考慮する場合には性能パラメータを推定することができるが、推定できる性能パラメータの数はセンサの数以下に限られていた。

【先行技術文献】

【特許文献】

【0004】

【特許文献1】特開平８−１２３５０６号公報

【発明の概要】

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

【0005】

しかしながら、例えばガスタービンエンジンの如き動的システムにおいては、エンジンの回転数、温度、圧力などセンサにより検出できる観測値が限定的であり、推定できる性能パラメータの数も限られていた。したがって、外部から観察できないエンジン内部にある各要素の経年劣化や異常を十分に特定することができなかった。

【0006】

このため、動的システムにおいては、センサの数より多い性能パラメータについても経時的な変化を考慮して推定できるようにすることが求められている。

【0007】

この出願に係る発明は、上述の実情に鑑みて提案されるものであって、センサの数より多い性能パラメータについても経時的な変化を考慮して推定することができるような動的システムの推定装置及び方法を提供することを目的とする。

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

【0008】

この発明に係る動的システムの推定装置は、動的システムの状態及び性能パラメータを推定する動的システムの推定装置である。動的システムから所定の観測値を検出する所定の数のセンサと、前記センサにて検出した観測値に基づいて、前記動的システムの動的モデルを用いて前記動的システムの状態を推定するカルマンフィルタとを含んでいる。前記動的システムは、性能パラメータにより経時的な性能の変化が記述される要素を含み、前記カルマンフィルタは、前記要素に替えて一次遅れ要素を導入することにより、前記所定の数より多い性能パラメータの推定も可能とする。

【0009】

前記性能パラメータは、前記一次遅れ要素の出力値の入力値に対する比として与えられてもよい。前記動的システムはガスタービンエンジンであって、前記観測値は、ガスタービンエンジンの回転数、温度及び圧力を含んでもよい。

【0010】

この発明に係る動的システムの推定方法は、動的システムの状態及び性能パラメータを推定する動的システムの推定方法である。所定の数のセンサにて動的システムから所定の観測値を検出するステップと、前記センサにて検出した観測値に基づいて、前記動的システムの動的モデルを用いたカルマンフィルタにて前記動的システムの状態を推定するステップとを含んでいる。前記動的システムは、性能パラメータにより経時的な性能の変化が記述される要素を含み、前記カルマンフィルタは、前記要素に替えて一次遅れ要素を導入することにより、前記所定の数より多い性能パラメータの推定も可能とする。

【発明の効果】

【0011】

この発明によると、センサの数より多い性能パラメータについても経時的な変化を考慮して推定することができるようになる。したがって、動的システムにおいてセンサの数が限定されている場合にも、センサの数より多い要素の経時的な性能の変化を知ることができる。

【図面の簡単な説明】

【0012】

【図1】ガスタービンエンジン推定装置の概略的な構成を示す図である。

【図2】ガスタービンエンジンの実機モデルを示す図である。

【図3】積分要素の構成を示す図である。

【図4】ファン・圧縮機要素を説明する図である。

【図5】ファン・圧縮機要素と要素間容積の関係を説明する図である。

【図6】カルマンフィルタの構成を示す図である。

【図7】一次遅れ要素の構成を示す図である。

【図8】一次遅れを説明する図である。

【図9】図５において乗算器を一次遅れに置き換えた図である。

【図10】実施例１の結果を示すグラフである。

【図11】実施例２の結果を示すグラフである。

【図12】実施例３の結果を示すグラフである。

【図13】実施例４の結果を示すグラフである。

【図14】従来のカルマンフィルタを含むガスタービンエンジン推定装置の構成を示す図である。

【図15】カルマンフィルタによる補正を説明する図である。

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

【0013】

以下、動的システムの推定装置及び推定方法の実施の形態として、ガスタービンエンジン推定装置について図面を参照して詳細に説明する。

【0014】

図１は、ガスタービンエンジン推定装置の概略的な構成を示す図である。このガスタービンエンジン推定装置において、図中の上側にはガスタービンエンジンの実機モデルが配置されている。この実機モデルは、後述するように、ガスタービンエンジンの構成を簡易化して動作を模擬するように、ガスタービンエンジンから高圧圧縮機、燃焼器を含む構成を取り出してモデル化したものである。また、燃焼器の後段には高圧タービンの一部の構成も含まれている。図中では、高圧圧縮機を単に圧縮機と記載している。

【0015】

ガスタービンエンジン推定装置において、図中の下側にはカルマンフィルタが配置されている。このカルマンフィルタは、ガスタービンエンジンの実機モデルのダイナミクスを再現する動的モデルとして実機モデルと同様の構成を含んでいる。なお、図中において、カルマンフィルタによる推定値にはハットをつけて実機モデルの変数と区別することにする。

【0016】

このカルマンフィルタの圧縮機においては、実機モデルに性能パラメータ（流量変化特性係数）ｑｗｃｉを導入する乗算器に対応する部分には、この性能パラメータｑｗｃｉを推定するための一次遅れ要素が含まれている。本実施の形態では、この一次遅れ要素を導入することによって性能パラメータを状態変数に変換し、センサの数より多い性能パラメータについても時間変化を推定することができるようにしている。このような一次遅れ要素についてはさらに後述する。

【0017】

図２は、図１に示したガスタービンエンジン推定装置における実機モデルの構成を示す図である。この実機モデルは、実際のガスタービンエンジンの動作を模擬するように簡略化して構成したものである。図中の圧縮機と燃焼器のブロックは、それぞれガスタービンエンジンの高圧圧縮機及び燃焼器に対応している。また、燃焼器のブロックの右側の部分は高圧タービンの一部に相当している。

【0018】

図３は、積分要素の構成を示す図である。本実施の形態では、離散時間でシミュレーションを行うため、離散時間の増分ｄｔを掛けてから積分を行っている。ペナント（三角旗）型のブロックは、入力値を定数倍するゲイン要素を意味している。なお、カルマンフィルタにおけるゲイン要素については、対応する符号に“´”を服することで実機モデルのゲイン要素との対応関係を示すことにする。

【0019】

図４は、ガスタービンエンジンのファン・圧縮機要素のモデル化について説明する図である。このファン・圧縮機要素は、図２に示した実機モデルにおける圧縮器に相当している。また、燃焼室の後段の高圧タービンの一部に対応する部分にも相当している。

【0020】

図４（ａ）に示すように、ファン・圧縮機要素は、断熱効率η、ロータ機械の回転数Ｎであり、要素入口空気圧力Ｐｃｉ、要素入口空気流量Ｗｃｉ、要素出口空気圧力Ｐｏ、要素出口空気流量Ｗｏ、抽気流量Ｗｂである。

【0021】

図４（ｂ）は、これらの変数について、ファン・圧縮機要素の要素性能と要素間容積における関係を示している。要素間容積は、ファン・圧縮機要素の出口容積に相当している。

【0022】

要素性能は、要素出口空気圧力Ｐｏの要素入口空気圧力Ｐｃｉに対する比Ｐｏ／Ｐｃｉと回転数Ｎに応じて要素入口空気流量Ｗｃｉと断熱効率ηが決まる静特性となっている。要素間容積は、容積出口空気流量Ｗｏと容積入口空気流量Ｗｉの差Ｗｏ−Ｗｉに応じて容積の圧力Ｐｏと質量Ｍが決まる動特性となっている。ここで、要素入口空気流量Ｗｉは要素入口空気流量Ｗｃｉから抽気流量Ｗｂを差し引いたＷｃｉ−Ｗｂである。

【0023】

図５は、図４（ｂ）で示したファン・圧縮機要素と要素間容積における関係を説明する図である。ファン・圧縮機要素には、要素入口圧力Ｐｃｉ、回転数Ｎ、要素入口比エンタルピｈｃｉが入力されている。また、後段の要素間要素から容積出口圧力Ｐ０が入力されている。

【0024】

ファン・圧縮機要素は、空気流量Ｗｃｉｚと容積出口圧力Ｐ０の容積入口圧力Ｐｃｉに対する圧力比Π＝Ｐ０／Ｐｃｉとの関係を示す特性曲線のマップに従い、当該圧力比Πに対応する空気流量Ｗｃｉ ｍａｐを決める。この値に性能パラメータ（流量特性変化係数）ｑｗｃｉを乗じた値を要素入口空気流量Ｗｃｉとして出力する。

【0025】

また、ファン・圧縮機要素は、空気流量Ｗｃｉと断熱効率ηの関係を示す特性曲線のマップに従い、要素入口空気流Ｗｃｉに対応する断熱効率ηを決める。そして、要素入口比エンタルピｈｃｉ、圧力比Π、断熱効率ηに基づいて計算した容積入口比エンタルピｈｉを出力する。

【0026】

前記の空気流量Ｗｃｉと前記圧力比Πの関係を示すマップと空気流量Ｗｃｉと断熱効率ηの関係を示すマップにおいては、複数の特性曲線が描かれている。これらの特性曲線は、例えばファン・圧縮機要素の効率によって選択される。図中では９０％の効率に相当する特性曲線が選択されている。

【0027】

要素間容積には、ファン・圧縮機要素からの要素入口空気流量Ｗｃｉから抽気流量Ｗｂを差し引いた要素入口空気流量Ｗｉが入力され、ファン・圧縮機要素から容積入口比エンタルピｈｉが入力されている。また、後段の要素から容積出口空気流量Ｗｏが入力されている。

【0028】

要素間容積においては、容積入口空気流量Ｗｉから容積出口容器流量Ｗｏを差引いた空気流量が質量の時間微分となっている。要素間容積は、この質量の時間微分を積分した質量Ｍを出力している。また、容積入口空気流量Ｗｉ、容積出口空気流量Ｗｏ、質量Ｍ等に基づいて容積出口エンタルピｈｏ、容積出口圧力ｐｏを計算して出力している。

【0029】

図２に示した実機モデルは、このようなファン・圧縮器要素と要素間容積の関係に基づいて次のように構成されている。

【0030】

圧縮機における空気流量Ｗ２５ｍａｐは、高圧圧縮機の特性曲線のマップにより定められる空気流量を模擬している。このため、Ｗ２５ｍａｐは、第１ゲイン要素１１によりロータ機械の回転数ＮＧに比例する値と、第２ゲイン要素１２により圧縮機の出口圧力Ｐ３に比例する値の和とされている。

【0031】

圧縮器では、図５のファン・圧縮機要素に示したように、空気流量Ｗ２５ｍａｐに性能パラメータ（流量特性変化係数）ｑｗｃｉを乗算して要素入口空気流量Ｗ２５としている。要素間容積においては、流入空気流量となるＷ２５から流出空気流量となるＷ３を差し引いたものが質量の時間微分となることから、Ｗ２５−Ｗ３を積分することによりＭ３を得ている。また、要素間容積における圧力Ｐ３は、第３ゲイン要素１３により質量Ｍ３に比例するものとしている。

【0032】

要素間容積の出口空気流量Ｗ３は、このＷ３の時間微分が入口圧より出口圧を差し引いた値に比例することから、入口圧となるＰ３から出口圧となるＰ４を差し引いたＰ３−Ｐ４を積分した値としている。

【0033】

燃焼器では、要素間容積において、質量Ｍ４の時間微分は、高圧圧縮機の出口空気流量Ｗ３に燃料流量ＷＦを加えた流入ガス流量から流出ガス流量となるＷ４を差し引いた値となっている。したがって、質量Ｍ４はこれらの値を積分することにより求められる。また、燃焼器のガス圧力Ｐ４は、第６ゲイン要素１６により質量Ｍ４に比例するものとされる。

【0034】

燃焼器の後段の高圧タービンに相当する部分においては、高圧タービンにおける特性曲線のマップにより定められるガス流量を模擬している。このため、高圧タービンの出口ガス流量Ｗ４は、第９ゲイン要素２１によりロータ機械の回転数ＮＧに比例するようにされた値と、第８ゲイン要素２０により燃焼器のガス圧力Ｐ４に比例するようにされた値の和とれている。

【0035】

ロータ機械の回転数ＮＧの時間微分は、高圧タービンのトルクから高圧圧縮機のトルクを差し引いた値に比例する。したがって、ロータ機械の回転数ＮＧはこの値を積分することにより求められる。高圧タービンのトルクと高圧圧縮機のトルクは、第５ゲイン要素１５により燃焼器のガス圧力Ｐ４に比例するようにされた値と、第７ゲイン要素１９により高圧圧縮機のガス圧力Ｐ３に比例するようにされた値とされている。

【0036】

図６は、図１に示したガスタービンエンジン推定装置におけるカルマンフィルタの構成を示す図である。本実施の形態のカルマンフィルタは、図２で示した実機モデルと対比すると、実機モデルにおいて性能パラメータ（流量特性変化係数）ｑｗｃｉを乗算した乗算器が一次遅れ要素に置き換わっている点において相違している。

【0037】

図７は、一次遅れ要素の構成を示す図である。一次遅れ要素は、積分要素の構成を含み、所定の時定数τを有している。本実施の形態では、離散時間でシミュレーションするために離散時間の増分ｄｔを含んでいる。

【0038】

カルマンフィルタは、状態ｘ、性能パラメータｑ、センサパラメータｙ、入力ｕについて、次のような式によって表すことができる。ｔは離散時間であり、離散時間の増分はΔｔである。ｘはｎ行、ｑはｐ行、ｙはｍ行、ｕはｌ行の列ベクトルである。

【数1】

【0039】

一般に、状態ｘはロータ機械の回転数、各要素の内部エネルギー、エンタルピなど動的システムにおける状態を記述する。性能パラメータは、流量・効率等の動的システムを構成する要素の性能を規定する特性変化係数である。センサパラメータは、ロータ機械の回転数、温度、圧力等の動的システムの外部からセンサにて検出された観測値である。

【0040】

具体的に、図６に示したカルマンフィルタは、状態ｘとして次のようなパラメータを含んでいる。

【表1】

【0041】

性能パラメータｑとして次のようなパラメータを含んでいる。

【表2】

【0042】

センサｙとして次のようなパラメータを含んでいる。

【表3】

【0043】

入力ｕとして次のようなパラメータを含んでいる。

【表4】

【0044】

式（１）は、次のように線形化することができる。

【数2】

【0045】

式（２）を行列表示すると、次のようになる。この式（３）に基づいて定常カルマンフィルタを設計する。

【数3】

【0046】

ただし、次の（ｎ＋ｐ）行の列ベクトルを規定した。

【数4】

【0047】

さらに、離散形の式（１）を連続形で記述すると、式（５）のようになる。ここで、１行目と２行目の式を式（４）で規定されるベクトルを用いて一つの式とすることができる。

【数5】

【0048】

カルマンフィルタは次のようになる。

【数6】

【0049】

線形化すると次のようになる。

【数7】

【0050】

ここで、性能パラメータを含む状態ｘ_ｑを推定するためには（Ａ、Ｃ）が可観測である必要があることから、次のような制約がある。すなわち、下記の式（８）のように、推定することができる性能パラメータの数は、センサの数以下でなければならない。この式（８）は、本明細書の末尾の付録１において導出する。例として示した図６のシステムでは、表２と表３のとおり、（８）式の両辺はそれぞれ１であり、成立している。
（性能パラメータｑの数）≦（センサｙの数） ・・・（８）

【0051】

ただし、状態ｘについては、センサ数の制約はなく、単に（Ａ、Ｃ）が可観測であればよい（このことについても付録１において示す）。このことから、推定したい生の性能パラメータｑのうち一部を状態ｘとして含ませることにより、式（７）を成立させることができる。換言すると、推定する性能パラメータの数をセンサの数より多くすることが可能になる。性能パラメータｑのうち一部を状態ｘとすることにより、（７）式の第１式で明らかなように状態ｘは動特性をもつこととなる。

【0052】

本実施の形態では、状態ｘにできるのは、動特性を持った積分要素の出力値であることを利用している。換言すると、積分値が状態ｘとなるので、性能パラメータｑに該当する要素を積分要素に置き換えることにより、性能パラメータｑを単なる係数から状態ｘの一部とすることができる。このことは、例えば性能パラメータｑを導入する乗算器を一次遅れ要素に置き換えることにより実現することができる。

【0053】

図８は、一次遅れ要素を説明する図である。図８（ａ）に示すように、所定の時定数τについて、一次遅れ要素への入力値Ｗは１／τ倍され、−１／τ倍された出力Ｗ´が加えられる。これらの和は積分されて出力値Ｗ´となる。

【0054】

図８（ｂ）において、図中の実線で示すようなステップ関数Ｗが一次遅れ要素に入力されると、図中の破線に示すようにステップ入力に次第に追随するような形状の出力Ｗ´が得られる。

【0055】

この出力Ｗ´は積分計算が必要であるので、状態ｘの一部とすることができる。一次遅れ要素の出力Ｗ´はカルマンフィルタにより調整され、一次遅れ要素への入力と差異が生じることになる。性能パラメータｑに相当する性能の経時的な変化は、一次遅れの出力値の入力値に対する比Ｗ´／Ｗにより求めることができる。

【0056】

図９は、図５の構成において性能パラメータ（流量特性変化係数）ｑｗｃｉを導入した乗算器を一次遅れ要素に置き換えたものである。このような一次遅れ要素により容積入口空気流量Ｗｃｉは積分値となり、状態ｘの一部とすることができる。

【0057】

この一次遅れ要素においては、ファン・圧縮機要素からの要素入口空気流量Ｗｃｉ ｍａｐを１／τ倍したものにこの一次遅れ要素の出力である容積流入空気質量Ｗｃｉを（−１／τ）倍したものを加え、さらに要素入口流量Ｗｃｉの変化量ΔＷｃｉを加えている。そして、これらの和を積分して出力となる容積流入空気質量Ｗｃｉを得ている。

【0058】

図９の構成において性能パラメータ（流量特性変化係数）ｑｗｃｉは使用されていないが、これに相当する値は要素入口空気流量Ｗｃｉ ｍａｐの容積流入空気質量Ｗｃｉに対する比として与えられる。換言すると、一次遅れ要素を導入することによって、本来は定数であった性能パラメータｑｗｃｉの経時的な変化を得ることができた。

【0059】

次に、このような構成を有するガスタービンエンジン推定装置においてシミュレーションを実施した。このシミュレーションは、離散時間で行った。このため、所定時間τの一次遅れ要素

【数8】

について、次のように離散化して適用した。この式の導出については、本明細書の末尾に付録２として記載した。

【数9】

【実施例1】

【0060】

ガスタービンエンジン推定装置の実施例１として、離散時間ｄｔと１次遅れの時間τをともに０．００２５とし、性能パラメータ（流量特性変化係数）のｑｗｃｉにおけるシステム雑音の大きさ（分散）を１×１０^−５としてシミュレーションを行った。

【0061】

このガスタービンエンジン推定装置の各ゲイン要素のゲインについては、次のように設定した。

【表5】

【0062】

図１０は、シミュレーションにより得られた容積入口空気流量Ｗ２５の時間変化を示している。図中において、実線は実機モデルで得られた真の値であり、破線はカルマンフィルタの一次遅れ要素により得られた推定値である。以下でも同様である。

【0063】

図中において、時間の経過により真の値と推定値の間に偏差が生じているが、この偏差の値は小さい。したがって、この実施例１において、カルマンフィルタはタービン実機の動作を良好に推定しているということができる。性能パラメータ（特性変化変数）を推定することなく、１次遅れを介して得られた状態ｗ２５ハットが推定できていることがいえる。

【実施例2】

【0064】

実施例２では、実施例１に対して性能パラメータ（流量特性変化係数）のｑｗｃｉにおけるシステム雑音大きさ（分散）を実施例１（１０^−５）に対して１×１０^１０と大きくしてシミュレーションを行った。他の構成は実施例１と同様である。図１１は、実施例２の容積入口空気流量Ｗ２５の時間変化を示している。時間の経過により真の値と推定値の間に偏差が生じているが、この偏差の値は実施例１と同様である。

【実施例3】

【0065】

実施例３では、実施例１に対して一次遅れの時間τを０．０２５と離散時間の１０倍に大きくしてシミュレーションを行った。他の構成は実施例１と同様である。図１２は、実施例３の容積入口空気流量Ｗ２５の時間変化を示している。時間の経過による真の値と推定値の間の偏差は実施例１と同様である。

【実施例4】

【0066】

実施例４では、実施例１に対して一次遅れの時間τを０．０２５と離散時間の１０倍に大きくし、さらに性能パラメータ（流量特性変化係数）のｑｗｃｉにおけるシステム雑音の大きさ（分散）を１×１０^１０と大きくしてシミュレーションを行った。図１３は、実施例４による容積入口空気流量Ｗ２５の時間変化を示している。真の値と推定値の偏差は、この実施例１〜４の中では最も小さくなった。

【比較例】

【0067】

図１４は、比較のために従来のカルマンフィルタを含むガスタービンエンジン推定装置の構成を示す図である。このガスタービンエンジン推定装置においては、本実施の形態とは異なり容積入口空気流量Ｗ２５の推定のために一次遅れ要素が設けられていない点において相違している。

【0068】

このような従来のカルマンフィルタでは、推定することができる性能パラメータの数は、センサの数以下でなければならなかった。なお、カルマンフィルタにおける推定値にはハットを付した。

【0069】

図１５は、図５に示したファン・圧縮要素と要素間容積の構成においてカルマンフィルタによる補正を施した図である。実機モデルの断熱効率ηは、カルマンフィルタによる補正により特性曲線のマップから求めた断熱効率ηｍａｐにさらに性能パラメータ（断熱効率特性変化係数）ｑηを乗じたηｍａｐ×ｑηとして与えられている。性能パラメータ（流量特性変化係数）ｑｗｃｉは、変化量Δｑｗｃｉを積分要素で積分して与えられている。

【0070】

要素間容積の質量Ｍの時間微分にはその変化量

【数10】

が加えられる。質量Ｍの時間微分は次の式で与えられ、性能パラメータｑｗｃｉによる補正を質量Ｍの時間微分の変化量によって代替するようにしている。

【数11】

【0071】

なお、本実施の形態においては、動的システムをガスタービンエンジンに適用したが、本発明はこれに限定されない。本発明は、各種エンジン、プラント等において、センサの数より多い性能パラメータを推定するために適用することができる。

【0072】

（付録１）
式（３）を参照すると、本実施の形態の可観測行列は、

【数12】

によって与えられる。この可観測行列がフルランクであれば可観測になる。

【0073】

ここで、

【数13】

【数14】

【数15】

である。

【0074】

行列Ｆ´、Ｂ´、Ｈ´は、行方向及び列方向にそれぞれ

【数16】

の要素を含んでいる。

【0075】

可観測行列は、次のように分解することができる。

【数17】

【0076】

ここで、行列Ｍ_Ｏ，Ｄ、Ｉ^（Ｄ）、Ｔ、Ｓ_Ｄは、行方向及び列方向にそれぞれ

【数18】

の要素を含んでいる。

【0077】

ここで、一般に行列のランクについて、次の関係式が成立する。

【数19】

この関係式を可観測行列が分解された式に適用すると、次のような関係式が得られる。

【数20】

【0078】

可観測行列のランク

【数21】

でフルランクとなるには、各分解要素のＩ^（Ｄ）、Ｔ、Ｓ_Ｄのランクがｎ＋ｐ以上となる必要がある。行列Ｓ_Ｄは、（ｎ＋ｍ）×（ｎ＋ｐ）行列であるため、ｎ＋ｐ≦ｎ＋ｍすなわちｐ≦ｍでないと、ランクがｎ＋ｐ以上になり得ない。

【0079】

したがって、ランクがｎ＋ｐとなるには、「性能パラメータｑの数ｐがセンサｙの数ｍ以下とならなくてはいけない」という制約があることになる。この結果を式（７）に適用すると、（Ａ、Ｃ）が可観測という必要から、性能パラメータの数がセンサの数以下でなければならないという式（８）が得られることになる。

【0080】

次に、行列Ｔには、性能パラメータを含まないエンジンモデルのシステム（Ｆ，Ｈ）の可観測行列Ｍ^（ｘ）_０が含まれている。

【数22】

【0081】

（Ｆ，Ｈ）が可観測のとき、Ｍ^（ｘ）_０はフルランクで

【数23】

となる。

【0082】

このとき、行列Ｔの構成から、明らかに

【数24】

となり、行列Ｔもフルランクとなる。したがって、性能パラメータを含まない状態ｘにセンサ数の制約はない。この結果を式（７）に適用すると、単に（Ａ、Ｃ）が可観測であればよいことになる。

【0083】

（付録２）
離散時間でシミュレーションするために１次遅れ要素を離散時間で記述する。１次遅れ要素は、次のように表される。

【数25】

この式は、次のように変形することができる。

【数26】

【0084】

この式の両辺に

【数27】

を掛けると次のようになる。

【数28】

さらに、次のようにまとめることができる。

【数29】

両辺を積分する。

【数30】

【0085】

ここで、ｔ−Δｔからｔの間にわたってＸ＝Ｘｎ（サンプルホールド）とすると、

【数31】

となり、さらに

【数32】

から次の式が得られる。なお、式中の斜線は両辺で相殺する数を示す。

【数33】

ここで、

【数34】

とすると、離散形の一次遅れの式は

【数35】

となる。

【符号の説明】

【0086】

１１、１１´ 第１ゲイン要素
１２、１２´ 第２ゲイン要素
１３、１３´ 第３ゲイン要素
１４、１４´ 第４ゲイン要素
１５、１５´ 第５ゲイン要素
１６、１６´ 第６ゲイン要素
１９、１９´ 第７ゲイン要素
２０、２０´ 第８ゲイン要素
２１、２１´ 第９ゲイン要素

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】