JP2004279338A

JP2004279338A - 原子炉の温度反応度係数計測方法

Info

Publication number: JP2004279338A
Application number: JP2003074055A
Authority: JP
Inventors: Harutsugu Mori; 治嗣森; Yuichi Kagami; 雄一加賀見; Shigeru Kanemoto; 茂兼本; Mitsuhiro Enomoto; 光広榎本; Shigeo Ehata; 茂男江畑; Tetsuo Tamaoki; 哲男玉置
Original assignee: Toshiba Corp; Tokyo Electric Power Co Inc
Current assignee: Toshiba Corp; Tokyo Electric Power Co Holdings Inc
Priority date: 2003-03-18
Filing date: 2003-03-18
Publication date: 2004-10-07

Abstract

【課題】中性子束と冷却材温度の観測信号から温度反応度係数を従来方法より、より短時間で計測することにある。
【解決手段】原子炉の中性子束の一点近似動特性モデル１ａと冷却材温度の動特性モデル１ｂで構成され前記原子炉の中性子束と冷却材温度の観測信号が入力される状態方程式変換部１と、この状態方程式変換部１より得られる状態方程式にカルマンフィルタを適用して対数尤度を計算する尤度計算部３と、この尤度計算ステップで求められた対数尤度を最大化することで温度反応度係数を推定する温度反応度係数推定部４とを備えて、前記原子炉の中性子束と冷却材温度の観測信号から温度反応度係数を測定する。
【選択図】図１

Description

【０００１】
【発明の属する技術分野】
本発明は、カルマンフィルタによる最尤法を用いて原子炉の温度反応度係数を計測する方法に関する。
【０００２】
【従来の技術】
従来、原子炉における反応度や減速材温度係数は、原子炉の制御、安全確保および設計の妥当性を検討する上で必ず必要なものであり、次のような方法により求められていた。
【０００３】
即ち、温度反応度係数は、次の一点炉近似動特性モデルの反応度と温度の比例定数Ｋ_Ｔで定義される。
【０００４】
ｄｎ／ｄｔ＝（ρ−β）ｎ／ｌ＋λＣ＋Ｖ_１
ｄＣ／ｄｔ＝βｎ／ｌ−λＣ＋Ｖ_２
ｄＴ／ｄｔ＝ｕ（ｔ）＋Ｖ_３ …… （１）
ρ＝Ｋ_ＴＴ＋ｃｏｎｓｔ
ここで、ｎ：中性子束（定常値で規格化）、Ｃ：遅発中性子先行核密度（定常値で規格化）、Ｔ：冷却材温度（℃）、ｕ：温度外乱、ρ：温度反応度フィードバック、Ｋ_Ｔ：温度反応度係数（０．４５ｃ／℃）、ｌ：中性子寿命（０．００００４ｓｅｃ）、β：遅発中性子割合（０．００５６）、λ：遅発中性子先行核崩壊定数（０．０８ｓｅｃ^−１）、Ｖ_ｉ：システムノイズである。
【０００５】
上記（１）式には、中性子束や温度のランダムな挙動を模擬するためのシステムノイズを加えているが、併せて下記のように計測にあたって必然的に加わる観測ノイズを考慮した観測モデルを仮定する必要がある。
【０００６】
従来、この観測ノイズは、システムノイズと明確に区別して表現されないことがあるが、最適なパラメータ推定法を考えるには、この両者の区別が必要なことから、これらを下式のように定義しておく。
【０００７】
Ｙ_１＝ｎ＋Ｗ_１
Ｙ_２＝Ｔ＋Ｗ_２ …… （２）
従来の温度反応度係数の推定では、まず炉内乃至炉外のセンサーで計測される中性子束信号から、下記の逆時間公式を用いて反応度を計算している。
【０００８】
【数１】

【０００９】
さらに、この反応度と温度計測値に下記の一次回帰式を仮定し、最小２乗法で係数を求める。ここで、εは観測誤差を、δは定常値からの偏差を表している。
【００１０】
δρ＝Ｋ_Ｔδ_Ｔ＋ε ……（４）
このとき、最小２乗法による係数は、観測信号の相関係数から次のように求められる。
【００１１】
Ｋ_Ｔ＝＜δρδＴ＞／＜δＴδＴ＞ ……（５）
ただし、＜＞は時間平均を意味するものとする。これは、従来用いられている典型的な推定手法である。
【００１２】
一方、最近では、原子炉の反応度測定方法および減速材温度係数測定方法として、計測信号で必ず存在する観測ノイズを考慮した因子分析法による推定方法が提案されている（特許文献１）。上記（４）式において、観測ノイズを考慮すると、下記式のように書き換えられる。
【００１３】
δρ＝Ｋ_Ｔ（δＴ＋ε_Ｔ）−ε_ρ ……（６）
このとき、（５）式と同じ最小２乗法で温度反応度係数を求めると、入力信号の観測ノイズが分母に入ってくるため、下式のように＜ε_Ｔε_Ｔ＞だけ推定値がバイアスを持っていることが分かる。
【００１４】
Ｋ_Ｔ＝＜δρδＴ＞／（＜δＴδＴ＞＋＜ε_Ｔε_Ｔ＞）……（７）
ここで、観測ノイズは、白色で状態変数と無相関、すなわち、＜δＴε_Ｔ＞＝０と仮定した。
【００１５】
このような場合、本来は下記の観測ノイズを明確に考慮したモデル、すなわち
【数２】

【００１６】
に従って、係数を推定する必要がある。これは、反応度と温度に比例関係があるという過程のもとでは、共通因子である真の温度変化を共通因子とした因子分析モデルになっていることが分かる。
【００１７】
従って、因子分析アルゴリズムを適用することにより、入力信号側が観測ノイズを持っていてもバイパスを持たない温度反応度係数を推定することが可能になる。
【００１８】
ただし、実用的には温度をランプ状に変化させて温度反応度係数を求める場合、入力の温度信号を移動平均処理することで、観測ノイズの影響を減らして、バイアスのない係数を推定することが一般的に行なわれている。また、ランプ状外乱の場合、（８）式で反応度と温度を下式のようにそれぞれ時間の関数としてフィッティングし、その傾きの比として温度反応度係数を求めることでバイアスのない推定値を得ることもできる。
【００１９】
δρ＝ａ_ρｔ＋ε_ρ
δＴ＝ａ_Ｔｔ＋ε_Ｔ
Ｋ_Ｔ＝ａ_ρ／ρ_Ｔ …… （９）
また、もう一つの原子炉の減速材温度係数測定方法として、前記した（６）式を周波数変換してその比をとることで温度係数を求めるノイズ分析手法による温度係数測定法がある（非特許文献１）。
【００２０】
この減速材温度係数測定法は、（６）式を周波数領域で表現した式と考え、両辺に温度信号の共役複素数をかけて平均を取ることで、下式のようなコヒーレンス関数と伝達関数（周波数応答関数）を得ることができるものである。
【００２１】
【数３】

【００２２】
ここで、＜＞は周波数領域での平均化を意味している。ただし、（１０）式を求める際に、観測ノイズは、状態変数と無相関という仮定を用いた。このとき、入力である温度信号の観測ノイズ＜ε_Ｔε_Ｔ＞が無視できる場合、伝達関数は、すべての周波数領域で一定値となり、これが温度反応度係数となる。また、温度外乱が十分大きければ、反応度の観測ノイズも無視できるため、コヒーレンス関数は１．０になる。
【００２３】
しかしながら、通常はこのような条件は満たさないため、例えば、ランプ状の外乱の場合、低周波数領域ではコヒーレンス関数が１．０に近くなるが、高周波領域では小さな値となる。この高周波数領域では、温度の観測ノイズも無視できなくなるため、すべての周波数領域の平均として求める温度反応度係数は、バイパスを持ってくることになる。
【００２４】
このように上記の周波数領域を用いた温度係数推定法では、コヒーレンス値の大きいところの伝達関数ゲインの平均で温度反応度係数を求めているが、これは下式のようにコヒーレンス関数の大きい部分だけの伝達関数の平均をとることで、観測ノイズ＜ε_Ｔε_Ｔ＞の影響の小さい部分を選択して平均することで、バイアスを小さくしていることになる。
【００２５】
【数４】

【００２６】
【特許文献１】
特開２００１−２５５３９１号公報
【００２７】
【非特許文献１】
日本原子力学会１９９８年秋の大会、Ｆ１０、上田他
「ノイズ分析手法による減速材温度係数測定法の実証」
【００２８】
【発明が解決しようとする課題】
このように従来の原子炉の温度反応度係数計測方法は、時間領域と周波数領域における反応度係数の推定法であるが、いずれも反応度の推定を介した後のフィッティングであり、直接観測される中性子束信号と温度信号に含まれる本来の情報を十分に活用しきれていない。このため、従来の方法では、温度反応度係数の推定値が真値からバイアスしたり、ばらつきが大きく計測に時間がかかるという問題があった。
【００２９】
本発明は上記のような事情に鑑みてなされたもので、中性子束と冷却材温度の観測信号から温度反応度係数を従来方法より、より短時間で計測することができる原子炉温度反応度係数計測方法を提供することを目的とする。
【００３０】
【課題を解決するための手段】
本発明は上記の目的を達成するため、次のような方法により原子炉の温度反応度係数を計測するものである。
【００３１】
請求項１に対応する発明は、原子炉の中性子束の一点近似動特性モデルと冷却材温度の動特性モデルで構成され前記原子炉の中性子束と冷却材温度の観測信号が入力される状態方程式変換ステップと、この状態方程式変換ステップより得られる状態方程式にカルマンフィルタを適用して対数尤度を計算する尤度計算ステップと、この尤度計算ステップで求められた対数尤度を最大化することで温度反応度係数を推定する温度反応度係数推定ステップとを備えて、前記原子炉の中性子束と冷却材温度の観測信号から温度反応度係数を測定する。
【００３２】
請求項２に対応する発明は、請求項１に対応する発明の原子炉の温度反応度係数計測方法において、前記尤度計算ステップで対数尤度を最大化するにあたっては、非線形計画法を適用する。
【００３３】
請求項３に対応する発明は、請求項１に対応する発明の原子炉の温度反応度係数計測方法において、尤度計算ステップは、複数の温度反応度係数を与えて対数尤度を計算し、温度反応度係数推定ステップは、その多項式フィッティングで最大値を見出して温度反応度係数を推定する。
【００３４】
請求項４に対応する発明は、請求項１乃至請求項３の何れかに対応する発明の原子炉の温度反応度係数計測方法において、前記冷却材温度の動特性モデルは、温度信号の回帰モデルで表し、その係数を温度観測信号のフィッティングで求めたものである。
【００３５】
請求項５に対応する発明は、請求項１乃至請求項３の何れかに対応する発明の原子炉の温度反応度係数計測方法において、前記冷却材温度の動特性モデルは、温度信号および中性子束信号の回帰モデルで表し、その係数を観測信号のフィッティングで求めたものである。
【００３６】
請求項６に対応する発明は、請求項１乃至請求項３の何れかに対応する発明の原子炉の温度反応度係数計測方法において、前記冷却材温度の動特性モデルは、設計で用いられているパラメータを用いたモデルである。
【００３７】
【発明の実施の形態】
以下本発明の実施の形態を図面を参照して説明する。
【００３８】
図１は本発明による原子炉の温度反応度係数計測方法を説明するための第１の実施形態を示す機能ブロック図である。
【００３９】
従来の方法では、（３）式による反応度推定を一旦行なった上で、温度反応度係数を推定していたが、本発明の第１の実施形態では、図１に示すように原子炉の中性子束の一点近似動特性モデル１ａと冷却材温度の動特性モデル１ｂで構成される状態方程式変換部１に原子炉の中性子束と冷却材温度の観測信号を直接入力し、その状態方程式をカルマンフィルタ予測計算部２によりカルマンフィルタによる状態を推定した後、尤度関数計算部３にてモデルの尤度関数を計算し、さらに温度反応度係数推定部４によりこの尤度関数を最小化するように温度反応度係数を推定するものである。なお、状態方程式変換部１には設計パラメータが与えられている。
【００４０】
ここで、本発明の第１の実施形態について詳細に説明する。
【００４１】
前述の（１），（２）式で表現したシステム状態空間モデルは、このままの非線形での扱いも可能であるが、数値計算の簡易さと実用性を考え、本実施形態では、線形近似と即発跳躍近似をすることで、下式のような線形の状態方程式に変換して用いる。このとき、状態変数Ｘは、先行核密度Ｃと、冷却材温度Ｔの二つの信号になり、観測変数に中性子束が現れる。
【００４２】
【数５】

【００４３】
次に、文献（赤池、北川偏、「時系列解析の実際ＩＩ」、１９９５年、朝倉書店発行）の手法に沿って、この離散系状態方程式に関してカルマンフィルタを適用すると、下記の逐次式が得られる。ここでは、同時に得られる対数尤度関数も併せて示している。
【００４４】
【数６】

【００４５】
ここで、Ｘ_ｎ _｜ _ｎ−１は、ｎ−１時点までの観測値を用いた状態量Ｘ_ｎの推定値で、Ｖ_ｎ _｜ _ｎ−１は、その分散である。状態量と分散の初期値Ｘ_０１０、Ｖ_０１０を与えることで、逐次、状態量、カルマンゲイン、対数尤度を求めることができる。
【００４６】
このとき、システムモデルの未知パラメータである温度反応度係数Ｋ_Ｔは、対数尤度を最大化する値として、非線形計画法で求めることができる。この非線形計画法は、市販されているアルゴリズム（例えばＭＡＴＡＢ^ＴＭ）を用いることで簡単に実現できる。
【００４７】
システムノイズ分散Ｑと観測ノイズ分散Ｒは、同時に推定することも可能であるが、ここでは反応度外乱を印加することを前提としているため、システムノイズの影響は小さいものとみなすことができると考え、Ｑ〜０とし、観測ノイズについては実信号の観測値を考慮して、外乱による温度変動幅の３０％程度の観測ノイズになるように仮定している。反応度外乱幅は、システムモデルの係数ａにより調整できるが、この係数も既知の値として仮定している。
【００４８】
このような仮定のもとでは、（１４）式の不確定パラメータは、温度反応度係数のみとなるため、（１４）式を繰返して解く数値計算法は短時間で行なうことができる。
【００４９】
図２は、本実施形態での中性子束、温度、さらに中性子束から計算した反応度の観測時間変化を示したものである。ここでは、本発明の性能を確かめるため、（１）式のモデルを用いたシミュレーションによりこの観測時系列信号を求めているが、（１）式のシステムノイズと（２）式の観測ノイズを変えて、５０回同じシミュレーションで観測時系列波形を求め、それぞれに本発明の方法を適用することで、推定性能の検証を行なった。
【００５０】
ここで、推定性能とは、推定値のバイアスのなさと、ばらつきの小ささを意味している。
【００５１】
図３は、この結果を示すグラフである。温度反応度係数の真値０．４５¢／℃に対して、カルマンフィルタを用いた最尤推定法ではバイアスもなく、かつばらつきも小さい推定値が得られることが分かる。なお、最尤法による最適値の推定は、ＭＡＴＬＡＢ^ＴＲによる非線形計画法の関数を利用した。
【００５２】
これに対して、従来法の（５）式を用いた反応度と温度の単純な一次回帰フィッティングによる値は、係数が小さい値にバイアスしているほか、ばらつきも大きい。バイアスの理由は、（７）式で述べたように入力信号である温度に含まれる観測ノイズを正しく処理していないためである。
【００５３】
図４には反応度と温度の観測信号に５点の移動平均をとって平滑化した後に、（５）式の最小２乗法フィッティングをして温度反応度係数を求めた結果を示すグラフであるが、やはり推定値のバイアスがなくなっていることが分かる。ただし、推定値のばらつきについては大きいままである。
【００５４】
図５は、（１１）式に従って、周波数解析による伝達関数の平均ゲインで、コヒーレンス値の大きい低周波領域の平均値から求めたものであるが、推定値のばらつきの大きさは同様である。
【００５５】
これらの結果をまとめると表１のようになる。
【００５６】
【表１】

【００５７】
ここで、表１において、
ＫａｌｍａｎＦｉｌｔｅｒ：本発明のカルマンフィルターによる最尤推定法による温度反応度係数の推定結果
ＳｉｍｐｌｅＦｉｔｔｉｎｇ：（５）式を用いた反応度と温度の単純な一次回帰フィティングによる温度反応度係数の推定結果
ＡｖｅｒａｇｅＦｉｔｔｉｎｇ：５点の移動平均をとって平滑化した後に、（５）式を用いた反応度と温度の一次回帰フィティングをして温度反応度係数を推定した結果
ＴｆＧａｉｎ：（１１）式に従って周波数解析による伝達関数から温度反応度係数を推定した結果
ＩｎｄｉｖｉｄｕａｌＦｉｔｔｉｎｇ：（９）式に従って反応度と温度をそれぞれ時間の関数で一次回帰フィティングをして温度反応度係数を推定した結果
である。
【００５８】
係数推定値のバイアスについては、入力側の観測ノイズの処理を前処理で行なうことで解決できるが、反応度と温度に混入する観測ノイズのために大きくなっていると考えられる推定値のばらつきは、従来法では解決できない。ＩｎｄｉｖｉｄｕａｌＦｉｔｔｉｎｇは、ランプ状外乱に対して適用できる方法で、反応度と温度をそれぞれ時間の一次式でフィッティングして、その傾きの比から係数を求める方法であり、因子分析法に似たやり方となる。これも、バイアスをなくすという観点では有効であるが、ばらつきを小さくするという観点では改善されない。
【００５９】
これに対して、本発明のカルマンフィルタによる最尤推定法は、状態方程式で規定した制約の中で最適パラメータを求めるため、このばらつきが少なくなっている。ばらつきの少なさは、計測時間の短縮を意味しており、前記した表１から１／７倍程度短縮できることになると言える。
【００６０】
このように本発明の第１の実施形態では、状態方程式にカルマンフィルタを用いた状態推定を適用することで、観測ノイズを明確に取扱い、さらにパラメータ設定のための評価関数を対数尤度として求めているため、推定性能を従来方法に比べてより向上させることができる。特に、推定値のばらつきの小ささは、短時間で高精度の推定が可能であり、実用的な価値も大きい。
【００６１】
図６は、本発明による原子炉の温度反応度係数計測方法を説明するための第２の実施形態を示す機能ブロック図で、図１と同一部分には同一符号を付してその説明を省略し、ここでは異なる点について述べる。
【００６２】
本発明の第２の実施形態では、原子炉の冷却材温度の動特性モデル１ｂを、温度信号および中性子束信号の回帰モデルで表し、その係数を温度動特性モデル同定手段５を用いて原子炉の中性子束信号および冷却材温度信号のフィティングで求めるものである。
【００６３】
前述した第１の実施形態で用いた（１）式では、核加熱状態を想定して、中性子束から温度へのフィードバックを無視し、また、ランプ状の人為外乱を想定したモデルとしたが、出力状態で人為外乱のない定常ゆらぎの場合、
ｄＴ／ｄｔ＝ａ_１ｎ−ａ_２Ｔ＋Ｖ_３ ……（１５）
といったモデルで表現することが必要である。ここでは、温度に関して一次遅れのモデルとしたが、ゆらぎの構造に依存して
【数７】

【００６４】
といった１入力１出力の自己回帰（ＡＲＸ）モデルで表現することも可能である。この係数は、温度と中性子束の観測値から（１４）式のカルマンフィルタアルゴリズムとは独立に求めることができるため、（１４）式の温度反応度係数推定アルゴリズムにそのまま組込むことが可能である。
【００６５】
これにより、本発明は中性子束から温度へのフィードバックがある出力状態での温度反応度係数の推定にも適用できる。
【００６６】
また、図７は第１の実施形態と同じ５０ケースのシミュレーションをした際の対数尤度関数値（符号を逆転して表示）を、温度反動時係数の関数として示したものである。この最小値が、最適な温度反応度係数推定値になるが、２次形の関数形状をしているため、数点の温度反応度係数における対数尤度関数値を求めて、これを二次関数でフィッティングしてその最小値を求めることで、温度反応度係数を推定することも可能である。これは第１の実施形態で用いた非線形計画法と比べて、推定精度が落ちる可能性もあるが、推定に際しての計算量が一定になるため、実用上の利点になることがある。
【００６７】
なお、上記実施形態では、冷却材温度の動特性モデルを、原子炉の冷却材の温度信号および中性子束信号の回帰モデルで表したが、冷却材の温度信号の回帰モデルで表しても良い。
【００６８】
【発明の効果】
以上述べたように本発明による原子炉の温度反応度係数測定方法によれば、原子炉の中性子束と冷却材の温度信号の動特性を表現した状態方程式にカルマンフィルタを適用し、さらに最尤法を用いて温度反応度係数を推定することにより、従来法である逆時間方程式により推定する反応度を介することがないので、よりばらつきの少ない推定が可能となる。これは、計測時間の短縮を意味し、計測に要する時間の節約から原子炉の稼働率の向上に役立てることができる。
【図面の簡単な説明】
【図１】本発明による原子炉の温度反応度係数計測方法を説明するための第１の実施形態を示す機能ブロック図。
【図２】同実施形態における計測信号の例で、観測ノイズを伴う中性子束と温度の時間波形、中性子束から計算される反応度の時間波形を示すグラフ。
【図３】同実施形態における温度反応度係数の推定結果を、従来の一次回帰フィッティング法と比較して示すグラフ。
【図４】同実施形態における温度反応度係数の推定結果を、従来の一次回帰フィッティング法と比較して示すグラフ。
【図５】同実施形態における温度反応度係数の推定結果を、従来の周波数解析法による伝達関数により求めた値と比較して示すグラフ。
【図６】本発明による原子炉の温度反応度係数計測方法を説明するための第２の実施形態を示す機能ブロック図。
【図７】同実施形態における対数尤度関数値を温度反応度係数の関数として示すグラフ。
【符号の説明】
１ａ…一点炉近似動特性モデル
１ｂ…温度動特性モデル
１…状態方程式変換部
２…カルマンフィルタ予測計算部
３…尤度関数計算部
４…温度反応度係数推定部
５…温度動特性モデル同定手段

Claims

原子炉の中性子束の一点近似動特性モデルと冷却材温度の動特性モデルで構成され前記原子炉の中性子束と冷却材温度の観測信号が入力される状態方程式変換ステップと、
この状態方程式変換ステップより得られる状態方程式にカルマンフィルタを適用して対数尤度を計算する尤度計算ステップと、
この尤度計算ステップで求められた対数尤度を最大化することで温度反応度係数を推定する温度反応度係数推定ステップとを備えて、
前記原子炉の中性子束と冷却材温度の観測信号から温度反応度係数を測定することを特徴とする原子炉の温度反応度係数計測方法。
請求項１記載の原子炉の温度反応度係数計測方法において、前記尤度計算ステップで対数尤度を最大化するにあたっては、非線形計画法を適用することを特徴とする原子の炉温度反応度係数計測方法。
請求項１記載の原子炉の温度反応度係数計測方法において、尤度計算ステップは、複数の温度反応度係数を与えて対数尤度を計算し、温度反応度係数推定ステップは、その多項式フィッティングで最大値を見出して温度反応度係数を推定することを特徴とする原子炉の温度反応度係数計測方法。
請求項１乃至請求項３の何れかに記載の原子炉の温度反応度係数計測方法において、前記冷却材温度の動特性モデルは、温度信号の回帰モデルで表し、その係数を温度観測信号のフィッティングで求めたものであることを特徴とする原子炉の温度反応度係数計測方法。
請求項１乃至請求項３の何れかに記載の原子炉の温度反応度係数計測方法において、前記冷却材温度の動特性モデルは、温度信号および中性子束信号の回帰モデルで表し、その係数を観測信号のフィッティングで求めたものであることを特徴とする原子炉の温度反応度係数計測方法。
請求項１乃至請求項３の何れかに記載の原子炉の温度反応度係数計測方法において、前記冷却材温度の動特性モデルは、設計で用いられているパラメータを用いたモデルであることを特徴とする原子炉の温度反応度係数計測方法。