(57)【要約】      （修正有）【課題】メソスケール領域での複数サイトにおける自然エネルギー発電機の発電出力の出力変動を分析して一般式で表現する。【解決手段】風力発電機等がＮ箇所に設けられ、Ｎ箇所で実測した風力発電機等の出力を総計してある周波数ｆにおける実測型総出力変動Ｓ（ｆ）を求め、Ｎ箇所で実測した風力発電機等の出力をＮ箇所それぞれについて個別に実測型個別出力変動Ｓｉ（ｆ）として求め、ある変動周期Ｔｘから同期が崩れランダムに移行すると仮定して数式で示す遷移型総出力変動Ｓｔｒａ（ｆ）により実測型総出力変動Ｓ（ｆ）を表現する。【選択図】図１４
2010239782.tif 000002
【特許請求の範囲】
【請求項１】
自然エネルギー発電機（発電設備も含む。）がＮ箇所に設けられた場合であって、Ｎ箇所で取得した時系列の実測データを時間成分のまま総計してから周波数成分に分けることによって、ある周波数ｆにおける実測型総出力変動Ｓ（ｆ）を求め、一方、Ｎ箇所で取得した時系列の実測データをＮ箇所それぞれについて個別に周波数成分に分け、ある周波数ｆにおける各箇所の実測データを実測型個別出力変動Ｓｉ（ｆ）として求め（ｉ＝１〜Ｎ）、
ある変動周期Ｔｘから同期が崩れランダムに移行すると仮定して下記数式（１）と（２）を引用する下記の数式（３）で示す遷移型総出力変動Ｓｔｒａ（ｆ）により実測型総出力変動Ｓ（ｆ）を表現するために、数式（３）で用いる変動周期Ｔｘを変化させた遷移型総出力変動Ｓｔｒａ（ｆ）と、実測型総出力変動Ｓ（ｆ）との誤差が小さくなる変動周期Ｔｘを、遷移時定数として求める自然エネルギー発電出力変動の分析方法。
【数１】
2010239782.tif 000016
【数２】
2010239782.tif 000017
【数３】
2010239782.tif 000018
【請求項２】
自然エネルギー発電機（発電設備も含む。）がＮ箇所に設けられた場合であって、
Ｎ箇所で取得した時系列の各実測データと、Ｎ箇所の実測データを時間成分のまま総計した時系列の総計データについて、出力変動幅算出処理を行い、
出力変動幅算出処理では、出力変動幅を評価する時間窓としての基準時間間隔の始点前後における短時間間隔での出力平均値と、基準時間間隔の終点前後における短時間間隔での出力平均値との差からなる出力変動幅を算出処理し、この算出処理を、時間窓を各実測データの時系列の先頭から末尾へ移動させながら繰り返すことによって、その時間窓における時系列の出力変動幅を求め、その時間窓における時系列の出力変動幅の確率分布に基づいてその時間窓における出力変動幅を決定するまでの一連の処理を行い、
出力変動幅算出処理を、時間窓の基準時間間隔を変化させながら繰り返すことによって、各基準時間間隔についての時間窓における出力変動幅を求め、
基準時間間隔の逆数をｆとし、時系列の総計データに対する各基準時間間隔についての時間窓における出力変動幅を、ｆを変数とする実測型総出力変動Ｓ（ｆ）として求め、
一方、Ｎ箇所の実測データに対する各基準時間間隔についての時間窓における出力変動を、ｆを変数とする実測型個別出力変動Ｓｉ（ｆ）としてそれぞれ求め（ｉ＝１〜Ｎ）、
ある変動周期Ｔｘから同期が崩れランダムに移行すると仮定して下記の数式（１）と（２）を引用する下記の数式（３）で示す遷移型総出力変動Ｓｔｒａ（ｆ）により実測型総出力変動Ｓ（ｆ）を表現するために、数式（３）で用いる変動周期Ｔｘを変化させた遷移型総出力変動Ｓｔｒａ（ｆ）と、実測型総出力変動Ｓ（ｆ）との誤差が小さくなる変動周期Ｔｘを、遷移時定数として求める自然エネルギー発電出力変動の分析方法。
【数４】
2010239782.tif 000019
【数５】
2010239782.tif 000020
【数６】
2010239782.tif 000021
【請求項３】
自然エネルギー発電機（発電設備も含む。）が現在設けられている場合であって、自然エネルギー発電機が将来導入されるときの発電出力の出力変動を推定するにあたって、
請求項１もしくは請求項２に記載の自然エネルギー発電出力変動の分析方法により求められた遷移時定数Ｔｘを用いる遷移型総出力変動Ｓｔｒａ（ｆ）又は現在までの測定により得られた実測型総出力変動Ｓ（ｆ）を表す総出力変動Ｓｐｒｅ（ｆ）と、下記数式（５）のゲイン関数Ｇ（ｆ）とに基づく下記数式（４）により、将来導入されるときの想定変動Ｓｅｓｔ（ｆ）を求める自然エネルギー発電出力変動の推定方法。但し、ゲイン関数Ｇ（ｆ）の遷移時定数Ｔｙおよびαを現在の自然エネルギー発電機の分布密度と将来のそれとの状況に応じて変化させる係数とする。
【数７】
2010239782.tif 000022
【数８】
2010239782.tif 000023
【請求項４】
請求項１又は請求項２に記載の自然エネルギー発電出力変動の分析方法により求められた遷移時定数Ｔｘと、その遷移時定数Ｔｘを用いる遷移型総出力変動Ｓｔｒａ（ｆ）又は現在までの測定により得られた実測型総出力変動Ｓ（ｆ）の何れかを表す総出力変動Ｓｐｒｅ（ｆ）と、下記数式（６）とにより、将来導入されるときの想定変動Ｓｅｓｔ（ｆ）を求める自然エネルギー発電出力変動の推定方法。但し、βを現在の自然エネルギー発電機の分布密度と将来のそれとの状況に応じて変化させる係数とする。
【数９】
2010239782.tif 000024
【請求項５】
Ｎ箇所で取得した自然エネルギー発電機の出力を分析するためにコンピュータに、
Ｎ個所の風力発電所における時系列の実測データを記憶装置に保存するデータ入力ステップと、
Ｎ箇所の実測データを時間成分のまま総計してから周波数成分に分けることによって、分析対象となる時間範囲に基づく周波数ｆにおける実測型総出力変動Ｓ（ｆ）を求め、一方、Ｎ箇所の実測データそれぞれについて個別に周波数成分に分け、周波数ｆにおける各箇所の自然エネルギー発電機の出力変動を実測型個別出力変動Ｓｉ（ｆ）として求める（ｉ＝１〜Ｎ）周波数分析ステップと、
実測型個別出力変動Ｓｉ（ｆ）を下記の数式（１）の通りに計算処理を行って、その計算処理結果を同期型出力変動Ｓcoh（ｆ）として記憶装置に保存すると共に、ランダム値計算ステップでは、実測型個別出力変動Ｓｉ（ｆ）を下記の数式（２）の通りに計算処理を行って、その計算処理結果をランダム型出力変動Ｓran（ｆ）として記憶装置に保存するコヒーレント・ランダム値計算ステップと、
変動周期Ｔｘを変数とし且つ下記数式（１）と（２）を引用する下記数式（３）で示す遷移型総出力変動Ｓｔｒａ（ｆ）と、実測型総出力変動Ｓ（ｆ）との誤差が小さくなる変動周期Ｔｘを遷移時定数として求める遷移時定数計算ステップ
を実行させるための自然エネルギー発電出力変動の分析プログラム。
【数１０】
2010239782.tif 000025
【数１１】
2010239782.tif 000026
【数１２】
2010239782.tif 000027
【請求項６】
Ｎ箇所で取得した自然エネルギー発電機の出力を分析するためにコンピュータに、
Ｎ個所の風力発電所における時系列の実測データを記憶装置に保存するデータ入力ステップと、
Ｎ箇所の実測データを時間成分のまま総計して総計データを記憶装置に保存する総計ステップと、
出力変動幅を評価する時間窓としての基準時間間隔の始点前後における短時間間隔での出力平均値と、基準時間間隔の終点前後における短時間間隔での出力平均値との差からなる出力変動幅を算出処理し、この算出処理を、時間窓を各実測データの時系列の先頭から末尾へ移動させながら繰り返すことによって、その時間窓における時系列の出力変動幅を求め、その時間窓における時系列の出力変動幅の確率分布に基づいてその時間窓における出力変動幅を決定する一連の処理を行う出力変動幅算出処理ステップと、
出力変動幅算出処理ステップを、時間窓の基準時間間隔を変化させながら繰り返すことによって、各基準時間間隔についての時間窓における出力変動幅を求めるループ処理ステップと、
基準時間間隔の逆数をｆとし、時系列の総計データに対する各基準時間間隔についての時間窓における出力変動幅を、ｆを変数とする実測型総出力変動Ｓ（ｆ）として求め、一方、Ｎ箇所の実測データに対する各基準時間間隔についての時間窓における出力変動を、ｆを変数とする実測型個別出力変動Ｓｉ（ｆ）としてそれぞれ求める（ｉ＝１〜Ｎ）周波数相当分析ステップと、
Ｎ箇所の実測データを時間成分のまま総計してから周波数成分に分けることによって、分析対象となる時間範囲に基づく周波数ｆにおける実測型総出力変動Ｓ（ｆ）を求め、一方、Ｎ箇所の実測データそれぞれについて個別に周波数成分に分け、周波数ｆにおける各箇所の自然エネルギー発電機の出力変動を実測型個別出力変動Ｓｉ（ｆ）として求める（ｉ＝１〜Ｎ）周波数相当分析ステップと、
実測型個別出力変動Ｓｉ（ｆ）を下記の数式（１）の通りに計算処理を行って、その計算処理結果を同期型出力変動Ｓcoh（ｆ）として記憶装置に保存すると共に、ランダム値計算ステップでは、実測型個別出力変動Ｓｉ（ｆ）を下記の数式（２）の通りに計算処理を行って、その計算処理結果をランダム型出力変動Ｓran（ｆ）として記憶装置に保存するコヒーレント・ランダム値計算ステップと、
変動周期Ｔｘを変数とし且つ下記数式（１）と（２）を引用する下記数式（３）で示す遷移型総出力変動Ｓｔｒａ（ｆ）と、実測型総出力変動Ｓ（ｆ）との誤差が小さくなる変動周期Ｔｘを遷移時定数として求める遷移時定数計算ステップ
を実行させるための自然エネルギー発電出力変動の分析プログラム。
【数１３】
2010239782.tif 000028
【数１４】
2010239782.tif 000029
【数１５】
2010239782.tif 000030
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、電力系統に自然エネルギー発電機が現在複数連系されている場合において、自然現象によって変動する発電出力を分析する方法、及びプログラムに関する。また、電力系統に自然エネルギー発電機が現在複数連系されている場合において、将来さらに自然エネルギー発電機が連系されることによって変動する出力を推定する方法に関する。
【背景技術】
【０００２】
風力や太陽光のような自然エネルギーは、CO2排出抑制の有力な手段として期待されているが、出力が自然現象によって変動する。このため、マクロな問題として、出力変動が電力系統の周波数変動を引き起こす恐れがある。
【０００３】
マクロな問題への対策については、現在のところ風力発電所などに二次電池やキャパシタを併設して出力変動を平滑化する研究などが行われている（非特許文献１）。また、海外に目を向けると、風力発電の導入を積極的に進めているドイツ・デンマーク等においては、すでに総電力需要のうちの３割以上を風力発電により供給するエリア、最軽負荷時における風力発電の割合が５割を超えるエリアなどが存在し、その系統運用上の問題点や対策が報告されている（非特許文献２）。しかし、風況はその国および地域の地形や気候特性により大きく異なることが予想されるため、日本国内でも各地域における実測データの収集および蓄積と、実測データに基づく分析が重要と考えられる。
【０００４】
＜スペクトル解析による風力出力変動の分析方法＞
電力系統の周波数は、図１に示すように対象とする負荷および発電の変動周期によって給電調整、ＬＦＣ等の幾つかの異なる手段を組み合わせることにより制御されている。このことから、風力発電の影響を評価するにあたっては、その出力変動がどのような周波数成分を持つかを調べることが有効であると考えられる。
【０００５】
大規模ウィンドファームの導入などにより、出力変動の短周期成分の平滑化効果が指摘されており（非特許文献３）、国内においても出力変動の実測データの収集・検討が実施されている（非特許文献３、４）。
【０００６】
しかし、メソスケール（〜数百km）での平滑化現象の分析は今後の課題とされており、比較的広い地域にわたる実測データによる検討は十分に行われているとは言い難い。
【０００７】
＜周波数成分の分析＞
また、周波数成分の分析には、横軸に周期ないし周波数を、縦軸に風力発電出力のパワースペクトル密度（ｋＷ２/Ｈｚ）をとるグラフを用いる例が多い（非特許文献６）。仮にＴ秒間のデータについてFast Fourier Transform （以後ＦＦＴ）解析を実施した場合、基本周波数１／Ｔ（Hz）の整数倍ごとの周波数成分が得られ、これを用いてパワースペクトルを算出・表示すると図２のようになり、周波数の高い部分が見にくくなる。
【０００８】
このため、便宜的にＦＦＴで算出された各次数のパワースペクトルを周波数の対数目盛り(0.001, 0.01, 0.1, 1, 10)に対して１０等分刻みで(例えば対数目盛り0.1〜1の間を１０等分刻みで)合計し、さらにその平方根を算出することで、縦軸を当該の周波数帯域の理想的なバンドパスフィルタを通過した出力値の標準偏差（kW）に相当する値で表示する手法がある（非特許文献５）。
【０００９】
この手法を用いた例として、定格出力２７５kWの風車について、ある程度の風が吹いた１日を無作為に抽出し、周波数成分に分解した結果を図３に示す。０．５秒周期あたりに塔体影効果（Tower  Shadow  Effect）による鋭いピークが見られる。また、全周波数領域にわたる変動成分の合計値は９２kWで、これは定格出力の３３％にもなる。
【００１０】
＜ランダム及び同期した出力変動の特性＞
周波数ｆにおける、Ｎ箇所の風力発電機ないしサイトの出力変動Si(f)（ｉ=1〜N）の総計をＳ（ｆ）とすると、もし各出力変動が完全に同期（coherent）していれば、同期仮説の総出力変動Ｓcoh（ｆ）として下記の数式（１）が成り立つ。
【数１】
2010239782.tif 000003
一方、もし各出力変動が完全にランダム（random）であれば、ランダム仮説の総出力変動Ｓran(f)として下記の数式（２）が成り立つ。
【数２】
2010239782.tif 000004
【００１１】
一つの風力発電サイト内では、１〜１０分程度を境にして、これより長周期では同期、短周期ではランダムに近くなるとの報告がなされている(非特許文献６)。これから、複数サイト間の平滑化に関しても、次のような傾向があると考えられる。
・「風が吹く日、吹かない日」というような気象レベルの遅い変動は同期する傾向がある。
・１０秒周期のような速い変動はランダムである。
【００１２】
＜時系列の実測データによる風力出力変動の分析＞
上述したスペクトルによる解析は、変動をランダムな周期的変動として扱うＬＦＣに対しては有効であるが、合計出力の突変が問題となる予備力評価には向かない。そこで、実測された風力発電の出力変動について、時系列の実測データをそのまま用いる分析も考えられる。時系列の実測データを用いる手法としては、一定時間間隔での出力変化(time variation)や、一定時間幅内の最大出力変動幅(widest span)により分析を行っている（非特許文献３、７、８、９）。なお、この明細書での時系列の実測データとは、電力系統の観測点で計測した時系列の電圧、電流データから求められた時系列の有効電力データのことをいう。
【先行技術文献】
【非特許文献】
【００１３】
【非特許文献１】新エネルギー・産業技術総合開発機構、エネルギー総合工学研究所：「蓄電池併設風力発電導入可能性調査  平成12年度調査報告書」、NEDO-NP-0004(2002)
【非特許文献２】L. Soder、 L. Hofmann、 A. Orths、H. Holttinen、 Y. Wan and A. Tuohy : “Experience From Wind Integration in SomeHigh Penetration Areas”、 IEEE Trans Energy Conversion、 Vol.22、 No. 1、 pp. 4-12(2007)
【非特許文献３】総合エネルギー資源調査会新エネルギー部会：風力発電系統連系対策小委員会中間報告書 (2004)
【非特許文献４】K. Abe、 O. Ishioka、 Y. Ichikawa and S. Enomoto: “The Analysis ofThe Wind Power Fluctuation and A Consideration of Smoothing Effect.”、 T. IEEJapan、 Vol. 121-B、 No.12、 pp. 1681-1689 (2001)阿部公哉・石岡修・市川嘉則・榎本重朗：「風力発電所の出力変動実績の分析と平滑化効果に関する一考察」、電学論B、 121、12、 pp. 1681-1689 (2001)
【非特許文献５】Y. Sekine: Power SystemEngineering、 Denki-Shoin Inc.、 pp. 48 (1976) (in Japanese)関根泰次：電力系統工学、電気書院、pp.48  (1976)
【非特許文献６】七原俊也：「風力発電電力系統安定化等調査の概要について」、季報エネルギー総合工学、Vol. 25、No. 4、 pp. 42-57 (2003)
【非特許文献７】P. Gardner、 H. Snodin、 A.Higgins and S. McGoldrick、 “The Impacts of Increased Levels of Wind Penetrationon the Electricity Systems of the Republic of Ireland and Northern Ireland: Final Report” (2003)
【非特許文献８】ISET e. V.、 “Wind EnergyReport Germany 2006” (2006)
【非特許文献９】Y. Wan、 “Wind Power PlantBehaviors: Analyses of Long-Term Wind Power Data”、 NREL Technical Report (2004)
【発明の概要】
【発明が解決しようとする課題】
【００１４】
本発明は、メソスケール領域での複数サイトにおける自然エネルギー発電機の発電出力の出力変動を分析して一般式で表現すること、また、自然エネルギー発電機を将来導入した時における発電出力の出力変動を推定するに当たって、平滑化効果を一般式で表現することを目的とする。
【課題を解決するための手段】
【００１５】
請求項１の発明は、自然エネルギー発電機（発電設備も含む。）がＮ箇所に設けられた場合であって、Ｎ箇所で取得した時系列の実測データを時間成分のまま総計してから周波数成分に分けることによって、ある周波数ｆにおける実測型総出力変動Ｓ（ｆ）を求め、一方、Ｎ箇所で取得した時系列の実測データをＮ箇所それぞれについて個別に周波数成分に分け、ある周波数ｆにおける各箇所の実測データを実測型個別出力変動Ｓｉ（ｆ）として求め（ｉ＝１〜Ｎ）、ある変動周期Ｔｘから同期が崩れランダムに移行すると仮定して前述の数式（１）と（２）を引用する下記の数式（３）で示す遷移型総出力変動Ｓｔｒａ（ｆ）により実測型総出力変動Ｓ（ｆ）を表現するために、数式（３）で用いる変動周期Ｔｘを変化させた遷移型総出力変動Ｓｔｒａ（ｆ）と、実測型総出力変動Ｓ（ｆ）との誤差が小さくなる変動周期Ｔｘを、遷移時定数として求める自然エネルギー発電出力変動の分析方法である。
【数３】
2010239782.tif 000005
【００１６】
請求項２の発明は、自然エネルギー発電機（発電設備も含む。）がＮ箇所に設けられた場合であって、Ｎ箇所で取得した時系列の各実測データと、Ｎ箇所の実測データを時間成分のまま総計した時系列の総計データについて、出力変動幅算出処理を行い、出力変動幅算出処理では、出力変動幅を評価する時間窓としての基準時間間隔の始点前後における短時間間隔での出力平均値と、基準時間間隔の終点前後における短時間間隔での出力平均値との差からなる出力変動幅を算出処理し、この算出処理を、時間窓を各実測データの時系列の先頭から末尾へ移動させながら繰り返すことによって、その時間窓における時系列の出力変動幅を求め、その時間窓における時系列の出力変動幅の確率分布に基づいてその時間窓における出力変動幅を決定するまでの一連の処理を行い、出力変動幅算出処理を、時間窓の基準時間間隔を変化させながら繰り返すことによって、各基準時間間隔についての時間窓における出力変動幅を求め、基準時間間隔の逆数をｆとし、時系列の総計データに対する各基準時間間隔についての時間窓における出力変動幅を、ｆを変数とする実測型総出力変動Ｓ（ｆ）として求め、一方、Ｎ箇所の実測データに対する各基準時間間隔についての時間窓における出力変動を、ｆを変数とする実測型個別出力変動Ｓｉ（ｆ）としてそれぞれ求め（ｉ＝１〜Ｎ）、ある変動周期Ｔｘから同期が崩れランダムに移行すると仮定して上述の数式（１）と（２）を引用する上述の数式（３）で示す遷移型総出力変動Ｓｔｒａ（ｆ）により実測型総出力変動Ｓ（ｆ）を表現するために、数式（３）で用いる変動周期Ｔｘを変化させた遷移型総出力変動Ｓｔｒａ（ｆ）と、実測型総出力変動Ｓ（ｆ）との誤差が小さくなる変動周期Ｔｘを、遷移時定数として求める自然エネルギー発電出力変動の分析方法である。
【００１７】
請求項３の発明は、自然エネルギー発電機（発電設備も含む。）が現在設けられている場合であって、自然エネルギー発電機が将来導入されるときの発電出力の出力変動を推定するにあたって、請求項１もしくは請求項２に記載の自然エネルギー発電出力変動の分析方法により求められた遷移時定数Ｔｘを用いる遷移型総出力変動Ｓｔｒａ（ｆ）又は現在までの測定により得られた実測型総出力変動Ｓ（ｆ）を表す総出力変動Ｓｐｒｅ（ｆ）と、下記数式（５）のゲイン関数Ｇ（ｆ）とに基づく下記数式（４）により、将来導入されるときの想定変動Ｓｅｓｔ（ｆ）を求める自然エネルギー発電出力変動の推定方法である。但し、ゲイン関数Ｇ（ｆ）の遷移時定数Ｔｙおよびαを現在の自然エネルギー発電機の分布密度と将来のそれとの状況に応じて変化させる係数とする。
【数４】
2010239782.tif 000006
【数５】
2010239782.tif 000007
【００１８】
請求項４の発明は、請求項１又は請求項２に記載の自然エネルギー発電出力変動の分析方法により求められた遷移時定数Ｔｘと、その遷移時定数Ｔｘを用いる遷移型総出力変動Ｓｔｒａ（ｆ）又は現在までの測定により得られた実測型総出力変動Ｓ（ｆ）の何れかを表す総出力変動Ｓｐｒｅ（ｆ）と、下記数式（６）とにより、将来導入されるときの想定変動Ｓｅｓｔ（ｆ）を求める自然エネルギー発電出力変動の推定方法である。但し、βを現在の自然エネルギー発電機の分布密度と将来のそれとの状況に応じて変化させる係数とする。
【数６】
2010239782.tif 000008
【００１９】
請求項５の発明は、Ｎ箇所で取得した自然エネルギー発電機の出力を分析するためにコンピュータに、Ｎ個所の風力発電所における時系列の実測データを記憶装置に保存するデータ入力ステップと、Ｎ箇所の実測データを時間成分のまま総計してから周波数成分に分けることによって、分析対象となる時間範囲に基づく周波数ｆにおける実測型総出力変動Ｓ（ｆ）を求め、一方、Ｎ箇所の実測データそれぞれについて個別に周波数成分に分け、周波数ｆにおける各箇所の自然エネルギー発電機の出力変動を実測型個別出力変動Ｓｉ（ｆ）として求める（ｉ＝１〜Ｎ）周波数分析ステップと、実測型個別出力変動Ｓｉ（ｆ）を前述の数式（１）の通りに計算処理を行って、その計算処理結果を同期型出力変動Ｓcoh（ｆ）として記憶装置に保存すると共に、ランダム値計算ステップでは、実測型個別出力変動Ｓｉ（ｆ）を前述の数式（２）の通りに計算処理を行って、その計算処理結果をランダム型出力変動Ｓran（ｆ）として記憶装置に保存するコヒーレント・ランダム値計算ステップと、変動周期Ｔｘを変数とし且つ前述の数式（１）と（２）を引用する前述の数式（３）で示す遷移型総出力変動Ｓｔｒａ（ｆ）と、実測型総出力変動Ｓ（ｆ）との誤差が小さくなる変動周期Ｔｘを遷移時定数として求める遷移時定数計算ステップを実行させるための自然エネルギー発電出力変動の分析プログラムである。
【００２０】
請求項６の発明は、Ｎ箇所で取得した自然エネルギー発電機の出力を分析するためにコンピュータに、Ｎ個所の風力発電所における時系列の実測データを記憶装置に保存するデータ入力ステップと、Ｎ箇所の実測データを時間成分のまま総計して総計データを記憶装置に保存する総計ステップと、出力変動幅を評価する時間窓としての基準時間間隔の始点前後における短時間間隔での出力平均値と、基準時間間隔の終点前後における短時間間隔での出力平均値との差からなる出力変動幅を算出処理し、この算出処理を、時間窓を各実測データの時系列の先頭から末尾へ移動させながら繰り返すことによって、その時間窓における時系列の出力変動幅を求め、その時間窓における時系列の出力変動幅の確率分布に基づいてその時間窓における出力変動幅を決定する一連の処理を行う出力変動幅算出処理ステップと、出力変動幅算出処理ステップを、時間窓の基準時間間隔を変化させながら繰り返すことによって、各基準時間間隔についての時間窓における出力変動幅を求めるループ処理ステップと、基準時間間隔の逆数をｆとし、時系列の総計データに対する各基準時間間隔についての時間窓における出力変動幅を、ｆを変数とする実測型総出力変動Ｓ（ｆ）として求め、一方、Ｎ箇所の実測データに対する各基準時間間隔についての時間窓における出力変動を、ｆを変数とする実測型個別出力変動Ｓｉ（ｆ）としてそれぞれ求める（ｉ＝１〜Ｎ）周波数相当分析ステップと、Ｎ箇所の実測データを時間成分のまま総計してから周波数成分に分けることによって、分析対象となる時間範囲に基づく周波数ｆにおける実測型総出力変動Ｓ（ｆ）を求め、一方、Ｎ箇所の実測データそれぞれについて個別に周波数成分に分け、周波数ｆにおける各箇所の自然エネルギー発電機の出力変動を実測型個別出力変動Ｓｉ（ｆ）として求める（ｉ＝１〜Ｎ）周波数相当分析ステップと、実測型個別出力変動Ｓｉ（ｆ）を前述の数式（１）の通りに計算処理を行って、その計算処理結果を同期型出力変動Ｓcoh（ｆ）として記憶装置に保存すると共に、ランダム値計算ステップでは、実測型個別出力変動Ｓｉ（ｆ）を前述の数式（２）の通りに計算処理を行って、その計算処理結果をランダム型出力変動Ｓran（ｆ）として記憶装置に保存するコヒーレント・ランダム値計算ステップと、変動周期Ｔｘを変数とし且つ前述の数式（１）と（２）を引用する前述の数式（３）で示す遷移型総出力変動Ｓｔｒａ（ｆ）と、実測型総出力変動Ｓ（ｆ）との誤差が小さくなる変動周期Ｔｘを遷移時定数として求める遷移時定数計算ステップを実行させるための自然エネルギー発電出力変動の分析プログラムである。
【発明の効果】
【００２１】
本発明の分析方法、分析プログラム及び推定方法によれば、いずれも遷移時定数を用いた一般式で容易に出力変動の影響の検討を実現することができる。
【図面の簡単な説明】
【００２２】
【図１】電力系統の周波数制御方法を示すグラフである。
【図２】ＦＦＴ解析の例を示すグラフである。
【図３】図２のグラフを簡略化したグラフである。
【図４】短周期の出力変動を分析したグラフである。
【図５】長周期の出力変動を分析したグラフである。
【図６】出力変動の分析手法のうち移動平均方式を示すグラフである。
【図７】各種分析手法による比較結果を示すグラフである。
【図８（Ａ）（Ｂ）】移動平均方式による時間窓の一定時間間隔を２０分、２０時間とした場合の出力変動を示すグラフである。
【図９】移動平均方式による時間窓の一定時間間隔を変化させた場合の出力変動を示すグラフである。
【図１０】最大出力変動幅を想定するグラフである。
【図１１】ゲイン関数を想定するグラフである。
【図１２】風力大量導入時の出力変動の予測を示すグラフである。
【図１３】風力発電の増設イメージを示す図である。
【図１４】分析プログラムの第一例を示すフローチャートである。
【図１５】分析プログラムの第二例を示すフローチャートである。
【発明を実施するための形態】
【００２３】
本発明は、背景技術で記載したように、一つの風力発電サイト内での傾向から、複数サイト間の風力出力変動の傾向を推測したことに基づくものである。つまり、「風が吹く日、吹かない日」というような気象レベルの遅い変動は同期するであろう。１０秒周期のような速い変動はランダムであろう。という傾向を推測したものである。そして、推測した傾向が正しいことを、比較的短周期の場合と、長周期の場合で次のように確認した。
【００２４】
＜比較的短周期の出力変動＞
そこで、このような気象レベルの遅い出力変動と速い出力変動との中間部分において、数十秒から数十分周期（ガバナフリーやＬＦＣが対処すべき周期）の出力変動について時系列の実測データを用いて傾向を調査した。対象とした風力発電サイトは、３箇所（Ａ〜Ｃとする）で、ほぼ一直線上に並んでおり、サイトＡＢ間は約９０km、サイトＢＣ間は約８０km離れている。３サイトの合計出力は6.3ＭＷである。各サイトの発電機出力合計値は、ＧＰＳ時計により時刻同期をとって１秒ごとに計測されている。約１日間のデータについて、背景技術の＜周波数成分の分析＞で説明した方法（縦軸を当該の周波数帯域の理想的なバンドパスフィルタを通過した出力値の標準偏差（ｋＷ）に相当する値で表示する手法）によりスペクトラムを作成した結果を図４に示す。図４において、同期（coherent）およびランダム（random）で示す曲線は、前述した数式（１）および（２）により求めたものであり、measuredで示す曲線は、全サイトの出力（時系列の実測データ）を時間成分のまま合計した総計データに対してスペクトル解析した結果である。今回の時系列の総計データでは比較的短周期の出力変動については、ランダム仮説に近い結果が得られた。
【００２５】
＜長周期の出力変動＞
次に、さらに長期の出力変動について確認する。同じ風力3サイトについて1秒サンプリングした時系列の実測データの1分平均をとり、2006年8月30日から2007年2月末までの約182日間（218分間）の連続データをFFT解析してスペクトラムを作成した。結果を図5に示す。同図から遅い出力変動は同期する一方、速い出力変動はランダムとなる傾向があることがわかる。
【００２６】
本発明の第一実施形態である風力出力変動の分析方法は、スペクトラムのこのような傾向を確認した上で創作された。つまり、ある変動周期Txから同期が崩れランダムに移行すると仮定すると、スペクトラム（遷移型出力変動）Stra(f)は前述した下記の数式(3)のように表現できる。式中、fは周波数（周期の逆数）、jは虚数単位である。また、以後Txを「遷移時定数」と呼ぶこととする。
【数７】
2010239782.tif 000009
【００２７】
数式(3)において、遷移時定数Txを変化させ、合成スペクトラムと遷移スペクトラムS(f)の誤差が最小となる遷移時定数を求めたところ、約17万秒（48時間）となった。この時間は、「風の吹く日、吹かない日」のような気象状況による影響が支配的な時間間隔である。給電調整では1週間を一つの単位とすることから、遅い同期した風力出力変動の考慮が重要になることが予想される。
【００２８】
次に、合計出力の突変が問題となる予備力評価に適した本発明の第二実施形態である、時系列データによる風力出力変動の分析方法について述べる。
【００２９】
複数サイト間の時系列の実測データによる風力出力変動の分析方法を下記の表１に示す。手法A〜Dは、従来から行われているものであり、一定時間間隔（基準時間間隔）での出力変化(time  variation)や、一定時間間隔内の最大出力変動幅(widest  span)により評価を実施している。このうち通常の平均方式を採用したtime  variationとは次の通りである。まず、時系列の実測データに対して基準時間間隔（Time  span：例えば１時間）ごとに、短時間間隔（基準時間間隔よりも短い、例えば１０分間）における出力変動の平均値を求める。その後に、時系列順に前後する両者の出力変化（平均値の差）を求めるものである。
【００３０】
【表１】
2010239782.tif 000010
【００３１】
Eの手法は今回採用した出力変動幅算出処理手法で、移動平均方式を採用したtime  variationである。これは、出力変動を評価する時間窓としての基準時間間隔、図6中のTime  spanを利用するものである。そして、時系列の実測データに対して時間窓の基準時間間隔（Time  span：例えば１時間）の始点前後と終点前後の短時間間隔（例えば、図６ではTime  spanの半分、３０分間）における出力変動の平均値を求めてから両者の差を求め、この時間窓を短時間間隔よりも短い時間、例えば１分ごとに移動させながら同様の計算処理を繰り返して、各時間窓における時系列の出力変動幅を求めるものである。
【００３２】
これらの手法の比較に用いた時系列の実測データは、上述の＜長周期の出力変動＞で用いた３サイトの１分サンプリングデータ182日分と同一のものであり、出力変動の評価時間を1時間とした。最大出力変動幅については、変化の符号を無視する場合もあるが、ここでは他の手法との整合性をとるため、符号を考慮した。
【００３３】
表２は各手法による分析結果を示したものである。各サイトの定格発電機出力の合計値を１p.u.とした。最大および最小の行は、各手法における正方向および負方向への変化幅が最大であったデータを示しており、その絶対値はA、Dを除けば０．５p.u.程度となっている。
【表２】
2010239782.tif 000011
【００３４】
±3σおよび±2σの行は、求めた変動幅の確率分布の標準偏差を算出し、これを3倍および2倍した値である。また、Top、Bot.0.13%およびTop、Bot.2.3%の行は、求めた変動幅の分布の両端から0.13%および2.3%の位置（それぞれ正規分布における±3σおよび±2σの位置）における変動幅の値である。
【００３５】
例えばTop 2.3%の値は+2σの値より大きくなっており、Bot. 2.3%の値は-2σの値より小さくなっている。これは、±2σ範囲の外側に、確率95.4%の範囲が来ることを意味している。A〜Eいずれの手法を用いても、今回測定した変動幅の分布は正規分布と比べ端の部分の確率が高いことがわかる。
【００３６】
図7では、出力変動幅１％刻みでその確率をプロットした（いわゆるアンブレラカーブ）。Aは平滑化していないため短周期成分を含む変動を評価することとなり、他の方式に比べて変動幅が大きい部分の確率が高くなる傾向にある。また、手法B〜Dにおいては、変動幅の大きい存在確率の少ない部分についてはデータが欠落しやすい。それに比べて、手法Eは、手法B〜Dと類似の傾向を示すが、手法B〜Dにおいてはデータの欠落しやすい部分（変動幅の大きい存在確率の少ない部分）についても、欠落ないカーブが得られることがわかる。以降はこのEの方法を用いることにする。
【００３７】
＜複数サイトでの出力変動幅の平滑化＞
ここでは複数のサイトがある場合、各サイト出力の変動幅と全サイト合計出力の変動幅の間にはどのような関係があるかを分析した。
【００３８】
対象とした時系列の実測データは、＜長周期の出力変動＞と同様の３サイトの182日間の1分サンプリングデータである。図6におけるTime  span:  Tを20分、および20時間としてプロットした例を図8に示す。Tが短い20分の場合には、個別サイトのアンブレラカーブに比べてサイト合計のアンブレラカーブが顕著に狭くなっているのに対し、Tが長い20時間の場合には、すべてのアンブレラカーブが重なってくることがわかる。これは、長周期の各サイト出力変動が同時に発生する傾向を示しており、＜長周期の出力変動＞の分析結果とも定性的に一致する。
【００３９】
そこで、出力変動幅を評価する時間窓を、10分から約10日まで変化させて、各時間窓における時系列の出力変動幅の確率分布に基づいてその時間窓における出力変動幅を求めた上で、各サイトの変動幅の単純合計、２乗和の平方根および全サイト合計値の変動幅を求めた。ここでは、確率分布に基づいて時間窓における出力変動幅はTop0.13%とBot0.13%の絶対値の平均値とした。ただし、予備力評価という観点からすれば、時間窓における出力変動幅には、どのような値を用いても良く、確率分布とは無関係にMax、Minを採用しても良い。なお、単純合計等を求める際には、周波数に相当する値として、基準時間間隔（単位：秒）の逆数ｆを採用する。ｘ軸を時間窓、ｙ軸を変動幅としてプロットした結果を図9に示す。以下、前述の検証結果にならって、単純合計を「同期(coherent)」、2乗和の平方根を「ランダム(random)」、また全サイト合計の出力変動幅を便宜的にmeasured（実測型出力変動）と呼ぶこととする。
【００４０】
時間窓が2時間程度より短い場合には、measuredはランダムとほぼ一致ないしそれより小さいが、それより長い周期においては、周期が長くなるにつれて、「ランダム」に比べて大きくなり、「同期」のカーブに近づいていく傾向にあることがわかる。
【００４１】
前述した数式（３）を用いて「ランダム」と「同期」の遷移を表現すると、遷移時定数の値は、約15000秒（約4.2時間）程度となり、先の48時間と比較すると1/10程度の値となっている。これは、図１においてはLFCよりも長周期のEDCの領域に相当することから、大量導入時にはEDC領域に、より大きな変動が生じることが想定される。
【００４２】
＜風力大量導入時の検討＞
これまでの検討から得られた知見を用いて、風力発電が大量に導入された場合を検討する。これが、本発明の第三実施形態である、風力出力変動の推定方法である。
【００４３】
＜海外データによる検討＞
風力発電大量導入時の想定にあたり、風力発電機出力の規模に応じて、最大出力変動幅がどの程度になるかを海外における時系列の実測データを用いて検討した。非特許文献８の15分平均の時系列の実測データを基に作成された風力発電機1台、1箇所のウィンドファーム、ドイツ全体の出力変動のアンブレラカーブを使用し、表１のCの方法に相当すると考えられる最大出力変動幅を求めた。図10に、発電機定格合計出力(MW)をx軸に、1時間最大出力変動幅(MW)をy軸にプロットした。図中の点Aは風力発電機1台、点BはAの発電機が含まれる1箇所のウィンドファーム、点Cはドイツ全体を示す。出力変動が1時間の領域で完全にランダムであれば、ｘの0.5乗に比例するはずだが、これよりも大きい0.88程度に比例している。
【００４４】
この理由としては、次のような説明が考えられる。
・非常に長周期の出力変動では、サイト内、サイト間とも変動が同期する傾向にあるため、変動幅は発電機定格合計出力のほぼ1乗に比例する。
・一方、非常に短周期の変動では、サイト内、サイト間とも変動はランダムであるため、変動幅は発電機定格合計出力のほぼ0.5乗に比例する。
・その中間的な1時間程度の領域では、サイト間ではランダムだが、サイト内は同期しているため、先ほどの例においては0.88乗に比例する結果となった。
【００４５】
以上の推測をもとに、現時点の測定により得られた変動幅が大量導入時にどの程度になるかについて検討する。
現在までの測定により得られた変動幅としては、図9において横軸の1/Time
span をfとし、measuredで示す曲線（実測型総出力変動S(f)）をfの関数とした総出力変動Spre(f)を用いる。なお、遷移型総出力変動Stra(f)も誤差が小さければ実測型総出力変動S(f)とほぼ同じとなるので、Stra(f)をSpre(f)としても良い。また、大量導入時の想定変動をSest(f)とする。Sest(f)とSpre(f)の間に前述した下記の数式(４)の関係式があるものと仮定する。
【数８】
2010239782.tif 000012
【００４６】
ここで、Sest(f)を算出するためには、G(f)を想定する必要がある。ここでは図11に例示するような、４つのケースを想定しながら、前述した下記の数式（５）を採用した。この式は風力出力がＭ倍になると同期変動はＭ倍に、ランダム変動は√Ｍ倍となる特性を考慮した式である。ゲイン関数Ｇ（ｆ）の遷移時定数Ｔｙおよびαを現在の自然エネルギー発電機の分布密度と将来のそれとの状況に応じて変化させる係数とする。
【数９】
2010239782.tif 000013
【００４７】
・ケース１は現在のサイト分布状況と同じように、サイト間距離が比較的離れた状況でサイト数が増加する場合を想定したものである。上記の数式（５）のG(f)としては、α＝０、Ｔｙ＝Ｔｘを採用する。
【００４８】
・ケース２は、風力適地が地域的に偏在しており、既存サイトの増設に近いような状況で風力が導入されていく場合を想定したものである。これらの近接した発電機群の出力は、比較的短周期から同期が始まることとなり、G(f)は、ケース１に比較して短周期側のゲインが高くなると想定される。そこで、数式（５）の遷移時定数Tyには、Ｔｘとは別にサイト近傍での遷移時定数T’を用いることとし、その値としては非特許文献文６を参考にT’=600秒を仮定した。さらに、同期からランダムへの移行について、次の3パターンを想定した。
【００４９】
a)ケース2aは同期からランダムに移行する角周波数を1/T’とした。すなわち数式（５）のα=0に相当する。この場合、ランダムから同期に移行する周波数は高くなる。
【００５０】
b)ケース2bはランダムから同期に移行する角周波数を1/T’とした。すなわち数式（５）のα=0.5に相当する。この場合、同期からランダムに移行する周波数が低くなる。
【００５１】
c)ケース2cはケース2aおよびケース2bの中間的な特性を考えた。すなわち数式（５）のα=0.25に相当する。
なお図11においては、Ｍ=100、Ty=15000、 T’=600としている。
【００５２】
非特許文献８中のデータによると、ドイツ全体で16.8GWの風力発電機が存在し、1時間幅の最大変動が4.4GWであったとの報告がある。そこで、当社（北陸電力株式会社）で測定された合計最大出力6.3MWにおける1時間変動幅をベースとして、仮に16.8GW導入された状況における変動幅を各想定ケースについて算出したところ、表3に示すようにケース2cがドイツでの実態（4.4GW）に近い値となった。
【表３】
2010239782.tif 000014
【００５３】
＜大量導入時の検討＞
直前の分析手法を用いて、図9の結果をもとに、風力大量導入時の出力変動幅について検討した。風力発電の定格出力合計を最大負荷の40%と仮定した例を図12に示す。変動幅は3σ相当（上下0.13%位置）の値とし、あわせて負荷の変動幅もプロットした。先の検討において比較的実態に近いと考えられたケース2cに基づく想定では、1時間程度の幅の部分で風力の変動幅が負荷の変動幅と同程度となり、また数時間〜数十分程度の部分でも両者はごく近接している。すなわち、風力大量導入時には数時間〜数10分の変動に及ぼす影響が大きくなる可能性が高いことが予想される。
【００５４】
負荷の変動幅は、日間変動のため約半日の部分に顕著なピークが出るが、その変動はかなりの精度で予測できる反面、現状では出力変動の予測精度はそれほど高くないため、この領域についても出力変動は無視できないと考えられる。よって、今後風力出力変動の予測精度向上が重要になってくると考えられる。また、前述したように風力の変動が、正規分布と比べ端の部分の存在確率が高い分布であることから、図示した負荷の変動をより大きく上回るような出力変動が発生する可能性も無視できないことに留意する必要があると考えられる。
【００５５】
・ケース３Ａは、極端な例として、全ての新設風車が既存サイトの増設により導入されていく場合を考えた（図１３ ケース3A）。この場合、数式（５）の遷移時定数Ty＝０として、大量導入時の想定変動Sest(f)を表現できる。
【００５６】
・ケース３Ｂは、全ての風車が満遍なく分布し、各サイトの規模は現状と同程度の場合を考えた（図１３ ケース3B）。遅い変動の振幅は、ケース3Aと同じくM・Stra(f)となろう。また、速い変動の振幅は、Sran(f)を現状のランダム仮説による遷移型総変動変動とすると、M0.5
・Sran(f)となろう。すなわち、数式（５）の代わりに前述した下記の数式（６）において、β=0.5とすることにより、大量導入時の想定変動Sest(f)を表現できる。
【数１０】
2010239782.tif 000015
【００５７】
・実態はケース３Ａと３Ｂの中庸に近いと考えられるので、ケース３ＡＢとして、上記の数式(６)において、例えばβ=0.75とすることにより、大量導入時の想定変動Sest(f)を表現できる。
【００５８】
上述した自然エネルギー発電出力変動の分析方法は、キーボードやマウス等からなる入力装置、ディスプレイ等からなる出力装置、分析プログラムの命令を順番に実行するＣＰＵ、分析プログラムや分析プログラムの実行に必要なデータ及び計算結果等を保存する記憶装置を構成要素とする標準的なコンピュータを用いて行われる。
【００５９】
図１４のフローチャートは前述した第一実施形態をコンピュータで実行させるための分析プログラムの第一例である。分析プログラムをコンピュータに実行させることにより、コンピュータが各種の手段として機能し、各種のステップが順番に行われる。
【００６０】
まず、データ入力ステップが行われる。データ入力ステップでは、複数個所の風力発電所における同一時刻範囲内の発電機出力データ（時系列の有効電力データ）を入力する行列状の入力フォームを記憶装置から読み込んで、出力装置に表示する。入力フォームは、見出し欄の項目の１列目に通し番号、２列目に日時（年月日及び時分秒）、３列目以降に各風力発電所の実測データ（有効電力データ）を並列して表示したものである。これら見出し欄の下に各内容の入力欄が時系列順（例えば１秒経過毎）に設けてある。この入力欄に、作業者が実測データ等を、入力装置から手作業で入力し、実測データ等の入力が完了したことを入力装置から入力すると、実測データ等が記憶装置に保存される。また、入力作業を簡略化するために、このコンピュータに別の記憶装置を接続した上で、この別の記憶装置に記憶されているファイルを入力装置から指定し、分析プログラムの入力フォームにファイルの内容を取り込む操作をすることによって、そのファイルの内容（入力欄に実測データが登録されたのもの）をこのコンピュータの記憶装置に保存させるものでも良い。続いて、対象指定ステップが実行される。
【００６１】
対象指定ステップでは、先のデータ入力ステップで入力された実測データのうち、分析対象となる時間範囲及び風力発電所並びに遷移時定数の探索範囲を指定する指定フォームを記憶装置から読み込んで、出力装置に表示する。指定フォームは、見出し欄の項目として時間範囲及び風力発電所並びに遷移時定数の探索範囲を表示し、各項目の下に指定欄を表示するものである。時間範囲の指定欄は、日時（何年何月何日何時何分何秒〜何年何月何日何時何分何秒）を指定（入力）するものである。また、風力発電所の指定欄は、風力発電所（例えばＡ発電所、Ｂ発電所）を指定（入力）するものである。遷移時定数の探索範囲の指定欄は、後述の遷移時定数計算ステップでのループ処理で用いられる数値としての初期値、増分、終値を、例えば１万、１０００、１００万と指定（入力）するものである。この値は、経験値である。そして、これらデータの指定が完了したこと入力手段から入力すると、これら指定データが記憶装置に保存されると共に、周波数分析ステップが実行される。また、これらデータの指定をしなくとも、予め経験値等のデフォルト条件が初期設定として記憶装置に保存されている場合には、対象指定ステップを通り越して、周波数分析ステップが実行される。
【００６２】
周波数分析ステップでは、データ入力ステップで入力された実測データのうち指定された対象について、以下の第一〜第三の処理を行う。
【００６３】
第一の処理。まず、先のステップで指定された時間範囲を記憶装置から読み込んで時間長（秒数）に換算し、時間長を記憶装置に保存する。便宜上、時間長がＴ秒に換算されたものとする。
【００６４】
第二の処理。次に、指定された対象の全ての風力発電所について、指定された時間範囲の実測データを記憶装置から読み込んで、同一時刻毎に実測データを総計（時間成分のまま総計）して総計データを求め、全風力発電所の総計データを記憶装置に保存する。また、全風力発電所の総計データについてＦＦＴ解析を行い、基本周波数１／Ｔ（Ｈｚ）、即ち周波数ｆの整数倍ごとの周波数成分を計算し、その計算結果を実測型総出力変動Ｓ（ｆ）として記憶装置に保存する。
【００６５】
第三の処理。次に、指定された対象について、風力発電所毎に、指定された時間範囲の実測データを記憶装置から読み込んで、ＦＦＴ解析を行い、周波数ｆの整数倍ごとの周波数成分を計算し、その計算結果を実測型個別出力変動Ｓｉ（ｆ）として記憶装置に保存する。指定された風力発電所の全てについて、これらの処理が行われる。なお、ここでのｉには、実測データが読み込まれた風力発電所の順番通りの数が付与される。以上の３つの処理が終わった後、コヒーレント、ランダム値計算ステップが実行される。
【００６６】
コヒーレント、ランダム値計算ステップでは、まず、コヒーレント値計算ステップが実行される。コヒーレント値の計算ステップでは、指定された風力発電所の全てについて先の第三の処理で計算された実測型個別出力変動Ｓｉ（ｆ）を記憶装置から読み込んで、前述した数式（１）の通りに計算処理を行って、その計算処理結果を同期型出力変動Ｓcoh（ｆ）として記憶装置に保存する。
【００６７】
次にランダム値計算ステップが実行される。ランダム値計算ステップでは、指定された風力発電所の各々について、先の第三の処理で計算された実測型個別出力変動Ｓｉ（ｆ）を記憶装置から読み込んで、前述した数式（２）の通りに計算処理を行って、その計算処理結果をランダム型出力変動Ｓran（ｆ）として記憶装置に保存する。
【００６８】
続いて、平滑化ステップが実行される。平滑化ステップでは、実測型総出力変動Ｓｉ（ｆ）、同期型出力変動Ｓcoh（ｆ）、及びランダム型出力変動Ｓran（ｆ）のそれぞれについて順番に同じ処理を行う。まず、実測型総出力変動Ｓｉ（ｆ）を記憶装置から読み込んで、周波数の0.001〜0.01の間、0.01〜0.1の間、0.1〜1の間、1〜10の間のそれぞれについて、対数スケールの均等な割合で等分（例えば１０等分）した範囲のデータを合計し、さらにその合計値の平方根を算出し、この平方根の結果を平滑化した実測型総出力変動Ｓｉ（ｆ）として記憶装置に保存する。同期型出力変動Ｓcoh（ｆ）、及びランダム型出力変動Ｓran（ｆ）についても同処理を行って、それぞれ平滑化した同期型出力変動Ｓcoh（ｆ）、及びランダム型出力変動Ｓran（ｆ）として、記憶装置に保存する。
【００６９】
続いて遷移時定数計算ステップが実行される。遷移時定数計算ステップでは、記憶装置に保存されている数値としての初期値、増分、終値（例えば１万、１０００、１００万）を変数Tx(i)とする遷移時定数の探索範囲においてループ処理を行う。ループ処理では、最初に、遷移時定数の探索範囲のうち、次の遷移時定数候補である変数Tx(i)、ここでは初期値を選択して、記憶装置から読み込む。その上で、平滑化した同期型出力変動Ｓcoh（ｆ）、及びランダム型出力変動Ｓran（ｆ）を記憶装置から読み込んで、前述した数式（３）により遷移型総出力変動Stra（ｆ）を計算し、変数iと関連付けた形でその計算結果を記憶装置に保存する。次に、この計算結果と平滑化した実測型総出力変動Ｓｉ（ｆ）とを記憶装置から読み込んで、最小二乗法によって誤差二乗和を算出する。今回の誤差二乗和Ｅrr（ｉ）を前回までの誤差二乗和Ｅrr（ｉ）と比較する。最初のループ処理のみ、初期設定された非常に大きな条件値を、前回までの誤差二乗和Ｅrr（ｉ）として用いる。そして、今回の誤差二乗和Ｅrr（ｉ）が小さければ、変数Tx(i)を遷移時定数Txとして記憶装置に保存し、最後の処理に移って、遷移時定数がその探索範囲の終値に達したか否かを判定する。一方、今回の誤差二乗和Ｅrr（ｉ）が小さくなければ、そのまま最後の処理に移って、遷移時定数候補である変数Tx(i)がその探索範囲の終値に達したか否かを判定する。ここでは、変数Tx(i)に初期値が使われているので、終値には達していない場合となり、このステップの最初の処理に戻り、遷移時定数の探索範囲のうち次の遷移時定数候補を、初期値に増分を加算した値として記憶装置に保存し、以後同様の処理を繰り返す。一方、処理を繰り返しているうちに、最後の処理で、遷移時定数候補である変数Tx(i)が終値に達し、その処理が終了した場合には、遷移時定数Txが確定し、計算結果出力ステップが実行される。
【００７０】
最後の計算結果出力ステップとして、出力装置の画面に遷移時定数Txを表示すると共に、行列状の出力フォームを記憶装置から読み込んで、出力装置に表示する。出力フォームは、見出し欄の項目の１列目にｆ、２列目以降に遷移時定数Txと同じ値の変数Tx(i)に関連付けた形の平滑化した同期型出力変動Ｓcoh（ｆ）、及びランダム型出力変動Ｓran（ｆ）、並びに平滑化した実測型総出力変動Ｓｉ（ｆ）の順で並列して表示したものである。これら見出し欄の下に各内容の出力欄が時系列順に設けてある。また、行列状の出力フォームを出力装置に表示する代わりに、遷移時定数Txと同じ値の変数Tx(i)に関連付けた形の平滑化した同期型出力変動Ｓcoh（ｆ）、及びランダム型出力変動Ｓran（ｆ）、並びに平滑化した実測型総出力変動Ｓｉ（ｆ）を、ｘ軸に出力（kＷ）を、ｙ軸にｆ（ｓ）をそれぞれ取った折れ線グラフの中で同時に出力装置に表示してもよい。これにより、分析プログラムの第一例の全てのステップが終了する。
【００７１】
図１５のフローチャートは前述した第二実施形態をコンピュータで実行させるための分析プログラムの第二例である。最初に第一例の分析プログラムと同様にデータ入力ステップが実行される。続いて、実行される対象指定ステップでは、先例の分析プログラムと同様の分析対象となる時間範囲、風力発電所、及び遷移時定数の探索範囲に加え、後述の出力変動幅算出処理ステップで用いる時間窓の基準時間間隔を指定する指定フォームを記憶装置から読み込んで、出力装置に表示する。指定フォームは、見出し欄の項目として時間範囲及び風力発電所並びに遷移時定数の探索範囲の他に、時間窓の基準時間間隔、時間窓の移動時間を表示し、各項目の下に指定欄を表示するものである。また、時間窓の基準時間間隔の指定欄は、ループ処理で用いられる数値としての初期値、増分、終値を、例えば初期値１０分、増分１０分、終値１０日と指定（入力）するものである。時間窓の移動時間の指定欄は、同じくループ処理ステップで用いられる数値を、例えば１分と指定（入力）するものである。そのほかの時間範囲及び風力発電所並びに遷移時定数の探索範囲の各指定欄は、先例の指定フォームと同様である。そして、これらデータの指定が完了したことを入力手段から入力すると、これら指定データが記憶装置に保存されると共に、総計ステップが実行される。
【００７２】
総計ステップでは、指定された対象の全ての風力発電所について、指定された時間範囲の実測データを記憶装置から読み込んで、同一時刻毎に実測データを総計して（時間成分のまま総計）、その総計データを記憶装置に保存する。
【００７３】
続いて、ループ処理ステップが実行される。ループ処理ステップの中では、出力変動幅算出処理ステップが実行される。出力変動幅算出処理ステップでは、時間窓の基準時間間隔の初期値１０分が記憶装置から読み込まれ、その１／２の値を計算し、計算値（５分）を短時間間隔として記憶装置に保存する。その後に、全風力発電所の総計データを記憶装置から読み込み、基準時間間隔１０分の始点前後と終点前後の短時間間隔５分における出力変動の平均値をそれぞれ算出し、両者の平均値の差を算出処理することを、総計データの先頭から末尾まで時間窓を移動時間１分ごとに移動させながら繰り返し、その時間窓における時系列の出力変動幅を求め、求めた結果を記憶装置に保存する。続いて、指定された対象の風力発電所の各々についてその指定された時間範囲の実測データについても、同様に、その時間窓における時系列の出力変動幅を記憶装置に保存する。
【００７４】
続いて、出力変動幅算出処理ステップでは、全風力発電所の総計データについて求められた時間窓における時系列の出力変動幅に対して確率分布をそれぞれ適用する。そして、時系列の出力変動幅に対する全体の個数を１００％として取り扱い、確率分布の上位Ｘ（例えば0.13）％、下位Ｘ％の位置における出力変動幅を抽出し、これら出力変動幅の絶対値の平均値を求め、その時間窓における出力変動幅として記憶装置に保存する。なお、0＜Ｘ≦100とする。各風力発電所の実測データついて求められた時間窓における時系列の出力変動幅に対しても、同様に、その時間窓における時系列の出力変動幅を求める。そして、ループ処理ステップの終端では、基準時間間隔がその終値に達しているか否かを判定する。そして、達していなければ、ループ処理ステップの始端に戻り、増分１０分を基準時間間隔に加算して時間窓を変更し、以後同様に、全風力発電所の総計データ及び各風力発電所の実測データについて、その時間窓における時系列の出力変動幅をそれぞれ求める。最終的にループ処理ステップの終端で基準時間間隔がその終値に達した場合には、周波数相当分析ステップが実行される。
【００７５】
周波数相当分析ステップでは、ループ処理ステップで求めた、総計データに対する各基準時間間隔についての時間窓における出力変動幅を、基準時間間隔の逆数ｆを用い、ｆを変数とする実測型総出力変動Ｓ（ｆ）として求める。また、ループ処理ステップで求めた、各風力発電所の発電機出力データに対する各基準時間間隔についての時間窓における出力変動幅を、ｆを変数とする実測型個別出力変動Ｓｉ（ｆ）としてそれぞれ求める。
【００７６】
以後、第一例と同様に、コヒーレント、ランダム値計算ステップ、遷移時定数計算ステップ、計算結果出力ステップが順次実行され、遷移時定数Tx、出力フォームが出力装置にそれぞれ表示され、これにより、分析プログラムの第二例の全てのステップが終了する。
【図１】
2010239782.tif 000031
【図１４】
2010239782.tif 000032
【図１５】
2010239782.tif 000033
【図２】
2010239782.tif 000034
【図３】
2010239782.tif 000035
【図４】
2010239782.tif 000036
【図５】
2010239782.tif 000037
【図６】
2010239782.tif 000038
【図７】
2010239782.tif 000039
【図８（Ａ）（Ｂ）】
2010239782.tif 000040
【図９】
2010239782.tif 000041
【図１０】
2010239782.tif 000042
【図１１】
2010239782.tif 000043
【図１２】
2010239782.tif 000044
【図１３】
2010239782.tif 000045