(57)【要約】【課題】  系統連系される風力発電機の出力電力振動の持続時間が短い場合であっても風力発電機の振動周波数を正確に推定でき風力発電機の運転状態判別方法を提供することである。【解決手段】  測定された振動分を含む風力発電機の出力電力信号を所定の周期でサンプリングし、理論振動周波数を中心とする移動平均法による帯域フィルタを通して出力電力信号の時系列データを入力して集積し、集積した出力電力信号の時系列データを予め用意した減衰振動の自己回帰モデルに適用し、時系列データを適用した自己回帰モデルの係数から出力電力信号の振動周波数を算出し、予め定めた条件を満たす振動周波数を振動周波数推定値として抽出し、振動周波数推定値に基づいて風力発電機の回転数を特定する。【選択図】  図１
2006246604.tif 000002
【特許請求の範囲】
【請求項１】
測定された振動分を含む極数切換型の風力発電機の出力電力信号を所定の周期でサンプリングして出力電力信号の時系列データを入力し、前記出力電力信号の振動として理論的に推定される理論振動周波数を中心とする移動平均法による帯域フィルタに前記時系列データを通し、前記帯域フィルタにより低周波分や高周波分が除去された出力電力信号の時系列データを集積し、集積した出力電力信号の時系列データを予め用意した減衰振動の自己回帰モデルに適用し、時系列データを適用した自己回帰モデルの係数から出力電力信号の振動周波数を算出し、予め定めた条件に基づいて得られた振動周波数の妥当性を判定し、予め定めた条件を満たす振動周波数を振動周波数推定値として抽出し、振動周波数推定値に基づいて風力発電機の回転数を特定することを特徴とする風力発電機の運転状態判別方法。
【請求項２】
測定された振動分を含む可変速型の風力発電機の出力電力信号を所定の周期でサンプリングして出力電力信号の時系列データを入力し、前記出力電力信号の振動が存在すると推定される存在領域をｎ等分した各領域の中心周波数を前記出力電力信号の振動として理論的に推定される理論振動数とし、その各理論振動数を中心とする移動平均法による帯域フィルタに前記時系列データを通し、各理論振動数毎に前記帯域フィルタにより低周波分や高周波分が除去された出力電力信号の時系列データを集積し、集積した各理論振動数毎の出力電力信号の時系列データを予め用意した減衰振動の自己回帰モデルに適用し、各理論振動数毎の時系列データを適用した自己回帰モデルの係数から出力電力信号の振動周波数を算出し、予め定めた条件に基づいて得られた各振動周波数の妥当性を判定し、予め定めた条件を満たす各振動周波数を各振動周波数推定値として抽出し、各振動周波数推定値に基づいて風力発電機の回転数を特定することを特徴とする風力発電機の運転状態判別方法。
【請求項３】
前記測定された振動分を含む風力発電機の出力電力信号を、前記帯域フィルタに２回通すことを特徴とする請求項１または２記載の風力発電機の運転状態判別方法。
【請求項４】
前記帯域フィルタを通過した風力発電機の出力電力信号のうち、所定の閾値以下の区間の出力電力信号は、除外することを特徴とする請求項１ないし３のいずれか一記載の風力発電機の運転状態判別方法。
【請求項５】
前記自己回帰モデルの係数のうち、振動を示す正弦項を含む係数の前記正弦項の絶対値が１以下の条件で、振動周波数の妥当性を判定することを特徴とする請求項１ないし４のいずれか一記載の風力発電機の運転状態判別方法。
【請求項６】
前記自己回帰モデルの係数のうち、減衰振動を示す時定数を含む係数の前記時定数の値が正である条件で、振動周波数の妥当性を判定することを特徴とする請求項１ないし５のいずれか一記載の風力発電機の運転状態判別方法。
【請求項７】
前記自己回帰モデルの係数のうち、減衰振動を示す時定数を含む係数の前記時定数の範囲は前記時系列データに基づいて推定した値より大なる条件で、振動周波数の妥当性を判定することを特徴とする請求項６記載の風力発電機の運転状態判別方法。
【請求項８】
前記自己回帰モデルの係数のうち、減衰振動を示す時定数を含む係数の範囲は、前記時系列データに基づいて推定した値以上で、かつ前記時定数を無限大とした場合の値に余裕値を加味した値以下の条件で、振動周波数の妥当性を判定することを特徴とする請求項６記載の風力発電機の運転状態判別方法。
【請求項９】
前記自己回帰モデルの係数のうち、減衰振動を示す時定数を含む係数の変動分が予め定めた閾値より小なる条件で、振動周波数の妥当性を判定することを特徴とする請求項１ないし８のいずれか一記載の風力発電機の運転状態判別方法。
【請求項１０】
振動周波数の妥当性の判定として、推定領域での振動周波数の平均値の±１５％以内の振動周波数を対象とすることを特徴とする請求項１ないし請求項９のいずれか一記載の風力発電機の運転状態判別方法。
【請求項１１】
振動周波数の妥当性の判定として、推定領域で振動周波数の一定継続している時間長が所定の時間長以下である振動周波数は除外することを特徴とする請求項１ないし請求項１０のいずれか一記載の風力発電機の運転状態判別方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、風力発電機が系統連系される系統連系システムにおける風力発電機の運転状態判別方法に関する。
【背景技術】
【０００２】
風力発電機として広く使用される交流発電機は、極数切換型誘導発電機と可変速型発電機との2つのタイプに大別される。極数切換型誘導発電機は極数を切り替えて回転数を切り替えるものであり、通常2種類の回転数を持っている。極数切換型誘導発電機は、極数切換時の突入電流による急激な電圧低下などの電力品質問題などがあり、電圧低下量は切換後の極数により異なるため、極数状態を示す回転数の把握が必要となる。
【０００３】
一方、可変速型発電機は、風速にあわせ回転数を変える超同期セルビウス誘導発電機や、系統と周波数変換器を通して系統連系され交流励磁が与えられて同期発電機を可変速運転するものである。可変速型発電機では、風速に会わせて回転数を変化させる複雑な制御を行うため、予めシミュレーションを行い電力品質上の問題の有無などを検討する必要がある。
【０００４】
すなわち、極数切換型誘導発電機と可変速型発電機とのいずれのタイプであっても、その基礎となるモデル化には発電機実測データが必要となり、その基礎データの一つに回転数がある。回転数は、風力発電機が連系される電力会社では直接測定できず、得られるデータは系統連系システムとの責任分界点の電力量計での二相電力計法での電圧及び電流しかない。
【０００５】
風力発電機の回転数を得るための手段として、風車のブレード回転による風力発電機の出力電力Ｐ（ｔ）の脈動を利用できる。すなわち、風車の塔付近は、塔体の影響で風速が弱まる傾向があり、それに起因して、風力発電機の出力が低下するため、風力発電機の出力電力Ｐ（ｔ）に脈動が生じるというタワーシャドウ効果が生じる。
【０００６】
例えば、極数切換型の風力発電機での２つの脈動の周波数(理論振動周波数)ｆt1、ｆt2は、ブレード回転数をｎi（ｒｐｍ）、ブレード枚数ＮBとすると、次式で与えられる。
【０００７】
［数１］
ｆti＝ｎiＮB／６０  …（１）
（ｉ＝１，２）
これより、風力発電機の出力電力Ｐ（ｔ）の脈動の周波数ｆｔｉが分かれば、風力発電機の回転数を知ることができる。なお、以後の説明を簡潔にするため、理論振動周波数を単にｆｔと表記する。
【０００８】
次に、理論振動周波数ｆｔの推定法を説明する。いま、風力発電機の出力電力Ｐ（ｔ）は（２）式に示すように正弦波であるとする。Ｐｍは出力電力Ｐ（ｔ）の波高値、θ０は位相である。
【０００９】
［数２］
Ｐ（ｔ）＝Ｐmcos（２πｆtｔ＋θ0）…（２）
上式の風力発電機出力Ｐ（ｔ）に対して、時刻ｔを中心に時間領域ｔ＋Ｔ／２〜ｔ−Ｔ／２の両端で方形波窓を適用しフーリェ変換を施す。この場合、積分時間Ｔは（３）式に示す関係にあるとする。
【００１０】
［数３］
２πｆtＴ＝２πＮ    …（３）
（Ｎは振動周期数）
そうすると、フーリェ変換値の絶対値｜χ（ｔ，ω）｜は、（４）式で示される。
【数４】
2006246604.tif 000003
【００１１】
ここで、ω＝ωｔ＋δω、ωｔ＝２πｆｔ、θ（ｔ）＝ωｔ＋θ0、δω：角周波数誤差である。
【００１２】
次に、（４）式のcos(δωＴ)項をマクローリン展開し４次項まで近似し、さらにδωは十分小とし、｛δω／（２ωｔ＋δω）｝2項は無視し、２ωt＋δω≒２ωtと近似すればフーリェ変換値の絶対値｜χ（ｔ，ω）｜は次の（５）式のように近似される。
【数５】
2006246604.tif 000004
【００１３】
さらに、（５）式の右辺の第４項を十分小として無視すれば、次の（６）式に示
す近似式が得られる。
【数６】
2006246604.tif 000005
【００１４】
この（６）式より根号内の第２項までを考察すると、δω＝０の場合にフーリェ変換値の絶対値｜χ（ｔ，ω）｜は最大となるので、このωを求めるものとする。なお、第３項までを考慮する場合には角周波数誤差δωに関する極値条件より、フーリェ変換値の絶対値｜χ（ｔ，ω）｜を最大とする角周波数誤差δωは次の（７）式の関係で示される。
【数７】
2006246604.tif 000006
【００１５】
（７）式から分かるように、角周波数誤差δωは時間関数であり、次の特性を持つことが分かる。
【００１６】
（ａ）角周波数誤差δωは振動周期数Ｎの２乗Ｎ2に反比例するので、フーリェ変換時に窓長である積分時間Ｔを大きくすれば急激に減少する。
【００１７】
（ｂ）cos(２ωtｔ＋２θ0)に比例するので、フーリェ変換すると、角周波数誤差δωのトレンドはｆtの２倍調波で振動する。
【００１８】
（７）式より振動周期数Ｎによる振動誤差ε（ε＝δω／ωt）の最大値は、cos(２ωtｔ＋２θ0)＝１のときであり、下記の表１のように与えられる。
【表１】
2006246604.tif 000007
【００１９】
表１に示すように、振動周期数ＮがＮ=１では振動誤差εが１５．２％と大きく、実用性に乏しいことが分かる。このため、振動誤差εの抑制を考慮した振動周期数ＮがＮ＝３が実用面から採用されてきた。
【発明の開示】
【発明が解決しようとする課題】
【００２０】
しかしながら、風車のタワーシャドウ効果による風力発電機の出力電力Ｐ（ｔ）の振動は３周期も持続せず、振動周波数の推定ができない場合もある。このため、振動周期数Ｎが２以下の場合でも推定が可能な方法が望まれる。また、振動周期数ＮがＮ＝３の場合には、振動周波数が３周期に亘って平均化される。
【００２１】
そのため、データ上、振動が持続しない区間(以後、非良質区間と呼ぶ)であっても、一見して振動が発生しているような結果を与える場合もある。従って、データの非良質区間を自動的に判断可能な推定方法も望まれている。
【００２２】
そこで、本出願人の発明者らは、系統連系される風力発電機の出力電力の振動周波数を短い振動周期数で正確に推定できる風力発電機の運転状態判別方法を発明し、特願２００４−１６８０２１号により特許出願した。
【００２３】
すなわち、風力発電機の出力電力信号の理論振動周波数ｆｔを中心とする移動平均法による帯域フィルタに通し、積分時間領域(１／ｆt)秒のフーリェ変換を施してそのフーリェ変換値の絶対値が最大となる周波数を推定振動周波数信号ｆ1として求め、順次２ｊ倍調波の誤差振動を含む推定振動周波数ｆk（ｋ＝１／２i）を移動平均法による高域フィルタに通し、積分時間領域(１／２iｆt)秒のフーリェ変換を施し推定振動周波数信号ｆk+1（ｋ＋１＝１／２i+1）を求め、さらに、推定振動周波数信号ｆk+1に対して積分時間領域(２／２i+1ｆｔ)秒のフーリェ変換を施して推定振動周波数信号ｆm（ｍ＝２／２i+1）を求め、数値デジタル低域フィルタを通して最終の周波数推定値ｆTを決定する。
【００２４】
これにより、風車発電機の出力電力の振動をより正確に抽出し、風力発電機の出力電力の振動が持続する区間をより正確に抽出するとともに、風力発電機の出力電力の振動の周波数推定値をより正確に推定できるようにした。
【００２５】
この特願２００４−１６８０２１号のものでは、風力発電機の出力電力振動がある程度継続した場合には、その振動の周波数推定を正確に推定できるが、出力電力振動の持続時間が短い場合には、振動が持続しない非良質区間と判定してしまい、良質区間が殆どないと判定してしまうことがあることが分かった。
【００２６】
本発明の目的は、系統連系される風力発電機の出力電力振動の持続時間が短い場合であっても風力発電機の振動周波数を正確に推定でき、風力発電機の回転数と出力電力との関係などの運転状態判別を行うことができる風力発電機の運転状態判別方法を提供することである。
【課題を解決するための手段】
【００２７】
請求項１の発明に係わる風力発電機の運転状態判別方法は、測定された振動分を含む極数切換型の風力発電機の出力電力信号を所定の周期でサンプリングして出力電力信号の時系列データを入力し、前記出力電力信号の振動として理論的に推定される理論振動周波数を中心とする移動平均法による帯域フィルタに前記時系列データを通し、前記帯域フィルタにより低周波分や高周波分が除去された出力電力信号の時系列データを集積し、集積した出力電力信号の時系列データを予め用意した減衰振動の自己回帰モデルに適用し、時系列データを適用した自己回帰モデルの係数から出力電力信号の振動周波数を算出し、予め定めた条件に基づいて得られた振動周波数の妥当性を判定し、予め定めた条件を満たす振動周波数を振動周波数推定値として抽出し、振動周波数推定値に基づいて風力発電機の回転数を特定することを特徴とする。
【００２８】
請求項２の発明に係わる風力発電機の運転状態判別方法は、測定された振動分を含む可変速型の風力発電機の出力電力信号を所定の周期でサンプリングして出力電力信号の時系列データを入力し、前記出力電力信号の振動が存在すると推定される存在領域をｎ等分した各領域の中心周波数を前記出力電力信号の振動として理論的に推定される理論振動数とし、その各理論振動数を中心とする移動平均法による帯域フィルタに前記時系列データを通し、各理論振動数毎に前記帯域フィルタにより低周波分や高周波分が除去された出力電力信号の時系列データを集積し、集積した各理論振動数毎の出力電力信号の時系列データを予め用意した減衰振動の自己回帰モデルに適用し、各理論振動数毎の時系列データを適用した自己回帰モデルの係数から出力電力信号の振動周波数を算出し、予め定めた条件に基づいて得られた各振動周波数の妥当性を判定し、予め定めた条件を満たす各振動周波数を各振動周波数推定値として抽出し、各振動周波数推定値に基づいて風力発電機の回転数を特定することを特徴とする。
【００２９】
請求項３の発明に係わる風力発電機の運転状態判別方法は、請求項１または２の発明において、前記測定された振動分を含む風力発電機の出力電力信号を、前記帯域フィルタに２回通すことを特徴とする。
【００３０】
請求項４の発明に係わる風力発電機の運転状態判別方法は、請求項１ないし３のいずれか一の発明において、  前記帯域フィルタを通過した風力発電機の出力電力信号のうち、所定の閾値以下の区間の出力電力信号は、除外することを特徴とする。
【００３１】
請求項５の発明に係わる風力発電機の運転状態判別方法は、請求項１ないし４のいずれか１の発明において、前記自己回帰モデルの係数のうち、振動を示す正弦項を含む係数の前記正弦項の絶対値が１以下の条件で、振動周波数の妥当性を判定することを特徴とする。
【００３２】
請求項６の発明に係わる風力発電機の運転状態判別方法は、請求項１ないし５のいずれか１の発明において、前記自己回帰モデルの係数のうち、減衰振動を示す時定数を含む係数の前記時定数の値が正である条件で、振動周波数の妥当性を判定することを特徴とする。
【００３３】
請求項７の発明に係わる風力発電機の運転状態判別方法は、請求項６の発明において、前記自己回帰モデルの係数のうち、減衰振動を示す時定数を含む係数の前記時定数の範囲は前記時系列データに基づいて推定した値より大なる条件で、振動周波数の妥当性を判定することを特徴とする。
【００３４】
請求項８の発明に係わる風力発電機の運転状態判別方法は、請求項６の発明において、前記自己回帰モデルの係数のうち、減衰振動を示す時定数を含む係数の範囲は、前記時系列データに基づいて推定した値以上で、かつ前記時定数を無限大とした場合の値に余裕値を加味した値以下の条件で、振動周波数の妥当性を判定することを特徴とする。
【００３５】
請求項９の発明に係わる風力発電機の運転状態判別方法は、請求項１ないし８のいずれか１の発明において、前記自己回帰モデルの係数のうち、減衰振動を示す時定数を含む係数の変動分が予め定めた閾値より小なる条件で、振動周波数の妥当性を判定することを特徴とする。
【００３６】
請求項１０の発明に係わる風力発電機の運転状態判別方法は、請求項１ないし９のいずれか１の発明において、振動周波数の妥当性の判定として、推定領域での振動周波数の平均値の±１５％以内の振動周波数を対象とすることを特徴とする。
【００３７】
請求項１１の発明に係わる風力発電機の運転状態判別方法は、請求項１ないし１０のいずれか１の発明において、振動周波数の妥当性の判定として、推定領域で振動周波数の一定継続している時間長が所定の時間長以下である振動周波数は除外することを特徴とする。
【発明の効果】
【００３８】
本発明によれば、出力電力信号の時系列データに基づいて出力電力信号の減衰振動の自己回帰モデルを作成し、自己回帰モデルの係数から出力電力信号の振動周波数を算出するので、短い振動期間の時系列データでも出力電力信号の振動周波数を算出できる。また、自己回帰モデルとして先験的に減衰振動波形モデルを採用するので、分解能は時系列データを得るサンプリング間隔で規定されない。また、振動周波数を算出するための計算量が少なくて済む。さらに、減衰振動波形モデルを仮定するので、実データから得られる風力発電機の時定数との対応から、推定結果である振動周波数の妥当性の判定が容易に行える。
【発明を実施するための最良の形態】
【００３９】
以下、本発明の実施の形態を説明する。図１は、本発明の実施の形態に係わる風力発電機の運転状態判別方法のフローチャートである。この実施の形態では極数切換型の風力発電機を対象としている。極数切換型の風力発電機の出力電力は、極数切換型の風力発電機が連系される系統連系システムの電力量計で測定される。測定された極数切換型の風力発電機の出力電力信号にはタワーシャドウ効果による脈動分が含まれている。
【００４０】
まず、測定された振動分を含む極数切換型の風力発電機の出力電力信号Ｐ（ｔ）を所定のサンプリング周期で読み込み（Ｓ１）、出力電力信号の振動として理論的に推定される理論振動周波数ｆｔを中心とする移動平均法による帯域フィルタに通し、低周波分や高周波分が除去された出力電力信号ＰB（ｔ）を求める（Ｓ２）。そして、所定のサンプリング数の時系列データＰB,n（ｔ）（ｎ＝１、２、…Ｎ）を集積する（Ｓ３）。
【００４１】
次に、出力電力信号の時系列データＰB,n（ｔ）に基づいて出力電力信号の減衰振動の自己回帰モデルを作成する（Ｓ４）。帯域フィルタ適用後の電力ＰB（ｔ）は、一般に振動領域は減衰振動であるので、自己回帰モデルＡＲＭ（Auto Regressive Model）として２次モデルを採用する。すなわち、後述するように、振動を示す正弦項を含む係数と減衰振動を示す時定数を含む係数との２個の係数を用いて、複数個の時系列データＰB,n（ｔ）の関係式を導き自己回帰モデルを作成する。
【００４２】
そして、自己回帰モデルの関係式の２個の係数から出力電力信号ＰＢ（ｔ）の振動周波数ｆを算出し（Ｓ５）、算出された振動周波数ｆの妥当性を予め定めた条件に基づいて判定する（Ｓ６）。振動周波数ｆの妥当性の判定の仕方については後述する。この判定により、予め定めた条件を満たす振動周波数を振動周波数推定値ｆTとして抽出し（Ｓ７）、振動周波数推定値ｆTに基づいて風力発電機の回転数を特定する（Ｓ８）。（１）式のｆtiに振動周波数推定値ｆTを代入して風力発電機の回転数ｎiを特定する。
【００４３】
以下、各ステップＳ１〜Ｓ７の処理内容につき詳細に説明する。まず、帯域フィルタについて説明する。風力発電機の回転数の推定を行うに際し、風力発電機の出力信号Ｐ（ｔ）に含まれるタワーシャドウ効果による振動を抽出するわけであるが、タワーシャドウ効果による振動に低周波振動が重畳することがあり、低周波振動が重畳すると基本波の振動誤差が発生する。そこで、これを抑制するとともに、高周波振動を除去するために帯域フィルタを適用する。すなわち、極数切換型の風力発電機では（１）式で与えられる理論振動周波数は２種類しかないので、これらを通過率が最大の周波数とする帯域フィルタを適用する。
【００４４】
この場合の帯域フィルタとしては、位相が対周波数特性を持たない移動平均法の低域フィルタと高域フィルタとを組み合わせたものを使用する。帯域フィルタはフィルタ中心の周波数を理論振動周波数ｆtとし、かつ、その通過率ＫＢの対周波数（ｆ）特性が次式で与えられるものを採用する。
【数８】
2006246604.tif 000008
【００４５】
また、帯域フィルタは風力発電機の出力電力信号Ｐ（ｔ）に対し、必要に応じて連続2回適用する。これにより、風力発電機の出力電力信号Ｐ（ｔ）に含まれる低周波振動及び高周波振動を除去する。
【００４６】
図２は、本発明の実施の形態で使用する帯域フィルタの通過率ＫＢの対周波数特性図である。図２に示すように、帯域フィルタは中心の周波数を理論振動周波数ｆtとし、通過率ＫBが０．４であるフィルタである。
【００４７】
図３は、測定された風力発電機の出力電力信号Ｐ（ｔ）に連続2回に亘って帯域フィルタを適用した場合の風力発電機の出力電力信号ＰB（ｔ）の対周波数特性図である。図３では、測定された風力発電機の出力電力信号Ｐ（ｔ）と帯域フィルタを適用した後の風力発電機の出力電力信号ＰB（ｔ）とを示しており、時刻４秒の近傍で２周期程度の振動のみを持つＰ（ｔ）とＰB（ｔ）とを比較図示しているが、帯域フィルタにより十分平滑化されていることが分かる。
【００４８】
次に、帯域フィルタを通過した風力発電機の出力電力信号ＰB（ｔ）のうち、所定の閾値ＰB,C以下の区間の出力電力信号ＰB（ｔ）は、必要に応じて除外する。これは、タワーシャドウ効果による振動振幅が小の場合には推定結果の信頼性は低く、出力電力信号ＰB（ｔ）の脈動が小の領域での推定結果を除外することが望ましいからである。このため、除外の指標としての出力電力信号ＰＢ（ｔ）の振幅値に対する閾値ＰB,Cは、以下のように設定する。
【００４９】
まず、発電機は定格電圧ＶNで一定、かつ力率も１であるとすると、電流脈動分により電力脈動分が決定される。従って、Ａ／Ｄ変換器の１ビットに対応した有効電力の跳び幅Δｐは、ＭビットのＡ／Ｄ変換器の入力レンジＶinと、定格電流ＩNに対し電圧換算された入力信号レベルＶIとを用いて表すと、次の（９）式で与えられる。
【数９】
2006246604.tif 000009
【００５０】
ここで、β＝２M-1（ＶI／Ｖin）である。
【００５１】
さらに、帯域フィルタを２回通し帯域フィルタの通過率ＫＢが０．４であることを考慮すれば、出力電力信号ＰB（ｔ）の飛び幅ΔｐBは０．４Δｐになる。これより、図４に示すように、出力電力信号ＰB（ｔ）の閾値ＰB,CはＡ／Ｄ変換器の最小分解能（ΔｐＢ／２）が臨界値となる。実際のトレンドは正弦波ではないので余裕度α（α0≧α≧１）を用いて、閾値ＰB,Cは（１０）式に示すように設定する。
【数１０】
2006246604.tif 000010
【００５２】
次に、出力電力信号ＰＢ（ｔ）の時系列データＰB,n（ｔ）は、サンプリング周期Δｔ毎に得られるので、順次、サンプリング周期Δｔ毎の時系列データＰB,n（ｔ）（ｎ＝１、２、…Ｎ）を集積する。
【００５３】
自己回帰モデルは、帯域フィルタ適用後の電力ＰB（ｔ）の振動領域が減衰振動であることから、（１１）式に示すような２次モデルを採用する。
【数１１】
2006246604.tif 000011
【００５４】
ここで、Ｐは出力電力信号ＰＢ（ｔ）の波高値、ωは出力電力信号ＰB（ｔ）の角周波数、θは出力電力信号ＰB（ｔ）の位相、Ｔは減衰の時定数である。そして、（１１）式を書き換えると（１２）式が得られる。
【数１２】
2006246604.tif 000012
【００５５】
この（１２）式のｔに各々のサンプリング時点の時刻ｎΔｔ（ｎ＝１，２，…Ｎ）を代入すると、各々の時系列データＰB,n（ｔ）（ｎ＝１、２、…Ｎ）となる。いま、（１２）式のｔに時刻（ｎ−２）Δｔを代入すると、時系列データＰB,n-2は（１３）式で示される。
【数１３】
2006246604.tif 000013
【００５６】
そして、時間を進める演算子Ｚ（Ｚ＝ελΔｔ）を導入すると、時系列データＰB,n-1と時系列データＰB,n-2との関係は（１４）式で示される。
【数１４】
2006246604.tif 000014
【００５７】
（１３）式及び（１４）式より、Ｐhελ(n-2)Δt及び（Ｐhελ(n-2)Δt）*を求めると、（１５）式が得られる。
【数１５】
2006246604.tif 000015
【００５８】
一方、時系列データＰB,nと時系列データＰB,n-2との関係は、時間を進める演算子Ｚ（Ｚ＝ελΔｔ）を用いて表すと、（１６）式で示される。
【数１６】
2006246604.tif 000016
【００５９】
（１６）式に（１５）式で求めたＰhελ(n-2)Δt及び（Ｐhελ(n-2)Δt）*を代入すると、（１７）式が得られる。
【数１７】
2006246604.tif 000017
【００６０】
時間を進める演算子Ｚとその複素共役Ｚ*との和（Ｚ＋Ｚ*）をａ1、時間を進める演算子Ｚとその複素共役Ｚ*との積（ＺＺ*）をａ2とおくと、（１８）式が得られる。
【数１８】
2006246604.tif 000018
【００６１】
（１８）式から分かるように、時系列データＰB,nは、１サンプリング前の時系列データＰB,n-1に係数ａ1を乗算したものと、２サンプリング前の時系列データＰB,n-2に係数ａ2 を乗算したものとの和で示され、係数ａ1及び係数ａ2は（１９）式で示される。
【数１９】
2006246604.tif 000019
【００６２】
（１９）式から分かるように、係数ａ1は振動を示す正弦項を含む係数であり、係数ａ2は減衰振動を示す時定数Ｔを含む係数である。この係数ａ1、ａ2は、４個の時系列データから求めることができる。すなわち、３個の時系列データＰB,n、ＰB,n-1、ＰB,n-2からなる（１８）式と、１サンプリング前の３個の時系列データＰB,n-1、ＰB,n-2、ＰB,n-3からなる（１８）式との連立方程式により、係数ａ1、ａ2を求めることができる。つまり、４個の時系列データＰB,n、ＰB,n-1、ＰB,n-2、ＰB,n-3から（１８）式を用いて連立方程式を立てて係数ａ1、ａ2を求めることができる。
【００６３】
本発明の実施の形態では、測定精度を考慮しＮ個の連続データ列を用い最小自乗法により係数ａ1、ａ2を求める。そして、得られた係数ａ1、ａ2より、（１９）式を用いて振動周波数を求める。まず、減衰時定数Ｔは（１９）式の第２式より（２０）式で求められる。
【数２０】
2006246604.tif 000020
【００６４】
また、振動周波数ｆは（１９）式の第１式に第２式を代入し、ω＝２πｆであることから、振動周波数ｆは（２１）式で求められる。
【数２１】
2006246604.tif 000021
【００６５】
図５は、風力発電機の出力電力信号ＰB、振動周波数ｆ及び減衰振動を示す時定数Ｔを含む係数ａ2のトレンド図である。図５に示すように、風力発電機の出力電力信号ＰBは１４６ｓ〜１４９ｓの間で脈動を生じている。一方、算出された振動周波数ｆは、１４５．８５ｓ付近で立ち上がり、１４５．９ｓ〜１４７．８５ｓの間で、０．９Ｈｚ〜１．０Ｈｚの範囲でほぼ一定の振動が継続している。そして、１４７．８５ｓ付近で立ち下がり、１４７．９ｓ〜１４７．４ｓの間で、再び、０．８Ｈｚ〜０．９Ｈｚの範囲でほぼ一定の振動が継続している。さらに、１４７．４ｓ付近で立ち下がり、１４７．６ｓ〜１４８．１５ｓの間で、０．６５Ｈｚ〜０．７５Ｈｚの範囲でもほぼ一定の振動が継続している。１４８．３ｓ以降においては振動周波数ｆは大きく変動しており、ほぼ一定の振動が継続している期間は短い。
【００６６】
この振動周波数ｆの特性から視察で判断すると、１４５．９ｓ〜１４８．１５ｓの間の振動周波数ｆが風力発電機の出力電力信号ＰBのタワーシャドウ効果による脈動の振動周波数であると推定できる。
【００６７】
本発明の実施の形態では、振動周波数ｆを精度良く得るために、測定結果である振動周波数ｆの妥当性の判断基準を設け妥当な振動周波数ｆを判定する。すなわち、自己回帰モデルによりタワーシャドウ効果による振動を減衰振動として取扱うため、実波形との対応が容易であり、得られた振動周波数ｆを妥当性の判定は、良否基準の設定が明確になることが期待できる。
【００６８】
まず、第１の基準として、（２０）式のcos項は絶対値は１以下なので、これを満足しないものは除外する。すなわち、（２２）式の条件を満足しない振動周波数ｆは除外する。
【数２２】
2006246604.tif 000022
【００６９】
次に、タワーシャドウ効果の振動は、風力発電機の制御系が働き制動するため、減衰振動になると思われる。そこで、減衰振動の条件より時定数Ｔは正との基準を導入する。
【００７０】
さらに、実データからの推定時定数を用いての第２の基準として、実データの推定減衰時定数Ｔはブレード回転が機械的なので、通常はＴ≧１秒の範囲に設定される。また、サンプリング周期Δｔは０．０２ｓ程度と設定される。これにより、係数ａ2は（１９）式の第２式より、時定数を無限大とした場合の値は１となる。
【００７１】
ただし、うなり状の波形もあるため係数ａ2は１をはるかに超える場合もある。また、出力急変時には制御系が即応するためＴ＜１の場合もある。これにより、係数ａ2の範囲による基準を次のように設定する。
【００７２】
［数２３］
ａ2,max≧ａ2≧ａ2,min     …（２３）
例えば、ａ2,maxとしては時定数Ｔを無限大とした場合の値１に余裕値、例えば０．０２〜０．０１を加味した値１．０２〜１．０１程度、ａ2,minとしては０．９６〜０．９９程度の値を用いる。これにより、例えば、１４５ｓ〜１４６ｓの間及び１４９ｓ〜１５０ｓの間の振動周波数ｆを排除できる。
【００７３】
次に、振動の減衰が滑らかな場合は係数ａ2はほぼ一定となるので、第３の基準として、係数ａ2の変化分の絶対値が所定の閾値よりも小である領域の振動周波数ｆを採用する。この場合、その前処理として移動平均法により平滑化された係数[ａ2]を求め、その平滑化された係数[ａ2]を用いて、平滑化された係数[ａ2]の変化分の絶対値｜Δ[ａ2,k]｜が所定の閾値[ａ2,c]よりも小であることで判定する。（２４）式はその判定式である。
【００７４】
［数２４］
Δ[ａ2,k]＝｜[ａ2,k]−[ａ2,k-1]｜＜[ａ2,c]    …（２４）
ここで、kは時間ステップ、[ａ2,c]は[ａ2,k]の閾値である。
【００７５】
図６は移動平均法により平滑化された係数[ａ2]及びその平滑化された係数[ａ2]の変化分の絶対値｜Δ[ａ2,k]｜のトレンド図である。係数[ａ2]の変化分の絶対値｜Δ[ａ2,k]｜の閾値[ａ2,k]の値を適切に設定することで、例えば、１４５ｓ〜１４６ｓの間及び１４９ｓ〜１５０ｓの間の振動周波数ｆを排除できる。
【００７６】
さらに、これらの基準が厳しく、得られる推定値の領域が少ないため、視察では妥当と思われる領域も除外されることがある。そこで、その防止策として、推定領域に隣接した領域の振動周波数ｆの値が適切であるか否かを判定する次なる第４の基準を導入する。すなわち、推定領域の振動周波数ｆの平均値ｆavに対し、隣接領域での振動周波数ｆの値が（２５）式の条件を満たす領域まで推定領域を拡大する。これにより、振動周波数の妥当性の判定は、推定領域での振動周波数の平均値の±１５％以内の振動周波数を対象とする。
【００７７】
［数２５］
１．１５ｆav≧ｆ≧０．８５ｆav  …（２５）
また、振動周波数の妥当性の判定として、推定領域で振動周波数ｆの一定継続している時間長が所定の時間長以下である振動周波数ｆは除外する。
【００７８】
図７は、得られた振動周波数推定値ＴTの分布図である。測定結果である振動周波数ｆに対し、妥当性の判断基準を適用して得られた振動周波数推定値ｆTの分布示している。振動周波数推定値ｆTは、１４６．２ｓ〜１４７．０ｓの間に０．８Ｈｚ〜１．０Ｈｚの範囲で分布している。振動周波数ｆの特性から視察で判断した場合とほぼ同じ範囲に分布している。図７の場合には、振動周波数推定値ｆTはほぼ０．９Ｈｚ程度であることが分かる。振動周波数推定値ｆTを（１）式のｆtiに代入して風力発電機の回転数ｎiを算出し特定する。
【００７９】
以上の説明では、極数切換型の風力発電機の場合について説明したが、次に、可変速型の風力発電機について適用する場合について説明する。可変速型の発電機では、そのタワ−シャドウ効果による理論振動周波数ｆtは極数切換型の発電機とは異なり、一般に同期速度に対応する周波数ｆ0を中心に約±３０％の領域で変動する。そこで、出力電力信号の振動が存在すると推定される存在領域、すなわち同期速度に対応する周波数ｆ0を中心に約±３０％の領域をｎ等分し、そのｎ等分した各領域の中心周波数を出力電力信号Ｐ（ｔ）の振動として理論的に推定される理論振動数ｆti（ｉ＝１，２，３，…ｎ）とする。実用的には、例えば３分割し、分割した各領域において、その各理論振動数ｆt1、ｆt2、ｆt3を中心とする図１に示した処理を行う。
【００８０】
図８は極数切換型の風力発電機の出力電力Ｐ（ｔ）の特性及び可変速型の風力発電機の出力電力Ｐ（ｔ）の特性の特性図である。図８（ａ）に示すように、極数切換型の風力発電機の出力電力Ｐ（ｔ）の特性は理論振動数ｆtを中心としてタワ−シャドウ効果による振動を示す特性となるが、可変速型の風力発電機の出力電力Ｐ（ｔ）の特性は、図８（ｂ）に示すように、極数切換型に比しタワ−シャドウ効果による出力脈動量が極めて小である。これは、風力発電機の制御系が出力電力Ｐ（ｔ）の変動を抑制するための制御を行うためである。
【００８１】
このため、可変速型の風力発電機の出力電力Ｐ（ｔ）に対して、図１に示した処理をそのまま適用したとき、例えば、同期速度に対応する周波数ｆ0を中心とする帯域フィルタを２回連続適用したときには、周波数の存在範囲端部で約３５％の減衰があるため、この近傍での振動が見落とされる懸念がある。この端部での減衰に対応するために、同期速度に対応する周波数ｆ0を中心に約±３０％の領域をｎ等分し、そのｎ等分した各領域の中心周波数を出力電力信号Ｐ（ｔ）の振動として理論的に推定される理論振動数ｆti（ｉ＝１，２，３，…ｎ）とする。
【００８２】
例えば、周波数の存在範囲を３等分することで端部減衰が約１０％と抑制される。さらに、各領域の中心周波数を理論振動数に模擬することで、それらの各領域に対して、図１の処理にて周波数推定を行えば、可変速型の風力発電機でも周波数推定が可能である。この場合、擬似的な理論振動数理論振動数ｆti（ｉ＝１，２，３）は次のように定義される。
【数２６】
2006246604.tif 000023
【００８３】
ここで、ｆtmaxはタワーシャドウ効果による周波数の最大値、ｆtminはタワーシャドウ効果による周波数の最小値である。
【００８４】
次に、可変速型の風力発電機の周波数推定においても、帯域フィルタを通過した風力発電機の出力電力信号ＰB（ｔ）のうち、所定の閾値ＰB,C以下の区間の出力電力信号ＰB（ｔ）は、必要に応じて除外する。
【００８５】
図９は所定の閾値ＰB,Cを適用しない場合の周波数推定値ｆTのトレンド図である。図９から分かるように、所定の閾値ＰB,Cを適用しない場合には、１５０秒以降ではほぼ同時刻に大きく異なる推定周波数が存在し、物理的に不合理なことが分かる。その理由は、振動振幅が小の場合には推定結果の信頼性は低く、出力電力信号ＰB（ｔ）の脈動が小の領域での推定結果を除外することが望ましいからである。
【００８６】
そこで、図１０に示すように、帯域フィルタを通過した出力電力信号ＰB（ｔ）に所定の閾値ＰB,Cを適用し、振動振幅が小である時間領域での振動を除外する。この結果、図１１に示すように、物理的に不合理なものが除外され、風力発電機の周波数推定値ｆTのトレンドが得られる。
【００８７】
以上述べたように、本発明の実施の形態によれば、出力電力信号の時系列データに基づいて出力電力信号の減衰振動の自己回帰モデルを作成し、自己回帰モデルの係数から出力電力信号の振動周波数を算出するので、短い振動期間の時系列データでも出力電力信号の振動周波数を算出できる。
【００８８】
すなわち、フーリェ変換による場合は一定振動モデルであるため、例えば、風力発電機の制御系が風速の変化などで発生した振動の抑制を開始し減衰させた場合、特に振動の持続時間が短い場合には、その減衰振動に対応しきれないが、本発明の実施の形態では、自己回帰モデルとして先験的に減衰振動波形モデルを採用するので、短い振動期間の時系列データでも出力電力信号の振動周波数を算出できる。
【００８９】
一方、自己回帰モデルの係数ａ1、ａ2を最小自乗法で求めるので、分解能は時系列データを得るサンプリング間隔で規定されないし、フーリェ変換による場合に比較して振動周波数を算出するための計算量が少なくて済む。さらに、減衰振動波形モデルを仮定するので、実データから得られる風力発電機の時定数との対応から、推定結果である振動周波数の妥当性の判定が容易に行える。
【図面の簡単な説明】
【００９０】
【図１】本発明の実施の形態に係わる風力発電機の運転状態判別方法のフローチャート。
【図２】本発明の実施の形態で使用する帯域フィルタの通過率ＫＢの対周波数特性図。
【図３】本発明の実施の形態において、測定された風力発電機の出力電力信号Ｐ（ｔ）に連続2回に亘って帯域フィルタを適用した場合の風力発電機の出力電力信号ＰB（ｔ）の対周波数トレンド図。
【図４】本発明の実施の形態における所定の閾値ＰB,Cと帯域フィルタを通過した出力電力信号ＰB（ｔ）との関係を示すトレンド図。
【図５】本発明の実施の形態において、風力発電機の出力電力信号ＰB、振動周波数ｆ及び減衰振動を示す時定数Ｔを含む係数ａ2のトレンド図。
【図６】本発明の実施の形態において、移動平均法により平滑化された係数[ａ2]及びその平滑化された係数[ａ2]の変化分の絶対値｜Δ[ａ2,k]｜のトレンド図。
【図７】本発明の実施の形態において、得られた振動周波数推定値ＴTの分布図。
【図８】極数切換型の風力発電機の出力電力Ｐ（ｔ）の特性及び可変速型の風力発電機の出力電力Ｐ（ｔ）の特性の特性図。
【図９】帯域フィルタを通過した風力発電機の出力電力信号ＰB（ｔ）に所定の閾値ＰB,Cを適用しない場合の周波数推定値ｆTのトレンド図。
【図１０】帯域フィルタを通過した出力電力信号ＰB（ｔ）と所定の閾値ＰB,Cを示すトレンド図。
【図１１】帯域フィルタを通過した風力発電機の出力電力信号ＰB（ｔ）に所定の閾値ＰB,Cを適用した場合の周波数推定値ｆTのトレンド図。
【符号の説明】
【００９１】
Ｓ１…出力電力読み込み処理、Ｓ２…帯域フィルタリング処理、Ｓ３…時系列データ集積処理、Ｓ４…自己回帰モデル作成処理、Ｓ５…振動周波数算出処理、Ｓ６…振動周波数妥当性判定処理、Ｓ７…振動周波数推定値抽出処理、Ｓ８…回転数特性処理
【図１】
2006246604.tif 000024
【図２】
2006246604.tif 000025
【図３】
2006246604.tif 000026
【図４】
2006246604.tif 000027
【図５】
2006246604.tif 000028
【図６】
2006246604.tif 000029
【図７】
2006246604.tif 000030
【図８】
2006246604.tif 000031
【図９】
2006246604.tif 000032
【図１０】
2006246604.tif 000033
【図１１】
2006246604.tif 000034