(57)【要約】      （修正有）【課題】既存のカーネル回帰法に比して予測値計算の正確度を向上させることができる発電所計測器性能監視予測方法を提供する。【解決手段】計測器信号を行列表示し、正規化し、訓練用、最適化用及びテスト用に分離し、主成分を抽出し、最適化用データを用いてＳＶＲモデルの最適定数を反応分析表面法により求め、最適定数を用いて訓練モデルを生成し、与えられた正規化出力を元の範囲に逆正規化して変数の予測値を求める逆正規化を行い、既存のカーネル回帰法に比べて予測値計算の正確度を向上させる。【選択図】図２
2010256328.tif 000002
【特許請求の範囲】
【請求項１】
全体データを行列の形に表示する行列表示ステップ；
前記全体データをデータセットとして正規化する正規化ステップ；
前記正規化されたデータセットを訓練用、最適化用及びテスト用に分離する分離ステップ；
正規化された各データセットの主成分を抽出する主成分抽出ステップ；
反応表面分析法を用いて最適化用データの予測値誤差を最適化するＳＶＲモデルの最適定数を求める最適定数算出ステップ；
最適定数を用いてＳＶＲモデルを生成する訓練モデル生成ステップ；
正規化されたテストデータを入力としてカーネル関数を求めＳＶＲモデルの出力値を予測する訓練モデル出力予測ステップ；及び
与えられた正規化出力を元の範囲に逆正規化して変数の予測値を求める逆正規化ステップを含むことを特徴とする発電所計測器性能監視予測方法。
【請求項２】
前記行列表示ステップでは、全体データを次の数式のように行列の形に表示することを特徴とする請求項１に記載の発電所計測器性能監視予測方法。
【数１】
2010256328.tif 000062
ここで、Ｘは全体データセットであり、Ｘｔｒは訓練用データセットであり、Ｘｏｐは最適化用データセットであり、Ｘｔｓはテスト用データセットであり、３ｎは全体データの個数であり、ｍは計測器の番号である。
【請求項３】
前記正規化ステップでは、次の数式を用いて全体データを正規化することを特徴とする請求項１に記載の発電所計測器性能監視予測方法。
【数２】
2010256328.tif 000063
ここで、ｉ＝１、２、．．．、３ｎである。
【請求項４】
前記分離ステップでは、データセットＺを訓練用、最適化用及びテスト用の三つに分離し、前記分離されたデータは各々Ｚｔｒ、Ｚｏｐ、Ｚｔｓとし、次の数式のようにｎ×ｍの大きさを有することを特徴とする請求項１に記載の発電所計測器性能監視予測方法。
【数３】
2010256328.tif 000064
ここで、ｉ＝０、１、２、．．．、ｎ−１である。
【請求項５】
前記主成分抽出ステップでは、主成分の分散（すなわち、共分散マトリックスのＥｉｇｅｎｖａｌｕｅ）を大きさ順に並べ、百分率分散値が最も大きい主成分から、その累積和が９９．５％以上になるまでのＺｔｒ、Ｚｏｐ、Ｚｔｓに対する主成分（Ｐｔｒ、Ｐｏｐ、Ｐｔｓ）を選択することで正規化された各データセットＺｔｒ、Ｚｏｐ、Ｚｔｓの主成分を抽出することを特徴とする請求項１に記載の発電所計測器性能監視予測方法。
【請求項６】
各データセットＺｔｒ、Ｚｏｐ、Ｚｔｓから各変数の平均値を引き、これをＡマトリックスとし、次の数式で表すことを特徴とする請求項５に記載の発電所計測器性能監視予測方法。
【数４】
2010256328.tif 000065
【請求項７】
次の数式を用いて、ＡＴＡの根（ｅｉｇｅｎｖａｌｕｅ）λを求め、Ａのｓｉｎｇｕｌａｒ  ｖａｌｕｅ  Ｓを求め、０を除いた根（ｅｉｇｅｎｖａｌｕｅ）λを降順でまとめ、これをλ１、λ２、．．．、λｍと設定することを特徴とする請求項６に記載の発電所計測器性能監視予測方法。
【数５】
2010256328.tif 000066
【請求項８】
ｎ×ｎマトリックスであるＡＡＴのｅｉｇｅｎｖｅｃｔｏｒを求め、ｕｎｉｔａｒｙ  ｍａｔｒｉｘ  Ｕを求め、次の数式からｅｉｇｅｎｖａｌｕｅ  λを求め、これを代入して各ｅｉｇｅｎｖａｌｕｅ  λに対するｎ×１であるｅｉｇｅｎｖｅｃｔｏｒ  ｅ１、ｅ２、．．．、ｅｍを求めることを特徴とする請求項７に記載の発電所計測器性能監視予測方法。
【数６】
2010256328.tif 000067
【請求項９】
次の数式を用いて各主成分の分散を求めることを特徴とする請求項８に記載の発電所計測器性能監視予測方法。
【数７】
2010256328.tif 000068
【請求項１０】
次の数式を用いて各主成分の分散を全体主成分の分散の和で割って百分率を求めることを特徴とする請求項９に記載の発電所計測器性能監視予測方法。
【数８】
2010256328.tif 000069
【請求項１１】
百分率分散（％σp）が最も大きいものから累積計算して、所望の百分率分散（例えば９９．９８％）までの主成分ｐ個を選択することを特徴とする請求項１０に記載の発電所計測器性能監視予測方法。
【請求項１２】
主成分を次の数式で計算して抽出することを特徴とする請求項１１に記載の発電所計測器性能監視予測方法。
【数９】
2010256328.tif 000070
【請求項１３】
次の数式でｍ次元の入力変数ｘ１、ｘ２、．．． 、ｘｍをｐ次元の主成分θ１、θ２、．．．、θｐで圧縮することを特徴とする請求項１２に記載の発電所計測器性能監視予測方法。
【数１０】
2010256328.tif 000071
ここで、ｐはｍ以下の整数である。
【請求項１４】
第ｋの出力に対するＳＶＲ（Ｓｕｐｐｏｒｔ  Ｖｅｃｔｏｒ  Ｒｅｇｒｅｓｓｉｏｎ）で求める最適回帰式（ＯＲＬ：Ｏｐｔｉｍｕｍ  Ｒｅｇｒｅｓｓｉｏｎ  Ｌｉｎｅ）を次の数式で表すことを特徴とする請求項１３に記載の発電所計測器性能監視予測方法。
【数１１】
2010256328.tif 000072
ここで、ｋは１、２、．．．、ｍである。
【請求項１５】
第ｋの出力変数ｙ（ｋ）に対してε‐ｉｎｓｅｎｓｉｔｉｖｅ  Ｌｏｓｓ  Ｆｕｎｃｔｉｏｎを次の数式のように定義すれば、ｙ（ｋ）に対するＯＲＬを求めるための最適化問題はその次の数式のように表すことを特徴とする請求項１４に記載の発電所計測器性能監視予測方法。
【数１２】
2010256328.tif 000073
ここで、ｋは１、２、．．．、ｍであり、ξｋｊとξ＊ｋｊはスラック変数（Ｓｌａｃｋ  Ｖａｒｉａｂｌｅ）を意味する。
【請求項１６】
前記の最適化問題は次の数式のような双対問題に変えて表すことを特徴とする請求項１５に記載の発電所計測器性能監視予測方法。
【数１３】
2010256328.tif 000074
ここで、ｋは１、２、．．．、ｍである
【請求項１７】
ラグランジュ乗数λｋ，ｊとλ＊ｋ，ｊを次の数式に代入して、ＡＡＳＶＲの第ｋの出力変数に対するＯＲＬを決定することを特徴とする請求項１６に記載の発電所計測器性能監視予測方法。
【数１４】
2010256328.tif 000075
【請求項１８】
元のデータを高次元空間に非線形的に写像した結果をベクトルФ （・）とすれば、非線形写像の結果であるФ（・）の内積として定義されるカーネル（Ｋｅｒｎｅｌ）関数を用いて次の数式のような最適の非線形回帰式を求めることを特徴とする請求項１７に記載の発電所計測器性能監視予測方法。
【数１５】
2010256328.tif 000076
【請求項１９】
バイアス項は、次の数式のように任意のＳｕｐｐｏｒｔ  Ｖｅｃｔｏｒであるθｒとθｓを用いて計算されることを特徴とする請求項１７に記載の発電所計測器性能監視予測方法。
【数１６】
2010256328.tif 000077
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、発電所計測器性能監視予測方法に関し、特には、原子力発電所が運転中の状態で、発電所の安全監視用計測器の性能を常時監視するために必要な技術に関する。
【背景技術】
【０００２】
一般的に、全ての発電設備には、運転性向上と安全性確保を目的として多数の計測器を設け、前記計測器からリアルタイムに信号を取得して、発電所監視系統と保護系統に用いている。特に、原子力発電所の安全系統に関連する計測チャネルには、計測信号の正確度と信頼度を保障するために、多重計測器概念を採用しており、運営技術指針書上に記載のように核周期（約１８ヶ月）ごとに点検と矯正を実行している。世界中の原子力発電所は、無駄に行われる計測器矯正業務を改善するために、条件に基づく監視（ＣＢＭ、Ｃｏｎｄｉｔｉｏｎ  Ｂａｓｅｄ  Ｍｏｎｉｔｏｒｉｎｇ）方法論を開発して点検と矯正周期を長くする技術を開発している。
【０００３】
図１は、一般的な計測器性能常時監視システムのブロック構成図であって、Ａｕｔｏ−Ａｓｓｏｃｉａｔｉｖｅモデルと呼ばれるシステムである。図１を参照すれば、一般的な計測器性能常時監視システムは、予測モデル（ｐｒｅｄｉｃｔｉｏｎ  ｍｏｄｅｌ）、比較モジュール（ｃｏｍｐａｒｉｓｏｎ  ｍｏｄｕｌｅ）、判断ロジック（ｄｅｃｉｓｉｏｎ  ｌｏｇｉｃ）で構成されているが、計測信号を予測モデルに入力すると入力測定値に対するモデルの予測値を出力し、比較モジュールにより測定値と予測値との差分を判断ロジックに入力し連続的に監視すれば、計測器のドリフトと故障を感知することができる。
【０００４】
一方、Ａｒｇｏｎｎｅ  Ｎａｔｉｏｎａｌ  ＬａｂｏｒａｔｏｒｙはＭＳＥＴ（Ｍｕｌｔｉｖａｒｉａｔｅ  Ｓｔａｔｅ  Ｅｓｔｉｍａｔｉｏｎ  Ｔｅｃｈｎｉｑｕｅ）を開発して米国特許を取得し、ＳｍａｒｔＳｉｇｎａｌ  Ｃｏｒｐｏｒａｔｉｏｎ及びＥｘｐｅｒｔ  Ｍｉｃｒｏｓｙｓｔｅｍｓ社は同特許を商業的に使用できるよう製品化した。Ｅｘｐｅｒｔ  Ｍｉｃｒｏｓｙｓｔｅｍｓ社は米国内Ｐａｌｏ  Ｖｅｒｄｅ、Ｌｉｍｅｒｉｃｋ  １／２、ＴＭＩ、Ｖ．Ｃ．Ｓｕｍｍｅｒ、Ｓｅｑｕｏｙａｈ１、Ｓａｌｅｍ  １号機でＭＳＥＴを用いる製品を現場に設置して、計測チャネルオンライン監視を行っている。また、ＳｍａｒｔＳｉｇｎａｌ  Ｃｏｒｐｏｒａｔｉｏｎはその後ＭＳＥＴに対する特許を使用できなくなり、カーネル回帰法に基づく計測器性能常時監視技術を開発した。
【０００５】
計測器の予測値を計算するための方法としては、線形回帰分析法（Ｌｉｎｅａｒ  Ｒｅｇｒｅｓｓｉｏｎ  Ｍｅｔｈｏｄ）が一般的に使用される。この方法は、予測しようとする計測器信号と線形的な相関関係の高い他の計測器信号を選択し、予測値と測定値との誤差二乗和が最小になるよう回帰係数を求める方法であって、次の数１により表現できる。
【０００６】
【数１】
2010256328.tif 000003
【０００７】
前記の線形回帰分析法は、既知の従属変数と独立変数で回帰係数が決まると、未知の従属変数に対する独立変数を予測できる。
【０００８】
しかしながら、既存の線形回帰分析方法論は、従属変数同士で互いに線形的連関性の大きい場合、多重共線性の問題が生じ、従属変数に含まれる小さいノイズに対して独立変数に大きい誤差が発生する。
【０００９】
カーネル回帰法は、既存の線形回帰分析方法論や神経回路網のように入力と出力との相関関係を最適化する回帰係数や加重値のような媒介変数（Ｐａｒａｍｅｔｅｒ）を使用せず、選別された測定データをメモリベクトルとして格納し、測定信号セットに対するメモリベクトル内の訓練データセットのユークリッド距離からカーネルの加重値を求め、これをメモリベクトルに適用して計測器の予測値を求める非媒介変数回帰法（Ｎｏｎｐａｒａｍｅｔｒｉｃ  ｒｅｇｒｅｓｓｉｏｎ  ｍｅｔｈｏｄ）である。カーネル回帰法のような非媒介変数回帰法は、入出力関係が非線形状態のモデルと信号雑音に強いという長所を有する。
【００１０】
既存のカーネル回帰法の計算手順は、下記の５段階の過程からなる。
【００１１】
１段階：訓練データを次の数２の行列の形に表す。
【００１２】
【数２】
2010256328.tif 000004
【００１３】
ここで、Ｘはメモリベクトルに格納される訓練データの行列であって、ｎは訓練データ個数であり、ｍは計測器の番号である。
【００１４】
２段階：次の数３を用いて第１の計測器信号セットに対する訓練データのユークリッド距離の和を求める。
【００１５】
【数３】
2010256328.tif 000005
【００１６】
ここで、ｘは訓練データ、ｑはテストデータまたは問い合わせデータ（Ｑｕｅｒｙ  ｄａｔａ）、ｔｒｎは訓練データの番号、ｊは計測器の番号である。
【００１７】
３段階：次のカーネル関数を含める数４を用いて各々の訓練データセットと与えられたテストデータセットに対する加重値を求める。
【００１８】
【数４】
2010256328.tif 000006
【００１９】
前記の数４においては、加重関数として用いられたガウスカーネルは次のように定義される。
【００２０】
【数５】
2010256328.tif 000007
【００２１】
４段階：テストデータの予測値は、次の数６のように各々の訓練データに加重値を掛けた後、加重値の和で割って求める。
【００２２】
【数６】
2010256328.tif 000008
【００２３】
５段階：全体のテストデータに対する予測値を求めるために、２段階から４段階の過程を繰り返す。
このような従来のカーネル回帰法は、非線形状態のモデルと信号雑音に強い長所を有している。
【００２４】
しかし、前記の従来のカーネル回帰法は、選別された測定データをメモリベクトルとして格納し、測定信号セットに対するメモリベクトル内の訓練データセットのユークリッド距離からカーネルの加重値を求め、これをメモリベクトルに適用して計測器の予測値を求めるので、出力予測値の分散が大きくなり、これにより線形回帰分析法に比べて正確度の低いという問題点がある。
【発明の概要】
【発明が解決しようとする課題】
【００２５】
本発明の目的は、上記のような従来のカーネル回帰法の予測正確度の低い問題点を解消するために提案されたものであり、発電所データの正規化、主成分抽出、反応表面分析法を用いるＳＶＲ（Ｓｕｐｐｏｒｔ  Ｖｅｃｔｏｒ  Ｒｅｇｒｅｓｓｉｏｎ）モデル回帰式のパラメータ（カーネル帯域幅σ、損失関数ε−、ペナルティＣ）最適化、これを用いるＳＶＲモデル具現、出力予測値の逆正規化方法等を用いて発電所システムをモデリングした後、計測器信号を監視することにより、既存のカーネル回帰法に比して予測値計算の正確度を向上させることができる発電所計測器性能監視予測方法を提供することにある。
【課題を解決するための手段】
【００２６】
上述の目的を達成するために、本発明の発電所計測器性能監視予測方法は、発電所データの主成分を抽出しＳＶＲ（Ｓｕｐｐｏｒｔ  Ｖｅｃｔｏｒ  Ｒｅｇｒｅｓｓｉｏｎ）方法を用いてシステムを各種ケースに対してモデリングした後、回帰式の三つのパラメータを反応表面分析法を用いて最適化し、これを用いて発電所システムをモデリングした後、計測器信号を監視することにより、既存の広く使用されるカーネル回帰法に比して予測値計算の正確度を向上させることができる。
【発明の効果】
【００２７】
本発明の発電所計測器性能監視用予測方法は、原子力発電所の安全監視チャネルに用いられる計測器の性能を発電所が運転中の状態でオンラインで監視できるようにすることで、計測器の誤作動をリアルタイムに監視して計測器信頼度を向上させ、原子力発電所の計測器矯正周期を現在燃料交替周期の１８ヶ月から最大８年に延長することにより、矯正費用と放射線区域における矯正作業従事者の放射線被爆を低減し、無駄な矯正回数を減らすことで、誤矯正による発電所不時停止を予防し、発電所予防停止期間を短縮して発電所利用率を増加させることができる効果を有する。
【図面の簡単な説明】
【００２８】
【図１】従来の発電所計測器の性能常時監視システムのブロック構成図である。
【図２】本発明の一実施例にかかる発電所計測器性能監視予測方法を説明するための手順を示すフローチャートである。
【図３】図２の発電所計測器性能監視予測方法に従って生成されたＳＶＲモデルである。
【図４】ＳＶＲによる最適回帰式の一般的な概念図である。
【図５】モデルパラメータが三つである場合の中心合成計画における実験点を示す。
【図６Ａ】正確度テストのための原子力発電所原子炉炉心出力データである。
【図６Ｂ】正確度テストのための原子力発電所原子炉炉心出力データである。
【図７Ａ】正確度テストのための原子力発電所加圧器水位データである。
【図７Ｂ】正確度テストのための原子力発電所加圧器水位データである。
【図８Ａ】正確度テストのための原子力発電所蒸気発生器蒸気流量データである。
【図８Ｂ】正確度テストのための原子力発電所蒸気発生器蒸気流量データである。
【図９Ａ】正確度テストのための原子力発電所蒸気発生器狭域水位データである。
【図９Ｂ】正確度テストのための原子力発電所蒸気発生器狭域水位データである。
【図１０Ａ】正確度テストのための原子力発電所蒸気発生器圧力データである。
【図１０Ｂ】正確度テストのための原子力発電所蒸気発生器圧力データである。
【図１１Ａ】正確度テストのための原子力発電所蒸気発生器広域水位データである。
【図１１Ｂ】正確度テストのための原子力発電所蒸気発生器広域水位データである。
【図１２Ａ】正確度テストのための原子力発電所蒸気発生器主給水流量データである。
【図１２Ｂ】正確度テストのための原子力発電所蒸気発生器主給水流量データである。
【図１３Ａ】正確度テストのための原子力発電所タービン出力データである。
【図１３Ｂ】正確度テストのための原子力発電所タービン出力データである。
【図１４Ａ】正確度テストのための原子力発電所１次側充電流量データである。
【図１４Ｂ】正確度テストのための原子力発電所１次側充電流量データである。
【図１５Ａ】正確度テストのための原子力発電所残熱除去流量データである。
【図１５Ｂ】正確度テストのための原子力発電所残熱除去流量データである。
【図１６Ａ】正確度テストのための原子力発電所原子炉上部冷却材温度データである。
【図１６Ｂ】正確度テストのための原子力発電所原子炉上部冷却材温度データである。
【発明を実施するための形態】
【００２９】
以下、本発明が属する技術分野で通常の知識を有する者が本発明を容易に実施できる程度に詳細の説明するために、本発明の最も好ましい実施例を添付の図面を参照して詳しく説明する。本発明の目的、作用、効果を含めてその他の目的、特徴、そして動作上の利点が好ましい実施例の説明によってさらに明確になるであろう。
【実施例】
【００３０】
図２は、本発明の一実施例にかかる発電所計測器性能監視予測方法を説明するための手順を示すフローチャートである。図２に示すように、本発明の一実施例にかかる発電所計測器性能監視予測方法は、全体データを行列の形に表示する行列表示ステップＳ１００と、前記全体データをデータセットとして正規化する正規化ステップＳ２００と、データセットを訓練用、最適化用、テスト用に分離する分離ステップＳ３００と、正規化された各データセットの主成分を抽出する主成分抽出ステップＳ４００と、反応表面分析法を用いて最適化用データの予測値誤差を最適化するＳＶＲモデルの最適定数を求める最適定数算出ステップＳ５００と、最適定数を用いてＳＶＲモデルを生成する訓練モデル生成ステップＳ６００と、正規化されたテストデータを入力としてカーネル関数を求めＳＶＲモデルの出力値を予測する訓練モデル出力予測ステップＳ７００と、与えられた正規化出力を元の範囲に逆正規化して変数の予測値を求める逆正規化ステップＳ８００とを含む。
【００３１】
前記の構成による本発明の一実施例にかかる発電所計測器性能監視予測方法の作用は次の通りである。
【００３２】
まず、行列表示ステップＳ１００では、全体データセット（Ｘ）を次の数７のように行列の形に表す。
【００３３】
【数７】
2010256328.tif 000009
【００３４】
ここで、Ｘは全体データセットであり、Ｘｔｒは訓練用データセットであり、Ｘｏｐは最適化用データセットであり、Ｘｔｓはテスト用データセットであり、３ｎは全体データの個数であり、ｍは計測器の番号である。
【００３５】
次に、正規化ステップＳ２００では，次の数８を用いて全体データを正規化する。
【００３６】
【数８】
2010256328.tif 000010
【００３７】
ここで、ｉ＝１、２、．．．、３ｎである。
【００３８】
正規化された全体データセット（ｓｅｔ）Ｚは次の数９のように表すことができる。
【００３９】
【数９】
2010256328.tif 000011
【００４０】
前記分離ステップＳ３００では、データセットＺを訓練（Ｔｒａｉｎｉｎｇ）、最適化（Ｏｐｔｉｍｉｚａｔｉｏｎ）、テスト（Ｔｅｓｔ）用の三つに分離する。本発明の実施例において、三つに分離されたデータは各々Ｚｔｒ、Ｚｏｐ、Ｚｔｓとし、次の数１０のようにｎ×ｍの大きさを有することができる。
【００４１】
【数１０】
2010256328.tif 000012
【００４２】
ここで、ｉ＝０、１、２、．．．、ｎ−１である。
【００４３】
前記主成分抽出ステップＳ４００では、正規化された各データセットＺｔｒ、Ｚｏｐ、Ｚｔｓの主成分を抽出する。主成分の分散（すなわち、共分散マトリックスのＥｉｇｅｎｖａｌｕｅ）を大きさ順に並べ、百分率分散値が最も大きい主成分から、その累積和が使用者が定めた９９．５％以上になるまでのＺｔｒ、Ｚｏｐ、Ｚｔｓに対する主成分（Ｐｔｒ、Ｐｏｐ、Ｐｔｓ）を選択する。
【００４４】
-主成分を求める方法
主成分分析（ＰＣＡ）は、線形変換による多数の入力変数の少数の変数への圧縮に有効な方法である。このとき、圧縮された変数を主成分（Ｐｒｉｎｃｉｐａｌ  Ｃｏｍｐｏｎｅｎｔ）と呼ぶ。ＰＣＡは変数間相関関係を用いて元の次元のデータを二乗和（Ｓｕｍ  ｏｆ  Ｓｑｕａｒｅｓ）が最大化する低次元の超平面に放射する。
以下では、主成分を抽出する過程を説明する。
ｍ次元の入力変数をｘ１、ｘ２、．．．、ｘｍとし、これらの線形結合で生成される新しい変数をθ１、θ２、．．．、θｍとすれば、これらの関係は次の数１１及び数１２のように表現できる。
【００４５】
【数１１】
2010256328.tif 000013
【００４６】
【数１２】
2010256328.tif 000014
【００４７】
このとき、θ１、θ２、．．．、θｍを線形システムの主成分と言い、便宜上第１の主成分を最も重要な主成分と決める。最も重要な主成分とは、入力変数の変化を最も多く説明する、すなわち、分散の最も大きい主成分を意味する。線形変換Ｑは、次の数１３及び数１４の条件を満たすように決定する。
【００４８】
【数１３】
2010256328.tif 000015
【００４９】
【数１４】
2010256328.tif 000016
【００５０】
前記数１３の条件は、変換後のスケールをそのまま維持するようにし、前記数１４の条件は、変換後には変数間に互いに相関関係が無いようにする。
【００５１】
主成分は次のような手順を通じて求められる。
イ．各データセットＺｔｒ、Ｚｏｐ、Ｚｔｓから各変数の平均値を引き、これをＡマトリックスとする。本発明の実施例では、Ｚｔｒを例示して説明し、これを数式で表すと次の数１５のようになる。
【００５２】
【数１５】
2010256328.tif 000017
【００５３】
ロ．次の数１６乃至数１８を用いて、ＡＴＡの根（ｅｉｇｅｎｖａｌｕｅ）λを求め、Ａのｓｉｎｇｕｌａｒ  ｖａｌｕｅ  Ｓを求める。数１６から求めた０を除いた根（ｅｉｇｅｎｖａｌｕｅ）λを降順でまとめ、これをλ１、λ２、．．．、λｍとする。
【００５４】
【数１６】
2010256328.tif 000018
【００５５】
【数１７】
2010256328.tif 000019
【００５６】
【数１８】
2010256328.tif 000020
【００５７】
ハ．ｎ×ｎマトリックスであるＡＡＴのｅｉｇｅｎｖｅｃｔｏｒを求め、ｕｎｉｔａｒｙ  ｍａｔｒｉｘ  Ｕを求める。次の数１９からｅｉｇｅｎｖａｌｕｅ  λを求め、これを数２０に代入して各ｅｉｇｅｎｖａｌｕｅ  λに対するｎ×１であるｅｉｇｅｎｖｅｃｔｏｒ  ｅ１、ｅ２、．．．、ｅｍを求める。
【００５８】
【数１９】
2010256328.tif 000021
【００５９】
【数２０】
2010256328.tif 000022
【００６０】
ニ．次の数２１を用いて各主成分の分散を求める。
【００６１】
【数２１】
2010256328.tif 000023
【００６２】
ホ．次の数２２及び数２３を用いて各主成分の分散を全体主成分の分散の和で割って百分率を求める。
【００６３】
【数２２】
2010256328.tif 000024
【００６４】
【数２３】
2010256328.tif 000025
【００６５】
ヘ．百分率分散（％σp）が最も大きいものから累積計算して、所望の百分率分散（例えば９９．９８％）までの主成分ｐ個を選択する。
【００６６】
ト．主成分を次の数２４で計算して抽出する。
【００６７】
【数２４】
2010256328.tif 000026
【００６８】
チ．Ｚｏｐ、Ｚｔｓに対しても上記のような手順により主成分を抽出することができる。
【００６９】
次の表１は、抽出された主成分の分散を示す。
【００７０】
【表１】
2010256328.tif 000027
【００７１】
本発明の実施例では７個の主成分を使用するが、前記表１を参照すれば、７個の主成分を使用する場合、全体分散の９９．９８％を説明でき、残りの主成分を放棄することで発生する情報の損失は０．０２％に過ぎないことが分かる。
【００７２】
‐ＳＶＲ（Ｓｕｐｐｏｒｔ  Ｖｅｃｔｏｒ  Ｒｅｇｒｅｓｓｉｏｎ）モデリング
ｍ次元の入力変数ｘ１、ｘ２、．．． 、ｘｍをｐ次元の主成分θ１、θ２、．．．、θｐで圧縮するのは、次の数２５のように表すことができる。
【００７３】
【数２５】
2010256328.tif 000028
【００７４】
ここで、ｐはｍ以下の整数である。
【００７５】
第ｋの出力に対するＳＶＲで求める最適回帰式（ＯＲＬ：Ｏｐｔｉｍｕｍ  Ｒｅｇｒｅｓｓｉｏｎ  Ｌｉｎｅ）は次の数２６のように表すことができる。
【００７６】
【数２６】
2010256328.tif 000029
【００７７】
ここで、ｋは１、２、．．．、ｍである。
【００７８】
第ｋの出力変数ｙ（ｋ）に対してε‐ｉｎｓｅｎｓｉｔｉｖｅ  Ｌｏｓｓ  Ｆｕｎｃｔｉｏｎを次の数２７のように定義すれば、ｙ（ｋ）に対する最適回帰式（ＯＲＬ）を求めるための最適化問題は次の数２８のように表すことができる。
【００７９】
【数２７】
2010256328.tif 000030
【００８０】
【数２８】
2010256328.tif 000031
【００８１】
前記数２７及び数２８において、ｋは１、２、．．．、ｍである。そして、ξｋ，ｉとξ＊ｋ，ｉは、図４に示したスラック変数（Ｓｌａｃｋ  Ｖａｒｉａｂｌｅ）を意味する。図４は、ＳＶＲによる最適回帰線の概念図である。ここで、θｉは、ベクトルθの第ｉの要素でなく、ｘに対する第ｉの観測値ベクトルに対応する主成分ベクトルを意味する。
【００８２】
前記の最適化問題は次の数２９のような双対問題に変えて表すことができる。
【００８３】
【数２９】
2010256328.tif 000032
【００８４】
ここで、ｋは１、２、．．．、ｍである。
【００８５】
これにより、ラグランジュ乗数λｋ，ｉとλ＊ｋ，ｉを次の数３０に代入して、ＡＡＳＶＲの第ｋの出力変数に対する最適回帰式（ＯＲＬ）を決定することができる。
【００８６】
【数３０】
2010256328.tif 000033
【００８７】
上述した内容は最適の線形回帰式を求める過程である。元のデータを高次元空間に非線形的に写像した結果をベクトルФ（・）とすれば、非線形写像の結果であるФ（・）の内積として次の数３１のように定義される関数をカーネル（Ｋｅｒｎｅｌ）と呼ぶ。
【００８８】
【数３１】
2010256328.tif 000034
【００８９】
高次元空間で最適回帰式（ＯＲＬ）を探そうとするとき、個別のФ（ｘｉ）とФ（ｘｊ）を全て知るべき理由が無く、ただカーネル関数Ｋ（ｘｉ、ｘｊ）さえ知っていれば十分である。本発明の実施例においては、Ｇａｕｓｓｉａｎ  Ｒａｄｉａｌ  Ｂａｓｉｓ  Ｆｕｎｃｔｉｏｎを使用した。
【００９０】
カーネル関数Ｋ（θｉ、θ）＝Ф（θｉ）ＴФ（θ）を用いれば、次の数３２のような最適の非線形回帰式を求めることができる。
【００９１】
【数３２】
2010256328.tif 000035
【００９２】
ここで、バイアス項は、次の数３３のように任意のＳｕｐｐｏｒｔ  Ｖｅｃｔｏｒであるθｒとθｓを用いて計算されることができる。
【００９３】
【数３３】
2010256328.tif 000036
【００９４】
以上から得られた最適の非線形回帰式はｌｏｓｓ  ｆｕｎｃｔｉｏｎの定数ε（ｅｐｓｉｌｏｎ）、双対目的関数Ф（ｗｋ、ξｋ）にペナルティＣ、ＲＢＦ（Ｒａｄｉａｌ  Ｂａｓｉｓ  Ｆｕｎｃｔｉｏｎ）をカーネルとして用いる場合、カーネル帯域幅σが予め与えられるのが好ましい。この過程を繰り返して各々の出力に対する総ｍ個のＳＶＲを得て、図３に示したＡＡＳＶＲを構築することができる。
【００９５】
前記最適定数算出ステップＳ５００では、反応表面分析法を用いて最適化データ（Ｚｏｐ）の予測値誤差を最小化するＳＶＲモデルの最適定数（σ、ε、Ｃ）を求める。
【００９６】
イ．ＳＶＲモデルパラメータであるｓｉｇｍａ（σ）、ｅｐｓｉｌｏｎ（ε）、Ｃを各々ｖ１、ｖ２、ｖ３とする。
ロ．ｖ１、ｖ２、ｖ３に対する探索範囲を各々決める。適切な探索範囲は事前経験や小規模の予備実験により把握できる。本発明の実施例においては、各々ｖ１は０．５６４８６−１．６３５１４、ｖ２は０．０１０５３４−０．０３９９６６、ｖ３は２．１０７６−７．９９３３とする。
【００９７】
ハ．探索範囲の上限を各々Ｕ１、Ｕ２、Ｕ３とし、下限をＬ１、Ｌ２、Ｌ３とし、次の数３４のようにモデルパラメータを標準化する。
【００９８】
【数３４】
2010256328.tif 000037
【００９９】
ニ．標準化されたモデルパラメータｘ１、ｘ２、ｘ３の探索範囲を考慮して実験点、すなわちモデル性能の評価地点を定める。このために、統計的な実験計画の一つの中心合成計画を用いることができる。中心合成計画により定められる実験点を３次元空間で表現すれば、図５のようになる。
【０１００】
ホ．中心合成計画による実験点は８個の頂点、１個の中心点、そして６個の軸点で構成される。実験誤差の大きさを推定するためには、５回内外繰り返して実験を施す。軸点の座標は予測分散に関する統計的な性質を考慮してα＝２３／４＝１．６８１７９に定める。
【０１０１】
α＝［要因実験点の数］１／４
【０１０２】
中心点で５回繰り返した場合に、中心合成計画によるｘ１、ｘ２、ｘ３の実験点は次の表２のようになる。
【０１０３】
【表２】
2010256328.tif 000038
【０１０４】
ヘ．前記表２に指示するようにモデルパラメータｖ１、ｖ２、ｖ３の値を決めて表３（中心合成計画による（ｖ１、ｖ２、ｖ３）の実験点）を得られるが、これがモデルを得るために直接使用されるモデルパラメータの値となる。
【０１０５】
【表３】
2010256328.tif 000039
【０１０６】
ト．前記表３の各実験点でＺｔｒとＰｔｒを用いてＳＶＲモデルのｂｅｔａベクトルとｂｉａｓ定数を各々得る。そのためにはデータセットＺｔｒを用いる。実際、中心点に該当するＮｏ．１５からＮｏ．１９までは同じモデルが得られる。
【０１０７】
チ．各モデルの正確度を評価するためにデータセットＰｏｐをｍ個のＡＡＳＶＲに入力して最適化データの正規化された予測値数３５を求める。これから次の数３６を用いて出力モデルの正確度、すなわちＭＳＥを計算する。５個の中心点は同じモデルであるので、Ｐｏｐを５個に分割し、分割された下部データセットから別のＭＳＥ計算結果を得る。
【０１０８】
【数３５】
2010256328.tif 000040
【０１０９】
【数３６】
2010256328.tif 000041
【０１１０】
ここで、ｚｉｊはＰｏｐのうちセンサーｉの第ｊの入力データを表し、Ｚｉｊはモデルによる推定値である。モデルはＰｔｒを用いて得られたことを想起する。本発明の実施例において、実験結果である実験によるＭＳＥ計算結果は、次の表４に表すようになる。
【０１１１】
【表４】
2010256328.tif 000042
【０１１２】
リ．反応表面を求めるときには、ＭＳＥの代わりにｌｏｇ（ＭＳＥ）を使用する。このことを考慮して、モデルパラメータｘ１、ｘ２、ｘ３とｌｏｇ（ＭＳＥ）との間の反応表面を推定する。反応表面は次の数３７のような２次模型を仮定する。
【０１１３】
【数３７】
2010256328.tif 000043
【０１１４】
ここで、ｅはランダム誤差を意味する。
【０１１５】
本発明の実施例で推定された反応表面は下記のようになる。
ｌｏｇ（ＭＳＥ）＝−８．３４９２，−０．２１３１ｘ１，＋０．７７１６ｘ２，−０．０９５２ｘ３，＋０．２０１０ｘ１２，−０．０７５３ｘ２２＋０．０７９９ｘ３２
【０１１６】
ヌ．推定された反応表面式を用いてｌｏｇ（ＭＳＥ）を最小化するｘ１、ｘ２、ｘ３の最適条件を求める。２次反応表面を仮定したので、最適条件は偏微分を通して確認する。すなわち、次の数３８を同時に満たすｘ１、ｘ２、ｘ３を求める。
【０１１７】
【数３８】
2010256328.tif 000044
【０１１８】
本発明の実施例で得られた反応表面の場合、最適条件は（ｘ＊１、ｘ＊２、ｘ＊３）＝（１．５４３８、−０．６８１８、１．５９２９）である。
ル．最適条件（ｘ＊１、ｘ＊２、ｘ＊３）を次の数３９を用いて元の単位に換算する。
【０１１９】
【数３９】
2010256328.tif 000045
【０１２０】
本発明の実施例において、ｖ＊１＝１．３９１０＝σ＊、ｖ＊２＝０．０００５＝ε＊、ｖ＊３＝６．７９５１＝Ｃ＊になり、この条件で予測されたｌｏｇ（ＭＳＥ）は−９．９４４６なので、ここに指数を利用してＭＳＥに換算すると０．００００４８になる。
【０１２１】
前記訓練モデル生成ステップＳ６００では、ステップＳ５００で求めたＳＶＲモデルの三つの最適定数（ε＊、Ｃ＊、σ＊）と訓練データの主成分（Ｐｔｒ）、第１の信号の訓練データ（Ｚｔｒの第１列）を入力としてｑｕａｄｒａｔｉｃ  ｐｒｏｇｒａｍを用いて最適化問題を解き、ラグランジュ乗数の差分であるβ１（ｎ×１）とバイアス定数ｂ１を求めて、図３のＳＶＲ１のモデルを生成する。同じ方法で第２から第ｍのセンサーに対してこれを繰り返して行い、β２、β３、．．．、βｍとｂ２、ｂ３、．．．、ｂｍを求めＳＶＲ２からＳＶＲｍのモデルを生成して、図３のように全センサーに対するＳＶＲモデルを構築する。
【０１２２】
前記訓練モデルの出力予測ステップＳ７００では、訓練データの主成分（Ｐｔｒ）、テストデータの主成分（Ｐｔｓ）を用いてＧａｕｓｓｉａｎ  Ｒａｄｉａｌ  Ｂａｓｉｓ  Ｆｕｎｃｔｉｏｎのカーネル関数（Ｋ（ｎ×ｎ））を求め、ステップＳ６００で求めたＳＶＲモデルのラグランジュ乗数の差分β１、バイアス定数ｂ１を用いてＳＶＲ１の出力を求める。同じ方法で第２から第ｍのセンサーに対してこれを繰り返して行い、ＳＶＲ２からＳＶＲｍの出力のモデル予測値を求めることができる。
【０１２３】
これを数式で表すと、次の数４０のようになる。
【０１２４】
【数４０】
2010256328.tif 000046
【０１２５】
前記逆正規化ステップＳ８００では、ステップＳ７００で得られた正規化されたテストデータの予測値を元の範囲に逆正規化して元のスケールの各センサーに対する予測値を求める。これを数式で表すと、次の数４１のようになる。
【０１２６】
【数４１】
2010256328.tif 000047
【０１２７】
本発明の一実施例にかかる発電所計測器性能監視予測方法の優秀性を確認するために、実際原子力発電所の出力を０％から１００％に上昇させる途中１、２次系統で測定した計測器信号データを用いて既存の方法論と比較した。分析に使用されたデータは総１１個のセンサーにて測定された値である。
次の表５は、従来のカーネル回帰法と本発明の一実施例にかかる発電所計測器監視方法による計測器予測値の正確度を比較した表である。
【０１２８】
【表５】
2010256328.tif 000048
【０１２９】
前記表５で分析に使用されたデータは、次のような総１１個のセンサーにて測定された値である。
１．原子炉出力（％）
２．加圧器水位（％）
３．蒸気発生器蒸気流量（Ｍｋｇ／ｈｒ）
４．蒸気発生器狭域水位データ（％）
５．蒸気発生器圧力データ（Ｋｇ／ｃｍ２）
６．蒸気発生器広域水位データ（％）
７．蒸気発生器主給水流量データ（Ｍｋｇ／ｈｒ）
８．タービン出力データ（ＭＷｅ）
９．原子炉冷却材充電流量データ（ｍ３／ｈｒ）
１０．残熱除去流量データ（ｍ３／ｈｒ）
１１．原子炉上部冷却材温度データ（℃）
【０１３０】
上記のデータのモデル入力及び予測値を時間の関数にしてグラフで描くと、図６乃至図１６のようになる。
【０１３１】
図６Ａ及び図６Ｂは、原子力発電所原子炉炉心出力データに関し、図６Ａは、数７のテスト用入力データＸｔｓ１に該当し、図６Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ１に対する推定データ数４２を示す。
【０１３２】
【数４２】
2010256328.tif 000049
【０１３３】
図７Ａ及び図７Ｂは、原子力発電所加圧器水位データに関し、図７Ａは、数７のテスト用入力データＸｔｓ２に該当し、図７Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ２に対する推定データ数４３を示す。
【０１３４】
【数４３】
2010256328.tif 000050
【０１３５】
図８Ａ及び図８Ｂは、原子力発電所蒸気発生器蒸気流量データに関し、図８Ａは、数７のテスト用入力データＸｔｓ３に該当し、図８Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ３に対する推定データ数４４を示す。
【０１３６】
【数４４】
2010256328.tif 000051
【０１３７】
図９Ａ及び図９Ｂは、原子力発電所蒸気発生器狭域水位データに関し、図９Ａは、数７のテスト用入力データＸｔｓ４に該当し、図９Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ４に対する推定データ数４５を示す。
【０１３８】
【数４５】
2010256328.tif 000052
【０１３９】
図１０Ａ及び図１０Ｂは、原子力発電所蒸気発生器蒸気圧力データに関し、図１０Ａは、数７のテスト用入力データＸｔｓ５に該当し、図１０Ｂは，本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ５に対する推定データ数４６を示す。
【０１４０】
【数４６】
2010256328.tif 000053
【０１４１】
図１１Ａ及び図１１Ｂは、原子力発電所蒸気発生器広域水位データに関し、図１１Ａは、数７のテスト用入力データＸｔｓ６に該当し、図１１Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ６に対する推定データ数４７を示す。
【０１４２】
【数４７】
2010256328.tif 000054
【０１４３】
図１２Ａ及び図１２Ｂは、原子力発電所蒸気発生器主給水流量データに関し、図１２Ａは、数７のテスト用入力データＸｔｓ７に該当し、図１２Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ７に対する推定データ数４８を示す。
【０１４４】
【数４８】
2010256328.tif 000055
【０１４５】
図１３Ａ及び図１３Ｂは、原子力発電所タービン出力データに関し、図１３Ａは、数７のテスト用入力データＸｔｓ８に該当し、図１３Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ８に対する推定データ数４９を示す。
【０１４６】
【数４９】
2010256328.tif 000056
【０１４７】
図１４Ａ及び図１４Ｂは、原子力発電所１次側充電流量データに関し、図１４Ａは、数７のテスト用入力データＸｔｓ９に該当し、図１４Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ９に対する推定データ数５０を示す。
【０１４８】
【数５０】
2010256328.tif 000057
【０１４９】
図１５Ａ及び図１５Ｂは、原子力発電所残熱除去流量データに関し、図１５Ａは、数７のテスト用入力データＸｔｓ１０に該当し、図１５Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ１０に対する推定データ数５１を示す。
【０１５０】
【数５１】
2010256328.tif 000058
【０１５１】
図１６Ａ及び図１６Ｂは、原子力発電所原子炉上部冷却材温度データに関し、図１６Ａは、数７のテスト用入力データＸｔｓ１１に該当し、図１６Ｂは、本発明のアルゴリズムを用いて予測した数４１のテスト用入力Ｘｔｓ１１に対する推定データ数５２を示す。
【０１５２】
【数５２】
2010256328.tif 000059
【０１５３】
表５において、正確度は予測モデルを運転監視に適用するにおいて最も基本的な尺度となる。殆ど正確度はモデル予測値と実際測定値との平均自乗誤差で表す。
【０１５４】
次の数５３は、一つの計測器に対する正確度を表す数式である。
【０１５５】
【数５３】
2010256328.tif 000060
【０１５６】
ここで、Ｎはテストデータの数、数５４は第ｉのテストデータに対するモデルの推定値、ｘｉは第ｉのテストデータの測定値を表す。
【０１５７】
【数５４】
2010256328.tif 000061
【０１５８】
本発明の一実施例にかかる発電所計測器性能監視予測方法は、計測器信号の主成分を抽出し、最適化用データを用いてＳＶＲモデルの最適定数を反応分析表面法により求め、またモデルを訓練データを用いて訓練させて、既存のカーネル回帰法に比して予測値計算の正確度を向上させる。
【０１５９】
上述した本発明の詳細な説明では本発明の好ましい実施例を参照して説明したが、本発明の属する技術分野における通常の知識を有する者であれば、本発明の技術的思想及び技術領域を逸脱しない範囲内で、本発明の様々な修正及び変更が可能であることが理解できるはずである。したがって、本発明の技術的範囲は、明細書の詳細な説明に記載の内容に限定されるのではない。
【図１】
2010256328.tif 000078
【図２】
2010256328.tif 000079
【図３】
2010256328.tif 000080
【図４】
2010256328.tif 000081
【図５】
2010256328.tif 000082
【図６Ａ】
2010256328.tif 000083
【図６Ｂ】
2010256328.tif 000084
【図７Ａ】
2010256328.tif 000085
【図７Ｂ】
2010256328.tif 000086
【図８Ａ】
2010256328.tif 000087
【図８Ｂ】
2010256328.tif 000088
【図９Ａ】
2010256328.tif 000089
【図９Ｂ】
2010256328.tif 000090
【図１０Ａ】
2010256328.tif 000091
【図１０Ｂ】
2010256328.tif 000092
【図１１Ａ】
2010256328.tif 000093
【図１１Ｂ】
2010256328.tif 000094
【図１２Ａ】
2010256328.tif 000095
【図１２Ｂ】
2010256328.tif 000096
【図１３Ａ】
2010256328.tif 000097
【図１３Ｂ】
2010256328.tif 000098
【図１４Ａ】
2010256328.tif 000099
【図１４Ｂ】
2010256328.tif 000100
【図１５Ａ】
2010256328.tif 000101
【図１５Ｂ】
2010256328.tif 000102
【図１６Ａ】
2010256328.tif 000103
【図１６Ｂ】
2010256328.tif 000104