JP5300575B2 - ボイラ動特性試験解析装置 - Google Patents

Description

本発明は、高負荷変化率における発電出力と蒸気温度の追従性能を同時に向上させるのに好適な事業用ボイラ装置に係わり、特に、追従性能の最適化のための前処理としてのボイラ運転動特性試験解析装置に関する。

ボイラにおける発電出力（ＭＷ）制御に着目した水燃比制御又は給水制御の最適化のため、前処理としての動特性試験の解析方法に関する技術は、例えば、非特許文献１に紹介されていて公知である。この公知の解析方法では、ステップ応答の出力からムダ時間τと遅れ時間ＴおよびゲインＫを後述する方法にて算定し、出力を入力の（１次遅れ＋ムダ時間）近似で表現している（後述する図７を参照）。

種々の制御入力から種々の出力への信号伝達について、上述したような（１次遅れ＋ムダ時間）のシステムとして近似し、この近似に基づき出力から制御入力へのフィードバックの係数を決めている。しかし、この非特許文献１によるステップ応答解析方法では２つの課題がある。

１つ目は、先ず、（１次遅れ＋ムダ時間）系で近似できない入出力ループが取り扱えないことである。この方法では応答の初期特性のみを見てこれを（１次遅れ＋ムダ時間）系としてフィットしてしまうため、近似の度合いが非常に悪く、従って制御の最適化もできない。

２つ目は、この方法では、ムダ時間が実際の値より大きく算出されるという課題がある。これはボイラのような制御ループ間の干渉が強いときに好適な、ムダ時間値の異なる制御ループ間にバイアスを作用させるような制御方式、において課題となる。このようなバイアスではムダ時間の差だけの遅れあるいはリードが必要となるが、いずれもムダ時間差が大きいほどその制御回路による表現の誤差が大きくなる。

すなわち、非特許文献１のような従来技術による動特性試験結果の解析方法にあっては、（１次遅れ＋ムダ時間）で近似できない系を表現できないという課題と、ムダ時間が大きく算出され制御ループ間のバイアスを生成する際の遅れあるいはリード要素の制御回路による表現が大きな誤差を含む、という２つの課題が存在している。

「プロセス制御」、松原正一著、養賢堂出版、１９８６年発行

まず、上記非特許文献１に記載のような従来技術のボイラ動特性試験の解析方法の課題について図面を用いて説明する。初めに、従来技術によるボイラ動特性試験の解析方法と解析装置の構成を図１３に示す。以下の説明は、主として従来技術とその課題を示すが、この説明は、本発明における基礎、基盤となる技術をも示すものである。

図１３に示すボイラ動特性試験の解析装置は、応答波形記憶装置（データロガーに収納されたデジタル値のデータ）１０１と、応答解析装置１１２と、ＰＩＤ定数最適化装置１１４と、から構成され、応答解析装置１１２で得られたＰＩＤ最適値がボイラ制御装置１０６に入力される。

応答波形記憶装置１０１は、制御入力個数だけの行数と出力個数だけの列数を有する行列の形式で構成されるが、ここでは、制御入力がタービン・マスタｕ１、ボイラ・マスタｕ２と水燃比ｕ３の３個であり、出力が発電出力ｙ１、主蒸気圧力ｙ２及び１ＳＨ（１次過熱器）入口エンタルピｙ３の３個である場合を説明する。この場合には、図１３の行列の中身は図３の通りであり（入出力はいずれも無次元化されている）、行はα＝１（タービン・マスタ）、α＝２（ボイラ・マスタ）とα＝３（水燃比）、列はｊ＝１（発電出力）、ｊ＝２（主蒸気圧力）とｊ＝３（１次過熱器入口エンタルピ）である。ここで、応答波形記憶装置１０１に例示された９個の波形は、図７に示すようなステップ状の波形を入力（タービンマスタ（タービン出力指令値：例．後述するガバナの閉アクション）、ボイラマスタ（ボイラの蒸気流量指令値：例．後述する給水・燃料の同時減少）、水燃比（給水流量と燃料量との比）としたときに、ボイラの試運転時にその動特性試験の結果得られた出力（発電出力、主蒸気圧力、１次ＳＨ入口エンタルピ）、すなわち、応答トレンドである。

具体的には応答波形記憶装置には３×３＝９通りの（時間、応答値）の組５１１，５１２，５１３，５２１，５２２，５２３，５３１，５３２と５３３が、データ収録装置の時間刻みごとに格納されている。この９通りが順不同であるが一つずつ、くまなくスイッチを介して応答解析装置１１２にてデータ解析される。なお、この応答解析装置１１２でのデータ解析は人力もしくは自動化されたアルゴリズムによるものである。

応答解析装置１１２は、次に詳述する図７の方法によりデータを解析しすべての応答トレンドを１次遅れ時間Ｔ、ムダ時間τと応答のゲインＫでパラメトライズする。すなわち、応答解析装置１１２は、入力が応答波形記憶装置１０１からの波形（応答トレンド）であり、その出力が１次遅れ時間Ｔ、ムダ時間τ、ゲインＫである。このパラメトライズにより図１３の行列の対角成分に対応する制御ループのＰＩＤ最適値が、例えばＺｉｅｇｌｅｒ−Ｎｉｃｈｏｌｓの方法によりＰＩＤ定数最適化装置１１４にて得られ、ボイラ制御装置１０６に入力される。ここでもＰＩＤ定数最適化装置１１４は人力（手計算）であることが多い。

次に、応答解析装置１１２でのトレンドの解析方法について、図７を用いて応答パラメータの算出方法を示す。入力ｕ（５０１）を時点Ｏ（５０８）にてステップ変化させたときの出力ｙ（５０２）（実線）の時間変化が実線で示されている。従来技術においては出力ｙの変曲点Ｉ（５０３）を応答結果のデジタル値もしくは目視により求め、それとゼロ点あるいは漸近値５０５との交点をそれぞれＡ（５０４）およびＢ（５０７）とし、Ｂからの垂線と直線ＯＡとの交点をＣ（５０６）とする。これにより、ムダ時間をτ＝ＯＡ、１次遅れ時間をＴ＝ＡＣ、さらにゲインをＫ＝ＢＣと算出する（入力ｕ、出力ｙともに無次元化されていて、ｕのステップ幅は１に規格化されている）。これらのＫ、Ｔとτを使って入出力応答の（１次遅れ＋ムダ時間）近似が得られるが、それを破線６０２で示す（破線６０２は、τ，Ｔ，Ｋを用いた入出力応答における実線５０２の近似特性曲線）。この方法によりボイラの入出力応答の一例とそれらを近似したトレンドを図８に示す。

図８は第１行がタービン・マスタのステップ、第２行がボイラ・マスタのステップさらに第３行が水燃比のステップを示しており、いずれもステップ幅はマイナス（減）である（図７の例示ではステップ入力は増であるが、図８ではステップ減の例示）。出力が列方向に示されており、第１列が発電出力、第２列が主蒸気圧力さらに第３列が１次ＳＨ入口エンタルピのそれぞれ応答トレンドである。すなわち、実線５１１はタービン・マスタのステップ減に対する発電出力の応答トレンド、５１２は同じく主蒸気圧力の応答トレンド、また５１３は１次過熱器入口エンタルピの応答トレンドである。また、実線５２１，５２２と５２３はそれぞれボイラ・マスタのステップ減に対する発電出力、主蒸気圧力と１次過熱器入口エンタルピの応答トレンドであり、実線５３１，５３２と５３３は水燃比ステップ減の下での発電出力、主蒸気圧力と１次過熱器入口エンタルピの応答トレンドである。

図８において、破線を併記している応答トレンドがあり、これはタービン・マスタのステップ応答に関しては６１２、ボイラ・マスタのステップ応答については６２１と６２２、さらに水燃比のステップ応答では３つの出力すべてで６３１，６３２と６３３であり、これらの破線は、実線の応答トレンドを近似する（１次遅れ＋ムダ時間）系の応答トレンドである。

ところが、図８から直ちに分かるように、実線５１１，５１３と５２３を近似できる単純な（１次遅れ＋ムダ時間）系は存在しない。すなわち、実線５１１は前述したようにタービン・マスタのステップ変化の下での発電出力の変化を示し、ガバナ（タービンの入り口側に設けられた調速機）の閉アクションに伴い初期一時的に発電出力が下がるが、同時に主蒸気圧力が高まるため発電出力が増方向に転じこれら両者の作用の重ね合わせによって漸近的に発電出力が元の値に戻っていくことが分かる。

また、実線５１３のタービン・マスタのステップ減に伴う１次過熱器入口エンタルピの応答では、初期のガバナ閉のアクションで蒸気流量が減じることによる入熱量の減少が初期のエンタルピの低下を引き起こすものの、漸近的には主蒸気圧力が高めにキープされることが効いてエンタルピは元の値より高めの値に漸近することを示している。さらに、実線５２３ではボイラ・マスタのステップ状の減少により給水も燃料も同時に減少するが、初期は給水減少の効果がより早く現われてエンタルピの上昇となっていること、また漸近的には逆に燃料の減少効果がより顕著にあらわれエンタルピが元の値より低めキープとなることを示している。

これらの３ケースではいずれも出力は初期にいったん上昇あるいは下降した後で漸近的にゼロあるいは元の値より高め・低めの点に収束していくが、このようなトレンドは単純な（１次遅れ＋ムダ時間）近似では表現できない。さらに、たとえ（１次遅れ＋ムダ時間）系として近似できている制御ループであっても、図８のステップ初期を拡大した図９に見るように、ムダ時間を大きく評価し過ぎている、という問題がある。

すなわち、図９の第３行第１列の水燃比のステップ減での発電出力の応答トレンド５３１とその図７の従来アルゴリズムになる（１次遅れ＋ムダ時間）近似６３１の比較を例とすると、実機応答のムダ時間真値τ３１０（５００）はたかだか２０秒程度と見積もられるが、図７の従来方法によればムダ時間がτ３１（６００）が１００秒程度と算定されている。すなわち、ボイラの試運転での動特性試験ではムダ時間はトレンド５３１からその真値τは符号５００で表されて２０秒であるのに対して、このトレンド５３１を基にした応答解析装置１１２では近似トレンド６３１となり、この近似トレンドのムダ時間は符号６００で表され、１００秒程度となる。このような状況が引き起こす問題について、図を用いて以下詳しく説明する。

図１０は干渉する制御ループ間（例．ｕ２−ｙ２ループとｕ３−ｙ３ループとの間）に互いにバイアス（例．バイアス項４５２と４５３）を作用することで干渉を無くする、あるいは抑制する制御方式を示すものであって、このようなバイアスを前置することによりボイラ制御のＰＩＤ制御パラメータの最適化がより容易になるものである。ここでは、ボイラ入力として、図３で説明したタービン・マスタｕ１、ボイラ・マスタｕ２と水燃比ｕ３の内で、ボイラ・マスタｕ２（４０２）と水燃比ｕ３（４０３）のバイアスについて示している。

図１０において、Ｇ２２（８２２），Ｇ３２（８３２），Ｇ２３（８２３）とＧ３３（８３３）は、それぞれ、ボイラ・マスタｕ２（４０２）から主蒸気圧力４９２へ、水燃比ｕ３（４０３）から主蒸気圧力４９２へ、ボイラ・マスタｕ２（４０２）から１ＳＨ入口エンタルピｙ３（４９３）へ、水燃比ｕ３（４０３）から１ＳＨ入口エンタルピｙ３（４９３）へ、の作用を示す記号であり、制御装置ではなくボイラ・タービンの動特性を表現するものである。ボイラ・マスタｕ２（４０２）は主として主蒸気圧力ｙ２（４９２）を制御するがその作用は同時に１次ＳＨ入口エンタルピｙ３（４９３）への外乱ともなる。

同様に、水燃比ｕ３（４０３）は主として１ＳＨ入口エンタルピｙ３（４９３）を制御するがその作用は同時に主蒸気圧力ｙ２（４９２）への外乱としても作用する。そこで、図１０の制御回路では、ボイラ・マスタｕ２（４０２）に水燃比ｕ３（４０３）からのバイアス項４５２を加算し、また逆に、水燃比ｕ３（４０３）にはボイラ・マスタｕ２（４０２）からのバイアス項４５３を加算している。これにより、バイアス項４５２あるいは４５３が外乱項とキャンセルするように、バイアス項４５２あるいは４５３の信号を補償器Ｄ３２（９３２）あるいはＤ２３（９２３）で生成する。

ところが、これらの制御ループ｛ｕ２，ｕ３｝→｛ｙ２，ｙ３｝では、ｕ２→ｙ２，ｙ３の伝達におけるムダ時間がｕ３→ｙ２，ｙ３の伝達におけるムダ時間より小さい。すなわち、ボイラ・マスタｕ２（４０２）の方が水燃比ｕ３（４０３）より早く主蒸気圧力ｙ２（４９２）を、また、１ＳＨ入口エンタルピｙ３（４９３）についても動かすことが出来る。たとえば、制御ループのムダ時間が図１１と図１２の記載値（図７の従来方法で算定したムダ時間の値）であったとすると、主蒸気圧力の制御ループｕ２→ｙ２において外乱となる水燃比ｕ３（４０３）の作用４３２はτ３２＝２０秒が経って初めてｙ２（４９２）に現われる。

一方、ボイラ・マスタｕ２（４０２）はτ２２＝５秒でｙ２（４９２）に作用が現われるので、水燃比ｕ３（４０３）からのバイアスは２０−５＝１５だけの遅れを持たせておけば同時にｙ２（４９２）への作用となって現れ、これをキャンセルすることが可能になる（ボイラマスタへのバイアスとして１５の遅れを持たせれば、入力のｕ２にはムダ５（４２２）が加算されて２０の遅れとなり、符号４３２からの２０の遅れと同時に現れてキャンセルされる）。同様に図１２に示すように、水燃比ｕ３（４０３）→１ＳＨ入口エンタルピｙ３（４９３）へのボイラ・マスタｕ２（４０２）からの外乱をキャンセルするには、１２−１８＝−６秒だけ、すなわち６秒のリード時間要素が必要である。

ところが、これらの遅れ時間要素９３２やリード時間要素９２３にあっては、その遅れもしくはリード時間の値が小さいほど実現性、従って制御性が良い。特に、リード時間要素は微分要素を少なくも１回以上適用することで表現され、微分要素は微分パラメータによるパラメータ依存性を持つのでこのリード時間は小さければ小さいほど良好である。

従って、ムダ時間の算定値は小さければ小さいほど良いが、図９の破線に見るように、実際のムダ時間５００よりも破線（図７に示す従来の解析方法）は大きなムダ時間６００を算定してしまっている（従来方法ではムダ時間を過大評価で算定している）、という問題がある。

本発明の目的は、試運転動特性試験で得られたトレンドデータを従来の単純な（一時遅れ＋ムダ時間）系ではなくて、時間微分特性を考慮したアルゴリズムで解析することで、（一時遅れ＋ムダ時間）系で近似できない制御ループに対しても適用可能とし、さらに、ムダ時間を実際の値に近い値として算定することを可能とするボイラ運転動特性試験解析装置を提供することにある。

前記課題を解決するために、本発明は主として次のような構成を採用する。
ガバナにより発電出力を制御し、給水又は燃料で蒸気圧力及び蒸気温度を制御するボイラ運転制御を最適化する前処理としてのボイラ運転動特性試験装置において、
水燃比を含む複数の制御入力を印加したときの発電出力を含む複数の出力の応答波形を収納する応答波形記憶装置と、前記応答波形記憶装置からの応答波形を解析し、一次遅れ、ムダ時間、ゲインを含むパラメータとして出力する応答解析装置と、前記応答解析装置からのパラメータを用いて、ボイラ制御ループのＰＩＤ最適値を算出するＰＩＤ定数最適化装置と、前記応答解析装置からのパラメータを用いて、ボイラ制御ループ間の非干渉化のための非干渉補償装置に付与する非干渉補償定数を計算する非干渉補償定数計算装置と、を有し、
前記応答解析装置は、前記応答波形に対して、１次遅れ＋ムダ時間で近似できる制御ループに対して当該制御ループに近似する一の伝達関数を用いて解析するとともに、１次遅れ＋ムダ時間で近似できない制御ループに対しても当該制御ループに近似する他の伝達関数を用いて解析し、さらに、前記一の伝達関数と前記他の伝達関数は、解析されるムダ時間が前記応答波形による実際のムダ時間に近い値として評価されるものであるボイラ運転動特性試験装置。

また、前記応答解析装置は、応答波形を、１次遅れとムダ時間で近似される波形と、１次遅れとムダ時間とそれらの時間微分との線形和で近似される波形との２種類に分類し、さらに、前記近似される波形から算出されるムダ時間を、ムダ時間真値と前記ムダ時間真値からの偏差とに分解して多項式近似で表現する機能を有するボイラ運転動特性試験解析装置。

本発明によれば、試運転動特性試験で得られるトレンド・データを従来のような単純な（１次遅れ＋ムダ時間）系ではなく、その時間微分特性を考慮してタービン・マスタやボイラ・マスタの操作に対する応答をより詳細に近似したパラメータをトレンド・データの解析から得ることができる。

また、本発明によれば、ムダ時間についても従来のような単純なアルゴリズムで得られる値よりもより実際に近い値を算出することができる。

これら二つの新規要素によって、互いに干渉するボイラ制御ループの非干渉化のためのバイアス入力のための補償要素をより精度よく計算でき、特にバイアスにおいて必要となる遅れ時間とリード時間要素については、その遅れ時間値とリード時間値をいずれも従来アルゴリズムによるムダ時間を用いるときに比べて小さな値に抑えることができ、制御回路における数式演算の近似精度を向上することができる。

本発明の実施形態に係るボイラ運転動特性試験解析装置の全体構成及び関連構成を示すブロック図である。 本実施形態に係るボイラ運転動特性試験解析装置の内部構成である応答解析装置の構成と機能を表す図である。 本実施形態に係るボイラ運転動特性試験解析装置の内部構成である応答波形記憶装置に収録されるトレンド・データの詳細波形を示す図である。 本実施形態に係るボイラ運転動特性試験解析装置の内部構成である応答解析装置で解析された波形に対して算出されたパラメータを用いて、図３の収録した応答波形に近似した近似トレンド・データの詳細波形を示す図である。 本実施形態による近似トレンド・データに対して、従来技術の解析による近似トレンド・データと比較したトレンド・データの初期特性を示す図である。 本実施形態に係るボイラ運転動特性試験解析装置の内部構成である非干渉補償定数計算装置からの出力を用いて非干渉補償装置を前置したボイラ系に対するステップ応答トレンド・データに対して、非干渉補償装置の前置が無いボイラ系に対するステップ応答トレンド・データと比較したトレンド・データの特性を示す図である。 従来技術に関する応答波形解析方法におけるアルゴリズムを明示する原理的な説明図である。 従来技術に関する応答波形解析方法による近似トレンド・データを併記したステップ応答のトレンド・データの詳細を示す図である。 図８に示すトレンド・データの初期特性を示す図である。 本実施形態に係るボイラ運転動特性試験解析装置の内部構成である非干渉補償定数計算装置の演算を説明するための図である。 図１０における一つのバイアスにおいてムダ時間差の遅れを介したムダ時間の補償における信号の流れを示す図である。 図１０におけるもう一つのバイアスにおいてムダ時間差の進み（リード時間）を介したムダ時間の補償における信号の流れを示す図である。 従来技術に関するボイラ運転動特性試験解析装置の全体構成を示すブロック図である。

本発明の実施形態に係るボイラ運転動特性試験解析装置について、図面を参照しながら以下詳細に説明する。図１は本発明の実施形態に係るボイラ運転動特性試験解析装置の全体構成及び関連構成を示すブロック図である。図２は本実施形態に係るボイラ運転動特性試験解析装置の内部構成である応答解析装置の機能を表す図である。図３は本実施形態に係るボイラ運転動特性試験解析装置の内部構成である応答波形記憶装置に収録されるトレンド・データの詳細波形を示す図である。図４は本実施形態に係るボイラ運転動特性試験解析装置の内部構成である応答解析装置で解析された波形に対して算出されたパラメータを用いて図３の収録した応答波形に近似した近似トレンド・データの詳細波形を示す図である。図５は本実施形態による近似トレンド・データに対して、従来技術の解析による近似トレンド・データと比較したトレンド・データの初期特性を示す図である。

また、図６は本実施形態に係るボイラ運転動特性試験解析装置の内部構成である非干渉補償定数計算装置からの出力を用いて非干渉補償器を前置したボイラ系に対するステップ応答トレンド・データに対して、非干渉補償器の前置が無いボイラ系に対するステップ応答トレンド・データと比較したトレンド・データの特性を示す図である。図７は従来技術に関する応答波形解析方法におけるアルゴリズムを明示する原理的な説明図である。図８は従来技術に関する応答波形解析方法による近似トレンド・データを併記したステップ応答のトレンド・データの詳細を示す図である。図９は図８に示すトレンド・データの初期特性を示す図である。

また、図１０は本実施形態に係るボイラ運転動特性試験解析装置の内部構成である非干渉補償定数計算装置の演算を説明するための図である。図１１は図１０における一つのバイアスにおいてムダ時間差の遅れを介したムダ時間の補償における信号の流れを示す図である。図１２は図１０におけるもう一つのバイアスにおいてムダ時間差の進み（リード時間）を介したムダ時間の補償における信号の流れを示す図である。図１３は従来技術に関するボイラ運転動特性試験解析装置の全体構成を示すブロック図である。

図１において、本発明の実施形態に係るボイラ運転動特性試験解析装置は、応答波形記憶装置１０１と、応答解析装置１０２と、非干渉補償定数計算装置１０３と、ＰＩＤ定数最適化装置１０４と、から構成される。ここで、非干渉補償定数計算装置１０３の計算結果は非干渉補償装置１０５に入力され、ＰＩＤ定数最適化装置１０４の計算結果はボイラ制御装置１０６に入力される。

図１に示す応答波形記憶装置１０１は、図１３に示す従来技術における応答波形記憶装置１０１と同様であり、上記の「発明が解決しようとする課題」の欄で説明したように、図３の例示のように、９通りのトレンド・データ（ボイラ試運転で得られた生のデータ）が順不同であるが一つずつくまなくスイッチを介して応答解析装置１０２に送られてデータ解析される。ここで、本実施形態における応答解析装置１０２は、図７に示す従来の解析方法（上記の「発明が解決しようとする課題」の欄で説明した方法）を拡張した解析アルゴリズムを具備した自動化装置である（図７に示す，単純な（１次遅れ＋ムダ時間）系による解析アルゴリズムを拡張したもの）。この拡張した解析アルゴリズムが本発明の特徴の１つであり、応答解析装置１０２の内部構成と機能を図２に示す。

本実施形態に関する応答解析装置１０２は、応答波形記憶装置１０１からの出力データである応答波形に対して、ステップ入力が行われた時間（ｔｉｍｅ＝０）におけるＩｃ（Ｉｎｉｔｉａｌ Ｃｏｎｄｉｔｉｏｎ；初期状態）やステップ入力後に応答が安定したＦｃ（Ｆｉｎａｌ Ｃｏｎｄｉｔｉｏｎ；最終状態）からの飛び出しの無い波形１と、ＩｃまたはＦｃからの飛び出しの有る波形２と、を分類する機能１０２０を有している。さらに、応答解析装置１０２は、波形１のパラメータを計算する機能１０２１を有していて、その波形１のパラメータとして、ゲインＫ、１次遅れ時間Ｔ、ムダ時間τ、実機応答のムダ時間真値τ０を出力する。また、応答解析装置１０２は、波形２のパラメータを計算する機能１０２２を有していて、その波形２のパラメータとして、ゲインＫ、１次遅れ時間Ｔ、定数ａ（後述する）、ムダ時間τ、実機応答のムダ時間真値τ０を出力する。

再度、定性的に説明すると、応答解析装置１０２においては、図３と図８に示すように、応答波形記憶装置１０１からの波形にＩｃ（Ｉｎｉｔｉａｌ Ｃｏｎｄｉｔｉｏｎ；初期状態）やＦｃ（Ｆｉｎａｌ Ｃｏｎｄｉｔｉｏｎ；最終状態）からの飛び出しがあるか無いかで分類付け又は仕分けし、波形のパラメータを出力するものである。ここで、Ｉｃからの飛び出しがあるのは、図３の例で、１行１列（５１１）、１行３列（５１３）、２行３列（５２３）であり、飛び出しがあれば、１次遅れの微分値を１次遅れに加算する。Ｆｃの値に対してオーバーシュートして飛び出すような波形例は、図３には例示されていない。

本実施形態に関する応答解析装置１０２では、図３の内で、第１行の第１列（５１１）と第３列（５１３）、および第２行の第３列（５２３）の３つのケース（ａ）と、それ以外のケース（ｂ）とで、それぞれ別のデータ解析（１０２１と１０２２）がなされる。

まず、（ｂ）のケース（（１次遅れ＋ムダ時間）系で近似できる入出力ループ）は、初めに従来技術における図６の方法にてムダ時間τ、１次遅れ時間ＴとゲインＫが採取される。すなわち、ゲインＫはステップ入力５０１の漸近値に対するステップ応答結果５０２の漸近値（垂線ＢＣの長さ）の比率として算定される。また、ムダ時間はＯＡ長として、遅れ時間はＡＣ長として算出される（図７を参照）。

次いで、ムダ時間τは真値τ０とそれからの偏差Δτ＝τ−τ０に分解される。従来技術の解析方法から得られるムダ時間τは、図９から明らかなように真のムダ時間（図９の５００）より大きい値（図９の６００）が算出されることが多い。そこで、より真実に近い図９の５００（収録したトレンドデータ自体のムダ時間）をτ０と取り、このτ０をムダ時間要素で表現し、残りのΔτ＝τ−τ０を高次遅れで近似する。この高次遅れとしてはＰａｄｅ多項式近似を採用する。

すなわち、上記ケース（ｂ）である、９−３＝６、すなわち６つのケースでは入出力ループは、
Ｇ＝Ｋ×ｅｘｐ（−τ０×ｓ）×Ｐ２（ｓ；Δτ）×１／（１＋Ｔｓ） （１）
で表現される。このＧは、図５に示す３行１列の特性６３１（従来方法）から特性７３１を求める関数である。ここで、Ｋはゲイン、τ０はムダ時間の真値（実際値）、ｓは微分要素、ｅｘｐ（−τ０×ｓ）は純粋ムダ時間要素、Ｐ２（ｓ；Δτ）はＰａｄｅ多項式近似、１／（１＋Ｔｓ）は１次遅れ要素、をそれぞれ表す。また、Ｐ２（ｓ；Δτ）は以下で与えられる。

Ｐ２（ｓ；Δτ）＝
１／（１＋Δτ×ｓ＋Δτ^２／２！×ｓ^２＋・・・＋Δτ^ｎ／ｎ！×ｓ^ｎ） （２）
このｎは例えばｎ＝２である。

次に、ケース（ａ）（（１次遅れ＋ムダ時間）系で近似できない入出力ループ）では、１次遅れとその時間微分との線形和がとられ、和の係数（ａ^−１）は漸近値がデータと合うように算出される。すなわち、３つのケースの入出力ループにおける伝達関数の数式表現は、
Ｇ＝Ｋ×ｅｘｐ（−τ０×ｓ）×Ｐ２（ｓ；Δτ）×（１＋ａ^−１ｓ）／（１＋Ｔｓ） （３）
であり、Ｐ２（ｓ；Δτ）は式（２）と同じ多項式である。式（３）に示す伝達関数によって、図５の１行１列に示す波形５１１から波形７１１が求まる。式（３）が式（２）と異なるのは、第４番目の項の分子が１から（１＋ａ^−１ｓ）となることである。

繰り返して説明すると、式（１）と式（２）は、図３に示す６つのケースに適用し、式（３）と式（２）は３つのケース（１−１，１−３，２−３）に適用するものであり、式（１）と式（３）は、図５の例えば、波形７３１と波形７１１をそれぞれ求めるための関数であって、波形５３１と波形５１１にそれぞれ可成りよい近似をあらわすものである。式（３）の線形和（１＋ａ^−１ｓ）の１がＴ，τ系であり、ａ^−１ｓが時間微分を表しており、式（１）と式（３）は、その基本的な機能上の差異としては、オーバーシュート／アンダーシュートの表現で云えば、式（１）は図３に示すすべてのｙがＩｃとＦｃの内側にあり、式（３）は図３に示す或るｙがＩｃとＦｃから飛び出すものを表す。図３に示す応答波形を図５に示す一点鎖線に変換する手段（関数）は、応答波形（実線）からムダ時間の真値τ０の観察と、求めたＴ，Ｋ，τ，ａとによって、Δτ＝τ−τ０とし、Ｐ２（ｓ；Δτ）を計算し、基本式のｅ（−τ×ｓ）をｅ（−τ０×ｓ）に置き換えて変換するものである。

ここで、上記式（３）において、Ｋ，Ｔ，τとａはそれぞれ以下のように算定される。はじめに、Ｋは図７と同様に出力漸近値の入力ステップ量に対する比率として算出される。そこで、ムダ時間τ、遅れ時間ＴでゲインがＫである系を、ステップ応答データから差し引くと、これは、ゲインがＫ×ａ^−１でムダ時間τ、遅れ時間Ｔの時間変化率を表す微分特性の系であることがわかる。この微分特性の系を時間積分すればムダ時間τ、遅れ時間Ｔの系になるので、このτとＴは図７の方法で算出する。ここで、ａは、図３に示す或るｙがＩｃとＦｃから飛び出す挙動を規定する定数（ｓｅｃ^−１）であり、ＫとＫａ^−１の２つを計測することができるので、ａ^−１はＫａ^−１／Ｋで計算することができる。

なお、ステップ応答データから差し引くべき系のムダ時間τも遅れ時間Ｔも予め分からないため、推定値から出発して推定と算定値が等しくなるまで収束計算を行う。実際には図３に示すような系にあっては２、３回の収束でムダ時間τと遅れ時間Ｔを計算でき、収束に伴う計算長時間化や発散などの困難はない。

上述したような本実施形態に関する応答解析装置１０２における解析アルゴリズムでの伝達関数の数式表現による応答トレンドの近似が、図４の一点鎖線７１１〜７３３として表されている。また、その初期特性を拡大したものを図５に示す。図５では一点鎖線（本実施形態の特徴の１つである解析アルゴリズムの数式表現によるトレンドデータ）とともに破線６１２，６２１，６２２，６３１，６３２と６３３も併記するが、これは図７の従来技術になる（１次遅れ＋ムダ時間）近似である。図４と図５で一点鎖線は９ケースすべてを近似している。図５に示すように、従来技術による破線は６ケースしか近似できていないことに注視すべきである。

再度繰り返して説明すると、図５において、水燃比入力−発電出力のトレンドを例とした場合、実線の符号５３１は応答波形記憶装置に収納された応答波形であってムダ時間は実線が下降傾向をし始めるまでの時間であってこれが実際のムダ時間τ０である。破線の符号６３１は図７に示す従来の解析アルゴリズムで解析された応答トレンドであり、これによると、ムダ時間τは下降開始の時間であってτ０の略数倍（略３倍）の時間である。そこで、本実施形態では、式（１）と式（２）に示す伝達関数を用いて実線５３１を解析すると一点鎖線の符号７３１となり、この一点鎖線７３１によると、一点鎖線が下降し始めるまでの時間であるムダ時間はτ０に相当近似することとなる（略１倍）。

さらに、タービンマスタ−発電出力（図３に示す１行１列）のトレンドを例とした場合、図７に示す従来の解析アルゴリズムでは近似することができないが（図８で破線の特性が示されていない）、式（３）と式（２）に示す伝達関数を用いて、収納された応答波形である実線５１１を解析すると、一点鎖線７１１となって実線に可成り近似させることができる（図４と図５に示す符号７１１を参照）。応答波形５１１に漸近性のある波形が得られ、この波形を基にして、実際のムダ時間τ０に近いムダ時間τ、応答波形に近似した一次遅れＴ、応答波形に近似したゲインＫが出力可能となる。

次に、非干渉補償定数計算装置１０３では制御ループ間の非干渉化のために、本実施形態にあっては次のような補償要素が設計され、その各パラメータが算定される。すなわち、ボイラ・マスタｕ２と水燃比ｕ３の非干渉化を例として（図１０を参照）、詳細に数式を用いて説明すると以下のとおりである。本実施形態では、非干渉化バイアスとして、ボイラ・マスタｕ２から水燃比ｕ３への補償要素Ｄ２３と、水燃比ｕ３からボイラ・マスタｕ２への補償要素Ｄ３２の２つが前置される。はじめに、補償要素Ｄ２３は次式によって計算される。なお、補償要素Ｄ３２とＤ２３は、図１に示す非干渉補償装置１０５の構成要素である。

Ｄ２３＝−Ｋ２３／Ｋ３３×（１＋Ｔ２３・ｓ）／（１＋Ｔ３３・ｓ）
×Ｐ２（ｓ；Δτ２３）／Ｐ２（ｓ；Δτ３３）
×（１−（τ２３０−τ３３０）・ｓ＋（τ２３０−τ３３０）^２・ｓ^２／２！
＋・・・＋（τ２３０−τ３３０）^ｍ・ｓ^ｍ／ｍ！）
×（１＋ａ２３^−１・ｓ） （４）
であり、通常ｍ＝２をとる。Ｋ２３はｕ２系からｕ３系への非干渉化バイアスのループにおけるゲイン、Ｋ３３はｕ３系からｙ３系へのループにおけるゲイン、Ｔ２３とτ２３は同様なループの一次遅れ時間とムダ時間を表す。また、τ２３０は、同様な非干渉化バイアスのループにおけるムダ時間の真値を表す。

また、Ｄ３２は以下で与えられる。
Ｄ３２＝−Ｋ３２／Ｋ２２×（１＋Ｔ３２・ｓ）／（１＋Ｔ２２・ｓ）
×Ｐ２（ｓ；Δτ３２）／Ｐ２（ｓ；Δτ２２）
×ｅｘｐ（−（τ３２０−τ２２０）・ｓ） （５）
以上のとおり、補償要素Ｄ２３とＤ３２は、微分要素ｓ，ｓ^２、遅れ要素１／（１＋Ａ・ｓ）、Ｐ２（ｓ；Δτ）および純粋ムダ時間要素ｅｘｐ（−τ０・ｓ）の組み合わせで構成されるが、非干渉補償定数計算装置１０３ではこれらの補償要素に対して上式に示すパラメータ値が計算・設定される。また、微分要素は微分係数γを含む不完全微分要素で表現され、このパラメータγも非干渉補償定数計算装置１０３で設定される。式（４）のＤ２３）、式（５）のＤ３２は非干渉化補償装置１０５の構成要素であり、図１０に示す符号９２３、符号９３２に示すとおりである。

また、ＰＩＤ定数最適化装置１０４は、式（４）および式（５）で計算され、非干渉補償定数計算装置１０３から出力されるパラメータを用いて設定された補償要素を前置したボイラ／タービン・プラントに対して、ＰＩＤによる制御を最適化するための比例定数、積分時間と微分時間を非干渉化された制御ループそれぞれに対し算出する。すなわち、図６に非干渉化補償を前置した系のステップ特性を示すが、上段図の実線２２２と２２３がボイラ・マスタｕ２のステップ入力による出力ｙ２（主蒸気圧力）とｙ３（１ＳＨ入口エンタルピ）トレンドであり、下段図の実線２３２と２３３が水燃比ｕ３のステップ入力による出力ｙ２（主蒸気圧力）とｙ３（１ＳＨ入口エンタルピ）トレンドであり、併記されている破線３２２，３２３，３３２と３３３は非干渉化補償を前置していないときの応答である。

図６に示す特性から分かるように、応答は「対角的」になっている。すなわち、ボイラ・マスタｕ２の作用はほぼ主蒸気圧力ｙ２にのみ、また、水燃比ｕ３の作用はほぼ１ＳＨ入口エンタルピｙ３にのみ作用し、ボイラ・マスタｕ２から１ＳＨ入口エンタルピｙ３への外乱作用と水燃比ｕ３から主蒸気圧力ｙ２への外乱作用は十分に抑制されている。

そこで、ＰＩＤ定数最適化装置１０４では、ｕ２→ｙ２すなわち実線２２２と、ｕ３→ｙ３すなわち実線２３３のそれぞれの応答を（１次遅れ＋ムダ時間）でフィットし、そのゲインＫ２もしくはＫ３、１次遅れ時間Ｔ２もしくはＴ３と、ムダ時間τ２とτ３を算定する。この算定は従来技術の図７の手法を適用する。そして、得られたゲイン、１次遅れ時間とムダ時間を使って次のようにＰＩＤ制御の最適パラメータを算出する。先ずｕ２→ｙ２の制御ループでは最適値は以下のとおりである。

Ｐ２＝α２×Ｔ２／（Ｋ２×τ２） （６）
ＴＩ２＝β２×τ２ （７）
ＴＤ２＝γ２×τ２ （８）
ここで、Ｐ２はｕ２→ｙ２ループにおける比例定数、ＴＩ２は積分定数、ＴＤ２は微分定数のそれぞれの最適値を表す。

次に、ｕ３→ｙ３の最適値は以下のとおりである。

Ｐ３＝α３×Ｔ３／（Ｋ３×τ３） （９）
ＴＩ３＝β３×τ３ （１０）
ＴＤ３＝γ３×τ３ （１１）
ここで、Ｐ３はｕ３→ｙ３ループにおける比例定数、ＴＩ３は積分定数、ＴＤ３は微分定数のそれぞれの最適値を表す。

また、ＰＩＤ制御装置の構成は、制御偏差を入力して操作量を出力する伝達関数が
Ｃ２（ｓ）＝Ｐ２×（１＋１／（ＴＩ２×ｓ）＋ＴＤ２×ｓ） （１２）
Ｃ３（ｓ）＝Ｐ３×（１＋１／（ＴＩ３×ｓ）＋ＴＤ３×ｓ） （１３）
で与えられるものである。以上においてパラメータは、例えばα２＝α３＝１．２、β２＝β３＝２、γ２＝γ３＝０．５である。

ここで、本発明の実施形態を取り纏めると、次のような構成を備え、機能乃至作用を奏することを特徴とするものである。すなわち、先ず、（１次遅れ＋ムダ時間）で近似できない制御ループに対して、当該制御ループに近似する伝達関数の数式表現（図４に示す一点鎖線の伝達関数）を求め、さらに、ムダ時間τを実際のムダ時間値τ０に出来る限り近い値として評価して伝達関数を構成する。すなわち、（１次遅れ＋ムダ時間）で近似できない制御ループについては、その時間微分との線形和を取ることで応答波形との漸近特性を表現する。また、ムダ時間については、（１次遅れ＋ムダ時間）あるいはその時間微分との線形和で表現される数式において、ムダ時間τを真値と考えられる値τ０とそれからの誤差Δτ＝τ―τ０に分解し、ΔτをＰａｄｅ近似で数式表現するものである。

そして、本実施形態の具体的な構成として、ガバナで出力を制御し、給水もしくは燃料で蒸気圧力及び蒸気温度を制御するボイラ制御を最適化するための前処理としての運転動特性試験装置において、応答波形記憶装置１０１、応答解析装置１０２、非干渉補償定数計算装置１０３とＰＩＤ定数最適化装置１０４を有するものである。

ここで、前記応答解析装置１０２は、応答波形（応答波形記憶装置１０１に収納された波形であり、図５で実線、例．符号５３１）について、１次遅れとムダ時間で近似される波形（図５で破線、例．６３１）と、１次遅れとムダ時間とその時間微分との線形和で近似される波形（図５で一点鎖線、例．７３１，７１１）との２種類に分類し、前記近似される波形６３１から算出されるムダ時間を真値（τ０）と真値からの偏差（Δτ＝τ−τ０）とに分解して多項式近似（式（２））で表現する機能を奏させるものである。

さらに、前記非干渉補償定数計算装置１０３は、ムダ時間真値の差で演算される時間値が正であるときは遅れ要素、また前記時間値が負であるときはリード要素をそれぞれ計算する機能（式（４）と式（５））を奏するものである。また、前記ＰＩＤ定数最適化装置１０４は、前記非干渉補償定数計算装置１０３にて算出され、出力されたパラメータ値を用いた非干渉補償要素を介したボイラ系のステップ応答に対して、ゲイン、遅れ時間とムダ時間の組み合わせとしてＰＩＤ定数を算出する機能（式（６）〜式（１３））を奏するものである。

そして、上述した構成を備えることにより、試運転動特性試験において得られた出力応答トレンドをすべて数式表現で近似できるため、制御パラメータの最適化が可能である。また、ムダ時間が実際の値にできるだけ近い値として算定されるため、バイアス前置において必要な遅れあるいはリード時間値が可能な限り小さく設定でき、このため遅れ時間あるいはリード時間要素の制御回路としての実現性能が良いものとなる。

１０１ 応答波形記憶装置
１０２ 応答解析装置
１０３ 非干渉補償定数計算装置
１０４ ＰＩＤ定数最適化装置
１０５ 非干渉補償装置
１０６ ボイラ制御装置
ＢＣ ゲインＫ（図７）
５００ 実機応答のムダ時間真値（図９）
６００ 従来の解析方法によるムダ時間（図９）
ＯＡ ムダ時間τ（図７）
ＡＣ 一次遅れ時間Ｔ（図７）

Claims (4)

1. ガバナにより発電出力を制御し、給水又は燃料により蒸気圧力及び蒸気温度を制御するボイラ運転制御を最適化する前処理としてのボイラ運転動特性試験装置において、
   水燃比を含む複数の制御入力を印加したときの発電出力を含む複数の出力の応答波形を収納する応答波形記憶装置と、
   前記応答波形記憶装置からの応答波形を解析し、一次遅れ、ムダ時間、ゲインを含むパラメータとして出力する応答解析装置と、
   前記応答解析装置からのパラメータを用いて、ボイラ制御ループのＰＩＤ最適値を算出するＰＩＤ定数最適化装置と、
   前記応答解析装置からのパラメータを用いて、ボイラ制御ループ間の非干渉化のための非干渉補償装置に付与する非干渉補償定数を計算する非干渉補償定数計算装置と、を有し、
   前記応答解析装置は、前記応答波形に対して、１次遅れ＋ムダ時間で近似できる制御ループに対して当該制御ループに近似する一の伝達関数を用いて解析するとともに、１次遅れ＋ムダ時間で近似できない制御ループに対しても当該制御ループに近似する他の伝達関数を用いて解析し、さらに、前記一の伝達関数と前記他の伝達関数は、解析されるムダ時間が前記応答波形による実際のムダ時間に近い値として評価されるものである
   ことを特徴とするボイラ運転動特性試験装置。
2. 請求項１において、
   前記応答解析装置は、応答波形を、１次遅れとムダ時間で近似される波形と、１次遅れとムダ時間とそれらの時間微分との線形和で近似される波形との２種類に分類し、
   さらに、前記近似される波形から算出されるムダ時間を、ムダ時間真値と前記ムダ時間真値からの偏差とに分解して多項式近似で表現する機能を有する
   ことを特徴とするボイラ運転動特性試験解析装置。
3. 請求項１において、
   前記非干渉補償定数計算装置は、ムダ時間真値の差で演算される時間値が正であるときは遅れ要素、また、前記時間値が負であるときはリード要素をそれぞれ計算する機能を有する
   ことを特徴とするボイラ運転動特性試験解析装置。
4. 請求項１において、
   前記ＰＩＤ定数最適化装置は、前記非干渉補償定数計算装置にて計算され出力されたパラメータを用いた非干渉補償要素を介したボイラ系のステップ応答に対して、ゲイン、遅れ時間とムダ時間の組み合わせとしてＰＩＤ定数を算出する機能を有する
   ことを特徴とするボイラ運転動特性試験解析装置。

Images