(57)【要約】      （修正有）【課題】高調波歪を有する平衡３相電力系統(３ＰＨ)及び非平衡３相電力系統における基本周波数及び位相を推定する方法を提供する。【解決手段】信号のパラメーターを推定する方法であって、前記信号は、高調波歪を有する正弦波の３相混合体であり、該方法は、前記信号の４分の１サイクル以下の間、前記信号のサンプルを取得するステップと、前記サンプルに分類器を適用して、雑音パラメーター及び信号パラメーターを得るステップと、初期パラメーターに推定器を適用して、最終パラメーターを得るステップと、を含む。前記分類器は、多重信号分類器(ＭＵＳＩＣ)であり、該多重信号分類器は、前記推定器用に前記信号から最も強い高調波を抽出するステップと、次の反復の前に、前記信号から１次高調波を取り除くステップと、を反復的に実行する。【選択図】図３
2014153354.tif 000002
【特許請求の範囲】
【請求項１】
信号のパラメーターを推定する方法であって、前記信号は、高調波歪を有する正弦波の３相混合体であり、該方法は、
前記信号の４分の１サイクル以下の間、前記信号のサンプルを取得するステップと、
前記サンプルに分類器を適用して、雑音パラメーター及び信号パラメーターを得るステップと、
初期パラメーターに推定器を適用して、最終パラメーターを得るステップと、
を含む、信号のパラメーターを推定する方法。
【請求項２】
前記分類器は、多重信号分類器(ＭＵＳＩＣ)であり、該多重信号分類器は、
前記推定器用に前記信号から最も強い高調波を抽出するステップと、
次の反復の前に、前記信号から１次高調波を取り除くステップと、
を反復的に実行する、請求項１に記載の方法。
【請求項３】
前記分類器は、
前記サンプルから共分散行列を構成するステップと、
前記共分散行列の固有値分解から雑音部分空間を求めるステップと、
前記雑音部分空間の多項式係数を求めるステップと、
前記多項式係数の累乗根を前記抽出するステップに渡すステップと、
を含む、請求項２に記載の方法。
【請求項４】
前記パラメーターは、前記信号の周波数、位相、及び振幅を含む、請求項１に記載の方法。
【請求項５】
前記信号にクラーク変換を適用するステップ、
を更に含む、請求項１に記載の方法。
【請求項６】
サンプルの数が、前記信号の基本周波数の４分の１サイクルよりも長い時間期間で得られるサンプルの数よりも少ない、請求項１に記載の方法。
【請求項７】
前記推定器は加重最小二乗(ＷＬＳ)推定器である、請求項１に記載の方法。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、包括的には、３相電力系統に関し、より詳細には、電力系統における周波数及び位相等のパラメーターを推定することに関する。
【背景技術】
【０００２】
電力系統において信頼性を確保することにとって主要な関心事は、電力の品質と、信号パラメーターを所定の範囲内に維持することとである。電力信号の主要なパラメーターは、電圧振幅、位相角、及び基本周波数を含む。
【０００３】
電力系統では、周波数及び位相の推定が、商用送電網と分散発電機との間の同期を維持するのに用いられ、また、発電及び負荷の整合を介した電力系統安定化にも用いられる。しかしながら、電力系統の周波数及び位相は、過渡事象、例えば、電力サージ及び電圧ディップの間に極めて高速に変動する可能性があり、周波数及び位相を十分な精度で追跡することが、非常に困難になる。したがって、高調波歪及び雑音の存在下での高速かつ正確な周波数及び位相の推定は、依然、課題として残されている。
【０００４】
電力系統における周波数及び位相を測定する様々な方法が知られている。ゼロクロス検出手順は、広く用いられ、実施するのが容易である。しかしながら、その手順は、正確な測定を得るのに適切なフィルタリング及び長時間のサンプリングウィンドウを必要とする。
【０００５】
離散的フーリエ変換(ＤＦＴ)及び再帰的ＤＦＴ手順を用いるフィルタリング方法も、長時間のウィンドウを必要とし、このウィンドウは、正確な結果を得るには基本サイクル(周期)の倍数でなければならない。
【０００６】
３相位相ロックループ(３ＰＬＬ)回路が、３相(３ＰＨ)電力信号の周波数を推定及び追跡するのに頻繁に用いられる。３ＰＬＬ回路を構成する数多くの技法が知られている。しかしながら、ほとんどの場合、３ＰＬＬは、基本サイクルの１サイクルを超える過渡応答、又は良くても基本電圧成分の約半分のサイクルの応答時間を示す。周波数を推定するほとんどの手順は、少なくともその系統の１つの相全体の測定に基づいている。したがって、それらの手順は、追跡されている相が１サイクルの僅かの非常に短い電圧ディップ又は過渡事象を受けるとき、十分に動作しない。
【０００７】
クラーク(Clarke)のα、β変換が、３ＰＨ量を１相電力系統における複素量に変換するのに広く用いられ、これは、従来の高調波信号モデルを有する。この１相系統に基づくパラメーター推定器は、全ての３相電圧からの情報を利用することによってより正確である。
【０００８】
３ＰＨ電圧からの全体的な情報は、周波数推定器を設計するときに考慮することができる。しかしながら、高調波を除去するのに基本サイクルの倍数を必要とする適応ＦＩＲフィルターが依然として必要とされる。多重信号分類器(ＭＵＳＩＣ：MUltiple Signal Classifier)及び回転不変性による信号パラメーターの推定技法(ＥＳＰＲＩＴ：Estimation of Signal Parameters via Rotational Invariance Technique)を高調波信号の基本周波数及び位相の高分解能の推定に用いることができる。これについては、非特許文献１及び非特許文献２を参照されたい。分類器についてのＭＵＳＩＣ及び推定技法についてのＥＳＰＲＩＴという用語が本明細書において用いられることに留意されたい。なぜならば、このようにして、これらの技法がこの文献において及び当業者に知られているからである。例えば、特許文献１(1988)を参照されたい。
【０００９】
図１は、電力系統に関する従来技術のパラメーター推定方法１００を示している。ＷＬＳ推定器１２０への入力１１０は、正弦波信号の３相混合体１１０であり、信号パラメーター１３０が得られる。
【００１０】
図２は、別の従来技術のパラメーター推定方法２００を示している。分類器２２０、すなわちＭＵＳＩＣへの入力２１０は、正弦波信号の１相混合体２１０であり、信号パラメーター２３０が得られる。
【先行技術文献】
【特許文献】
【００１１】
【特許文献１】米国特許第４，７５０，１４７号
【非特許文献】
【００１２】
【非特許文献１】Schmidt「Multiple emitter location and signal parameter estimation」(IEEE Trans. Antennas Propag., vol. 34, no. 3, pp. 276-280, Mar. 1986)
【非特許文献２】Roy及びKailath「Esprit-estimation of signal parameters via rotational invariance techniques」(IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 7, pp. 984-995, Jul. 1989)
【発明の概要】
【００１３】
この発明の実施の形態は、高調波歪を有する平衡３相電力系統(３ＰＨ)及び非平衡３相電力系統における基本周波数及び位相を推定する方法を提供する。利点として、本方法は、４分の１サイクルのみの信号サンプル、例えば、５０Ｈｚの公称周波数で動作する送電網の場合、４ｋＨｚのサンプリング周波数を用いると２０個のサンプルに対して作用する。
【００１４】
具体的には、３相電力系統のモデルが、クラーク変換を用いて、雑音によって劣化した単一相高調波信号モデルに変換される。
【００１５】
新規な加重最小二乗(ＷＬＳ)パラメーター推定器は、信号の高調波構造を用いることに基づいている。
【００１６】
実際の電力系統には高調波が存在するため、本方法は、高調波周波数の推定値を利用することによってＷＬＳ推定器用に初期推定を精緻化するように反復的なものにされる。
【００１７】
本方法は、４分の１サイクルのサンプルのみを用いて、従来の推定器、すなわち、ＭＵＳＩＣ、ＥＳＰＲＩＴ、又はＷＬＳ推定器単独のものよりも性能が優れている。
【図面の簡単な説明】
【００１８】
【図１】加重最小二乗を用いた従来技術の信号パラメーター推定方法の流れ図である。
【図２】ＭＵＳＩＣを用いた従来技術の信号パラメーター推定方法の流れ図である。
【図３】この発明の実施の形態による信号パラメーター推定方法の流れ図である。
【図４】この発明の実施の形態による信号パラメーター推定方法の詳細の流れ図である。
【図５】この発明の実施の形態による信号パラメーター推定方法の詳細の流れ図である。
【図６】この発明の実施の形態による信号パラメーター推定システムの概略図である。
【図７】信号対雑音比の関数としてのクラメールラオ(Cramer-Rao)の下限のグラフである。
【発明を実施するための形態】
【００１９】
図３は、３相電力系統の周波数及び位相等のパラメーターを推定する方法を示している。この方法への入力は、正弦波信号の３相(３ＰＨ)混合体のサンプル３１０である。利点として、サンプルは、信号の基本サイクル(周期)の４分の１よりも短い信号から取得される。サンプルをより多く取得するほど、計算複雑度及びレイテンシーの増加を犠牲にして、周波数推定の精度は確実に改善される。
【００２０】
分類器、例えば、低複雑度の高速累乗根(fast-root：高速ルート)多重信号分類器(ＦＲＭ)がサンプルに適用されて(３２０)、初期パラメーター３２１が得られる。次に、加重最小二乗(ＷＬＳ)推定器が初期パラメーターに適用されて(３３０)、パラメーターが精緻化され、最終パラメーター３４０が得られる。
【００２１】
この方法は、反復的であり、当該技術分野において知られているように、メモリ及び入出力インターフェースに接続されたプロセッサ３００において実行することができる。
【００２２】
図４は、上記方法ステップをより詳細に示している。分類器４１０が信号３１０に適用される。最大振幅を有する高調波である１次高調波が、ステップ３３０用に抽出される。この抽出された信号成分は、ＭＵＳＩＣの次の適用の前に信号から取り除かれる(４３０)。
【００２３】
図５は、ステップ４１０の詳細を示している。サンプルから共分散行列が求められる(５１０)。次に、雑音部分空間が、高速データ投影方法を用いて求められる(５２０)。この雑音部分空間の多項式係数が、シューア(Schur)アルゴリズムを用いて求められる(５３０)。これについては、Sayed他「A survey of spectral factorization method」(Numerical Linear Algebra with Applications, vol. 8, pp. 467-496, Jul. 2001)を参照されたい。多項式の累乗根５４０がステップ４２０に渡される。
【００２４】
図６は、パラメーター推定を実施する等価システムを示している。この図に示す変数は、以下で詳細に説明される。ここで、ｖ＝［ｖ(１)…ｖ(Ｎ)］Ｔは、入力サンプルベクトルであり、
【００２５】
【数１】
2014153354.tif 000003
【００２６】
は、出力パラメーター推定値である。推定されたｋ次高調波成分は、
【００２７】
【数２】
2014153354.tif 000004
【００２８】
であり、ここで、
【００２９】
【数３】
2014153354.tif 000005
【００３０】
である。
【００３１】
３ＰＨ電力系統における信号モデル
Ｋ−１個の高調波歪を有する平衡３相(３ＰＨ)電力系統では、信号は、以下のようにモデル化することができる。
【００３２】
【数４】
2014153354.tif 000006
【００３３】
式中、
【００３４】
【数５】
2014153354.tif 000007
【００３５】
であり、
【００３６】
【数６】
2014153354.tif 000008
【００３７】
であり、ｎは、サンプルの時刻であり、ｎ＝１，２，…，Ｎであり、ω０は、０＜|ω０|＜πを有する基本周波数であり、位相φ∈［０，２π)であり、Ｖｋ＞０は、ｋ次高調波の振幅である。雑音ベクトルεａｂｃ(ｎ)＝［εａ(ｎ)  εｂ(ｎ)  εｃ(ｎ)］Ｔは、共分散行列σ２Ｉを有するゼロ平均ガウス確率ベクトルである。
【００３８】
ｉ＝ａ，ｂ，ｃを有する信号サンプル
【００３９】
【数７】
2014153354.tif 000009
【００４０】
について、この方法は、未知のパラメーター、すなわち、信号の基本周波数ω０及び位相φを推定する。この発明は、非平衡３相系統にも適用可能である。非平衡系統における電圧ｖａ(ｎ)、ｖｂ(ｎ)及びｖｃ(ｎ)の処理は、平衡３相系統について与えられた例と何ら異なるところはない。
【００４１】
単一相系統における高調波信号モデル
クラーク変換を式(１)に適用することによって、以下のように、αβ固定基準フレーム(stationary reference frame)において対応する信号が得られる。
【００４２】
【数８】
2014153354.tif 000010
【００４３】
ここで、変換行列は、以下のとおりである。
【００４４】
【数９】
2014153354.tif 000011
【００４５】
次に、上記変換から、複素電圧を以下のように形成することができる。
【００４６】
【数１０】
2014153354.tif 000012
【００４７】
θ＝ω０ｎ＋φを３ＰＨ系統の位相角と定義した場合、電圧は、以下のように表すことができる。
【００４８】
【数１１】
2014153354.tif 000013
【００４９】
式中、ε(ｎ)＝εα(ｎ)＋ｊεβ(ｎ)であり、
【００５０】
【数１２】
2014153354.tif 000014
【００５１】
は、共分散行列
【００５２】
【数１３】
2014153354.tif 000015
【００５３】
を有するゼロ平均ガウス確率ベクトルである。
【００５４】
クラーク変換を適用する利益は以下のとおりである。
１)３つの高調波成分の倍数は打ち消され、信号は、
【００５５】
【数１４】
2014153354.tif 000016
【００５６】
個の高調波しか含まない。これによって、パラメーター推定に関する測定量を削減することができる。
２)各高調波の指数符号は、式(５)において交互に変わる。換言すれば、隣接する高調波からの２つの未知のパラメーター間の距離は増加する。したがって、推定精度を改善することができる。
３)任意の位相ディップ、過渡事象、又は小さな停電を受けるとき、各単一相からの測定値を合成して、より信頼できる電力系統を提供することができる。
【００５７】
したがって、電力系統における基本周波数及び位相の推定に関するこの等価高調波信号モデルは、以下の式となる。
【００５８】
【数１５】
2014153354.tif 000017
【００５９】
式中、
【００６０】
【数１６】
2014153354.tif 000018
【００６１】
であり、
【００６２】
【数１７】
2014153354.tif 000019
【００６３】
であり、
【００６４】
【数１８】
2014153354.tif 000020
【００６５】
である。合計Ｋ個の高調波がこのモデルに関与し、ｋ次高調波の次数はｌｋである。
【００６６】
ｖ＝［ｖ(１)…ｖ(Ｎ)］Ｔ及びε∈ＣＮ×１を
【００６７】
【数１９】
2014153354.tif 000021
【００６８】
から同様に形成することにする。その場合、式(６)は、以下のように書き換えることができる。
【００６９】
【数２０】
2014153354.tif 000022
【００７０】
式中、
【００７１】
【数２１】
2014153354.tif 000023
【００７２】
であり、Ｈ(ω０)∈ＣＮ×Ｋは、以下の式によって与えられる第ｋ列を有するバンデルモンド(Vandermonde)行列を示す。
【００７３】
【数２２】
2014153354.tif 000024
【００７４】
信号サンプル
【００７５】
【数２３】
2014153354.tif 000025
【００７６】
からω０、φ及びＡｋが推定される。
【００７７】
ＷＬＳ周波数位相推定器
従来のＷＬＳ推定器が、本発明者らの信号モデル用に拡張される。従来のＷＬＳ推定器については、Stoica他「Computationally efficient parameter estimation for harmonic sinusoidal signals」(Signal Process, vol. 80, pp. 1937-1944, 2000)を参照。ＷＬＳは、高調波の構造を無視し、ＭＵＳＩＣ及びＥＳＰＲＩＴ等の従来の正弦波パラメーター推定器を用いて、初期推定値３２１、すなわち、ωｋについての
【００７８】
【数２４】
2014153354.tif 000026
、φｋについての
【００７９】
【数２５】
2014153354.tif 000027
【００８０】
及びＡｋについての
【００８１】
【数２６】
2014153354.tif 000028
【００８２】
を得る。次に、ＷＬＳ技法は、電力系統における高調波の構造、すなわち、ωｋ＝ｌｋω０及びφｋ＝ｌｋφを用いて、初期推定を精緻化し(３３０)、最終パラメーターを生成する。初期推定器３２０からの振幅推定値とＷＬＳからの振幅推定値とは同等であることに留意されたい。なぜならば、これらの振幅は、高調波構造に従わず、ＷＬＳを適用することは、Ａｋの推定値を精緻化することにならないからである。
【００８３】
具体的には、
【００８４】
【数２７】
2014153354.tif 000029
【００８５】
及びη＝［φ，ω０］Ｔとする。その場合、ランク２の行列
【００８６】
【数２８】
2014153354.tif 000030
【００８７】
が存在する。ここで、
【００８８】
【数２９】
2014153354.tif 000031
【００８９】
である。その結果、
【００９０】
【数３０】
2014153354.tif 000032
【００９１】
となる。
【００９２】
【数３１】
2014153354.tif 000033
【００９３】
をＭＵＳＩＣからのζの対応する初期推定値とする。ηのＷＬＳ推定器の推定値
【００９４】
【数３２】
2014153354.tif 000034
【００９５】
は、
【００９６】
【数３３】
2014153354.tif 000035
【００９７】
によって与えられる。
【００９８】
高調波構造を考慮しないＣＲＢ行列から得られる加重行列Ｗは、以下のとおりである。
【００９９】
【数３４】
2014153354.tif 000036
【０１００】
次に、式(９)におけるコスト関数が、以下のように書き換えられる。
【０１０１】
【数３５】
2014153354.tif 000037
【０１０２】
式(９)の解は、次のように表すことができる。式(１１)をφに関して微分すると、以下の式が得られる。
【０１０３】
【数３６】
2014153354.tif 000038
【０１０４】
上記式を０に等しくすると、その結果、φの推定値は以下のようになる。
【０１０５】
【数３７】
2014153354.tif 000039
【０１０６】
ここで、式の右辺は、未知のパラメーターω０に依存する。上記式(１３)を式(１１)に代入し、ω０に関して微分すると、ω０の以下のＷＬＳ推定値が得られる。
【０１０７】
【数３８】
2014153354.tif 000040
【０１０８】
上記
【０１０９】
【数３９】
2014153354.tif 000041
【０１１０】
を用いて、式(１３)内のω０と置き換えると、以下のようにφのＷＬＳ推定値を得ることができる。
【０１１１】
【数４０】
2014153354.tif 000042
【０１１２】
したがって、式(１４)及び式(１５)は、周波数パラメーター及び位相パラメーターのＷＬＳ推定値の閉形式を提供する。
【０１１３】
【数４１】
2014153354.tif 000043
【０１１４】
及び
【０１１５】
【数４２】
2014153354.tif 000044
【０１１６】
は、高調波構造を用いることによる、それぞれ
【０１１７】
【数４３】
2014153354.tif 000045
【０１１８】
及び
【０１１９】
【数４４】
2014153354.tif 000046
【０１２０】
にわたる加重線形回帰であるので、ＷＬＳ推定器は、高調波歪から有用な情報を抽出することができる。性能は、多数のサンプルに対してＣＲＢに非常に近い。ＣＲＢの導出については付録を参照されたい。
【０１２１】
サンプルが限られた電力系統用の改良型ＷＬＳ推定器
先ず、ＷＬＳ推定器が電力信号における非常に少数のサンプルに適用されるとき、例えば、４分の１サイクルのサンプルのみが推定に利用可能であるときのＷＬＳ推定器のための本発明者らの方法を説明する。次に、反復ＭＵＳＩＣに基づく本発明者らの改良型ＷＬＳ(ＩＷＬＳ)推定器を説明する。
【０１２２】
３ＰＨ電力系統におけるＷＬＳの不都合点
従来のＷＬＳ推定器は、非主要高調波が１次高調波に匹敵する振幅を有するときは十分に動作する。あいにく、３ＰＨ電力系統における高調波歪が非常に小さい場合、従来のＷＬＳ推定器は、少数のサンプルについてその利点を失う。その理由は、サンプルが限られている場合、弱い電力を有する非主要高調波の未知のパラメーターについての推定精度が保証されないからである。したがって、従来のＷＬＳ推定器において用いられる信頼できない情報によって、推定性能は減少する。
【０１２３】
したがって、本発明者らは、サンプルが非常に少数であっても、ＷＬＳ推定器の信頼できる初期推定を提供することができる方法を提供する。
【０１２４】
ＩＷＬＳパラメーター推定
初期パラメーター３２１用の推定器として、ＭＵＳＩＣが用いられる。ＭＵＳＩＣは、１次高調波の未知のパラメーターを追跡する(４２０)。これは、本発明者らの反復的方法にとって所望の特性である。その上、ＭＵＳＩＣは、一般に、サンプル数Ｎに対するＭの選択の影響を受けやすい。ここで、Ｍは、ＭＵＳＩＣにおいて用いられるサブベクトルの長さである。これは、平均する際に多くのサブベクトルを有することと、その一方で、高調波の十分な次元を保持することとの間の固有のトレードオフである。Ｍ＝４Ｎ／５であるとき、ＭＵＳＩＣは、１次高調波の未知のパラメーターに関して良好な推定を常に提供することがわかっており、したがって、特にサンプル数が少ないときは、Ｍ＝４Ｎ／５が好ましい。
【０１２５】
本発明者らの推定方法は、ＭＵＳＩＣを用いた初期推定と、ＷＬＳを用いた精緻化推定との２ステップ手順を用いる。サンプルが少ない場合、小さな振幅を有する高調波の未知のパラメーターの初期推定は不正確であるので、反復ＷＬＳ(ＩＷＬＳ)推定器を提供して、初期推定が精緻化される。
【０１２６】
具体的には、サンプルベクトルｖに関して、ＭＵＳＩＣを用いて、Ｋ個の全ての高調波の初期パラメーター
【０１２７】
【数４５】
2014153354.tif 000047
【０１２８】
を推定することができる。しかしながら、最大振幅を有する高調波、すなわち、１次高調波に対応する最も正確で信頼できる推定
【０１２９】
【数４６】
2014153354.tif 000048
【０１３０】
のみが保持される。
【０１３１】
【数４７】
2014153354.tif 000049
【０１３２】
を用いると、１次高調波
【０１３３】
【数４８】
2014153354.tif 000050
【０１３４】
６０１を近似的に構成することができる。したがって、この推定された高調波成分が、サンプルベクトルから取り除かれ、次いで、ＭＵＳＩＣが、この時点でＫ−１個の高調波しか有しないこのベクトルに対して再び適用され、Ｋ−１個の残りの１次高調波に対応する推定を収集することによって２次高調波の推定
【０１３５】
【数４９】
2014153354.tif 000051
【０１３６】
が得られる。Ｋ個のステップを用いて反復が続けられる。精緻化推定フェーズ用にＫ個の全ての初期パラメーター
【０１３７】
【数５０】
2014153354.tif 000052
【０１３８】
を収集することができる。
【０１３９】
図７は、独立した周波数ω１、ω２、ω３を有する信号に関与する３つの高調波を有するＣＲＢを周波数の関数として示している。対応する振幅は、それぞれＡ１＝１、Ａ２＝０．０６、Ａ３＝０．０５である。サンプル数は２０であり、ＣＲＢは、
【０１４０】
【数５１】
2014153354.tif 000053
【０１４１】
に関して求められる。ここで、σ２は雑音の電力であり、加重行列Ｗは式(１０)に与えられている。より大きな振幅を有する高調波ほど、その周波数推定値についてより低いＣＲＢを有する。したがって、３つの周波数推定値の精度(ａｃｃ)又は信頼性は、以下の式を満たす。
【０１４２】
【数５２】
2014153354.tif 000054
【０１４３】
式中、
【０１４４】
【数５３】
2014153354.tif 000055
【０１４５】
は、
【０１４６】
【数５４】
2014153354.tif 000056
【０１４７】
の精度を示す。
【０１４８】
式(１４)における周波数推定器及び式(１５)における位相推定器は、それぞれ
【０１４９】
【数５５】
2014153354.tif 000057
【０１５０】
及び
【０１５１】
【数５６】
2014153354.tif 000058
【０１５２】
にわたる加重線形結合であることを想起されたい。しかしながら、４分の１サイクルのサンプルのみ、すなわち低ＳＮＲを用いて、Ｋ個のパラメーターの全てとは限らないが、それらのパラメーターを正確に推定することができる。したがって、最初のｋ個(ｋ＜Ｋ)の高調波の推定値のみを用いて、ω０及びφが推定される。
【０１５３】
換言すれば、Ｋ回の全ての反復を経由するのではなく、ｋ回の反復しか実行されない。サンプル数が少ないとき、初期推定においてのみ反復が実施されることに留意されたい。したがって、反復によって導入される計算複雑度は許容可能である。
【０１５４】
計算複雑度の解析及び高速実施
本発明者らの方法は、ｋ回反復するＭＵＳＩＣを用いた初期パラメーター推定と、ＷＬＳ推定器を用いた精緻化推定との２つの主なステップを含む。ＷＬＳ推定器は、基本的には、Ｏ(１)の計算複雑度を有する初期推定に対する線形演算であり、したがって、ＩＷＬＳ推定器の計算コストは、主として、ＭＵＳＩＣの計算コストによって決まる。
【０１５５】
サンプル共分散行列を求めるのに用いられるベクトルの長さは、Ｍ＝４Ｎ／５である。この計算負担は、主として、２３Ｍ３フロップス(flops)を有するサンプル共分散行列の固有値分解と、２Ｍ−２次多項式の累乗根を求めることとに起因している。この多項式の２Ｍ−２個の累乗根を見つけることは、
【０１５６】
【数５７】
2014153354.tif 000059
【０１５７】
フロップスを必要とするその対応するコンパニオン行列の２Ｍ−２個の固有値を見つけることと等価である。したがって、ＭＵＳＩＣの支配的な複雑度は、以下となる。
【０１５８】
【数５８】
2014153354.tif 000060
【０１５９】
適応的な状況では、提案した手順を再帰的に適用することによってモデルパラメーターを追跡することは、多くの時間を要する。したがって、ＩＷＬＳ推定器の高速実施の明らかな必要性が存在する。高速累乗根ＭＵＳＩＣ(ＦＲＭ)推定は、多項式の累乗根を見つける複雑度を大幅に削減することによって、計算負担がより少ない効率的なアルゴリズムを提供する。これについては、Zhuang他「Fast root-music for arbitrary arrays」(Elec. Letters, vol. 46, no. 2, pp. 174-176, Jan. 2010)を参照されたい。
【０１６０】
２Ｍ−２の次数を有する多項式ｆ(ｚ)は、ｃを正の定数として、ｆ(ｚ)＝ｃｆ１(ｚ)ｆ１＊(１／ｚ＊)に因数分解することができるので、ｆ(ｚ)の累乗根は、反数の共役対で現れる。半分の累乗根、すなわち、Ｍ−１の次数を有するｆ１(ｚ)の累乗根を見つけることで十分である。
【０１６１】
高速累乗根ＭＵＳＩＣに基づく推定器
高速累乗根ＭＵＳＩＣは、以下のように要約することができる。
１．)以下のサンプル共分散行列Ｒを求める。
【０１６２】
【数５９】
2014153354.tif 000061
【０１６３】
式中、ｘｊ∈ＣＭ×１は、ｖをトランケートしたものである。
２．)Ｒの固有値分解を行って雑音部分空間Ｅを取得し、Ａ(ｍ，ｎ)を第(ｍ，ｎ)エントリーとして有するＡ＝ＥＥＨを構成する。
３．)２Ｍ−２次多項式ｆ(ｚ)の部分係数を、
【０１６４】
【数６０】
2014153354.tif 000062
【０１６５】
によって求める。ここで、ｉ＝−(Ｍ−１)，−(Ｍ−２)，…，０である。
４．)以下のシューアアルゴリズムを用いて、ｆ１(ｚ)のＭ個の係数を見つける。
(ａ)ｂｉに関する初期行列
【０１６６】
【数６１】
2014153354.tif 000063
【０１６７】
を構成する。
(ｂ)ｉ＝１，２，…，Ｉについて、以下のステップを反復する。
ｉ．Ｂｉ＝Ｂｉ−１Ｕｉ、ここで、
【０１６８】
【数６２】
2014153354.tif 000064
【０１６９】
は、以下のように定義される。
【０１７０】
【数６３】
2014153354.tif 000065
【０１７１】
式中、γ＝Ｂｉ−１(１，２)／Ｂｉ−１(１，１)である。
ｉｉ．Ｂｉの第１列を変更せずに維持したまま、第２列を１要素分だけシフトアップする。
(ｃ)ｆ１(ｚ)の係数は、ｂ１，ｉ＊によって与えられる。ここで、ｂ１，ｉは、Ｂｉの第１列である。
５．)ｆ１(ｚ)のＭ−１個の累乗根を見つける。これは、その対応するコンパニオン行列ＭのＭ−１個の固有値を見つけることと等価である。
【０１７２】
上記手順において、ステップ２は、固有値分解によって雑音部分空間ベクトルを求めるのに２３Ｍ３フロップスを必要とする。この多くの時間を要する演算は、再帰的に更新されるＥによって回避することができる。
【０１７３】
高速投影方法
したがって、以下の高速投影方法が提供される。
【０１７４】
前の時刻：Ｅ(ｎ−１)、新たなデータ：ｖ(ｎ)、ｎ＞Ｎ
【０１７５】
適用：
【０１７６】
【数６４】
2014153354.tif 000066
【０１７７】
この発明の実施の形態は、４分の１サイクル以下のサンプルを用いて、３相(３ＰＨ)電力系統における周波数及び位相の推定問題を解決する。高調波歪を有する３ＰＨ電力系統のこのモデルは、雑音によって劣化した単一相高調波信号モデルに変換される。
【０１７８】
新規なＷＬＳ周波数推定器及びＷＬＳ位相推定器が提供される。非主要高調波の低電圧特性及び少数のサンプルに起因して、反復手順を用いてＷＬＳ推定器の初期パラメーターを精緻化する改良型ＷＬＳ推定器が提供される。
【０１７９】
本明細書において与えた例は、平衡３ＰＨ電力系統についてのものであるが、この発明は、この発明に従って説明したアルゴリズムを変更することなく、電圧及び位相が非平衡である３ＰＨ電力系統等の他の３ＰＨ電力系統にも用いることができる。
【０１８０】
付録
クラメールラオの下限
式(６)における複素データモデルについて提起されたパラメーター推定問題のＣＲＢを導出する。
【０１８１】
【数６５】
2014153354.tif 000067
【０１８２】
を定義するとともに、ａＴ＝［Ａ１，…，ＡＫ］とすることによって、式(７)は以下のように書き換えられる。
【０１８３】
【数６６】
2014153354.tif 000068
【０１８４】
式中、
【０１８５】
【数６７】
2014153354.tif 000069
【０１８６】
であり、その第ｋ列は、
【０１８７】
【数６８】
2014153354.tif 000070
【０１８８】
によって定義される。
【０１８９】
スレピアンバングズ(Slepian-Bangs)の公式を用いることによって、この問題のためのＣＲＢ行列は、以下の式によって与えられる。
【０１９０】
【数６９】
2014153354.tif 000071
【０１９１】
次に、以下のように偏導関数が評価される。
【０１９２】
【数７０】
2014153354.tif 000072
【０１９３】
式中、
【０１９４】
【数７１】
2014153354.tif 000073
【０１９５】
である。ここで、
【０１９６】
【数７２】
2014153354.tif 000074
【０１９７】
及び
【０１９８】
【数７３】
2014153354.tif 000075
【０１９９】
の第ｋ列は、それぞれｃｋ＝ｊｌｋｂｋ及び
【０２００】
【数７４】
2014153354.tif 000076
【０２０１】
によって与えられる。
【０２０２】
それゆえ、以下の
【０２０３】
【数７５】
2014153354.tif 000077
【０２０４】
の結果が得られる。
【図１】
2014153354.tif 000078
【図２】
2014153354.tif 000079
【図３】
2014153354.tif 000080
【図４】
2014153354.tif 000081
【図５】
2014153354.tif 000082
【図６】
2014153354.tif 000083
【図７】
2014153354.tif 000084
【外国語明細書】
2014153354.tif 000085
2014153354.tif 000086
2014153354.tif 000087
2014153354.tif 000088
2014153354.tif 000089
2014153354.tif 000090
2014153354.tif 000091
2014153354.tif 000092
2014153354.tif 000093
2014153354.tif 000094
2014153354.tif 000095
2014153354.tif 000096
2014153354.tif 000097
2014153354.tif 000098
2014153354.tif 000099
2014153354.tif 000100
2014153354.tif 000101
2014153354.tif 000102
2014153354.tif 000103
2014153354.tif 000104
2014153354.tif 000105
2014153354.tif 000106
2014153354.tif 000107
2014153354.tif 000108
2014153354.tif 000109
2014153354.tif 000110
2014153354.tif 000111
2014153354.tif 000112
2014153354.tif 000113
2014153354.tif 000114
2014153354.tif 000115