(57)【要約】【課題】  運転中のプラントに外乱を与えることなく、ボイド反応度係数などの、炉心動特性に影響する主要なパラメータ、および定常ゆらぎを引き起こす雑音源を随時に推定し、炉心の監視に利用し得る方法を提供すること。【解決手段】  沸騰水型原子炉の中性子束信号および炉心流量信号の変動成分から、周波数解析により統計的特徴量を求め、前記沸騰水型原子炉炉心部の動特性モデルに複数の雑音源を仮定することにより前記統計的特徴量の理論値を求め、前記統計的特徴量と前記統計的特徴量の理論値との差を最小にするように動特性モデル内のパラメータまたは雑音源の大きさを求め、前記パラメータまたは雑音源の大きさを監視に用いることを特徴とする沸騰水型原子炉の炉心監視方法。【選択図】  図１
2005134291.tif 000002
【特許請求の範囲】
【請求項１】
沸騰水型原子炉の中性子束信号および炉心流量信号の変動成分から、周波数解析により統計的特徴量を求め、
前記沸騰水型原子炉炉心部の動特性モデルに複数の雑音源を仮定することにより前記統計的特徴量の理論値を求め、
前記統計的特徴量と前記統計的特徴量の理論値との差を最小にするように前記動特性モデル内のパラメータまたは雑音源の大きさを求め、
前記パラメータまたは雑音源の大きさを監視に用いることを特徴とする、沸騰水型原子炉の炉心監視方法。
【請求項２】
請求項１記載の方法において、
前記統計的特徴量として、炉心流量信号から中性子束信号への伝達関数を用いて、前記動特性モデル内のパラメータであるボイド反応度係数を求めることを特徴とする、沸騰水型原子炉の炉心監視方法。
【請求項３】
請求項１記載の方法において、
前記統計的特徴量として、前記炉心流量信号および前記中性子束信号のパワースペクトル密度、ならびに前記炉心流量信号と前記中性子束信号との間の伝達関数およびコヒーレンス関数を用いて前記動特性モデル内のパラメータである前記ボイド反応度係数および前記雑音源の大きさを同時に求めることを特徴とする、沸騰水型原子炉の炉心監視方法。
【請求項４】
請求項１ないし３の何れかに記載の方法において、
前記動特性モデルおよび複数の前記雑音源から前記炉心流量信号および前記中性子束信号に関する統計的特徴量を求める際に、前記炉心流量信号および前記中性子束信号の観測ノイズをパラメータとして持たせることを特徴とする、沸騰水型原子炉の炉心監視方法。
【請求項５】
請求項４記載の方法において、
周波数領域の統計的特徴量の差を用いて前記パラメータを推定する方式に替えて時間領域の最尤推定法を用いることを特徴とする、沸騰水型原子炉の炉心監視方法。
【請求項６】
請求項１ないし５の何れかに記載の方法において、
推定された動特性パラメータと、炉心動特性モデルを用いて、前記動特性モデルの固有値から求まる炉心の安定性尺度を用いることを特徴とする、沸騰水型原子炉の炉心監視方法。
【請求項７】
請求項１ないし６の何れかに記載の方法において、
前記動特性モデルとして、原理から求められる物理モデルと、観測時系列データからフィッティングにより求められるフィッティングモデルとを組み合わせて用いることを特徴とする、沸騰水型原子炉の炉心監視方法。  
【請求項８】
請求項１ないし６の何れかに記載の方法において、
物理モデルから推定するパラメータを、ボイド反応度係数と、高周波帯域を持つボイドノイズと、低周波帯域を持つボイドノイズとすることを特徴とする、沸騰水型原子炉の炉心監視方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、沸騰水型原子炉の炉心監視方法に係り、とくにボイド反応度係数などの炉心動特性パラメータならびに炉心の出力変動を引き起こす雑音源の大きさを測定して炉心の監視に用いる方法に関する。
【背景技術】
【０００２】
沸騰水型原子炉は、反応度増加による出力上昇時には燃料温度上昇に伴うドップラー効果やボイド率の増加に伴う負のボイド反応度効果により、出力係数は常に負に維持できるように設計されている。このような沸騰水型原子炉の基本的特性は、実機プラントの運転により既に確認されているものであるが、運転中にボイド反応度係数を測定することができれば、原子力プラントの安全な運転に役立つのみならず、従来の運用に比べてより現実的な炉心運用の可能性が期待される。
【０００３】
従来、ボイド反応度係数を測定する方法としては、プラントに人為的外乱を加え、その応答から評価する方法が提案されている。この方法は、ボイド反応度係数に限らず、対象系の動特性を規定するパラメータの推定方法として一般的に知られるものである（例えば非特許文献１参照）。
【０００４】
以下、図９を参照して、このような動特性パラメータを推定する従来の方法について説明する。
【０００５】
図９において、従来の動特性パラメータ推定方法では、人為的な外乱である制御入力を実機プラントに与え、その制御入力ｕおよび出力信号ｙを時間Ｔだけ測定する。ボイド反応度係数を推定する場合の外乱は、例えば再循環流量設定点の変更であり、ｕとして炉心流量、ｙとして中性子束を測定する。
【０００６】
動特性シミュレーションコードに推定対象パラメータの初期推定値Ｋを与え、入力に測定されたｕを印加して出力推定値ｙｙｐを求める。下式に示す評価基準、 
Ｊ＝∫0T（ｙｙｐ−ｙ）２ dt
を最小とするようなパラメータ推定値の補正値δＫを算出し、（Ｋ＋δＫ）を新たな推定値として再度シミュレーション解析を行い、出力推定値ｙｐを求める。これを繰り返し行なって、収束したときの推定値をＫの測定値とするものである。
【０００７】
このようなボイド反応度係数の測定方法においては、外乱を印加することで中性子束に有意の変化を生じさせて比較的大きなＳ／Ｎ比で測定することが可能であり、精度よく測定できる利点を持っている。
【０００８】
しかしながら、この方法では、炉心流量信号に大きな外乱を与える必要があり、商用運転原子炉で日常的な計測に用いることは難しい。
【０００９】
一方、特許文献１では、炉心流量信号と中性子束信号の定常ゆらぎを用いて、ボイド反応度係数を求める方法を考案している。
【００１０】
図１０は、この方法を示しており、信号測定結果から伝達関数を求め、さらに、その伝達関数の特徴量として、平均ゲイン、減衰率、共鳴周期などを求めて、その特徴量から相関式を用いてボイド反応度係数を推定するものである。
【００１１】
この方法は、定常ゆらぎ信号からボイド反応度係数を求めることができる初めての方法として重要であるが、相関式を予めシミュレーションで求めておく必要がある。
【特許文献１】特開2001-042078号公報
【特許文献２】特願2003−074055号
【非特許文献１】C.L.Medler & Chih-Chi Hsu: An Algorithm for Nonlinear System Parameter Identi-fication, IEEE Trans. Automatic Control, pp.726-728, Dec. 1969
【発明の開示】
【発明が解決しようとする課題】
【００１２】
上述した従来のボイド反応度係数の測定方法においては、プラント状態に有意の過渡変化が現れるだけの外乱を印加することが必要である。
【００１３】
しかしながら、発電を目的とする原子炉の場合、このような外乱の印加は、電力の安定供給という面からみてできるだけ避けねばならず、通常は起動試験の中でのみ行われるものである。
【００１４】
したがって、従来のボイド反応度係数の測定方法は、発電開始後のサイクル運転中の原子炉には適用できないという課題があった。
【００１５】
また、定常ゆらぎを用いる方法では、予め相関式をシミュレーションで求めておく必要があり、実用に際しての使い難さが問題であった。
【００１６】
本発明は上述の点を考慮してなされたもので、運転中のプラントに外乱を与えることなく、ボイド反応度係数などの、炉心動特性に影響する主要なパラメータ、および定常ゆらぎを引き起こす雑音源を随時に推定し、炉心の監視に利用し得る方法を提供することを目的とする。
【課題を解決するための手段】
【００１７】
上記目的を達成するため、本発明では、
沸騰水型原子炉の中性子束信号および炉心流量信号の変動成分から、周波数解析により統計的特徴量を求め、前記沸騰水型原子炉炉心部の動特性モデルに複数の雑音源を仮定することにより前記統計的特徴量の理論値を求め、前記統計的特徴量と前記統計的特徴量の理論値との差を最小にするように動特性モデル内のパラメータまたは雑音源の大きさを求め、前記パラメータまたは雑音源の大きさを監視に用いることを特徴とする、沸騰水型原子炉の炉心監視方法、
を提供するものである。
図１は、本発明に係る沸騰水型原子炉の炉心監視方法の基本的流れを示す。ゆらぎ成分からなる時系列信号は、周波数変換を用いて、パワースペクトル密度、伝達関数、コヒーレンス関数といった、周波数の関数としての統計的特徴量に変換することができる。
【００１８】
また、炉心部の動特性モデルを用いると、雑音源および観測ノイズを仮定することで、統計的特徴量に対応させることのできる量を理論的に求めることができる。この理論的に求める量を、統計的特徴量の理論値と呼ぶ。そして、統計的特徴量と統計的特徴量の理論値とを一致させるように、動特性モデル内のパラメータ、ならびに雑音源の大きさを求めることができる。
【００１９】
本発明で用いる物理モデルは、沸騰水型原子炉の１点炉近似動特性モデルと呼ばれるもので、本来、非線形のモデルである。ただし、定常値からの小さな変動を扱うということで線形化近似を行い、さらに中性子信号に関して即発跳躍近似を行って下記(1)式の形の状態方程式とする。この(1)式内のＸが状態変数で、Ｙが観測変数である中性子束信号を意味する。
【００２０】
また、炉心流量Ｗ(ｔ)は入力信号として扱って、中性子束から炉心流量へのフィードバックは無視している。また、Ｑはシステム雑音源、Ｒは観測ノイズで、白色の周波数特性を仮定している。そして、状態変数の中で、ボイド発生に関わる雑音源として、ボイドのランダムな発生を駆動するＶ(ｔ)に加えて、チャンネル入口流量のゆらぎや蒸気流量のゆらぎに起因した圧力変動などを代表する雑音源として、ＷＬ(ｔ)を状態変数として仮定している。
【００２１】
後者に関しては、結果的には、ボイド変動を通して中性子に影響するため、モデル上のＶ(ｔ)と同じ場所に印加されることになる。ただし、物理的仮定として、前者が広帯域、すなわち時定数(ＴＶ)の短い雑音源を想定したのに対して、後者は低周波帯域、すなわち時定数(ＴＬ)の長い雑音源を想定して、このモデルを用いることにした。
【数１】
2005134291.tif 000003
【００２２】
図２は、本発明の炉心動特性パラメータ監視法で用いる各記号を列挙したものである。
【００２３】
図３は、上記線形モデルを周波数領域のブロック線図で表現したものである。ここでは、観測ノイズも含めて、上記(1)式と同等のモデルとなっている。
【００２４】
この周波数領域のブロック線図内の記号を用いて、炉心流量から中性子束への伝達関数を求めると、下記(2)式が得られる。炉心の動特性を支配する最も重要なパラメータは、ボイド反応度係数 Kαである。この(2)式では、ボイド反応度係数 Kαが分子に入っていることから伝達関数のゲインに比例し、しかも分母にも入っているため、安定度を支配する極の位置を支配することが分かる。
【００２５】
本発明では、この伝達関数が、ボイド反応度係数に応じて変化することを利用して、観測データから周波数解析で計算される伝達関数と、下記(2)式の伝達関数の理論値とを一致させるよう、非線形計画法を用いてボイド反応度係数 Kαを求めている。
【数２】
2005134291.tif 000004
【００２６】
また上記(1)式より、観測信号の自己・相互パワースペクトル密度（ＡＰＳＤ／ＣＰＳＤ）は下記(3)式、
【数３】
2005134291.tif 000005
【００２７】
により求めることができる。この(3)式中、ＰＵ(S) は炉心流量のＡＰＳＤを、ＰＱ は白色性を仮定したシステムノイズを、PR は白色性を仮定した観測ノイズをそれぞれ示し、＊は複素共役を示すものとする。
【００２８】
こうして求めた伝達関数やＡＰＳＤなどは、観測されるゆらぎ信号から求まる統計量に対応しているため、両者の差を適切に定義した評価関数を作ることで、動特性モデル内の諸パラメータを推定することができる。
【００２９】
しかしながら、モデル内には、動特性パラメータ以外に、雑音源を定義するパラメータまで含まれており、これらの全てを単一の評価関数で一括して推定することは困難であると考えられる。そのため、本発明では、周波数領域で定義する下記の評価関数を用いる。
【数４】
2005134291.tif 000006
【００３０】
ここで、AP(f) は、中性子束信号の自己パワースペクトル密度の物理モデルからの計算値と、ゆらぎ信号の周波数変換による計算値とを示す。また、CF(f) は炉心流量と中性子束のコヒーレンス関数を、さらにG(f) は炉心流量から中性子束への伝達関数のゲインをそれぞれ示す。また、Ｗ１（ｆ），Ｗ２（ｆ），Ｗ３（ｆ）は適当に設定可能な重み関数で、例えばＷｉ（ｆ）＝１とかＷｉ（ｆ）＝１／ｆ（ｉ＝１，２，３）といった形に設定するものである。
【００３１】
これら３種類の特徴量のうち、自己パワースペクトル密度AP(f) の形状は、雑音源の大きさおよびボイド反応度係数の値に依存して変化する。また、伝達関数ゲインG(f) は、ボイド反応度係数の値のみに依存して変化する。さらに、コヒーレンス関数CF(f) は、複数雑音源それぞれの大きさの比に応じて変化する。
【００３２】
このように、本発明によれば、異なる特徴量を組み合わせた評価関数にすることで、動特性モデル内のボイド反応度係数および雑音源の大きさを、評価関数への感度の違いを利用して同時に推定することができる。
【００３３】
また、上記(1)式の時間領域でのモデル表現を用いると、観測時系列信号にカルマンフィルタを適用して尤度関数を求めることができる。この尤度関数は、動特性パラメータや雑音源の大きさの関数となっているので、これを最大にするように、ボイド反応度係数のような動特性パラメータや、雑音源の大きさを推定することもできる。ただし、この場合、上記(1)式の入力として用いている炉心流量の動的挙動も、観測ノイズを含めてモデル化する必要がある。
【００３４】
そのため、本発明では、下記(5)式のように、事前に自己回帰モデルを用いて炉心流量のゆらぎを観測ノイズまで含めて求めておき、上記(1)式と組み合わせて用いることで、特許文献２と同じ手法を適用して、ボイド反応度係数や雑音源の大きさを推定することができる。
【００３５】
ここで、特許文献２の手法とは、原子炉の中性子束の一点近似モデルと冷却材温度の動特性モデルとを用い、原子炉の中性子束信号および冷却材温度の観測信号が入力されて状態方程式に変換し、この状態方程式にカルマンフィルタを適用して対数尤度を求め、この対数尤度を最大化することで温度反応度係数を推定するものである。
【００３６】
本発明では、中性子束信号を用いることは同様であるが、冷却材温度の観測信号が炉心流量信号の変動成分に置き換わり、温度反応度係数に相当するものがボイド反応度係数となる。
【数５】
2005134291.tif 000007
【発明の効果】
【００３７】
本発明によれば、定常運転中に観測されるプラント信号の変動成分のみを用いて測定するため、プラントに外乱を与える必要が無く、サイクル運転中にオンラインで随時あるいは定期的にボイド反応度係数を測定することができる。
【００３８】
また、事前に動特性シミュレーションコードなどで作成する相関式も必要なく、簡便にボイド反応度係数を推定できる。
【００３９】
さらには、炉心の出力ゆらぎを引き起こす雑音源を推定したり、炉心の出力変動の安定性の尺度である減幅比を推定することで、より詳細な炉内監視に役立てることができる。
【発明を実施するための最良の形態】
【００４０】
（実施例１）
図１は、本発明の第１の実施例を示すもので、入力のプラントデータとして、炉心流量信号と中性子束信号とを用いる。この両信号は、デジタル値に変換された時系列データであるものとすると、それを周波数変換して平均操作を加えることで、パワースペクトル密度（ＡＰＳＤ）、伝達関数（ＴＦゲイン、位相）、コヒーレンス関数（ＣＦ）という、周波数の関数として表現できる特徴量に変換できる。
【００４１】
一方、図３に示す周波数領域で定義した炉心動特性モデルにおいて、ボイド反応度係数 Kα、ならびに高周波数帯域ボイドノイズ（標準偏差＝σV）、低周波数帯域ボイドノイズ（標準偏差＝σL）の初期値を仮定し、さらに炉心流量信号のパワースペクトル密度を入力とすることで、上記(2),(3)式にしたがって、中性子束のパワースペクトル密度、ならびに炉心流量と中性子束との間のコヒーレンス関数および伝達関数を求めることができる。
【００４２】
さらに、先に求めたこれらの計測値との差を上記(4)式で求める。この(4)式は、下記のように、Kα，σV，σLの関数としてみなすことができる。ここで、fmin，fmaxは、ゆらぎ信号として注目すべき周波数帯域を示したもので、この第１の実施例では、０．０１Ｈz〜１．２５Ｈzの周波数帯域を選んでいる。
【数６】
2005134291.tif 000008
【００４３】
この初期値として与えた３つのパラメータを、非線形計画法によるアルゴリズムにしたがって順次変更して行くと、上記(6)式を最小にする値として、３つのパラメータの最適値を決めることができる。
【００４４】
図４は、このように最適パラメータを決めた後の、パワースペクトル密度、コヒーレンス関数、伝達関数（ゲインと位相）を示す。計測値と物理モデルによる周波数領域での特徴量とがよく一致しており、求められた３つのパラメータが適切であることが分かる。
【００４５】
図５は、本発明の第１の実施例における炉心動特性パラメータ推定後のモデルで評価した、中性子束のパワースペクトル密度の各雑音源からの寄与の分解表現図である。この図５において、Ｗは炉心流量寄与を、Ｑｖは高周波数帯域ボイドノイズ寄与を、ＱLは低周波数帯域ボイドノイズ寄与をそれぞれ示している。
【００４６】
そして、中性子束のパワースペクトル密度の計算値が、入力信号である炉心流量と、高周波数帯域雑音源と、低周波数帯域雑音源の３つの寄与から構成されていることを示すものである。各雑音源からの寄与をこのような形で表示することで、例えば何らかの異常により雑音源の大きさが変わったような場合に、その変化を捉えることができる。したがって、ボイド反応度係数の変化の監視と併せて、炉心の異常監視に役立てることができる。
【００４７】
（実施例２）
図６は、本発明の第２の実施例による雑音源の推定結果を示すものである。ここでは、図３の周波数領域のモデルではなく、上記(1)式と上記(5)式とを組み合わせたモデルを用い手いる。そして、中性子束信号および炉心流量信号を入力して、特許文献２に記載の手法で用いたカルマンフィルタを適用することにより尤度関数値を求め、さらにその尤度関数値を最大にするように、３つのパラメータ Kα，σV，σLを非線形計画法(最尤推定法)で求めた結果を示している。
【００４８】
この図６は、最適パラメータを決めた後の、パワースペクトル密度、コヒーレンス関数、伝達関数（ゲインと位相）を示す。そして、計測値と物理モデルによる周波数領域での特徴量とは、周波数帯域の特徴量から求めた図４と同様によく一致しており、求められた３つのパラメータが適切であることが分かる。
【００４９】
このように、本発明では、周波数領域だけではなく、時間領域のデータからも炉心監視に必要なパラメータを推定することができ、監視の多様性を図ることができる。
【００５０】
（実施例３）
図７は、本発明の第３の実施例における、ボイド反応度係数の経年変化を推定した結果を示している。この第３の実施例は、第１の実施例による監視を長期間にわたって定期的に行い、沸騰水型原子炉の運転の推移に伴うボイド反応度係数の変化を調べたものである。
【００５１】
この図７に併せて示した炉心性能計算コードによる予測計算値からも分かるように、運転の継続に伴い燃料が燃焼することで、ボイド反応度係数は、若干ではあるが増加傾向にあることが分かる。この傾向が、本発明による監視方法で的確に捉えられており、炉内監視に役立て得ることが明らかである。
【００５２】
また、この第３の実施例では、上記(1)式の動特性モデルが、線形の微分方程式で表現されていることを利用して、その固有値を求め、その固有値から炉心の出力変動の安定性を調べることも行っている。
【００５３】
図８は、本発明の第３の実施例における炉心安定性の減幅比が経年変化する様子の推定結果を示すものである。これは、安定性の尺度としてよく用いられる減幅比の結果を示したものであるが、ボイド反応度係数の増加に伴って、僅かではあるものの炉心の安定性が悪化していることが分かる。
【００５４】
この安定性は、従来の技術では、中性子束のゆらぎ信号に自己回帰モデルなどをフィッティングして求め、その自己回帰モデルの固有値ないしインパルス応答から、安定性の尺度である減幅比を求めている。図８にはその結果も比較して示してあるが、本発明の方法に比べてばらつきが大きく、本発明の有用性が分かる。
【図面の簡単な説明】
【００５５】
【図１】本発明による炉心動特性パラメータ監視法のブロック線図。
【図２】本発明の炉心動特性パラメータ監視法で用いる物理モデルのパラメータの説明図。
【図３】本発明の炉心動特性パラメータ監視法で用いる物理モデルの周波数領域で表したブロック線図。
【図４】本発明の第1の実施例における炉心動特性パラメータ推定後の統計的特徴量の計測値（破線）とモデル計算値（実線）の比較図。
【図５】本発明の第1の実施例における炉心動特性パラメータ推定後のモデルで評価した中性子束のパワースペクトル密度の各雑音源からの寄与の分解表現図。
【図６】本発明の第２の実施例における最尤推定法を用いた炉心動特性パラメータ推定後の統計的特徴量の計測値（破線）とモデル計算値（実線）との比較図。
【図７】本発明の第３の実施例におけるボイド反応度係数の、経年変化の推定結果の例を示す図。
【図８】本発明の第３の実施例における炉心安定性減幅比の、経年変化の推定結果の例を示す図。
【図９】従来の外乱印加によるボイド反応度係数推定法のブロック線図
【図１０】従来の定常ゆらぎ信号を用いたボイド反応度係数推定法のブロック線図。
【図１】
2005134291.tif 000009
【図２】
2005134291.tif 000010
【図３】
2005134291.tif 000011
【図４】
2005134291.tif 000012
【図５】
2005134291.tif 000013
【図６】
2005134291.tif 000014
【図７】
2005134291.tif 000015
【図８】
2005134291.tif 000016
【図９】
2005134291.tif 000017
【図１０】
2005134291.tif 000018