JP2007053836A - 臨界故障除去時間算出方法、臨界故障除去時間算出プログラム、及び記録媒体 - Google Patents

Abstract

【課題】制御対象が故障後に回復可能となる時間と、故障後に回復不可能となる時間との臨界となる臨界故障除去時間を精度良く迅速に求める。
【解決手段】制御対象の臨界状態を示す関数ｆに用いられる、臨界故障除去時間の関数である多次元状態変数ｘ^０と、関数ｆをゼロとする多次元状態変数ｘ^ｍ＋１と、多次元状態変数ｘ^０とｘ^ｍ＋１との間で離散される複数の多次元状態変数ｘ^ｋ（１≦ｋ≦ｍ：ｋは整数）と、ｘ^０乃至ｘ^ｍ＋１の中で相互に隣接する多次元状態変数ｘ^ｋ及びｘ^ｋ＋１の間のユークリッド距離εとを、

が最小となる値に決定し、多次元状態変数ｘ^０から臨界故障除去時間を求めてなる。
【選択図】図１

Description

本発明は、制御対象が故障後に回復可能となる時間と、この制御対象が故障後に回復不可能となる時間と、の臨界となる臨界故障除去時間を求める臨界故障除去時間算出方法、及びこの臨界故障除去時間算出方法を情報処理装置に実施させる臨界故障除去時間算出プログラム、並びにこの臨界故障除去時間算出プログラムを記録した記録媒体に関する。

過渡安定度解析は、一般に、制御対象の非線形現象の過渡安定度を解析するのに有用な周知の方法である。この制御対象が例えば電力系統の場合、過渡安定度解析は、電力系統の運用において安全性を維持するために重要な役割を果たす。

電力系統の故障に起因する擾乱に対してこの電力系統が如何なる動的な応答をするかを調べる方法として、例えば、電力系統の初期状態から該当の状態までを数値積分により逐次計算するシミュレーション法が知られている（例えば、非特許文献１参照。）。

一方、前述した非特許文献１に開示されるように、電力系統のいわゆる過渡エネルギーに基づいて、この電力系統の安定度を評価するエネルギー関数法が知られている。このエネルギー関数法により、多数存在する上記擾乱に対して安定度評価ができる。また、前述した非特許文献１に開示されるように、様々なエネルギー関数法の中で、上記安定度判別の基になる電力系統のいわゆる臨界エネルギーを算出する方法が提案されており、ＢＣＵ法と称されている。このＢＣＵ法により、後述する支配的不安定平衡点（ＣＵＥＰ：controlling unstable equilibrium point）における臨界エネルギーが算出される。

図６を参照しつつ、前述したシミュレーション法及びエネルギー関数法による過渡安定度解析の対象となり得る電力系統の動的挙動の一例について説明する。同図は、制動無しの１機無限大母線系統の非線形現象を位相平面における軌跡で表わした模式図である。同図に例示されるように、上記１機の発電機の位相角δ及び角速度ωの軌跡（電力系統の状態を示す軌跡）には、故障軌跡１における異なる３つの点ＰＢ、ＰＣ、ＰＤをそれぞれ起点とする３種類がある。ここで、故障軌跡１とは、例えば母線近傍に短絡事故が発生した場合の電力系統の状態を示す軌跡を意味し、電力系統が安定状態にある安定平衡点ＰＡを起点とする例えば放物線形状をなすものである。

故障が除去された後に安定状態に戻ることが可能な電力系統の状態を示す軌跡２は、安定平衡点ＰＡの周りで閉じる例えばループ形状をなし、これは発電機の位相角δ及び角速度ωがそれぞれ所定範囲内で振動しながらやがては安定点に落ち着くことに対応する。また、故障が除去された後に安定状態に戻ることが不可能な電力系統の状態を示す軌跡４は、発電機の位相角δ及び角速度ωがともに発散するような形状をなし、これは発電機の同期が保たれなくなり脱調することに対応する。つまり、電力系統の状態は、安定平衡点ＰＡから故障軌跡１に沿って時間変化している際に、その故障が点ＰＢで相対的に早く除去された場合、軌跡２に乗り移ることにより安定状態への回復が可能となる一方、上記故障が点ＰＤで相対的に遅く除去された場合、軌跡４に乗り移ることにより安定状態への回復が不可能となる。

軌跡２及び軌跡４の間の臨界軌跡３は、その数理論上の特異点としての支配的不安定平衡点ＰＥを有するものである。この支配的不安定平衡点ＰＥは数理論上の仮想的な点である。即ち、電力系統の状態は、安定平衡点ＰＡから故障軌跡１に沿って時間変化している際に、その故障が点ＰＣで除去された場合、臨界軌跡３に沿って数理論上無限大の時間をかけて支配的不安定平衡点ＰＥに到達するとされている。

もし前述した安定平衡点ＰＡから点ＰＣまでの時間より短い時間で故障が除去されれば、電力系統の状態は、安定状態に回復可能となる。この安定平衡点ＰＡから点ＰＣまでの時間は、臨界故障除去時間と称されており、前述したシミュレーション法及びエネルギー関数法を用いた過渡安定度解析の１つは、この臨界故障除去時間を電力系統のモデルに基づいて算出することにある。前述した非特許文献１に開示されるように、過渡安定度解析では、臨界軌跡３が数値的に求められて、支配的不安定平衡点ＰＥにおける発電機の位相角δ及び角速度ωとともに、臨界故障除去時間が算出される。

例えば、臨界故障除去時間の算出に前述したシミュレーション法を用いる場合、電力系統のモデルに対して予め臨界故障除去時間を設定し、これにより算出された電力系統の状態が最終的に安定状態に回復するか否かを判別し、この判別結果が前述した臨界状態に対応するものとなるまで上記設定、算出、及び判別を繰り返し実施する。一方、臨界故障除去時間の算出に前述したエネルギー関数法の中の例えばＢＣＵ法を用いる場合、電力系統のモデルから得られるエネルギー曲面上で、前述した点ＰＣ及び点ＰＥ間に対応する位相角δの経路を見いだすべくこの曲面の勾配を算出する。更に、このＢＣＵ法では、前述した非特許文献１に開示されるように、この経路の探索をより効率的に実施可能とするShadowing法（BCU-Shaadowing法）が提案されている。また更に、このＢＣＵ法では、前述した非特許文献１に開示されるように、前述した点ＰＣから点ＰＥまでの臨界軌跡における電力系統の離散的な状態を示す多次元状態変数を一括して求める方法が、本出願人により提案されている。
餘利野直人、斉藤岳志、亀井敬史、及び佐々木博司著、「過渡安定性評価のための新しい解析手法の提案」、電気学会論文誌Ｂ、１２４巻１２号、２００４年、ｐｐ．１４２５−１４３０

しかしながら、制御対象の臨界故障除去時間を算出する場合、前述したシミュレーション法は、制御対象の状態の時間変化を忠実に再現するが故に計算時間がかかるため、例えば故障時における実時間（リアルタイム）解析には向いていない。また、前述したBCU-Shadowing法を含めたエネルギー関数法においても、前述した臨界エネルギーを算出することや、前述したエネルギー曲面上での位相角δの経路を精度良く探索すること等は困難であるとされている。よって、シミュレーション並みの精度を保ちつつリアルタイム解析をすることはやはり困難である。

一方、前述した離散的な多次元状態変数を一括して求める方法は、例えばリミッタ等を制御対象のモデルに含めることができないために、その適用可能なモデルが限定されるという問題点がある。よって、例えば電力系統以外の制御対象の非線形現象の過渡安定度解析に対してこの方法を適用することはできず、これは制御対象の適用範囲を狭めることになる。

本発明は、かかる課題に鑑みてなされたものであり、その目的とするところは、制御対象の臨界故障除去時間を精度良く迅速に算出する故障除去時間算出方法を提供することにある。

前記課題を解決するための発明は、制御対象が故障後に回復可能となる時間と、前記制御対象が前記故障後に回復不可能となる時間と、の臨界となる臨界故障除去時間を求める臨界故障除去時間算出方法であって、前記故障後に回復可能となる前記制御対象の状態と、前記故障後に回復不可能となる前記制御対象の状態と、の臨界となる前記制御対象の臨界状態を示す関数ｆに用いられる、前記臨界故障除去時間の関数である多次元状態変数ｘ^０と、前記関数ｆをゼロとする多次元状態変数ｘ^ｍ＋１と、前記多次元状態変数ｘ^０とｘ^ｍ＋１との間で離散される複数の多次元状態変数ｘ^ｋ（１≦ｋ≦ｍ：ｋは整数）と、前記多次元状態変数ｘ^０乃至ｘ^ｍ＋１の中で相互に隣接する多次元状態変数ｘ^ｋ及びｘ^ｋ＋１の間のユークリッド距離εと、を、

が最小となる値に決定し、前記多次元状態変数ｘ^０から前記臨界故障除去時間を求めてなる。

上式（以下、２乗総和式と称する）における２乗総和の対象を誤差ベクトルμ^ｋと称する場合、多次元状態変数ｘ^０乃至ｘ^ｍ＋１及びユークリッド距離εを、２乗総和式が最小となる値に決定することは、μ^ｋを略ゼロとする（ｍ＋１）個の多元連立方程式の解を求めることに等しい。例えば数値計算を実施して、μ^ｋを厳密にゼロとする（ｍ＋１）個の多元連立方程式の解を求める場合、０からｍまでのｋについて、μ^ｋを厳密にゼロとする等式からｘ^ｋを逐次求めることになる。しかし、数値計算の場合、ｋが大きくなるにつれてその誤差が累積して大きくなり、臨界状態を特徴付ける多次元状態変数であるｘ^ｍ＋１において最大の計算誤差が生じる虞がある。一方、μ^ｋを略ゼロとする（ｍ＋１）個の多元連立方程式の解を求めるべく、２乗総和式が最小となる値にｘ^０乃至ｘ^ｍ＋１及びεを決定すれば、前述した計算誤差を、ｘ^０乃至ｘ^ｍ＋１にわたって分散させることができる。よって、μ^ｋを厳密にゼロとする（ｍ＋１）個の多元連立方程式を直接解く場合に比べて、２乗総和式が最小となる値にｘ^０乃至ｘ^ｍ＋１及びεを決定することは、臨界状態に対するより高精度の計算結果を与えることになり、結果的にこれはｘ^０から求められる臨界故障除去時間の精度の向上につながる。

また、２乗総和式を最小にする場合、以下述べる理由により、例えば制御対象の臨界状態に対する所定の制約条件を示す関数ｇ（ｘ^ｋ）がゼロより大きい等の不等式を考慮することが可能である一方、μ^ｋを厳密にゼロとする（ｍ＋１）個の多元連立方程式の解を求める場合にはこの不等式を考慮できない。つまり、上記多元連立方程式は冗長性を有するものであるため、その変数であるｘ^０乃至ｘ^ｍ＋１及びεの実質的な個数は、方程式の個数（ｍ＋１）よりも少ない。一般に、変数の個数が方程式の個数よりも少ない場合、多元連立方程式では唯一解が求まらない。そこで、上記μ^ｋを厳密にゼロとする（ｍ＋１）個の多元連立方程式に対する唯一解ではない解を求めるためには、最小二乗公式が適用される。しかし、この最小二乗公式においては、不等式を考慮できないとされている。よって、上記μ^ｋを厳密にゼロとする（ｍ＋１）個の多元連立方程式の解を求める場合にはこの不等式を考慮できない。

この所定の制約条件は、広く一般の制御対象の臨界状態に対して付加され得るものである。よって、２乗総和式が最小となる値にｘ^０乃至ｘ^ｍ＋１及びεを決定することにより、例えば電力系統以外の一般の制御対象の臨界故障除去時間を算出できる。

また、かかる臨界故障除去時間算出方法において、前記制御対象は、電力系統であり、前記多次元状態変数ｘ^ｋ（０≦ｋ≦ｍ＋１：ｋは整数）は、前記電力系統における発電機の位相角度を成分とする変数ベクトルであり、前記関数ｆは、前記多次元状態変数ｘ^ｋに対して、当該多次元状態変数の時間微分ｄｘ^ｋ／ｄｔを与える関数ベクトルである、こととしてもよい。

また、かかる臨界故障除去時間算出方法において、前記多次元状態変数ｘ^０は、前記臨界故障除去時間のｎ次式（ｎ≧１：ｎは整数）を含み、前記電力系統の臨界状態における故障除去時の状態を示す変数ベクトルである、こととしてもよい。

また、かかる臨界故障除去時間算出方法において、前記多次元状態変数ｘ^ｍ＋１は、前記関数ｆの支配的不安定平衡点である特異点を示す変数ベクトルである、こととしてもよい。

また、前記課題を解決するための発明は、制御対象が故障後に回復可能となる時間と、前記制御対象が前記故障後に回復不可能となる時間と、の臨界となる臨界故障除去時間を求める臨界故障除去時間算出方法を情報処理装置に実施させる臨界故障除去時間算出プログラムであって、前記故障後に回復可能となる前記制御対象の状態と、前記故障後に回復不可能となる前記制御対象の状態と、の臨界となる前記制御対象の臨界状態を示す関数ｆに用いられる、前記臨界故障除去時間の関数である多次元状態変数ｘ^０と、前記関数ｆをゼロとする多次元状態変数ｘ^ｍ＋１と、前記多次元状態変数ｘ^０とｘ^ｍ＋１との間で離散される複数の多次元状態変数ｘ^ｋ（１≦ｋ≦ｍ：ｋは整数）と、前記多次元状態変数ｘ^０乃至ｘ^ｍ＋１の中で相互に隣接する多次元状態変数ｘ^ｋ及びｘ^ｋ＋１の間のユークリッド距離εと、を、

が最小となる値に決定し、前記多次元状態変数ｘ^０から前記臨界故障除去時間を求めてなる。

また、前記課題を解決するための発明は、制御対象が故障後に回復可能となる時間と、前記制御対象が前記故障後に回復不可能となる時間と、の臨界となる臨界故障除去時間を求める臨界故障除去時間算出方法を実施するための臨界故障除去時間算出プログラムを記録した、情報処理装置により読み取り可能な記録媒体であって、前記臨界故障除去時間算出プログラムは、前記故障後に回復可能となる前記制御対象の状態と、前記故障後に回復不可能となる前記制御対象の状態と、の臨界となる前記制御対象の臨界状態を示す関数ｆに用いられる、前記臨界故障除去時間の関数である多次元状態変数ｘ^０と、前記関数ｆをゼロとする多次元状態変数ｘ^ｍ＋１と、前記多次元状態変数ｘ^０とｘ^ｍ＋１との間で離散される複数の多次元状態変数ｘ^ｋ（１≦ｋ≦ｍ：ｋは整数）と、前記多次元状態変数ｘ^０乃至ｘ^ｍ＋１の中で相互に隣接する多次元状態変数ｘ^ｋ及びｘ^ｋ＋１の間のユークリッド距離εと、を、

が最小となる値に決定し、前記多次元状態変数ｘ^０から前記臨界故障除去時間を求めるプログラムである。

制御対象の臨界故障除去時間を精度良く迅速に算出できる。

図１を参照しつつ、本実施の形態の臨界故障除去時間算出方法について説明する。同図は、本実施の形態の制御対象の状態を示す離散化された多次元状態変数を表示する模式図である。尚、本実施の形態の制御対象の非線形現象は、この制御対象の状態を示す多次元状態変数ｘ及び多次元非線形関数ｆに適用される後記の「式１」の非線型方程式により記述されるものとする。

つまり、本実施の形態では、制御対象の臨界軌跡を「式１」の非線形方程式の解として求め、この解に基づいて臨界故障除去時間を算出する。

「式１」の方程式は時刻に依存する方程式であるため、離散的な時刻ｔ^ｋにおける離散的な解をｘ^ｋ（図１）とすると、いわゆる台形公式の近似により、後記の「式２」の等式が成立する。但し、ｋは、上記離散的な時刻の経過とともに制御対象の状態が断続的に移行する際の移行番号を示す整数である。

ここで、「式１」の方程式に基づき、ｆ（ｘ^ｋ）はｘ^ｋの時間微分ｄｘ^ｋ／ｄｔに等しいものとする。本実施の形態の離散的なｘ^ｋは、「式１」の方程式における連続的な多次元状態変数ｘに対応する離散的なベクトルを表わし、離散的な関数ｆ（ｘ^ｋ）は、「式１」の方程式における連続的な多次元関数ｆ（ｘ）に対応する離散的な関数ベクトルを表わす。尚、以後、ｘ^ｋを変数ベクトルと称することとする。

図１に例示されるように、ｘ^０（ｋ＝０）を、臨界故障除去時間（ＣＣＴ：critical clearing time）に対応する変数ベクトルとし、ｘ^ｕ（＝ｘ^ｍ＋１、但しｍは整数）を支配的不安定平衡点ＣＵＥＰにおける変数ベクトルとすれば、「式２」の等式におけるｘ^ｋ（０≦ｋ≦ｍ＋１）は、臨界故障除去時から、制御対象の状態が支配的不安定平衡状態に至るまでの臨界軌跡を構成するものとなる。尚、上記ｘ^０は、例えば図６に例示された故障軌跡１と臨界軌跡３との交点ＰＣに対応するものである。

但し、支配的不安定平衡点ＣＵＥＰに近づくにつれて「式２」の等式の右辺の時間差分（ｔ^ｋ＋１−ｔ^ｋ）は無限に大きくなるため、本実施の形態では、上記無限大を回避するために、図１に例示された隣接する２つの変数ベクトル間のユークリッド距離εを後記の「式３」の等式に基づいて定義しておく。

「式３」の等式に基づいて、上記無限大となる時間の項の代わりに等間隔のユークリッド距離εを用いることにより、「式２」の等式は後記の「式４」の等式に変換される。このユークリッド距離εは、支配的不安定平衡点ＣＵＥＰに近づいても無限大になることなく常に一定である。

ここで、「式４」の等式の右辺の０はゼロベクトルを意味し、よってこの等式は（ｍ＋１）個の多次元（ベクトルの次元）の連立方程式を意味することになる。つまり、「式４」の多元連立方程式の解を求めることは、無限大の時間を直接取り扱うことなく制御対象の臨界軌跡を構成する等間隔の点（図１）を求めることと等価となる。

本実施の形態の臨界故障除去時間算出方法では、「式４」の多元連立方程式の解を求めるために、この「式４」における左辺を誤差ベクトルとして、この誤差ベクトルの大きさ（ノルム）の総和を最小とするｘ^ｋ（０≦ｋ≦ｍ＋１）及びεを決定する後記「式５」の操作を実施する。

但し、「式５」の操作は、変数ベクトルｘ^０が臨界故障除去時間τに基づくことを示す後記の「式６」の等式と、変数ベクトルｘ^ｍ＋１における関数ｆ（即ち、ｄｘ^ｍ＋１／ｄｔ）は支配的不安定平衡点ＣＵＥＰに対応するゼロベクトルであることを示す後記の「式７」の等式と、からなる束縛条件の下に実施される。

本実施の形態の臨界故障除去時間算出方法では、上記操作により得られたｘ^０に「式６」の等式を適用することにより、臨界故障除去時間τが算出される。尚、上記操作にあたって、関数ｆにはいわゆる臨界状態に対応するエネルギーが与えられているものとする（臨界状態を示す関数）。

「式４」の多元連立方程式に基づけば、（ｍ＋１）個の誤差ベクトルは全てゼロとなることが理想的であるが、実際には、前述した台形公式の近似に起因する数値誤差が生じ、これが「式５」に例示される誤差ベクトルを形成する。そこで、本実施の形態では、「式５」の操作、即ち誤差ベクトルの大きさ（ノルム）の総和を最小とするｘ^ｋ（０≦ｋ≦ｍ＋１）及びεを求めることにより、臨界故障除去時から、制御対象の状態が支配的不安定平衡状態に至るまでの臨界軌跡に沿って誤差を分散させている。一方、もし、「式４」の多元連立方程式を、臨界故障除去時から、制御対象の状態が支配的不安定平衡状態に至るまで逐次解く場合、前述した誤差が逐次累積して支配的不安定平衡点ＣＵＥＰにおいて最大の誤差が生じる虞がある。よって、「式４」の多元連立方程式を直接解く場合に比べて、「式５」の操作を実施することは、臨界軌跡に対するより高精度の計算結果を与えることになり、結果的にこれは「式６」の等式により求められる臨界故障除去時間τの精度の向上につながる。

また、「式５」の操作を実施する場合、後記の「式８」の不等式による制約条件を更なる束縛条件として付加することが可能であるのに対し、「式４」の多元連立方程式を直接解くことにより前述した離散的な多次元状態変数を一括して求める方法では、このような束縛条件を考慮できない。

ここで、上記ｆ（ｘ）が例えば電力系統のモデルに対応する関数である場合には、ｇ（ｘ）は、このモデルに付加される発電機のリミッタ等を意味する関数である。

以下、「式４」の多元連立方程式を直接解く場合、「式８」の不等式による束縛条件を考慮できない理由を述べる。「式４」の多元連立方程式は冗長性を有するものであるため、その変数であるｘ^０〜ｘ^ｕ及びεの実質的な個数は、「式４」の多元連立方程式の個数（ｍ＋１）よりも少ない。一般に、変数の個数が方程式の個数よりも少ない場合、多元連立方程式では唯一解が求まらない。そこで、「式４」の多元連立方程式に対する唯一解ではない解を求めるためには、最小二乗公式が適用される。しかし、この最小二乗公式においては、不等式を考慮できないとされている。よって、「式４」の多元連立方程式を直接解く場合、「式８」の不等式による束縛条件を考慮できない。

以上から、「式６」及び「式７」の等式及び「式８」の不等式による束縛条件下で「式５」の操作を実施することは、広く電力系統以外の制御対象の非線形現象の過渡安定度解析に対しても適用可能となる。

＝＝＝１機無限大母線系統＝＝＝
前述した臨界故障除去時間算出方法を、１機無限大母線系統のモデルの非線形現象へ適用する適用例について説明する。

図２に例示されるように、本実施の形態の１機無限大母線系統１０のモデルは、自動電圧調整器（ＡＶＲ：Automatic Voltage Regulator）１００、ガバナ２００、及び制動を伴った同期発電機３００をそれぞれ示すモデルを備えて構成されている。同図は、本実施の形態の１機無限大母線系統のモデルを示す概念図である。このモデルによれば、前述した「式１」の方程式における変数ベクトル（多次元状態変数）ｘの各成分は、後記の「式９」〜「式１２」における内部位相角δ、角速度ω、過渡リアクタンス背後電圧Ｅ、及び機械入力Ｐ_Ｍとなる。

ここで、「式９」〜「式１２」の各等式の右辺が前述した関数ｆに相当する。また、ω_ｓは基準角速度を示し、Ｍは慣性定数を示し、Ｖ_ｂは無限大母線電圧を示し、ｄは制動係数を示し、Ｔ（Ｔ_ＡＶＲ、Ｔ_ＧＯＶ）は時定数を示し、Ｋ（Ｋ_ＡＶＲ、Ｋ_ＧＯＶ）はゲインを示し、Ｖ_ｔは端子電圧を示し、Ｖ_ｒｅｆは基準端子電圧を示し、Ｐ_Ｍｒｅｆは基準機械入力を示し、Ｘは過渡リアクタンス及び線路リアクタンスの合計を示すものである。

本実施の形態のモデルでは、時刻ｔ＝０において３相地絡故障のモデル（３ＬＧ）４００を想定し、時刻ｔ＝τにおいて２つの送電線５００、６００のうちの１つを開放することにより故障を除去するものとする。また、本実施の形態のモデルに対して、前述したシミュレーション法等を適用することにより、前述した故障軌跡が予め求められて、「式６」の束縛条件における変数ベクトルｘ^０は、例えばτの２次関数（ｎ次式）で記述され得るものとする。更に、「式９」〜「式１２」の各等式における右辺を０とし、「式７」の束縛条件に基づいて、ｘ^ｕ（＝ｘ^ｍ＋１）に相当する内部位相角δ、角速度ω、過渡リアクタンス背後電圧Ｅ、及び機械入力Ｐ_Ｍが予め求められているものとする。

本実施の形態では、「式９」〜「式１２」に示される関数ｆ及び変数ベクトルｘがｘ^０〜ｘ^ｕに離散化され、予め求められたｘ^０及びｘ^ｕを代入した誤差ベクトルについて、「式５」の操作が実施された。前述したように、ｘ^０はτの関数であるため、最小化のために最適化されるパラメータはｘ^０の代わりにτとされる一方、ｘ^ｕは予め求められたベクトルに固定された。つまり、最小化に際して、τ、ｘ^１〜ｘ^ｍ、及びεが最適化された。尚、最小化には、いわゆるニュートン法が適用された。即ち、τ、ｘ^１〜ｘ^ｍ、及びεの初期値から前述した誤差ベクトルの大きさ（ノルム）の総和が算出され、この総和が所定値に収束するまで、所定の繰り返し操作が実施された。また、前述した整数ｍを１から１００まで変化させつつ、同様の最小化の操作が実施された。例えばこのｍをより大きくすることは、前述した臨界軌跡の分割数をより多くして計算精度を向上させることを意味する。但し、ｍが大きくなれば計算時間はより長くなる。

表１において、異なる整数ｍ毎に、上記最小化のために最適化されたτ（臨界故障除去時間）（秒）と、上記繰り返し回数と、後述する所定の情報処理装置による計算時間（秒）と、を示す。一方、表２において、前述したシミュレーション法により算出された臨界故障除去時間（秒）及び後述する所定の情報処理装置による計算時間（秒）を、このシミュレーション法における数値積分の刻み毎に示す。例えばこの刻みをより小さくすることは、シミュレーション法の計算精度を向上させることを意味する。

表１に示される臨界故障除去時間は、比較的小さなｍについても、表２に示される最も精度の高い臨界故障除去時間1.3795-1.3800秒を0.01秒程度の誤差で再現していることがわかった。しかも、その計算時間は、シミュレーション法による計算時間に対するおよそ１０分の１であることがわかった。つまり、１機無限大母線系統に適用する場合、本実施の形態の臨界故障除去時間算出方法は、シミュレーション法並みの精度を保ちつつリアルタイム解析が可能な程度まで高速であると言える。

＝＝＝ＡＦ９及びＩＥＥＥ３０＝＝＝
前述した臨界故障除去時間算出方法を、Anderson & Fouadの３機９母線系統（ＡＦ９）のモデル（図３参照）及びIEEE６機３０母線系統（ＩＥＥＥ３０）のモデル（図４参照）の非線形現象へ適用する適用例について説明する。尚、図３は本実施の形態のＡＦ９のモデルを示す概念図であり、図４は本実施の形態のＩＥＥＥ３０のモデルを示す概念図である。

一般に、ｎ機の発電機を備えた電力系統のモデルによれば、前述した「式１」の方程式における変数ベクトル（多次元状態変数）ｘの成分は、後記の「式１３」における内部位相角δ_ｉ（１≦ｉ≦ｎ）となる。尚、このｉは、ｎ機の発電機のうちの第ｉ番目にあたる第ｉ発電機を示す整数である。

ここで、「式１３」のｎ個の等式の右辺が前述した関数ｆに相当する。また、Ｐ_ｍ（ｉ）は第ｉ発電機の機械入力を示し、Ｐ_ｅ（ｉ）は第ｉ発電機の電気入力を示し、Ｍ_ｉは第ｉ発電機の慣性定数を示し、Ｍ_ＴはＭ_ｉの総和を示すものである。但し、Ｐ_ｅ（ｉ）は第ｉ発電機の位相角δ_ｉから後記の「式１４」の等式に基づいて一義的に決定されるものである。

ここで、Ｙ_ｉｊは第ｉ発電機及び第ｊ発電機のアドミタンス行列の各要素を示し、Ｅ_ｉ及びＥ_ｊはそれぞれ第ｉ発電機及び第ｊ発電機の端子電圧を示し、α_ｉｊは第ｉ発電機及び第ｊ発電機に特有の定数を示すものである。

ところで、「式１３」のｎ個の等式の右辺は、本実施の形態の電力系統のモデルに対するポテンシャルエネルギーＶｐ（δ_ｉ）を各発電機の位相角度δ_ｉにより偏微分して得たｎ次元ベクトルに相当するものである。

本実施の形態のモデルでは、時刻ｔ＝０において複数の母線のうちの所定の母線近傍の３相地絡故障を想定し、時刻ｔ＝τにおいて２つの送電線のうちの１つを開放することにより故障を除去するものとする。また、前述したポテンシャルエネルギーＶｐ（δ_ｉ）の曲面上の第１極大点が予め求められ、「式６」の束縛条件における変数ベクトルｘ^０は、この第１極大点から推定されるものとする。更に、「式１３」のｎ個の等式における右辺を０とし、「式７」の束縛条件は、ｘ^ｕ（＝ｘ^ｍ＋１）に相当する第ｉ発電機の内部位相角δ_ｉ ^ｕに関する後記の「式１５」の等式により表わされる。

本実施の形態では、「式１３」に示される関数ｆ及び変数ベクトルｘがｘ^０〜ｘ^ｕに離散化され、「式１５」の等式の条件下で、前述した誤差ベクトルについて「式５」の操作が実施された。つまり、最小化に際して、ｘ^０〜ｘ^ｕ及びεが最適化された。尚、最小化には、いわゆるニュートン法が適用された。即ち、前述したｘ^０及びｘ^ｕの初期値に加えて、ｘ^１〜ｘ^ｍ及びεの初期値から前述した誤差ベクトルの大きさ（ノルム）の総和が算出され、この総和が所定値に収束するまで、所定の繰り返し操作が実施された。

＜＜＜Anderson & Fouadの３機９母線系統（ＡＦ９）＞＞＞
表３において、Anderson & Fouadの３機９母線系統（ＡＦ９）のモデルにおける異なる故障点毎に、前述した従来のBCU-Shadowing法により算出された臨界故障除去時間（秒）と、上記最小化のために最適化されたｘ^０から求められたτ（臨界故障除去時間）（秒）と、上記繰り返し回数と、を示す。

表３によれば、例えば、第１母線のＡ点（図３）で３相地絡故障が発生した場合、第１〜第４母線に対応する（各２本のうちの１本の）送電線を開放することにより故障を除去する臨界故障除去時間は0.33秒と算出された。この結果は、従来のBCU-Shadowing法により得られた結果と同じである。同様に、故障地点がＢ点〜Ｆ点、Ｈ点、及びＩ点（図３）である場合にも、本実施の形態の臨界故障除去時間算出方法により算出された臨界故障除去時間は、従来のBCU-Shadowing法により算出された臨界故障除去時間と一致している。尚、以上の本実施の形態の臨界故障除去時間算出方法は、前述した臨界軌跡の分割数を与えるｍを２に設定して実施されたものである。

また、表３によれば、第７母線のＧ点（図３）で３相地絡故障が発生した場合、第７及び第８母線に対応する（各２本のうちの１本の）送電線を開放することにより故障を除去する臨界故障除去時間は0.23秒と算出される一方、従来のBCU-Shadowing法では、支配的不安定平衡点ＣＵＥＰの探索に失敗し、よって臨界故障除去時間が算出できなかった。尚、本実施の形態の臨界故障除去時間算出方法は、前述した臨界軌跡の分割数を与えるｍを１８に設定して実施されたものである。

以上から、Ｇ点における結果に関しては、本実施の形態の臨界故障除去時間算出方法は、計算精度の点では、従来のBCU-Shadowing法よりも優れていると言える。

図５（ａ）は、前述したＡＦ９のモデルにおいて第２母線のＢ点で３相地絡故障が発生した場合の前述したポテンシャルエネルギーＶｐ（δ_１、δ_２、δ_３）をδ_１−δ_２面上の等高線として示した模式図である。同図に例示されるように、従来のBCU-Shadowing法では、前述した第１極大点から支配的不安定平衡点ＣＵＥＰまでの経路（図中の蛇行線）の探索に成功している。つまり、本実施の形態の臨界故障除去時間算出方法で得られた経路（図中の曲線）とは途中経路は異なるが、その端点である支配的不安定平衡点ＣＵＥＰは一致している。

一方、図５（ｂ）は、前述したＡＦ９のモデルにおいて第７母線のＧ点で３相地絡故障が発生した場合の前述したポテンシャルエネルギーＶｐ（δ_１、δ_２、δ_３）をδ_１−δ_２面上の等高線として示した模式図である。同図に例示されるように、従来のBCU-Shadowing法では、前述した第１極大点から支配的不安定平衡点ＣＵＥＰまでの経路（図中の蛇行線）の探索に失敗している。

以上のように、支配的不安定平衡点ＣＵＥＰにおいてポテンシャルエネルギーＶｐ（δ_１、δ_２、δ_３）の等高線が明確に鞍点の形状をなしている場合（図５（ａ））には、従来のBCU-Shadowing法でも支配的不安定平衡点ＣＵＥＰが求められるが、この等高線が鞍点ではなく複雑な形状をなしている場合（図５（ｂ））には、従来のBCU-Shadowing法では支配的不安定平衡点ＣＵＥＰが求められないが、本実施の形態の臨界故障除去時間算出方法によれば支配的不安定平衡点ＣＵＥＰが求められることがわかった。つまり、本実施の形態の臨界故障除去時間算出方法は、複雑な形状をなすポテンシャルエネルギーＶｐ（δ_ｉ）を有する電力系統モデルに関しては、従来のBCU-Shadowing法よりも優れている（高精度である）と言える。

＜＜＜IEEE６機３０母線系統（ＩＥＥＥ３０）＞＞＞
表４において、IEEE６機３０母線系統（ＩＥＥＥ３０）のモデルにおける異なる故障点毎に、前述した従来のBCU-Shadowing法により算出された臨界故障除去時間（秒）と、上記最小化のために最適化されたｘ^０から求められたτ（臨界故障除去時間）（秒）と、上記繰り返し回数と、を示す。

表４によれば、例えば、第２母線の近傍（図４）で３相地絡故障が発生した場合、第２〜第４母線に対応する（各２本のうちの１本の）送電線を開放することにより故障を除去する臨界故障除去時間は0.75秒と算出された。この結果は、従来のBCU-Shadowing法により得られた結果と同じである。同様に、故障地点が第５、第１１、及び第１３母線の近傍（図４）である場合にも、本実施の形態の臨界故障除去時間算出方法により算出された臨界故障除去時間は、従来のBCU-Shadowing法により算出された臨界故障除去時間と一致している。尚、以上の本実施の形態の臨界故障除去時間算出方法は、前述した臨界軌跡の分割数を与えるｍを２に設定して実施されたものである。

また、表４によれば、第１及び第８母線の近傍（図３）で３相地絡故障が発生した場合、臨界故障除去時間は0.82秒及び0.84秒とそれぞれ算出される一方、従来のBCU-Shadowing法では、支配的不安定平衡点ＣＵＥＰの探索に失敗し、よって臨界故障除去時間が算出できなかった。尚、以上の本実施の形態の臨界故障除去時間算出方法は、前述した臨界軌跡の分割数を与えるｍを２に設定して実施されたものである。

以上から、本実施の形態の臨界故障除去時間算出方法は、前述した第１及び第８母線の近傍における結果に関しては、従来のBCU-Shadowing法よりも優れている（高精度である）と言える。

＝＝＝情報処理装置及びプログラム＝＝＝
前述した実施の形態の臨界故障除去時間算出方法は、所定のＣＰＵ及び記憶部を備える情報処理装置を用いて実施される。この記憶部には、例えば「式９」〜「式１２」、又は、「式１３」及び「式１４」に示される関数ｆの情報を、電力系統データとして記憶するための容量が必要である。また、この記憶部には、前述した変数ベクトルｘ^０〜ｘ^ｕや誤差ベクトル等の情報を、作業用データとして記憶するための容量も必要である。更に、この記憶部は、例えば、前述したニュートン法のプログラムや、このニュートン法を用いてｘ^０〜ｘ^ｕ及びεの最適化を実施するプログラム等からなるプログラム又はプログラム群（臨界故障除去時間算出プログラム）を記憶している必要がある。尚、この記憶部は、本実施の形態の情報処理装置の内部又は外部に備えられるものである。ここで、外部に備えるとは、この情報処理装置と記憶部とが接続可能であることを意味するものである。また、上記プログラム又はプログラム群は、例えば、ＤＶＤやＣＤ等の光ディスク、或いは、ＭＯやフロッピーディスク等の磁気ディスクといった、情報処理装置により読み取り可能な記録媒体に記録されたものであってもよい。

本実施の形態の情報処理装置は、例えば電力系統の運用に携わる作業者が前述した初期値等を入力するためのキーボード及びマウスや、図５に例示される等高線図等を閲覧するためのディスプレイ等を更に備えていてもよい。

また、本実施の形態の情報処理装置は、例えば電力系統を安定化する系統安定化装置に組み込まれたものであってもよいし、或いは、独立したパーソナルコンピュータやワークステーション等であってもよい。尚、系統安定化装置は、例えば特開平１０−２８３２６号公報に開示された周知の装置である。

前述した実施の形態は、本発明の理解を容易にするためのものであり、本発明を限定して解釈するためのものではない。本発明は、その趣旨を逸脱することなく変更、改良されるとともに、本発明にはその等価物も含まれる。

前述した実施の形態の臨界故障除去時間算出方法は、電力系統の過渡安定度解析として用いられるものであったが、これに限定されるものではない。この臨界故障除去時間算出方法は、電力系統に限らず、一般に制御対象の非線形現象の過渡安定度を解析するのにも有用である。本発明によれば、一般の制御対象の臨界故障除去時間についても、精度良く迅速に算出できる。

本実施の形態の離散化された多次元状態変数を表示する模式図である。 本実施の形態の１機無限大母線系統のモデルを示す概念図である。 本実施の形態のＡＦ９のモデルを示す概念図である。 本実施の形態のＩＥＥＥ３０のモデルを示す概念図である。 本実施の形態の電力系統のモデルにおいて故障が発生した場合のポテンシャルエネルギーＶｐ（δ_１、δ_２、δ_３）をδ_１−δ_２面上の等高線として示した模式図である。 制動無しの１機無限大母線系統の非線形現象を位相平面における軌跡で表わした模式図である。

符号の説明

１ 故障軌跡
２ 故障が除去された後に安定状態に戻ることが可能な電力系統の状態を示す軌跡
３ 臨界軌跡
４ 故障が除去された後に安定状態に戻ることが不可能な電力系統の状態を示す軌跡
１０ １機無限大母線系統
１００ 自動電圧調整器（ＡＶＲ）
２００ ガバナ
３００ 同期発電機
４００ ３相地絡故障のモデル
５００、６００ 送電線

Claims (12)

1. 制御対象が故障後に回復可能となる時間と、前記制御対象が前記故障後に回復不可能となる時間と、の臨界となる臨界故障除去時間を求める臨界故障除去時間算出方法であって、
   前記故障後に回復可能となる前記制御対象の状態と、前記故障後に回復不可能となる前記制御対象の状態と、の臨界となる前記制御対象の臨界状態を示す関数ｆに用いられる、
   前記臨界故障除去時間の関数である多次元状態変数ｘ^０と、
   前記関数ｆをゼロとする多次元状態変数ｘ^ｍ＋１と、
   前記多次元状態変数ｘ^０とｘ^ｍ＋１との間で離散される複数の多次元状態変数ｘ^ｋ（１≦ｋ≦ｍ：ｋは整数）と、
   前記多次元状態変数ｘ^０乃至ｘ^ｍ＋１の中で相互に隣接する多次元状態変数ｘ^ｋ及びｘ^ｋ＋１の間のユークリッド距離εと、を、
   

   

   が最小となる値に決定し、前記多次元状態変数ｘ^０から前記臨界故障除去時間を求めることを特徴とする臨界故障除去時間算出方法。
2. 前記制御対象は、電力系統であり、
   前記多次元状態変数ｘ^ｋ（０≦ｋ≦ｍ＋１：ｋは整数）は、前記電力系統における発電機の位相角度を成分とする変数ベクトルであり、
   前記関数ｆは、前記多次元状態変数ｘ^ｋに対して、当該多次元状態変数の時間微分ｄｘ^ｋ／ｄｔを与える関数ベクトルである、
   ことを特徴とする請求項１に記載の臨界故障除去時間算出方法。
3. 前記多次元状態変数ｘ^０は、前記臨界故障除去時間のｎ次式（ｎ≧１：ｎは整数）を含み、前記電力系統の臨界状態における故障除去時の状態を示す変数ベクトルである、ことを特徴とする請求項２に記載の臨界故障除去時間算出方法。
4. 前記多次元状態変数ｘ^ｍ＋１は、前記関数ｆの支配的不安定平衡点である特異点を示す変数ベクトルである、ことを特徴とする請求項２又は３に記載の臨界故障除去時間算出方法。
5. 制御対象が故障後に回復可能となる時間と、前記制御対象が前記故障後に回復不可能となる時間と、の臨界となる臨界故障除去時間を求める臨界故障除去時間算出方法を情報処理装置に実施させる臨界故障除去時間算出プログラムであって、
   前記故障後に回復可能となる前記制御対象の状態と、前記故障後に回復不可能となる前記制御対象の状態と、の臨界となる前記制御対象の臨界状態を示す関数ｆに用いられる、
   前記臨界故障除去時間の関数である多次元状態変数ｘ^０と、
   前記関数ｆをゼロとする多次元状態変数ｘ^ｍ＋１と、
   前記多次元状態変数ｘ^０とｘ^ｍ＋１との間で離散される複数の多次元状態変数ｘ^ｋ（１≦ｋ≦ｍ：ｋは整数）と、
   前記多次元状態変数ｘ^０乃至ｘ^ｍ＋１の中で相互に隣接する多次元状態変数ｘ^ｋ及びｘ^ｋ＋１の間のユークリッド距離εと、を、
   

   

   が最小となる値に決定し、前記多次元状態変数ｘ^０から前記臨界故障除去時間を求めることを特徴とする臨界故障除去時間算出プログラム。
6. 前記制御対象は、電力系統であり、
   前記多次元状態変数ｘ^ｋ（０≦ｋ≦ｍ＋１：ｋは整数）は、前記電力系統における発電機の位相角度を成分とする変数ベクトルであり、
   前記関数ｆは、前記多次元状態変数ｘ^ｋに対して、当該多次元状態変数の時間微分ｄｘ^ｋ／ｄｔを与える関数ベクトルである、
   ことを特徴とする請求項５に記載の臨界故障除去時間算出プログラム。
7. 前記多次元状態変数ｘ^０は、前記臨界故障除去時間のｎ次式（ｎ≧１：ｎは整数）を含み、前記電力系統の臨界状態における故障除去時の状態を示す変数ベクトルである、ことを特徴とする請求項６に記載の臨界故障除去時間算出プログラム。
8. 前記多次元状態変数ｘ^ｍ＋１は、前記関数ｆの支配的不安定平衡点である特異点を示す変数ベクトルである、ことを特徴とする請求項６又は７に記載の臨界故障除去時間算出プログラム。
9. 制御対象が故障後に回復可能となる時間と、前記制御対象が前記故障後に回復不可能となる時間と、の臨界となる臨界故障除去時間を求める臨界故障除去時間算出方法を実施するための臨界故障除去時間算出プログラムを記録した、情報処理装置により読み取り可能な記録媒体であって、
   前記臨界故障除去時間算出プログラムは、
   前記故障後に回復可能となる前記制御対象の状態と、前記故障後に回復不可能となる前記制御対象の状態と、の臨界となる前記制御対象の臨界状態を示す関数ｆに用いられる、
   前記臨界故障除去時間の関数である多次元状態変数ｘ^０と、
   前記関数ｆをゼロとする多次元状態変数ｘ^ｍ＋１と、
   前記多次元状態変数ｘ^０とｘ^ｍ＋１との間で離散される複数の多次元状態変数ｘ^ｋ（１≦ｋ≦ｍ：ｋは整数）と、
   前記多次元状態変数ｘ^０乃至ｘ^ｍ＋１の中で相互に隣接する多次元状態変数ｘ^ｋ及びｘ^ｋ＋１の間のユークリッド距離εと、を、
   

   

   が最小となる値に決定し、前記多次元状態変数ｘ^０から前記臨界故障除去時間を求めるプログラムである、
   ことを特徴とする記録媒体。
10. 前記制御対象は、電力系統であり、
    前記多次元状態変数ｘ^ｋ（０≦ｋ≦ｍ＋１：ｋは整数）は、前記電力系統における発電機の位相角度を成分とする変数ベクトルであり、
    前記関数ｆは、前記多次元状態変数ｘ^ｋに対して、当該多次元状態変数の時間微分ｄｘ^ｋ／ｄｔを与える関数ベクトルである、
    ことを特徴とする請求項９に記載の記録媒体。
11. 前記多次元状態変数ｘ^０は、前記臨界故障除去時間のｎ次式（ｎ≧１：ｎは整数）を含み、前記電力系統の臨界状態における故障除去時の状態を示す変数ベクトルである、ことを特徴とする請求項１０に記載の記録媒体。
12. 前記多次元状態変数ｘ^ｍ＋１は、前記関数ｆの支配的不安定平衡点である特異点を示す変数ベクトルである、ことを特徴とする請求項１０又は１１に記載の記録媒体。
    
    

Images