(57)【要約】      （修正有）【課題】複数の発電機が連系した電力系統の臨界故障除去時間を精度良く迅速に算出する。【解決手段】故障除去時間τの関数であり、故障を除去した時の電力系統の状態を表す多次元状態変数ｘ０と、ｘ０を始点とする軌跡の終点を示す多次元状態変数ｘｍ＋１と、ｘ０とｘｍ＋１との間で離散化され、2010246289.tif 000084として定義される複数の多次元状態変数ｘｋ（１≦ｋ≦ｍ）と、ｘ０乃至ｘｍ＋１の中で相互に隣接するｘｋ及びｘｋ＋１の間のユークリッド距離εと、を定義し、2010246289.tif 000085として定義される誤差ベクトルμｋのノルムの総和として定義される目的関数に対し、制約条件の下で、目的関数を最小化する多次元状態変数ｘ０に基づいて、臨界となる故障除去時間τを求める。【選択図】図２
2010246289.tif 000002
【特許請求の範囲】
【請求項１】
複数の発電機が連系した電力系統が故障した後に回復可能となる時間と、前記電力系統が故障した後に回復不可能となる時間と、の臨界となる故障除去時間を求める臨界故障除去時間算出方法であって、
故障除去時間τの関数であり、前記故障を除去した時の前記電力系統の状態を表す多次元状態変数ｘ０と、
前記電力系統の状態の時間的変化を表した前記多次元状態変数ｘ０を始点とする軌跡の終点を示す多次元状態変数ｘｍ＋１（ｍは整数）と、
前記多次元状態変数ｘ０とｘｍ＋１との間で離散化され、
2010246289.tif 000071
として定義される電力系統方程式に従う複数の多次元状態変数ｘｋ（１≦ｋ≦ｍ＋１：ｋ、ｍは整数）と、
前記多次元状態変数ｘ０乃至ｘｍ＋１の中で相互に隣接する多次元状態変数ｘｋ及びｘｋ＋１の間のユークリッド距離εと、
を定義し、前記電力系統方程式、前記多次元状態変数ｘｋ（０≦ｋ≦ｍ＋１）、前記ユークリッド距離εを用いて
2010246289.tif 000072
として定義される誤差ベクトルμｋを用いて
2010246289.tif 000073
として定義される目的関数に対し、
前記多次元状態変数ｘ０から前記多次元状態変数ｘｍ＋１に至る前記軌跡は前記電力系統方程式に基づく所定の関数の特異点を通過する、という制約条件の下で、前記目的関数を最小化し、
前記目的関数を最小化する前記多次元状態変数ｘ０に基づいて、前記臨界となる故障除去時間τを求めること、を特徴とする臨界故障除去時間算出方法。
【請求項２】
請求項１に記載の臨界故障除去時間算出方法であって、
前記特異点は、前記電力系統に連系される発電機ｉ（ｉ＝１〜ｎ）の同期化力係数行列Ｋに基づく特異点としたこと、
を特徴とする臨界故障除去時間算出方法。
【請求項３】
請求項２に記載の臨界故障除去時間算出方法であって、
前記同期化力係数行列Ｋは、発電機ｉの加速電力Ｐａｉと位相角δｉを用いて、
2010246289.tif 000074
として定義されるｎ次元正方行列としたこと、
を特徴とする臨界故障除去時間算出方法。
【請求項４】
請求項２に記載の臨界故障除去時間算出方法であって、
前記同期化力係数行列は、発電機ｉの加速トルクＴａｉと位相角δｉを用いて、
2010246289.tif 000075
として定義されるｎ次元正方行列としたこと、
を特徴とする臨界故障除去時間算出方法。
【請求項５】
請求項３又は４に記載の臨界故障除去時間算出方法であって、
前記特異点は、
前記同期化力係数行列Ｋが非正則行列である第１条件と、
且つ前記同期化力係数行列Ｋのゼロ固有値に相当するｎ次元固有ベクトルｖの方向が角周波数ωｉを配列したｎ次元ベクトルωの方向に一致する第２条件と、
の両条件が成立する特異点としたこと、
を特徴とする臨界故障除去時間算出方法。
【請求項６】
請求項５に記載の臨界故障除去時間算出方法であって、
前記特異点は、
前記第１条件として、ｎ次元正方行列Ｐ及びｎ次元固有ベクトルｖの大きさを表すｎ次元ベクトルｑを用いて、
2010246289.tif 000076
が成立し、且つ前記第２条件として、スカラー係数ｋｓを用いて、
2010246289.tif 000077
が成立する特異点としたこと、
を特徴とする臨界故障除去時間算出方法。
【請求項７】
請求項６に記載の臨界故障除去時間算出方法であって、
前記特異点を前記軌跡の終点である前記多次元状態変数ｘｍ＋１とし、
前記誤差ベクトルμｍ＋１は、前記第１条件及び前記第２条件に基づき
2010246289.tif 000078
として定義され、対角要素が正である正方の対角行列Ｗを用いて
2010246289.tif 000079
として定義される二乗誤差関数を前記目的関数に加算した、
2010246289.tif 000080
として定義される関数を最小化する前記多次元状態変数ｘ０を求めること、
を特徴とする臨界故障除去時間算出方法。
【請求項８】
複数の発電機が連系する電力系統が故障した後に回復可能となる時間と、前記電力系統が故障した後に回復不可能となる時間と、の臨界となる故障除去時間を求めるために情報処理装置に実現させるプログラムであって、
前記情報処理装置に、
故障除去時間τの関数であり、前記故障を除去した時の前記電力系統の状態を表す多次元状態変数ｘ０と、
前記電力系統の状態の時間的変化を表した前記多次元状態変数ｘ０を始点とする軌跡の終点を示す多次元状態変数ｘｍ＋１（ｍは整数）と、
前記多次元状態変数ｘ０とｘｍ＋１との間で離散化され、
2010246289.tif 000081
として定義される電力系統方程式に従う複数の多次元状態変数ｘｋ（１≦ｋ≦ｍ：ｋは整数）と、
前記多次元状態変数ｘ０乃至ｘｍ＋１の中で相互に隣接する多次元状態変数ｘｋ及びｘｋ＋１の間のユークリッド距離εと、
を定義し、前記電力系統方程式、前記多次元状態変数ｘｋ（0≦ｋ≦ｍ+1）、前記ユークリッド距離εを用いて
2010246289.tif 000082
として定義される誤差ベクトルμｋを用いて
2010246289.tif 000083
として定義される目的関数に対し、
前記多次元状態変数ｘ０から前記多次元状態変数ｘｍ＋１に至る前記軌跡は前記電力系統方程式に基づく所定の関数の特異点を通過する、という制約条件の下で、前記目的関数を最小化し、
前記目的関数を最小化する前記多次元状態変数ｘ０に基づいて、前記臨界となる故障除去時間τを求める機能と、
を実現させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、臨界故障除去時間算出方法、プログラムに関する。
【背景技術】
【０００２】
本出願人は、制御対象の非線形現象の過渡安定度を解析する新たな手法として、特許文献１に開示するとおり、従来のシミュレーション法やエネルギー関数法よりも臨界故障除去時間を精度良く迅速に算出する臨界故障除去時間算出方法等を提案した。以下、特許文献１に開示される手法の概要について説明する。
【０００３】
まず、図１を参照して、臨界故障除去時間について説明する。尚、図１は、制動無しの１機無限大母線系統（制御対象の一例）の状態を発電機の位相角δ及び角速度ωの軌跡で表わした模式図である。故障除去後に安定となる場合の軌跡２と不安定となる場合の軌跡４の間にある臨界軌跡３は、支配的不安定平衡点ＰＥ（ＣＵＥＰ：controlling unstable equilibrium point）を有している。電力系統の状態は、安定状態を示す安定平衡点ＰＡから故障軌跡１に沿って時間変化している際、その故障が故障軌跡１上の点ＰＣで除去された場合、臨界軌跡３に沿って数理論上無限大の時間をかけて上記の支配的不安定平衡点ＰＥに到達するとされている。そこで、もし安定平衡点ＰＡから点ＰＣまでの時間より短い時間で故障が除去されれば、電力系統の状態は、安定状態に回復可能となる。この安定平衡点ＰＡから点ＰＣまでの時間が、臨界故障除去時間（ＣＣＴ：critical clearing time）と称されている。
【０００４】
つぎに、図２を参照して、特許文献１に開示された臨界故障除去時間算出方法について説明する。  
まず、故障軌跡１上の点であり且つ故障除去時の制御対象の状態（図１に示す点ＰＣ）を多次元状態変数ｘ０と定義する。尚、多次元状態変数ｘ０は、故障軌跡１上の点であるため、次の（式１）に示される故障除去時間τの関数として表すことができる。  
2010246289.tif 000003
【０００５】
また、多次元状態変数ｘ０を臨界軌跡３の始点とした故障除去後の制御対象の状態を、離散的な時刻ｔｋ(１≦ｋ≦ｍ＋１)の順に多次元状態変数ｘ１，ｘ２，・・ｘｍ，ｘｍ＋１と定義する。尚、多次元状態変数ｘｋ（０≦ｋ≦ｍ＋１）は、それぞれ複数の成分から成る多次元変数（ベクトル）である。例えば、特許文献１の段落番号３５には、多次元状態変数ｘの各成分として、１機無限大母線系統モデルにおける位相角δ、角速度ω、過渡リアクタンス背後電圧Ｅ、機械入力ＰＭ等が例示されている。
【０００６】
ここで、多次元状態変数ｘ０（ｋ＝０）が臨界故障除去時間に対応するベクトルであるものとし、多次元状態変数ｘｍ＋１（ｋ＝ｍ＋１）を支配的不安定平衡点ＣＵＥＰにおけるベクトルであるものとすれば、多次元状態変数ｘｋ（０≦ｋ≦ｍ＋１）は、制御対象の状態が臨界故障除去時の状態（ｘ０）から支配的不安定平衡状態（ｘｍ＋１）に至るまでの臨界軌跡３を構成する。尚、多次元状態変数ｘｋは、制御対象の非線形状態を表現する次の（式２）の多次元非線形方程式（電力系統方程式）の解として捉えることができる。  
2010246289.tif 000004
（式２）に対して台形公式を適用することで、多次元状態変数ｘｋの中で相互に隣接する多次元状態変数ｘｋ及びｘｋ＋１の間は、次の（式３）を満たす。尚、εは、相互に隣接する多次元状態変数ｘｋ及びｘｋ＋１の間のユークリッド距離を表している。  
2010246289.tif 000005
【０００７】
ところで、臨界軌跡３の終点である多次元状態変数ｘｍ＋１は、一般的に支配的不安定平衡点ＣＵＥＰとなることが知られている。そこで、（式３）に関する制約条件としては、臨界軌跡３の終点ｘｍ+1を支配的不安定平衡点ＣＵＥＰの所定値を表したｘｕとして指定する場合には次の（式４−１）となり、臨界軌跡３の終点ｘｍ＋１を支配的不安定平衡点ＣＵＥＰの状態を表した平衡条件として指定する場合には次の（式４−２）となる。  
2010246289.tif 000006
従って、（式１）、（式２）、（式３）及び臨界軌跡３の終点ｘｍ＋１の制約条件である（式４−１）又は（式４−２）による多元連立方程式を解くことによって、臨界軌跡３の始点ｘ０、ひいては臨界故障除去時間τを求めることができる。
【０００８】
尚、上記の（式３）の多元連立方程式を直接的に解くシミュレーション法を採用すると、台形公式の近似に起因する数値誤差が累積して臨界軌跡３の終点ｘｍ＋１で最大化する虞がある。このため、（式３）における左辺を誤差ベクトルとみなして、次の（式５）のように、この誤差ベクトル（ノルム）の総和を最小化するときの未知変数τ，ε，ｘ１，ｘ２，・・ｘｍ，ｘｍ＋１を一括して求めるようにする。  
2010246289.tif 000007
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特開２００７−５３８３６号公報
【発明の概要】
【発明が解決しようとする課題】
【００１０】
ところで、本出願人は、上記の制御対象を複数の発電機が連系する多機電力系統とした場合に、前回提案した上記の手法に従って臨界故障除去時間を算出したところ、誤差の点で更なる改善の必要性を認識することができた。
【００１１】
本発明は、かかる課題に鑑みてなされたものであり、その目的とするところは、複数の発電機が連系した電力系統の臨界故障除去時間を精度良く迅速に算出する故障除去時間算出方法及びプログラムを提供することにある。
【課題を解決するための手段】
【００１２】
前記課題を解決するための主たる本発明は、複数の発電機が連系した電力系統が故障した後に回復可能となる時間と、前記電力系統が故障した後に回復不可能となる時間と、の臨界となる故障除去時間を求める臨界故障除去時間算出方法であって、故障除去時間τの関数であり、前記故障を除去した時の前記電力系統の状態を表す多次元状態変数ｘ０と、前記電力系統の状態の時間的変化を表した前記多次元状態変数ｘ０を始点とする軌跡の終点を示す多次元状態変数ｘｍ＋１（ｍは整数）と、前記多次元状態変数ｘ０とｘｍ＋１との間で離散化され、
2010246289.tif 000008
または
2010246289.tif 000009
として定義される電力系統方程式に従う複数の多次元状態変数ｘｋ（１≦ｋ≦ｍ＋１：ｋ、ｍは整数）と、前記多次元状態変数ｘ０乃至ｘｍ＋１の中で相互に隣接する多次元状態変数ｘｋ及びｘｋ＋１の間のユークリッド距離εと、を定義し、前記電力系統方程式、前記多次元状態変数ｘｋ（０≦ｋ≦ｍ＋１）、前記ユークリッド距離εを用いて
2010246289.tif 000010
として定義される誤差ベクトルμｋを用いて
2010246289.tif 000011
として定義される目的関数に対し、前記多次元常態変数ｘ０から前記多次元状態変数ｘｍ＋１に至る前記軌跡は前記電力系統方程式に基づく所定の関数の特異点を通過する、という制約条件の下で、前記目的関数を最小化し、前記目的関数を最小化する前記多次元状態変数ｘ０に基づいて、前記臨界となる故障除去時間τを求めること、を特徴とする。
【００１３】
また、上記の臨界故障除去時間算出方法であって、前記特異点は、前記電力系統に連系される発電機ｉ（ｉ＝１〜ｎ）の同期化力係数行列Ｋに基づく特異点としてもよい。
【００１４】
また、上記の臨界故障除去時間算出方法であって、
前記同期化力係数行列は、発電機ｉの加速電力Ｐａｉと位相角δｉを用いて、
2010246289.tif 000012
として定義されるｎ次元正方行列としてもよい。
【００１５】
また、上記の臨界故障除去時間算出方法であって、前記同期化力係数行列は、発電機ｉの加速トルクＴａｉと位相角δｉを用いて、
2010246289.tif 000013
として定義されるｎ次元正方行列としてもよい。
【００１６】
また、上記の臨界故障除去時間算出方法であって、前記特異点は、前記同期化力係数行列Ｋが非正則行列である第１条件と、且つ前記同期化力係数行列Ｋのゼロ固有値に相当するｎ次元固有ベクトルｖの方向が角周波数ωｉを配列したｎ次元ベクトルωの方向に一致する第２条件と、の両条件が成立する特異点としてもよい。
【００１７】
また、上記の臨界故障除去時間算出方法であって、前記特異点は、前記第１条件として、ｎ次元固有ベクトルを用いて、
2010246289.tif 000014
が成立し、且つ前記第２条件として、スカラー係数ｋｓを用いて、
2010246289.tif 000015
が成立する特異点としてもよい。但し、ｑはゼロではない任意の値である。
【００１８】
また、上記の臨界故障除去時間算出方法であって、前記特異点を前記軌跡の終点である前記多次元状態変数ｘｍ＋１とし、前記誤差ベクトルμｍ＋１は、前記第１条件及び前記第２条件に基づき
2010246289.tif 000016
として定義され、任意の正の要素を対角要素とする適当な重み行列Ｗを用いて
2010246289.tif 000017
として定義される二乗誤差関数を前記目的関数に加算した、
2010246289.tif 000018
として定義される関数を最小化する前記多次元状態変数ｘ０を求めてもよい（手法１）。
【００１９】
また、第２条件に示すｖを第１条件に代入して消去すると、第１条件と第２条件の連立条件として、
2010246289.tif 000019
が得られ、この第３条件に基づき、
2010246289.tif 000020
として定義され、上記の重み行列Ｗを用いて、
2010246289.tif 000021
として定義される２乗誤差関数を前記目的関数に加算した
2010246289.tif 000022
として定義される関数を最小化する前記多次元状態変数Ｘ０を求めてもよい（手法２）。
【００２０】
また、上記の臨界故障除去時間算出方法であって、前記目的関数を最小化する過程で前記目的関数と前記多次元状態変数ｘｋ（0≦ｋ≦ｍ+1）とを対応づけて記憶してもよい。
【００２１】
また、前記課題を解決する主たるその他の本発明は、複数の発電機が連系する電力系統が故障した後に回復可能となる時間と、前記電力系統が故障した後に回復不可能となる時間と、の臨界となる故障除去時間を求めるために情報処理装置に実現させるプログラムであって、前記情報処理装置に、故障除去時間τの関数であり、前記故障を除去した時の前記電力系統の状態を表す多次元状態変数ｘ０と、前記電力系統の状態の時間的変化を表した前記多次元状態変数ｘ０を始点とする軌跡の終点を示す多次元状態変数ｘｍ＋１（ｍは整数）と、前記多次元状態変数ｘ０とｘｍ＋１との間で離散化され、
2010246289.tif 000023
または
2010246289.tif 000024
として定義される電力系統方程式に従う複数の多次元状態変数ｘｋ（１≦ｋ≦ｍ：ｋは整数）と、前記多次元状態変数ｘ０乃至ｘｍ＋１の中で相互に隣接する多次元状態変数ｘｋ及びｘｋ＋１の間のユークリッド距離εと、を定義し、前記電力系統方程式、前記多次元状態変数ｘｋ（0≦ｋ≦ｍ+1）、前記ユークリッド距離εを用いて
2010246289.tif 000025
として定義される誤差ベクトルμｋを用いて
2010246289.tif 000026
として定義される目的関数に対し、前記多次元常態変数ｘ０から前記多次元状態変数ｘｍ＋１に至る前記軌跡は前記電力系統方程式に基づく所定の関数の特異点を通過する、という制約条件の下で、前記目的関数を最小化し、前記目的関数を最小化する前記多次元状態変数ｘ０に基づいて、前記臨界となる故障除去時間τを求める機能と、を実現させるためのプログラムである。
【発明の効果】
【００２２】
複数の発電機が連系した電力系統の臨界故障除去時間を精度良く迅速に算出できる。
【図面の簡単な説明】
【００２３】
【図１】離散化された多次元状態変数を表示する模式図である。
【図２】制動無しの１機無限大母線系統の非線形現象を位相平面における軌跡で表わした模式図である。
【図３】臨界故障除去時間算出方法を説明するフローチャートである。
【図４】ＡＦ９のモデルを示す概念図である。
【図５】ＩＥＥＥ３０のモデルを示す概念図である。
【図６】情報処理装置の構成を示した図である。
【発明を実施するための形態】
【００２４】
＝＝＝臨界故障除去時間算出方法の概要説明＝＝＝
図２を参照しつつ、本発明の実施の形態の臨界故障除去時間算出方法について説明する。図２では、本発明の実施の形態の多機電力系統の状態を、離散的な時刻ｔｋ(０≦ｋ≦ｍ＋１；但しｍは整数)により離散化された多次元状態変数ｘｋ（０≦ｋ≦ｍ＋１）として表現している。尚、ｘｋは、次の（式６）に示される要素で構成されるものとし、以下では状態変数ベクトルと称することとする。  
2010246289.tif 000027
尚、（式６）において、ｘｇｉｋは、次の（式７）に示される離散的な時刻ｔｋにおける発電機ユニットｉ（ｉ=１〜Ｎ）の状態変数ベクトルを表現している。  
2010246289.tif 000028
（式７）に示される要素に関して、δｉｋは、発電機ユニットｉの時刻ｔｋにおける位相角（スカラー量）を表現しており、ωｉｋは、発電機ユニットｉの時刻ｔｋにおける角周波数（スカラー量）を表現している。尚、δｉｋ並びωｉｋの座標は、慣性中心座標系であってもよいし、それ以外の座標系でも良い。また、ｘＲｉｋは、発電機ユニットｉに関する制御系変数を含んだその他の状態変数ベクトルを表現しており、ｘｓｙｓｋは、上記以外の状態変数を全て含むベクトルを表現している。
【００２５】
尚、発電機ユニットのモデルとして簡略モデル（例えば、Ｘｄ’モデル）を使用する場合には、その状態変数にδとωが必ず含まれることになる。また、一般に、発電機の回転子の運動方程式のことを動揺方程式と呼んでおり、δは発電機の回転子の位相角であり、ωはδの時間ｔによる微分値であるため、簡略モデル以外の詳細モデルの場合であっても、その状態変数の中にはδとωが必ず含まれることになる。
【００２６】
また、ＡＶＲ（Automatic Voltage Regulator）、ガバナ、ＰＳＳ（Power System Stabilizer）等の制御系のモデルに依存して、対応する状態変数（発電機と制御系のセット）が決定される。例えば、発電機ユニットｉに関してＡＶＲとガバナを含む一機無限大母線系統の発電機ユニットのモデルを使用した場合、ｘｇｉｋ、ｘＲｉｋはそれぞれ次の（式８−１）、（式８−２）のようになる。尚、Ｅｉｋは、過渡リアクタンス背後電圧を表しており、Ｐｍｉｋは、機械入力を表している。  
2010246289.tif 000029
また、ｘｓｙｓｋの一例としては、発電機ユニットや変電所ユニット等に設置されるＳＶＣ（Static Var Compensator）やＳＶＧ（Static Var Generator）等の制御器の状態変数が挙げられる。
【００２７】
さらに、本発明の実施の形態の多機電力系統の状態は、（式６）で示される状態変数ベクトルｘｋと多次元関数ｆを用いた次の（式９）の非線型方程式（電力系統方程式）によって表現される。
【００２８】
この点についてさらに述べれば、一般の電力系統方程式は、（式６）に示した状態変数ｘと、従属変数ｙと用いて、
2010246289.tif 000030
のように表すことができる。（式Ａ１）は、多機電力系統に連系される発電機ｉ（ｉ＝１〜ｎ）の位相角δｉとその時間微分である角周波数ωｉの動揺を記述した、動揺方程式を必ず含む。また、（式Ａ２）は、多機電力系統の母線の電流・電圧、あるいは有効電力・無効電力を表現する系統方程式と、構成機器の特性を記述する代数法的式群により構成される。（式Ａ１）はｘの次元に等しい次元の微分方程式であり、（式Ａ２）はｙの次元に等しい次元の代数方程式であるので、（式Ａ２）を変形して、
2010246289.tif 000031
として、ｙを（式Ａ１）に代入すると次式のように表すことができる。  
2010246289.tif 000032
あるいは、これを標準形式でのモデル表現として以下のように表す。  
2010246289.tif 000033
【００２９】
なお、
2010246289.tif 000034
の記号は（ｄｘ／ｄｔ）を表すものとしている。ここで、（式９）は従属変数を消去した（式Ａ１）（式Ａ２）に等価なシステム表現で、制御工学等で標準的に用いられる形式である。なお、（式Ａ１）（式Ａ２）の形は電力系統のモデルに依存して変化するので、全てのモデルに対して（式Ａ３）のように数式で変形できず、式を用いて明示的に即ち解析的に表現できない場合もある。しかし、この場合も、ニュートン法などの数値解析により代数方程式（式Ａ２）を解いてｙｋを数値的に求め、これを(式Ａ１)に代入することで（式９）の値を算出できるほか、システムのヤコビ行列あるいは線形化システム行列Ａは（式Ａ７）あるいは（式Ａ１）（式Ａ２）のいずれの表現でも、
2010246289.tif 000035
あるいは
2010246289.tif 000036
のように算出できることが知られている等、（式Ａ１）と（式Ａ２）で表現されたシステムを（式Ａ７）として取り扱う方法論は理論的にも確立され、広く知られている。
【００３０】
このように、状態変数ベクトルｘｋは離散的な時刻ｔｋ、（０≦ｋ≦ｍ＋１）における一般の状態方程式
2010246289.tif 000037
あるいは、（式Ａ５）および（式Ａ６）と等価な標準形式での状態方程式
2010246289.tif 000038
の離散的な解を表現している。
【００３１】
つまり、本発明の実施の形態の臨界故障除去時間算出方法は、多機電力系統の臨界軌跡３を、（式９）の非線形方程式の解（つまり、状態変数ベクトルｘｋ）として求め、この解に基づいて臨界故障除去時間を算出するものである。さらに、「臨界軌跡３は（式９）の電力系統方程式に基づく所定の関数の特異点Ｓ（当該関数に適用される所定の基準が適用できなくなる点）を通過する」という新たに判明した知見に基づき、（式９）の解を求めるに際して、「臨界軌跡３の終点ｘｍ＋１が上記の特異点Ｓである」という制約条件を課することにした。
【００３２】
具体的には、式（９）で表現される電力系統方程式のうち、発電機ｉ（ｉ＝１〜ｎ）の位相角δｉ及び角周波数ωｉを用いた動揺方程式に注目する。動揺方程式は、多機電力系統のモデルとしてＸｄ’モデル等の簡略モデルを用いる場合には次の（式１０−１）が用いられ、詳細モデルを用いる場合には（式１０−２）が用いられる。  
2010246289.tif 000039
尚、Ｍｉは発電機ｉの慣性定数、Ｐａｉは発電機ｉの加速電力、Ｔａｉは発電機ｉの加速トルク、Ｐｅｉは発電機ｉの出力電力、Ｔｅｉは発電機ｉの出力トルク、Ｐｄｉは発電機ｉの出力電力の制動項、Ｔｄｉは発電機ｉの出力トルクの制動項である。
【００３３】
次に、動揺方程式（式１０−１）又は（式１０−２）の右辺に出現する項を用いて構成される所謂同期化力係数行列Ｋ（位相角δの関数）に基づく特異点Ｓを定義する。尚、発電機の同期化力とは、一般に、位相角δに対する発電機の加速電力Ｐａの微分係数ｄＰａ／ｄδ、又は位相角δに対する発電機の加速トルクＴａの微分係数ｄＴａ／ｄδとして定義される。また、多機電力系統の同期化力係数行列Ｋとは、一般に、各発電機ｉ間の上記同期化力を配列したｎ次元正方行列として定義される。同期化力が大きい場合には電力系統に擾乱が発生しても速やかに平衡状態に復帰するが、同期化力が小さい場合には発電機の位相角δが次第に広がり脱調に至る虞があることが知られている。
【００３４】
そして、臨界軌跡３の終点を前述した前回提案した手法（特開２００７−５３８３６号公報を参照）における臨界軌跡３の終点ｘｍ＋１（ｋ＝ｍ＋１）を所定の定数ｘuに束縛させるという制約条件（式４−１）又は（式４−２）に代えて、終点ｘｍ＋１を発電機の同期化力係数行列Ｋに基づく特異点Ｓに束縛させるという新たな制約条件（式１１）を導入した。  
2010246289.tif 000040
このように、終点ｘｍ＋１の束縛条件として従来の定数ｘuよりも普遍的な特異点Ｓを導入したことにより、最適化問題の制約条件の精度を向上させるとともに解探索の際に用いる目的関数の更なる改善に繋げることができ、臨界故障除去時間τの算出をより高精度に且つより迅速に行うことを可能とした。
【００３５】
＝＝＝臨界故障除去時間算出方法の詳細説明＝＝＝
以下では、本発明の実施の形態に係る臨界故障除去時間算出方法を詳細に説明する。（式９）の非線形方程式に台形公式の近似を適用することによって、次の（式１２）の等式が成立する。  
2010246289.tif 000041
但し、ｋは、上記離散的な時刻ｔｋの経過とともに多機電力系統の状態が断続的に移行する際の移行番号を示す整数を表している。また、ｘｋは、連続的な多次元状態変数ｘの連続的な時間推移における時刻ｔｋに対応する離散的なベクトルを表わしている。  
2010246289.tif 000042
はｘｋの時間微分ｄｘｔ／ｄｔを表記する記号である。（式１２）はいかなる動的システムに対しても適用可能であり、その際
2010246289.tif 000043
はＡｎｄｅｒｓｏｎ＆Ｆｏｕａｄの３機９母線系統（ＡＦ９）、ＩＥＥＥの６機３０母線系統（ＩＥＥＥ３０）等の電力系統モデル毎に、(式１３−１)の標準形式ではベクトル関数ｆ（ｘｋ）を用いて、あるいは(式１３−２)の一般形式では従属変数ｙｋを含むベクトル関数ｈ（ｘｋ、ｙｋ）とベクトル代数方程式により、以下のように定義される。  
2010246289.tif 000044
あるいは
2010246289.tif 000045
(式１３−２)の一般形式で表現されたシステムを、(式１３−１)の標準形式に対応付けてシステム解析を実施する方法論は広く知られているので、以下の説明では、(式１３−２)の一般形式で表現されたシステムの場合を含めて、標準形式を用いて説明する。
【００３６】
図１に例示されるように、ｋ＝０とした場合の変数ベクトルｘ０を、臨界故障除去時間ＣＣＴに対応する変数ベクトルとし、ｋ＝ｍ＋１とした場合の変数ベクトルｘｍ＋１を支配的不安定平衡状態における変数ベクトルとすると、（式１３）の等式を成立させるｘｋは、多機電力系統の状態が臨界故障除去時間の状態（＝ｘ０）から支配的不安定平衡状態（＝ｘｍ＋１）に至るまでの臨界軌跡３を構成することになる。但し、支配的不安定平衡点ＣＵＥＰに近づくにつれて（式１２）の等式の右辺の時間差分（ｔｋ＋１−ｔｋ）は無限に大きくなるため、本発明の実施の形態では、次の（式１４）に示す等式に基づいて、図１に例示された隣接する２つの変数ベクトル間のユークリッド距離εを定義する。  
2010246289.tif 000046
そして、（式１２）の等式の右辺の時間差分（ｔｋ＋１−ｔｋ）の代わりに（式１４）に示されるユークリッド距離εを用いることにより、（式１２）の等式は（式１５）の等式に変換される。尚、このユークリッド距離εは、支配的不安定平衡点ＣＵＥＰに近づいても無限大になることなく常に一定である。  
2010246289.tif 000047
（式１５）の等式の右辺の０は、ゼロベクトルを表しており、（式１５）の等式は、（ｍ＋１）個の多次元連立方程式を意味している。つまり、（式１５）の多元連立方程式の解を求めることは、無限大の時間を直接取り扱うことなく、多機電力系統の臨界軌跡３を構成する等間隔（ユークリッド距離ε）の点ｘｋを求めることと等価となる。
【００３７】
但し、（式１５）に示される（ｍ＋１）個の多元連立方程式の右辺は全てゼロとなることが理想的であるが、実際には、前述した台形公式の近似に起因した数値誤差が生じてくる。そこで、（式１５）の多元連立方程式の解を求める手法として、まず、（式１５）における左辺を誤差ベクトルμｋと定義した上で、（式１６）に示すように、誤差ベクトルμｋの大きさ（ノルム）の総和を目的関数Ｏとして定義する。  
2010246289.tif 000048
尚、μｋ’のダッシュは転置を表す。
【００３８】
そして、（式１５）の多元連立方程式の解を求める問題を、（式１６）の目的関数Ｏを最小化させる変数ベクトルｘｋ、ユークリッド距離ε、故障除去時間τ（つまり臨界故障除去時間τ）を求める最適化問題として次の（式１７）により定式化することにした。尚、ｘ０はτの関数であるため、次の（式１７）により最適化されるパラメータは、未知変数τ、ｘ１〜ｘｍ、及びεである。  
2010246289.tif 000049
【００３９】
尚、（式１７）において、ｘｋはＮ次元の実数体の集合に属しており、ユークリッド距離ε並びに故障除去時間τは実数体の集合に属しているものとする。また、同式に示される誤差ベクトルμｋ（０≦ｋ≦ｍ）は、次の（式１８）のとおり、（式１５）の左辺を表したものである。ちなみに、（式１７）に示される誤差ベクトルμｋ’は、誤差ベクトルμｋを転置したベクトルを表している。  
2010246289.tif 000050
また、（式１７）の解探索を実行するに際し、次の（式１９）及び（式２０）により定義される制約条件を課することにした。  
2010246289.tif 000051
【００４０】
（式１９）は、臨界軌跡３の始点ｘ０を、初期値ｘｐｒｅを起点とした故障軌跡１を表す故障除去時間τの関数ＸＦ（τ；ｘｐｒｅ）に束縛させることを表している。（式２０）は、終点ｘｍ＋１を発電機の同機化力係数行列Ｋに基づく特異点Ｓに束縛させることを表している。尚、同期化力係数行列Ｋは、式（１０−１）の動揺方程式を採用する場合には（式２１−１）が用いられ、式（１０−２）の動揺方程式を採用する場合には（式２１−２）が用いられる。  
2010246289.tif 000052
2010246289.tif 000053
同期化力係数行列Ｋに基づく特異点Ｓは、同期化力係数行列Ｋが非正則行列（特異）である第１条件と、且つ同期化力係数行列Ｋのゼロ固有値に相当するｎ次元固有ベクトルｖの方向が発電機ｉの角周波数ωｉを配列したｎ次元ベクトルωの方向に一致する第２条件と、の両条件が成立する特異点とした。
【００４１】
第１条件が成立する必要十分条件は、数論理上、次の（式２２−１）又は（式２２−２）が成立することである。  
2010246289.tif 000054
しかし、（式２２−１）又は（式２２−２）は数値計算が困難であるため、同期化力係数行列Ｋのゼロ固有値に相当する非ゼロのｎ次元固有ベクトルｖ、ｎ次元固有ベクトルｖの大きさを用いて、同期化力係数行列Ｋが非正則行列である必要十分条件（第１条件）を、次の（式２３）及び（式２４−１）のように定式化した。  
2010246289.tif 000055
尚、ｎ次元固有ベクトルｖは、方向のみが意味を持ち、大きさは任意である。従って、スカラーｑは、任意の数値として設定ができる。以下では、スカラーｑは便宜上１とし、（式２４−１）を次の（式２４−２）として表す。  
2010246289.tif 000056
【００４２】
第２条件が成立する必要十分条件は、スカラー係数ｋｓを用いると、次の（式２５−１）により定式化される。または、慣性中心ωｏを考慮して、次の（式２５−２）により定式化してもよい。以下では、（式２５−１）を前提に説明する。  
2010246289.tif 000057
【００４３】
以上より、特異点Ｓの条件は、（式２３）、（式２４−２）、（式２５−１）全てが成立することである。ここで、未知数は、ｎ次元固有ベクトルｖ、ｎ次元ベクトルω、スカラー係数ｋｓである。そして、これらの特異点Ｓの条件を臨界軌跡３の終点ｘｍ＋１の束縛条件として使用する。このとき、終点ｘｍ＋１の誤差ベクトルμｍ＋１を、次の（式２６）のとおり定義する。  
2010246289.tif 000058
【００４４】
また、（式２６）を用いた二乗誤差関数は、次の（式２７）により定義される。尚、Ｗは対角要素が正である正方の対角行列（対角要素のみ正で非対角要素は全てゼロの行列）である。  
2010246289.tif 000059
さらに、（式２７）の二乗誤差関数を（式１６）に示される目的関数Ｏに加算した目的関数Ｅが、次の（式２８）により定義される。  
2010246289.tif 000060
【００４５】
そして、上記の（式１７）等による最適化問題を、（式２８）の目的関数Ｅを最小化させる未知変数τ、ｘ１〜ｘｍ＋１、ε、ｖ、ｋｓを求める最適化問題へと帰着させる。即ち、（式１７）は、次の（式２９）に変形される。  
2010246289.tif 000061
【００４６】
（式２９）の最適化問題では、特異点Ｓの条件として（式２３）、（式２４−２）、(式２５−１)を全て用いて、終点Ｘｍ＋１の誤差ベクトルμｍ＋１を（式２６）により定義したが、μｍ＋１は特異点Ｓの等価条件を用いても定義できる。この等価条件とは、(式２３)、(式２４−２)、（式２５−１）より変数ｖとｋｓを消去して得られる（式Ｂ１）の条件を指す。  
2010246289.tif 000062
（式Ｂ１）に基づいて、終点Ｘｍ＋１の誤差ベクトルμｍ＋１を（式Ｂ２）で定義する。  
2010246289.tif 000063
上記で誤差ベクトルを定義した後は、前述の手順と全く同様に（式２７）および（式２８）により、２乗誤差関数および目的関数Ｅが定義される。このとき、（式Ｂ２）による特異点Ｓの条件を用いて、（式２８）の目的関数を最小化させる未知変数τ、ｘ１〜Ｘｍ＋１、εを求める最適化問題に帰着できる。すなわち、この場合も（式２９）と全く同じ形の（式Ｂ３）で記述できる。  
2010246289.tif 000064
ここで、（式２９）と(式Ｂ３)との相違点は、誤差関数と未知変数のみである。前者の誤差関数は（式２６）、後者の(式Ｂ２)で定義される。前者の未知変数は、τ、ｘ１〜Ｘｍ＋１、ε、ｖ、ｋｓであり、一方後者では変数ｖ、ｋｓが消去されたため、未知変数はτ、ｘ１〜Ｘｍ＋１、εとなっている。
【００４７】
以下、図３に示すフローチャートを用いて、本発明の実施の形態の臨界故障除去時間算出方法の全体的な流れを説明する。
【００４８】
まず、臨界故障除去時間の算出を行う多機電力系統を表したモデルを特定する。これにより、多次元関数ｆや、多次元状態変数ｘｋ（０≦ｋ≦ｍ＋１）が特定される（Ｓ３００）。つぎに、臨界軌跡３の始点ｘ０の束縛条件として（式１９）を設定し、臨界軌跡３の終点ｘｍ+1の束縛条件として特異点Ｓの条件（式２１）、（式２２−２）、（式２３−１）を設定する。これにより、（式２９）に示される最適化問題が特定される（Ｓ３０１）。
【００４９】
つぎに、（式２９）に示される最適化問題の解探索、つまり（式２８）に示される目的関数Ｅの最小化を実行する（Ｓ３０２）。この結果、臨界故障除去時の状態から支配的不安定平衡状態に至るまでの臨界軌跡３に沿って誤差を分散させながら、臨界軌跡２の始点ｘ０等を得る。そして、この得られた始点ｘ０に対して（式１９）の等式を適用することで臨界故障除去時間τを算出する（Ｓ３０３）。
【００５０】
尚、（式２８）に示される目的関数Ｅを最小化させる解探索を実行する過程で、目的関数Ｅと多次元状態変数ｘｋ（０≦ｋ≦ｍ＋１）との推移を所定のメモリに記憶しておき、解探索後に当該メモリに記憶された目的関数Ｅと多次元状態変数ｘk（０≦ｋ≦ｍ＋１）を時系列に表示等することにより、探索経路の確認や局所最適解に陥っていないか否かの確認等を視覚的に行ってもよい。
【００５１】
また、（式２９）の解探索の際に、次の（式３０）による不等式による制約条件を更に付加することができる。ｇ（ｘ）は、このモデルに付加される発電機のリミッタ等を意味する関数を採用することができる。これにより、最適化問題の定式化の精度を向上させることができる。  
2010246289.tif 000065
【００５２】
＝＝＝Ａｎｄｅｒｓｏｎ＆Ｆｏｕａｄの例題系統モデルによる性能評価＝＝＝
本発明の実施の形態の臨界故障除去時間算出方法を、Anderson & Fouadの３機９母線系統（ＡＦ９）のモデルの非線形現象へ適用することで、本発明の実施の形態の臨界故障除去時間算出方法の性能評価を行った。図４は、ＡＦ９のモデルを示した概念図である。
【００５３】
ＡＦ９のモデルを用いた性能評価にあたり、時刻ｔ＝０において複数の母線のうちの所定の母線近傍の故障点で３相地絡故障が発生し、時刻ｔ＝τにおいて２つの送電線のうちの１つを開放することにより故障が除去される場合を想定し、この場合において、予め定めておいた臨界故障除去時間ＣＣＴの厳密値（ｓ）と本発明の実施の形態の臨界故障除去時間算出方法によって算出された臨界故障除去時間ＣＣＴ（ｓ）との比較を行った。表１にＡＦ９のモデルを用いた性能評価の結果を示す。
【表１】
2010246289.tif 000066
【００５４】
表１は、臨界軌跡３の分割数ｍを１０に設定した場合の各故障ケースの結果を示している。「故障点」は、図４中のアルファベットが表す故障発生箇所を示しており、「故障線路」は、故障発生箇所を含んだ図４の数字が表す線路番号（両端母線番号）を示している。「ＣＣＴ厳密値」は、従来のシミュレーション法を繰り返すことで得られる、安定となる場合の臨界故障除去時間と不安定となる場合の臨界故障除去時間との両方を示しており、両方の時間の間（例えば、第１行目の場合、０．３４−０．３５とは、０．３４（ｓ）と０．３５（ｓ）との間）に臨界故障除去時間ＣＣＴの厳密値が存在することを意味している。「繰り返し回数」は、最適化問題をニュートン法で解いた際の反復回数（ただし、収束判定基準は10-5）を示しており、「計算時間」は計算に要したＣＰＵ時間を示している。
【００５５】
表１によれば、例えば、第１母線のＡ点で３相地絡故障が発生した場合、第１〜第４母線に対応する（各２本のうちの１本の）送電線を開放することにより故障を除去する場合の臨界故障除去時間ＣＣＴは、手法１では０．３４４６（ｓ）、手法２では０．３４２０（ｓ）と算出された。この結果は、従来のシミュレーション法により得られる臨界故障除去時間ＣＣＴの厳密値の範囲内（０．３４（ｓ）−０．３５（ｓ））に整合している。また、他の故障点Ｂ〜Ｉの場合に算出された臨界故障除去時間ＣＣＴについても同様のことがいえる。従って、本発明の実施の形態の臨界故障除去時間算出方法は、ＡＦ９のモデルにおいて、厳密且つ高速であり、安定した計算性能を示すことが分かる。
【００５６】
＝＝＝ＩＥＥＥの例題系統による性能評価＝＝＝
本発明の実施の形態の臨界故障除去時間算出方法を、ＩＥＥＥ６機３０母線系統（ＩＥＥＥ３０）のモデルの非線形現象へ適用することで、本発明の実施の形態の臨界故障除去時間算出方法の性能評価を行った。図５は、ＩＥＥＥ３０のモデルを示した概念図である。
【００５７】
ＩＥＥＥ３０のモデルを用いた性能評価にあたり、上記のＡＦ９のモデルを用いた性能評価の場合と同様に、時刻ｔ＝０において複数の母線のうちの所定の母線近傍の故障点で３相地絡故障が発生し、時刻ｔ＝τにおいて２つの送電線のうちの１つを開放することにより故障が除去される場合を想定し、この場合において、予め定めておいた臨界故障除去時間ＣＣＴの厳密値（ｓ）と本発明の実施の形態の臨界故障除去時間算出方法によって算出された臨界故障除去時間ＣＣＴ（ｓ）との比較を行った。表２にＩＥＥＥ３０のモデルを用いた性能評価の結果を示す。
【表２】
2010246289.tif 000067
【００５８】
表２は、臨界軌跡３の分割数ｍを１０に設定した場合の各故障ケースの結果を示している。「故障ケース」は、各故障ケースの識別子を示しており、「故障線路」は、故障発生箇所を含んだ図５の数字が表す線路番号（両端母線番号）を示している。「ＣＣＴ厳密値」は、従来のシミュレーション法を繰り返すことで得られる臨界故障除去時間ＣＣＴの厳密値が存在する範囲を示す。「繰り返し回数」は、最適化問題をニュートン法で解いた際の反復回数（ただし、収束判定基準は10-5）を示しており、「計算時間」は計算に要したＣＰＵ時間を示している。
【００５９】
表２によれば、例えば、第１母線の近傍で３相地絡故障が発生した場合、第１〜第２母線に対応する（各２本のうちの１本の）送電線を開放することにより故障を除去する場合の臨界故障除去時間ＣＣＴは、手法１では０．７８６４（ｃ）、手法２では０．７８８０（ｓ）と算出された。この結果は、従来のシミュレーション法により得られる臨界故障除去時間ＣＣＴの厳密値の範囲内（０．７８（ｓ）−０．７９（ｓ））に整合している。また、他の故障点Ｂ〜Ｉの場合に算出された臨界故障除去時間ＣＣＴについても同様のことがいえる。従って、本発明の実施の形態の臨界故障除去時間算出方法は、ＩＥＥＥ３０のモデルにおいて、厳密且つ高速であり、安定した計算性能を示すことが分かる。
【００６０】
＝＝＝情報処理装置、プログラム＝＝＝
本発明の実施の形態の臨界故障除去時間算出方法は、図４に示されるＣＰＵ１１０、メモリ１２０、キーボード、マウス等の入力装置１３０、液晶ディスプレイ、ブラウン管等の表示装置１４０を備える情報処理装置１００によって実行される。尚、情報処理装置１００は、例えば、電力系統の運用に携わる作業者が操作するＰＳＳとして機能するコンピュータやワークステーション等である。ＰＳＳは、例えば特開平１０−２８３２６号公報に開示された周知の装置である。
【００６１】
尚、メモリ１２０には、前述した臨界故障除去時間の算出を行うためのプログラム１２５が格納されており、ＣＰＵ１１０は、メモリ１２０からプログラム１２５を呼び出して実行することになる。また、メモリ１２０は、関数ｆに関する情報や、状態変数ベクトルｘ０〜ｘｍ＋１や誤差ベクトル等の中間データ等も記憶される。更に、メモリ１２０は、例えば、前述したニュートン法のプログラムや、このニュートン法を用いて未知変数ｘ０〜ｘｍ〜ｘｕ、τ、ε等の最適化を実施するプログラム等からなるプログラム又はプログラム群（臨界故障除去時間算出プログラム）を記憶している必要がある。尚、メモリ１２０は、本発明の実施の形態の情報処理装置１００の外部に備えてもよい。また、上記プログラム又はプログラム群は、例えば、ＤＶＤやＣＤ等の光ディスク、或いは、ＭＯやフロッピーディスク等の磁気ディスクといった、情報処理装置１００により読み取り可能な記録媒体に記録されたものであってもよい。
【００６２】
以上、本発明の実施の形態によれば、（式１６）の多元連立方程式を、臨界故障除去時の状態から支配的不安定平衡状態に至るまで逐次解くのではなく、（式２９）の解探索により一括して解くことによって、臨界軌跡３に対するより高精度な計算結果を与えることになり、ひいては臨界故障除去時間τの精度を向上化できる。
【００６３】
また、（式２９）の解探索の際に、終点ｘｍ＋１の束縛条件として従来の試行錯誤で指定する必要があった定数ｘuよりも安定性、確定性且つ普遍性を具備した特異点Ｓを導入したことにより、最適化問題の制約条件の精度を向上させるとともにその解探索の際に用いる目的関数の更なる改善を図り、臨界故障除去時間τの算出をより高精度に且つより迅速に行うことが可能となった。
【００６４】
前述した実施の形態は、本発明の理解を容易にするためのものであり、本発明を限定して解釈するためのものではない。本発明は、その趣旨を逸脱することなく変更、改良されるとともに、本発明にはその等価物も含まれる。
【００６５】
例えば、（式９）の代わりに
2010246289.tif 000068
【００６６】
を用い、
2010246289.tif 000069
【００６７】
を拘束条件として追加する。すなわち、
2010246289.tif 000070
【００６８】
を目的関数に追加する。そして、元の最小化問題の未知数（決定変数）に、さらにｙｋ、（ｋ＝０，・・・，ｍ＋１）を加えて、目的関数の最小化を行うこととしてもよい（ｙ０も未知数）。
【符号の説明】
【００６９】
１  故障軌跡
２  故障が除去された後に安定状態に戻ることが可能な電力系統の状態を示す軌跡
３  臨界軌跡
４  故障が除去された後に安定状態に戻ることが不可能な電力系統の状態を示す軌跡
１００  情報処理装置
１１０  ＣＰＵ
１２０  メモリ
１２５  プログラム
１３０  入力装置
１４０  表示装置
【図１】
2010246289.tif 000086
【図２】
2010246289.tif 000087
【図３】
2010246289.tif 000088
【図４】
2010246289.tif 000089
【図５】
2010246289.tif 000090
【図６】
2010246289.tif 000091