JP2010256328A - 発電所計測器性能監視予測方法 - Google Patents

Abstract

【課題】既存のカーネル回帰法に比して予測値計算の正確度を向上させることができる発電所計測器性能監視予測方法を提供する。
【解決手段】計測器信号を行列表示し、正規化し、訓練用、最適化用及びテスト用に分離し、主成分を抽出し、最適化用データを用いてＳＶＲモデルの最適定数を反応分析表面法により求め、最適定数を用いて訓練モデルを生成し、与えられた正規化出力を元の範囲に逆正規化して変数の予測値を求める逆正規化を行い、既存のカーネル回帰法に比べて予測値計算の正確度を向上させる。
【選択図】図２

Description

本発明は、発電所計測器性能監視予測方法に関し、特には、原子力発電所が運転中の状態で、発電所の安全監視用計測器の性能を常時監視するために必要な技術に関する。

一般的に、全ての発電設備には、運転性向上と安全性確保を目的として多数の計測器を設け、前記計測器からリアルタイムに信号を取得して、発電所監視系統と保護系統に用いている。特に、原子力発電所の安全系統に関連する計測チャネルには、計測信号の正確度と信頼度を保障するために、多重計測器概念を採用しており、運営技術指針書上に記載のように核周期（約１８ヶ月）ごとに点検と矯正を実行している。世界中の原子力発電所は、無駄に行われる計測器矯正業務を改善するために、条件に基づく監視（ＣＢＭ、Ｃｏｎｄｉｔｉｏｎ Ｂａｓｅｄ Ｍｏｎｉｔｏｒｉｎｇ）方法論を開発して点検と矯正周期を長くする技術を開発している。

図１は、一般的な計測器性能常時監視システムのブロック構成図であって、Ａｕｔｏ−Ａｓｓｏｃｉａｔｉｖｅモデルと呼ばれるシステムである。図１を参照すれば、一般的な計測器性能常時監視システムは、予測モデル（ｐｒｅｄｉｃｔｉｏｎ ｍｏｄｅｌ）、比較モジュール（ｃｏｍｐａｒｉｓｏｎ ｍｏｄｕｌｅ）、判断ロジック（ｄｅｃｉｓｉｏｎ ｌｏｇｉｃ）で構成されているが、計測信号を予測モデルに入力すると入力測定値に対するモデルの予測値を出力し、比較モジュールにより測定値と予測値との差分を判断ロジックに入力し連続的に監視すれば、計測器のドリフトと故障を感知することができる。

一方、Ａｒｇｏｎｎｅ Ｎａｔｉｏｎａｌ ＬａｂｏｒａｔｏｒｙはＭＳＥＴ（Ｍｕｌｔｉｖａｒｉａｔｅ Ｓｔａｔｅ Ｅｓｔｉｍａｔｉｏｎ Ｔｅｃｈｎｉｑｕｅ）を開発して米国特許を取得し、ＳｍａｒｔＳｉｇｎａｌ Ｃｏｒｐｏｒａｔｉｏｎ及びＥｘｐｅｒｔ Ｍｉｃｒｏｓｙｓｔｅｍｓ社は同特許を商業的に使用できるよう製品化した。Ｅｘｐｅｒｔ Ｍｉｃｒｏｓｙｓｔｅｍｓ社は米国内Ｐａｌｏ Ｖｅｒｄｅ、Ｌｉｍｅｒｉｃｋ １／２、ＴＭＩ、Ｖ．Ｃ．Ｓｕｍｍｅｒ、Ｓｅｑｕｏｙａｈ１、Ｓａｌｅｍ １号機でＭＳＥＴを用いる製品を現場に設置して、計測チャネルオンライン監視を行っている。また、ＳｍａｒｔＳｉｇｎａｌ Ｃｏｒｐｏｒａｔｉｏｎはその後ＭＳＥＴに対する特許を使用できなくなり、カーネル回帰法に基づく計測器性能常時監視技術を開発した。

計測器の予測値を計算するための方法としては、線形回帰分析法（Ｌｉｎｅａｒ Ｒｅｇｒｅｓｓｉｏｎ Ｍｅｔｈｏｄ）が一般的に使用される。この方法は、予測しようとする計測器信号と線形的な相関関係の高い他の計測器信号を選択し、予測値と測定値との誤差二乗和が最小になるよう回帰係数を求める方法であって、次の数１により表現できる。

前記の線形回帰分析法は、既知の従属変数と独立変数で回帰係数が決まると、未知の従属変数に対する独立変数を予測できる。

しかしながら、既存の線形回帰分析方法論は、従属変数同士で互いに線形的連関性の大きい場合、多重共線性の問題が生じ、従属変数に含まれる小さいノイズに対して独立変数に大きい誤差が発生する。

カーネル回帰法は、既存の線形回帰分析方法論や神経回路網のように入力と出力との相関関係を最適化する回帰係数や加重値のような媒介変数（Ｐａｒａｍｅｔｅｒ）を使用せず、選別された測定データをメモリベクトルとして格納し、測定信号セットに対するメモリベクトル内の訓練データセットのユークリッド距離からカーネルの加重値を求め、これをメモリベクトルに適用して計測器の予測値を求める非媒介変数回帰法（Ｎｏｎｐａｒａｍｅｔｒｉｃ ｒｅｇｒｅｓｓｉｏｎ ｍｅｔｈｏｄ）である。カーネル回帰法のような非媒介変数回帰法は、入出力関係が非線形状態のモデルと信号雑音に強いという長所を有する。

既存のカーネル回帰法の計算手順は、下記の５段階の過程からなる。

１段階：訓練データを次の数２の行列の形に表す。

ここで、Ｘはメモリベクトルに格納される訓練データの行列であって、ｎは訓練データ個数であり、ｍは計測器の番号である。

２段階：次の数３を用いて第１の計測器信号セットに対する訓練データのユークリッド距離の和を求める。

ここで、ｘは訓練データ、ｑはテストデータまたは問い合わせデータ（Ｑｕｅｒｙ ｄａｔａ）、ｔｒｎは訓練データの番号、ｊは計測器の番号である。

３段階：次のカーネル関数を含める数４を用いて各々の訓練データセットと与えられたテストデータセットに対する加重値を求める。

前記の数４においては、加重関数として用いられたガウスカーネルは次のように定義される。

４段階：テストデータの予測値は、次の数６のように各々の訓練データに加重値を掛けた後、加重値の和で割って求める。

５段階：全体のテストデータに対する予測値を求めるために、２段階から４段階の過程を繰り返す。
このような従来のカーネル回帰法は、非線形状態のモデルと信号雑音に強い長所を有している。

しかし、前記の従来のカーネル回帰法は、選別された測定データをメモリベクトルとして格納し、測定信号セットに対するメモリベクトル内の訓練データセットのユークリッド距離からカーネルの加重値を求め、これをメモリベクトルに適用して計測器の予測値を求めるので、出力予測値の分散が大きくなり、これにより線形回帰分析法に比べて正確度の低いという問題点がある。

本発明の目的は、上記のような従来のカーネル回帰法の予測正確度の低い問題点を解消するために提案されたものであり、発電所データの正規化、主成分抽出、反応表面分析法を用いるＳＶＲ（Ｓｕｐｐｏｒｔ Ｖｅｃｔｏｒ Ｒｅｇｒｅｓｓｉｏｎ）モデル回帰式のパラメータ（カーネル帯域幅σ、損失関数ε−、ペナルティＣ）最適化、これを用いるＳＶＲモデル具現、出力予測値の逆正規化方法等を用いて発電所システムをモデリングした後、計測器信号を監視することにより、既存のカーネル回帰法に比して予測値計算の正確度を向上させることができる発電所計測器性能監視予測方法を提供することにある。

上述の目的を達成するために、本発明の発電所計測器性能監視予測方法は、発電所データの主成分を抽出しＳＶＲ（Ｓｕｐｐｏｒｔ Ｖｅｃｔｏｒ Ｒｅｇｒｅｓｓｉｏｎ）方法を用いてシステムを各種ケースに対してモデリングした後、回帰式の三つのパラメータを反応表面分析法を用いて最適化し、これを用いて発電所システムをモデリングした後、計測器信号を監視することにより、既存の広く使用されるカーネル回帰法に比して予測値計算の正確度を向上させることができる。

本発明の発電所計測器性能監視用予測方法は、原子力発電所の安全監視チャネルに用いられる計測器の性能を発電所が運転中の状態でオンラインで監視できるようにすることで、計測器の誤作動をリアルタイムに監視して計測器信頼度を向上させ、原子力発電所の計測器矯正周期を現在燃料交替周期の１８ヶ月から最大８年に延長することにより、矯正費用と放射線区域における矯正作業従事者の放射線被爆を低減し、無駄な矯正回数を減らすことで、誤矯正による発電所不時停止を予防し、発電所予防停止期間を短縮して発電所利用率を増加させることができる効果を有する。

従来の発電所計測器の性能常時監視システムのブロック構成図である。 本発明の一実施例にかかる発電所計測器性能監視予測方法を説明するための手順を示すフローチャートである。 図２の発電所計測器性能監視予測方法に従って生成されたＳＶＲモデルである。 ＳＶＲによる最適回帰式の一般的な概念図である。 モデルパラメータが三つである場合の中心合成計画における実験点を示す。 正確度テストのための原子力発電所原子炉炉心出力データである。 正確度テストのための原子力発電所原子炉炉心出力データである。 正確度テストのための原子力発電所加圧器水位データである。 正確度テストのための原子力発電所加圧器水位データである。 正確度テストのための原子力発電所蒸気発生器蒸気流量データである。 正確度テストのための原子力発電所蒸気発生器蒸気流量データである。 正確度テストのための原子力発電所蒸気発生器狭域水位データである。 正確度テストのための原子力発電所蒸気発生器狭域水位データである。 正確度テストのための原子力発電所蒸気発生器圧力データである。 正確度テストのための原子力発電所蒸気発生器圧力データである。 正確度テストのための原子力発電所蒸気発生器広域水位データである。 正確度テストのための原子力発電所蒸気発生器広域水位データである。 正確度テストのための原子力発電所蒸気発生器主給水流量データである。 正確度テストのための原子力発電所蒸気発生器主給水流量データである。 正確度テストのための原子力発電所タービン出力データである。 正確度テストのための原子力発電所タービン出力データである。 正確度テストのための原子力発電所１次側充電流量データである。 正確度テストのための原子力発電所１次側充電流量データである。 正確度テストのための原子力発電所残熱除去流量データである。 正確度テストのための原子力発電所残熱除去流量データである。 正確度テストのための原子力発電所原子炉上部冷却材温度データである。 正確度テストのための原子力発電所原子炉上部冷却材温度データである。

以下、本発明が属する技術分野で通常の知識を有する者が本発明を容易に実施できる程度に詳細の説明するために、本発明の最も好ましい実施例を添付の図面を参照して詳しく説明する。本発明の目的、作用、効果を含めてその他の目的、特徴、そして動作上の利点が好ましい実施例の説明によってさらに明確になるであろう。

図２は、本発明の一実施例にかかる発電所計測器性能監視予測方法を説明するための手順を示すフローチャートである。図２に示すように、本発明の一実施例にかかる発電所計測器性能監視予測方法は、全体データを行列の形に表示する行列表示ステップＳ１００と、前記全体データをデータセットとして正規化する正規化ステップＳ２００と、データセットを訓練用、最適化用、テスト用に分離する分離ステップＳ３００と、正規化された各データセットの主成分を抽出する主成分抽出ステップＳ４００と、反応表面分析法を用いて最適化用データの予測値誤差を最適化するＳＶＲモデルの最適定数を求める最適定数算出ステップＳ５００と、最適定数を用いてＳＶＲモデルを生成する訓練モデル生成ステップＳ６００と、正規化されたテストデータを入力としてカーネル関数を求めＳＶＲモデルの出力値を予測する訓練モデル出力予測ステップＳ７００と、与えられた正規化出力を元の範囲に逆正規化して変数の予測値を求める逆正規化ステップＳ８００とを含む。

前記の構成による本発明の一実施例にかかる発電所計測器性能監視予測方法の作用は次の通りである。

まず、行列表示ステップＳ１００では、全体データセット（Ｘ）を次の数７のように行列の形に表す。

ここで、Ｘは全体データセットであり、Ｘｔｒは訓練用データセットであり、Ｘｏｐは最適化用データセットであり、Ｘｔｓはテスト用データセットであり、３ｎは全体データの個数であり、ｍは計測器の番号である。

次に、正規化ステップＳ２００では，次の数８を用いて全体データを正規化する。

ここで、ｉ＝１、２、．．．、３ｎである。

正規化された全体データセット（ｓｅｔ）Ｚは次の数９のように表すことができる。

前記分離ステップＳ３００では、データセットＺを訓練（Ｔｒａｉｎｉｎｇ）、最適化（Ｏｐｔｉｍｉｚａｔｉｏｎ）、テスト（Ｔｅｓｔ）用の三つに分離する。本発明の実施例において、三つに分離されたデータは各々Ｚｔｒ、Ｚｏｐ、Ｚｔｓとし、次の数１０のようにｎ×ｍの大きさを有することができる。

ここで、ｉ＝０、１、２、．．．、ｎ−１である。

前記主成分抽出ステップＳ４００では、正規化された各データセットＺｔｒ、Ｚｏｐ、Ｚｔｓの主成分を抽出する。主成分の分散（すなわち、共分散マトリックスのＥｉｇｅｎｖａｌｕｅ）を大きさ順に並べ、百分率分散値が最も大きい主成分から、その累積和が使用者が定めた９９．５％以上になるまでのＺｔｒ、Ｚｏｐ、Ｚｔｓに対する主成分（Ｐｔｒ、Ｐｏｐ、Ｐｔｓ）を選択する。

-主成分を求める方法
主成分分析（ＰＣＡ）は、線形変換による多数の入力変数の少数の変数への圧縮に有効な方法である。このとき、圧縮された変数を主成分（Ｐｒｉｎｃｉｐａｌ Ｃｏｍｐｏｎｅｎｔ）と呼ぶ。ＰＣＡは変数間相関関係を用いて元の次元のデータを二乗和（Ｓｕｍ ｏｆ Ｓｑｕａｒｅｓ）が最大化する低次元の超平面に放射する。
以下では、主成分を抽出する過程を説明する。
ｍ次元の入力変数をｘ_１、ｘ_２、．．．、ｘ_ｍとし、これらの線形結合で生成される新しい変数をθ_１、θ_２、．．．、θ_ｍとすれば、これらの関係は次の数１１及び数１２のように表現できる。

このとき、θ_１、θ_２、．．．、θ_ｍを線形システムの主成分と言い、便宜上第１の主成分を最も重要な主成分と決める。最も重要な主成分とは、入力変数の変化を最も多く説明する、すなわち、分散の最も大きい主成分を意味する。線形変換Ｑは、次の数１３及び数１４の条件を満たすように決定する。

前記数１３の条件は、変換後のスケールをそのまま維持するようにし、前記数１４の条件は、変換後には変数間に互いに相関関係が無いようにする。

主成分は次のような手順を通じて求められる。
イ．各データセットＺｔｒ、Ｚｏｐ、Ｚｔｓから各変数の平均値を引き、これをＡマトリックスとする。本発明の実施例では、Ｚｔｒを例示して説明し、これを数式で表すと次の数１５のようになる。

ロ．次の数１６乃至数１８を用いて、Ａ^ＴＡの根（ｅｉｇｅｎｖａｌｕｅ）λを求め、Ａのｓｉｎｇｕｌａｒ ｖａｌｕｅ Ｓを求める。数１６から求めた０を除いた根（ｅｉｇｅｎｖａｌｕｅ）λを降順でまとめ、これをλ_１、λ_２、．．．、λ_ｍとする。

ハ．ｎ×ｎマトリックスであるＡＡ^Ｔのｅｉｇｅｎｖｅｃｔｏｒを求め、ｕｎｉｔａｒｙ ｍａｔｒｉｘ Ｕを求める。次の数１９からｅｉｇｅｎｖａｌｕｅ λを求め、これを数２０に代入して各ｅｉｇｅｎｖａｌｕｅ λに対するｎ×１であるｅｉｇｅｎｖｅｃｔｏｒ ｅ_１、ｅ_２、．．．、ｅ_ｍを求める。

ニ．次の数２１を用いて各主成分の分散を求める。

ホ．次の数２２及び数２３を用いて各主成分の分散を全体主成分の分散の和で割って百分率を求める。

ヘ．百分率分散（％σ_p）が最も大きいものから累積計算して、所望の百分率分散（例えば９９．９８％）までの主成分ｐ個を選択する。

ト．主成分を次の数２４で計算して抽出する。

チ．Ｚｏｐ、Ｚｔｓに対しても上記のような手順により主成分を抽出することができる。

次の表１は、抽出された主成分の分散を示す。

本発明の実施例では７個の主成分を使用するが、前記表１を参照すれば、７個の主成分を使用する場合、全体分散の９９．９８％を説明でき、残りの主成分を放棄することで発生する情報の損失は０．０２％に過ぎないことが分かる。

‐ＳＶＲ（Ｓｕｐｐｏｒｔ Ｖｅｃｔｏｒ Ｒｅｇｒｅｓｓｉｏｎ）モデリング
ｍ次元の入力変数ｘ_１、ｘ_２、．．． 、ｘ_ｍをｐ次元の主成分θ_１、θ_２、．．．、θ_ｐで圧縮するのは、次の数２５のように表すことができる。

ここで、ｐはｍ以下の整数である。

第ｋの出力に対するＳＶＲで求める最適回帰式（ＯＲＬ：Ｏｐｔｉｍｕｍ Ｒｅｇｒｅｓｓｉｏｎ Ｌｉｎｅ）は次の数２６のように表すことができる。

ここで、ｋは１、２、．．．、ｍである。

第ｋの出力変数ｙ^（ｋ）に対してε‐ｉｎｓｅｎｓｉｔｉｖｅ Ｌｏｓｓ Ｆｕｎｃｔｉｏｎを次の数２７のように定義すれば、ｙ^（ｋ）に対する最適回帰式（ＯＲＬ）を求めるための最適化問題は次の数２８のように表すことができる。

前記数２７及び数２８において、ｋは１、２、．．．、ｍである。そして、ξ_ｋ，ｉとξ^＊ _ｋ，ｉは、図４に示したスラック変数（Ｓｌａｃｋ Ｖａｒｉａｂｌｅ）を意味する。図４は、ＳＶＲによる最適回帰線の概念図である。ここで、θ_ｉは、ベクトルθの第ｉの要素でなく、ｘに対する第ｉの観測値ベクトルに対応する主成分ベクトルを意味する。

前記の最適化問題は次の数２９のような双対問題に変えて表すことができる。

ここで、ｋは１、２、．．．、ｍである。

これにより、ラグランジュ乗数λ_ｋ，ｉとλ^＊ _ｋ，ｉを次の数３０に代入して、ＡＡＳＶＲの第ｋの出力変数に対する最適回帰式（ＯＲＬ）を決定することができる。

上述した内容は最適の線形回帰式を求める過程である。元のデータを高次元空間に非線形的に写像した結果をベクトルФ（・）とすれば、非線形写像の結果であるФ（・）の内積として次の数３１のように定義される関数をカーネル（Ｋｅｒｎｅｌ）と呼ぶ。

高次元空間で最適回帰式（ＯＲＬ）を探そうとするとき、個別のФ（ｘ_ｉ）とФ（ｘ_ｊ）を全て知るべき理由が無く、ただカーネル関数Ｋ（ｘ_ｉ、ｘ_ｊ）さえ知っていれば十分である。本発明の実施例においては、Ｇａｕｓｓｉａｎ Ｒａｄｉａｌ Ｂａｓｉｓ Ｆｕｎｃｔｉｏｎを使用した。

カーネル関数Ｋ（θ_ｉ、θ）＝Ф（θ_ｉ）^ＴФ（θ）を用いれば、次の数３２のような最適の非線形回帰式を求めることができる。

ここで、バイアス項は、次の数３３のように任意のＳｕｐｐｏｒｔ Ｖｅｃｔｏｒであるθ_ｒとθ_ｓを用いて計算されることができる。

以上から得られた最適の非線形回帰式はｌｏｓｓ ｆｕｎｃｔｉｏｎの定数ε（ｅｐｓｉｌｏｎ）、双対目的関数Ф（ｗ_ｋ、ξ_ｋ）にペナルティＣ、ＲＢＦ（Ｒａｄｉａｌ Ｂａｓｉｓ Ｆｕｎｃｔｉｏｎ）をカーネルとして用いる場合、カーネル帯域幅σが予め与えられるのが好ましい。この過程を繰り返して各々の出力に対する総ｍ個のＳＶＲを得て、図３に示したＡＡＳＶＲを構築することができる。

前記最適定数算出ステップＳ５００では、反応表面分析法を用いて最適化データ（Ｚｏｐ）の予測値誤差を最小化するＳＶＲモデルの最適定数（σ、ε、Ｃ）を求める。

イ．ＳＶＲモデルパラメータであるｓｉｇｍａ（σ）、ｅｐｓｉｌｏｎ（ε）、Ｃを各々ｖ_１、ｖ_２、ｖ_３とする。
ロ．ｖ_１、ｖ_２、ｖ_３に対する探索範囲を各々決める。適切な探索範囲は事前経験や小規模の予備実験により把握できる。本発明の実施例においては、各々ｖ_１は０．５６４８６−１．６３５１４、ｖ_２は０．０１０５３４−０．０３９９６６、ｖ_３は２．１０７６−７．９９３３とする。

ハ．探索範囲の上限を各々Ｕ_１、Ｕ_２、Ｕ_３とし、下限をＬ_１、Ｌ_２、Ｌ_３とし、次の数３４のようにモデルパラメータを標準化する。

ニ．標準化されたモデルパラメータｘ_１、ｘ_２、ｘ_３の探索範囲を考慮して実験点、すなわちモデル性能の評価地点を定める。このために、統計的な実験計画の一つの中心合成計画を用いることができる。中心合成計画により定められる実験点を３次元空間で表現すれば、図５のようになる。

ホ．中心合成計画による実験点は８個の頂点、１個の中心点、そして６個の軸点で構成される。実験誤差の大きさを推定するためには、５回内外繰り返して実験を施す。軸点の座標は予測分散に関する統計的な性質を考慮してα＝２^３／４＝１．６８１７９に定める。

α＝［要因実験点の数］^１／４

中心点で５回繰り返した場合に、中心合成計画によるｘ_１、ｘ_２、ｘ_３の実験点は次の表２のようになる。

ヘ．前記表２に指示するようにモデルパラメータｖ_１、ｖ_２、ｖ_３の値を決めて表３（中心合成計画による（ｖ_１、ｖ_２、ｖ_３）の実験点）を得られるが、これがモデルを得るために直接使用されるモデルパラメータの値となる。

ト．前記表３の各実験点でＺｔｒとＰｔｒを用いてＳＶＲモデルのｂｅｔａベクトルとｂｉａｓ定数を各々得る。そのためにはデータセットＺｔｒを用いる。実際、中心点に該当するＮｏ．１５からＮｏ．１９までは同じモデルが得られる。

チ．各モデルの正確度を評価するためにデータセットＰｏｐをｍ個のＡＡＳＶＲに入力して最適化データの正規化された予測値数３５を求める。これから次の数３６を用いて出力モデルの正確度、すなわちＭＳＥを計算する。５個の中心点は同じモデルであるので、Ｐｏｐを５個に分割し、分割された下部データセットから別のＭＳＥ計算結果を得る。

ここで、ｚ_ｉｊはＰｏｐのうちセンサーｉの第ｊの入力データを表し、Ｚ_ｉｊはモデルによる推定値である。モデルはＰｔｒを用いて得られたことを想起する。本発明の実施例において、実験結果である実験によるＭＳＥ計算結果は、次の表４に表すようになる。

リ．反応表面を求めるときには、ＭＳＥの代わりにｌｏｇ（ＭＳＥ）を使用する。このことを考慮して、モデルパラメータｘ_１、ｘ_２、ｘ_３とｌｏｇ（ＭＳＥ）との間の反応表面を推定する。反応表面は次の数３７のような２次模型を仮定する。

ここで、ｅはランダム誤差を意味する。

本発明の実施例で推定された反応表面は下記のようになる。
ｌｏｇ（ＭＳＥ）＝−８．３４９２，−０．２１３１ｘ_１，＋０．７７１６ｘ_２，−０．０９５２ｘ_３，＋０．２０１０ｘ_１ ^２，−０．０７５３ｘ_２ ^２＋０．０７９９ｘ_３ ^２

ヌ．推定された反応表面式を用いてｌｏｇ（ＭＳＥ）を最小化するｘ_１、ｘ_２、ｘ_３の最適条件を求める。２次反応表面を仮定したので、最適条件は偏微分を通して確認する。すなわち、次の数３８を同時に満たすｘ_１、ｘ_２、ｘ_３を求める。

本発明の実施例で得られた反応表面の場合、最適条件は（ｘ^＊ _１、ｘ^＊ _２、ｘ^＊ _３）＝（１．５４３８、−０．６８１８、１．５９２９）である。
ル．最適条件（ｘ^＊ _１、ｘ^＊ _２、ｘ^＊ _３）を次の数３９を用いて元の単位に換算する。

本発明の実施例において、ｖ^＊ _１＝１．３９１０＝σ^＊、ｖ^＊ _２＝０．０００５＝ε^＊、ｖ^＊ _３＝６．７９５１＝Ｃ^＊になり、この条件で予測されたｌｏｇ（ＭＳＥ）は−９．９４４６なので、ここに指数を利用してＭＳＥに換算すると０．００００４８になる。

前記訓練モデル生成ステップＳ６００では、ステップＳ５００で求めたＳＶＲモデルの三つの最適定数（ε^＊、Ｃ^＊、σ^＊）と訓練データの主成分（Ｐｔｒ）、第１の信号の訓練データ（Ｚｔｒの第１列）を入力としてｑｕａｄｒａｔｉｃ ｐｒｏｇｒａｍを用いて最適化問題を解き、ラグランジュ乗数の差分であるβ_１（ｎ×１）とバイアス定数ｂ_１を求めて、図３のＳＶＲ_１のモデルを生成する。同じ方法で第２から第ｍのセンサーに対してこれを繰り返して行い、β_２、β_３、．．．、β_ｍとｂ_２、ｂ_３、．．．、ｂ_ｍを求めＳＶＲ_２からＳＶＲ_ｍのモデルを生成して、図３のように全センサーに対するＳＶＲモデルを構築する。

前記訓練モデルの出力予測ステップＳ７００では、訓練データの主成分（Ｐｔｒ）、テストデータの主成分（Ｐｔｓ）を用いてＧａｕｓｓｉａｎ Ｒａｄｉａｌ Ｂａｓｉｓ Ｆｕｎｃｔｉｏｎのカーネル関数（Ｋ（ｎ×ｎ））を求め、ステップＳ６００で求めたＳＶＲモデルのラグランジュ乗数の差分β_１、バイアス定数ｂ_１を用いてＳＶＲ_１の出力を求める。同じ方法で第２から第ｍのセンサーに対してこれを繰り返して行い、ＳＶＲ_２からＳＶＲ_ｍの出力のモデル予測値を求めることができる。

これを数式で表すと、次の数４０のようになる。

前記逆正規化ステップＳ８００では、ステップＳ７００で得られた正規化されたテストデータの予測値を元の範囲に逆正規化して元のスケールの各センサーに対する予測値を求める。これを数式で表すと、次の数４１のようになる。

本発明の一実施例にかかる発電所計測器性能監視予測方法の優秀性を確認するために、実際原子力発電所の出力を０％から１００％に上昇させる途中１、２次系統で測定した計測器信号データを用いて既存の方法論と比較した。分析に使用されたデータは総１１個のセンサーにて測定された値である。
次の表５は、従来のカーネル回帰法と本発明の一実施例にかかる発電所計測器監視方法による計測器予測値の正確度を比較した表である。

前記表５で分析に使用されたデータは、次のような総１１個のセンサーにて測定された値である。
１．原子炉出力（％）
２．加圧器水位（％）
３．蒸気発生器蒸気流量（Ｍｋｇ／ｈｒ）
４．蒸気発生器狭域水位データ（％）
５．蒸気発生器圧力データ（Ｋｇ／ｃｍ^２）
６．蒸気発生器広域水位データ（％）
７．蒸気発生器主給水流量データ（Ｍｋｇ／ｈｒ）
８．タービン出力データ（ＭＷｅ）
９．原子炉冷却材充電流量データ（ｍ^３／ｈｒ）
１０．残熱除去流量データ（ｍ^３／ｈｒ）
１１．原子炉上部冷却材温度データ（℃）

上記のデータのモデル入力及び予測値を時間の関数にしてグラフで描くと、図６乃至図１６のようになる。

図６Ａ及び図６Ｂは、原子力発電所原子炉炉心出力データに関し、図６Ａは、数７のテスト用入力データＸ_ｔｓ１に該当し、図６Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ１に対する推定データ数４２を示す。

図７Ａ及び図７Ｂは、原子力発電所加圧器水位データに関し、図７Ａは、数７のテスト用入力データＸ_ｔｓ２に該当し、図７Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ２に対する推定データ数４３を示す。

図８Ａ及び図８Ｂは、原子力発電所蒸気発生器蒸気流量データに関し、図８Ａは、数７のテスト用入力データＸ_ｔｓ３に該当し、図８Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ３に対する推定データ数４４を示す。

図９Ａ及び図９Ｂは、原子力発電所蒸気発生器狭域水位データに関し、図９Ａは、数７のテスト用入力データＸ_ｔｓ４に該当し、図９Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ４に対する推定データ数４５を示す。

図１０Ａ及び図１０Ｂは、原子力発電所蒸気発生器蒸気圧力データに関し、図１０Ａは、数７のテスト用入力データＸ_ｔｓ５に該当し、図１０Ｂは，本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ５に対する推定データ数４６を示す。

図１１Ａ及び図１１Ｂは、原子力発電所蒸気発生器広域水位データに関し、図１１Ａは、数７のテスト用入力データＸ_ｔｓ６に該当し、図１１Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ６に対する推定データ数４７を示す。

図１２Ａ及び図１２Ｂは、原子力発電所蒸気発生器主給水流量データに関し、図１２Ａは、数７のテスト用入力データＸ_ｔｓ７に該当し、図１２Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ７に対する推定データ数４８を示す。

図１３Ａ及び図１３Ｂは、原子力発電所タービン出力データに関し、図１３Ａは、数７のテスト用入力データＸ_ｔｓ８に該当し、図１３Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ８に対する推定データ数４９を示す。

図１４Ａ及び図１４Ｂは、原子力発電所１次側充電流量データに関し、図１４Ａは、数７のテスト用入力データＸ_ｔｓ９に該当し、図１４Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ９に対する推定データ数５０を示す。

図１５Ａ及び図１５Ｂは、原子力発電所残熱除去流量データに関し、図１５Ａは、数７のテスト用入力データＸ_ｔｓ１０に該当し、図１５Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ１０に対する推定データ数５１を示す。

図１６Ａ及び図１６Ｂは、原子力発電所原子炉上部冷却材温度データに関し、図１６Ａは、数７のテスト用入力データＸ_ｔｓ１１に該当し、図１６Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘ_ｔｓ１１に対する推定データ数５２を示す。

表５において、正確度は予測モデルを運転監視に適用するにおいて最も基本的な尺度となる。殆ど正確度はモデル予測値と実際測定値との平均自乗誤差で表す。

次の数５３は、一つの計測器に対する正確度を表す数式である。

ここで、Ｎはテストデータの数、数５４は第ｉのテストデータに対するモデルの推定値、ｘ_ｉは第ｉのテストデータの測定値を表す。

本発明の一実施例にかかる発電所計測器性能監視予測方法は、計測器信号の主成分を抽出し、最適化用データを用いてＳＶＲモデルの最適定数を反応分析表面法により求め、またモデルを訓練データを用いて訓練させて、既存のカーネル回帰法に比して予測値計算の正確度を向上させる。

上述した本発明の詳細な説明では本発明の好ましい実施例を参照して説明したが、本発明の属する技術分野における通常の知識を有する者であれば、本発明の技術的思想及び技術領域を逸脱しない範囲内で、本発明の様々な修正及び変更が可能であることが理解できるはずである。したがって、本発明の技術的範囲は、明細書の詳細な説明に記載の内容に限定されるのではない。

Claims (19)

1. 全体データを行列の形に表示する行列表示ステップ；
   前記全体データをデータセットとして正規化する正規化ステップ；
   前記正規化されたデータセットを訓練用、最適化用及びテスト用に分離する分離ステップ；
   正規化された各データセットの主成分を抽出する主成分抽出ステップ；
   反応表面分析法を用いて最適化用データの予測値誤差を最適化するＳＶＲモデルの最適定数を求める最適定数算出ステップ；
   最適定数を用いてＳＶＲモデルを生成する訓練モデル生成ステップ；
   正規化されたテストデータを入力としてカーネル関数を求めＳＶＲモデルの出力値を予測する訓練モデル出力予測ステップ；及び
   与えられた正規化出力を元の範囲に逆正規化して変数の予測値を求める逆正規化ステップを含むことを特徴とする発電所計測器性能監視予測方法。
2. 前記行列表示ステップでは、全体データを次の数式のように行列の形に表示することを特徴とする請求項１に記載の発電所計測器性能監視予測方法。
   

   

   ここで、Ｘは全体データセットであり、Ｘｔｒは訓練用データセットであり、Ｘｏｐは最適化用データセットであり、Ｘｔｓはテスト用データセットであり、３ｎは全体データの個数であり、ｍは計測器の番号である。
3. 前記正規化ステップでは、次の数式を用いて全体データを正規化することを特徴とする請求項１に記載の発電所計測器性能監視予測方法。
   

   

   ここで、ｉ＝１、２、．．．、３ｎである。
4. 前記分離ステップでは、データセットＺを訓練用、最適化用及びテスト用の三つに分離し、前記分離されたデータは各々Ｚｔｒ、Ｚｏｐ、Ｚｔｓとし、次の数式のようにｎ×ｍの大きさを有することを特徴とする請求項１に記載の発電所計測器性能監視予測方法。
   

   

   ここで、ｉ＝０、１、２、．．．、ｎ−１である。
5. 前記主成分抽出ステップでは、主成分の分散（すなわち、共分散マトリックスのＥｉｇｅｎｖａｌｕｅ）を大きさ順に並べ、百分率分散値が最も大きい主成分から、その累積和が９９．５％以上になるまでのＺｔｒ、Ｚｏｐ、Ｚｔｓに対する主成分（Ｐｔｒ、Ｐｏｐ、Ｐｔｓ）を選択することで正規化された各データセットＺｔｒ、Ｚｏｐ、Ｚｔｓの主成分を抽出することを特徴とする請求項１に記載の発電所計測器性能監視予測方法。
6. 各データセットＺｔｒ、Ｚｏｐ、Ｚｔｓから各変数の平均値を引き、これをＡマトリックスとし、次の数式で表すことを特徴とする請求項５に記載の発電所計測器性能監視予測方法。
   

   
7. 次の数式を用いて、Ａ^ＴＡの根（ｅｉｇｅｎｖａｌｕｅ）λを求め、Ａのｓｉｎｇｕｌａｒ ｖａｌｕｅ Ｓを求め、０を除いた根（ｅｉｇｅｎｖａｌｕｅ）λを降順でまとめ、これをλ_１、λ_２、．．．、λ_ｍと設定することを特徴とする請求項６に記載の発電所計測器性能監視予測方法。
   

   
8. ｎ×ｎマトリックスであるＡＡ^Ｔのｅｉｇｅｎｖｅｃｔｏｒを求め、ｕｎｉｔａｒｙ ｍａｔｒｉｘ Ｕを求め、次の数式からｅｉｇｅｎｖａｌｕｅ λを求め、これを代入して各ｅｉｇｅｎｖａｌｕｅ λに対するｎ×１であるｅｉｇｅｎｖｅｃｔｏｒ ｅ_１、ｅ_２、．．．、ｅ_ｍを求めることを特徴とする請求項７に記載の発電所計測器性能監視予測方法。
   

   
9. 次の数式を用いて各主成分の分散を求めることを特徴とする請求項８に記載の発電所計測器性能監視予測方法。
   

   
10. 次の数式を用いて各主成分の分散を全体主成分の分散の和で割って百分率を求めることを特徴とする請求項９に記載の発電所計測器性能監視予測方法。
    

    
11. 百分率分散（％σ_p）が最も大きいものから累積計算して、所望の百分率分散（例えば９９．９８％）までの主成分ｐ個を選択することを特徴とする請求項１０に記載の発電所計測器性能監視予測方法。
12. 主成分を次の数式で計算して抽出することを特徴とする請求項１１に記載の発電所計測器性能監視予測方法。
    

    
13. 次の数式でｍ次元の入力変数ｘ_１、ｘ_２、．．． 、ｘ_ｍをｐ次元の主成分θ_１、θ_２、．．．、θ_ｐで圧縮することを特徴とする請求項１２に記載の発電所計測器性能監視予測方法。
    

    

    ここで、ｐはｍ以下の整数である。
14. 第ｋの出力に対するＳＶＲ（Ｓｕｐｐｏｒｔ Ｖｅｃｔｏｒ Ｒｅｇｒｅｓｓｉｏｎ）で求める最適回帰式（ＯＲＬ：Ｏｐｔｉｍｕｍ Ｒｅｇｒｅｓｓｉｏｎ Ｌｉｎｅ）を次の数式で表すことを特徴とする請求項１３に記載の発電所計測器性能監視予測方法。
    

    

    ここで、ｋは１、２、．．．、ｍである。
15. 第ｋの出力変数ｙ^（ｋ）に対してε‐ｉｎｓｅｎｓｉｔｉｖｅ Ｌｏｓｓ Ｆｕｎｃｔｉｏｎを次の数式のように定義すれば、ｙ^（ｋ）に対するＯＲＬを求めるための最適化問題はその次の数式のように表すことを特徴とする請求項１４に記載の発電所計測器性能監視予測方法。
    

    

    ここで、ｋは１、２、．．．、ｍであり、ξ_ｋｊとξ^＊ _ｋｊはスラック変数（Ｓｌａｃｋ Ｖａｒｉａｂｌｅ）を意味する。
16. 前記の最適化問題は次の数式のような双対問題に変えて表すことを特徴とする請求項１５に記載の発電所計測器性能監視予測方法。
    

    

    ここで、ｋは１、２、．．．、ｍである
17. ラグランジュ乗数λ_ｋ，ｊとλ^＊ _ｋ，ｊを次の数式に代入して、ＡＡＳＶＲの第ｋの出力変数に対するＯＲＬを決定することを特徴とする請求項１６に記載の発電所計測器性能監視予測方法。
    

    
18. 元のデータを高次元空間に非線形的に写像した結果をベクトルФ （・）とすれば、非線形写像の結果であるФ（・）の内積として定義されるカーネル（Ｋｅｒｎｅｌ）関数を用いて次の数式のような最適の非線形回帰式を求めることを特徴とする請求項１７に記載の発電所計測器性能監視予測方法。
    

    
19. バイアス項は、次の数式のように任意のＳｕｐｐｏｒｔ Ｖｅｃｔｏｒであるθ_ｒとθ_ｓを用いて計算されることを特徴とする請求項１７に記載の発電所計測器性能監視予測方法。
    

    

Images